newmat3.cpp
Go to the documentation of this file.
1 
6 
7 // Copyright (C) 1991,2,3,4: R B Davies
8 
9 //#define WANT_STREAM
10 
11 #include "include.h"
12 
13 #include "newmat.h"
14 #include "newmatrc.h"
15 
16 #ifdef use_namespace
17 namespace NEWMAT {
18 #endif
19 
20 
21 #ifdef DO_REPORT
22 #define REPORT { static ExeCounter ExeCount(__LINE__,3); ++ExeCount; }
23 #else
24 #define REPORT {}
25 #endif
26 
27 //#define MONITOR(what,storage,store)
28 // { cout << what << " " << storage << " at " << (long)store << "\n"; }
29 
30 #define MONITOR(what,store,storage) {}
31 
32 
33 // Control bits codes for GetRow, GetCol, RestoreRow, RestoreCol
34 //
35 // LoadOnEntry:
36 // Load data into MatrixRow or Col dummy array under GetRow or GetCol
37 // StoreOnExit:
38 // Restore data to original matrix under RestoreRow or RestoreCol
39 // DirectPart:
40 // Load or restore only part directly stored; must be set with StoreOnExit
41 // Still have decide how to handle this with symmetric
42 // StoreHere:
43 // used in columns only - store data at supplied storage address;
44 // used for GetCol, NextCol & RestoreCol. No need to fill out zeros
45 // HaveStore:
46 // dummy array has been assigned (internal use only).
47 
48 // For symmetric matrices, treat columns as rows unless StoreHere is set;
49 // then stick to columns as this will give better performance for doing
50 // inverses
51 
52 // How components are used:
53 
54 // Use rows wherever possible in preference to columns
55 
56 // Columns without StoreHere are used in in-exact transpose, sum column,
57 // multiply a column vector, and maybe in future access to column,
58 // additional multiply functions, add transpose
59 
60 // Columns with StoreHere are used in exact transpose (not symmetric matrices
61 // or vectors, load only)
62 
63 // Columns with MatrixColX (Store to full column) are used in inverse and solve
64 
65 // Functions required for each matrix class
66 
67 // GetRow(MatrixRowCol& mrc)
68 // GetCol(MatrixRowCol& mrc)
69 // GetCol(MatrixColX& mrc)
70 // RestoreRow(MatrixRowCol& mrc)
71 // RestoreCol(MatrixRowCol& mrc)
72 // RestoreCol(MatrixColX& mrc)
73 // NextRow(MatrixRowCol& mrc)
74 // NextCol(MatrixRowCol& mrc)
75 // NextCol(MatrixColX& mrc)
76 
77 // The Restore routines assume StoreOnExit has already been checked
78 // Defaults for the Next routines are given below
79 // Assume cannot have both !DirectPart && StoreHere for MatrixRowCol routines
80 
81 
82 // Default NextRow and NextCol:
83 // will work as a default but need to override NextRow for efficiency
84 
86 {
87  REPORT
88  if (+(mrc.cw*StoreOnExit)) { REPORT this->RestoreRow(mrc); }
89  mrc.rowcol++;
90  if (mrc.rowcol<nrows_val) { REPORT this->GetRow(mrc); }
91  else { REPORT mrc.cw -= StoreOnExit; }
92 }
93 
95 {
96  REPORT // 423
97  if (+(mrc.cw*StoreOnExit)) { REPORT this->RestoreCol(mrc); }
98  mrc.rowcol++;
99  if (mrc.rowcol<ncols_val) { REPORT this->GetCol(mrc); }
100  else { REPORT mrc.cw -= StoreOnExit; }
101 }
102 
104 {
105  REPORT // 423
106  if (+(mrc.cw*StoreOnExit)) { REPORT this->RestoreCol(mrc); }
107  mrc.rowcol++;
108  if (mrc.rowcol<ncols_val) { REPORT this->GetCol(mrc); }
109  else { REPORT mrc.cw -= StoreOnExit; }
110 }
111 
112 
113 // routines for matrix
114 
116 {
117  REPORT
118  mrc.skip=0; mrc.storage=mrc.length=ncols_val;
119  mrc.data=store+mrc.rowcol*ncols_val;
120 }
121 
122 
124 {
125  REPORT
126  mrc.skip=0; mrc.storage=mrc.length=nrows_val;
127  if ( ncols_val==1 && !(mrc.cw*StoreHere) ) // ColumnVector
128  { REPORT mrc.data=store; }
129  else
130  {
131  Real* ColCopy;
132  if ( !(mrc.cw*(HaveStore+StoreHere)) )
133  {
134  REPORT
135  ColCopy = new Real [nrows_val]; MatrixErrorNoSpace(ColCopy);
136  MONITOR_REAL_NEW("Make (MatGetCol)",nrows_val,ColCopy)
137  mrc.data = ColCopy; mrc.cw += HaveStore;
138  }
139  else { REPORT ColCopy = mrc.data; }
140  if (+(mrc.cw*LoadOnEntry))
141  {
142  REPORT
143  Real* Mstore = store+mrc.rowcol; int i=nrows_val;
144  //while (i--) { *ColCopy++ = *Mstore; Mstore+=ncols_val; }
145  if (i) for (;;)
146  { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore+=ncols_val; }
147  }
148  }
149 }
150 
152 {
153  REPORT
154  mrc.skip=0; mrc.storage=nrows_val; mrc.length=nrows_val;
155  if (+(mrc.cw*LoadOnEntry))
156  {
157  REPORT Real* ColCopy = mrc.data;
158  Real* Mstore = store+mrc.rowcol; int i=nrows_val;
159  //while (i--) { *ColCopy++ = *Mstore; Mstore+=ncols_val; }
160  if (i) for (;;)
161  { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore+=ncols_val; }
162  }
163 }
164 
166 {
167  // always check StoreOnExit before calling RestoreCol
168  REPORT // 429
169  if (+(mrc.cw*HaveStore))
170  {
171  REPORT // 426
172  Real* Mstore = store+mrc.rowcol; int i=nrows_val;
173  Real* Cstore = mrc.data;
174  // while (i--) { *Mstore = *Cstore++; Mstore+=ncols_val; }
175  if (i) for (;;)
176  { *Mstore = *Cstore++; if (!(--i)) break; Mstore+=ncols_val; }
177  }
178 }
179 
181 {
182  REPORT
183  Real* Mstore = store+mrc.rowcol; int i=nrows_val; Real* Cstore = mrc.data;
184  // while (i--) { *Mstore = *Cstore++; Mstore+=ncols_val; }
185  if (i) for (;;)
186  { *Mstore = *Cstore++; if (!(--i)) break; Mstore+=ncols_val; }
187 }
188 
189 void Matrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrMat(); } // 1808
190 
192 {
193  REPORT // 632
194  if (+(mrc.cw*StoreOnExit)) { REPORT RestoreCol(mrc); }
195  mrc.rowcol++;
196  if (mrc.rowcol<ncols_val)
197  {
198  if (+(mrc.cw*LoadOnEntry))
199  {
200  REPORT
201  Real* ColCopy = mrc.data;
202  Real* Mstore = store+mrc.rowcol; int i=nrows_val;
203  //while (i--) { *ColCopy++ = *Mstore; Mstore+=ncols_val; }
204  if (i) for (;;)
205  { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore+=ncols_val; }
206  }
207  }
208  else { REPORT mrc.cw -= StoreOnExit; }
209 }
210 
212 {
213  REPORT
214  if (+(mrc.cw*StoreOnExit)) { REPORT RestoreCol(mrc); }
215  mrc.rowcol++;
216  if (mrc.rowcol<ncols_val)
217  {
218  if (+(mrc.cw*LoadOnEntry))
219  {
220  REPORT
221  Real* ColCopy = mrc.data;
222  Real* Mstore = store+mrc.rowcol; int i=nrows_val;
223  // while (i--) { *ColCopy++ = *Mstore; Mstore+=ncols_val; }
224  if (i) for (;;)
225  { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore+=ncols_val; }
226  }
227  }
228  else { REPORT mrc.cw -= StoreOnExit; }
229 }
230 
231 // routines for diagonal matrix
232 
234 {
235  REPORT
236  mrc.skip=mrc.rowcol; mrc.storage=1;
237  mrc.data=store+mrc.skip; mrc.length=ncols_val;
238 }
239 
241 {
242  REPORT
243  mrc.skip=mrc.rowcol; mrc.storage=1; mrc.length=nrows_val;
244  if (+(mrc.cw*StoreHere)) // should not happen
245  Throw(InternalException("DiagonalMatrix::GetCol(MatrixRowCol&)"));
246  else { REPORT mrc.data=store+mrc.skip; }
247  // not accessed
248 }
249 
251 {
252  REPORT
253  mrc.skip=mrc.rowcol; mrc.storage=1; mrc.length=nrows_val;
254  mrc.data = mrc.store+mrc.skip;
255  *(mrc.data)=*(store+mrc.skip);
256 }
257 
259  // 800
260 
262  // not accessed
263 
265 {
266  REPORT
267  if (+(mrc.cw*StoreOnExit))
268  { REPORT *(store+mrc.rowcol)=*(mrc.data); }
269  mrc.IncrDiag();
270  int t1 = +(mrc.cw*LoadOnEntry);
271  if (t1 && mrc.rowcol < ncols_val)
272  { REPORT *(mrc.data)=*(store+mrc.rowcol); }
273 }
274 
275 // routines for upper triangular matrix
276 
278 {
279  REPORT
280  int row = mrc.rowcol; mrc.skip=row; mrc.length=ncols_val;
281  mrc.storage=ncols_val-row; mrc.data=store+(row*(2*ncols_val-row+1))/2;
282 }
283 
284 
286 {
287  REPORT
288  mrc.skip=0; int i=mrc.rowcol+1; mrc.storage=i;
289  mrc.length=nrows_val; Real* ColCopy;
290  if ( !(mrc.cw*(StoreHere+HaveStore)) )
291  {
292  REPORT // not accessed
293  ColCopy = new Real [nrows_val]; MatrixErrorNoSpace(ColCopy);
294  MONITOR_REAL_NEW("Make (UT GetCol)",nrows_val,ColCopy)
295  mrc.data = ColCopy; mrc.cw += HaveStore;
296  }
297  else { REPORT ColCopy = mrc.data; }
298  if (+(mrc.cw*LoadOnEntry))
299  {
300  REPORT
301  Real* Mstore = store+mrc.rowcol; int j = ncols_val;
302  // while (i--) { *ColCopy++ = *Mstore; Mstore += --j; }
303  if (i) for (;;)
304  { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore += --j; }
305  }
306 }
307 
309 {
310  REPORT
311  mrc.skip=0; int i=mrc.rowcol+1; mrc.storage=i;
312  mrc.length=nrows_val;
313  if (+(mrc.cw*LoadOnEntry))
314  {
315  REPORT
316  Real* ColCopy = mrc.data;
317  Real* Mstore = store+mrc.rowcol; int j = ncols_val;
318  // while (i--) { *ColCopy++ = *Mstore; Mstore += --j; }
319  if (i) for (;;)
320  { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore += --j; }
321  }
322 }
323 
325 {
326  REPORT
327  Real* Mstore = store+mrc.rowcol; int i=mrc.rowcol+1; int j = ncols_val;
328  Real* Cstore = mrc.data;
329  // while (i--) { *Mstore = *Cstore++; Mstore += --j; }
330  if (i) for (;;)
331  { *Mstore = *Cstore++; if (!(--i)) break; Mstore += --j; }
332 }
333 
335  // 722
336 
337 // routines for lower triangular matrix
338 
340 {
341  REPORT
342  int row=mrc.rowcol; mrc.skip=0; mrc.storage=row+1; mrc.length=ncols_val;
343  mrc.data=store+(row*(row+1))/2;
344 }
345 
347 {
348  REPORT
349  int col=mrc.rowcol; mrc.skip=col; mrc.length=nrows_val;
350  int i=nrows_val-col; mrc.storage=i; Real* ColCopy;
351  if ( +(mrc.cw*(StoreHere+HaveStore)) )
352  { REPORT ColCopy = mrc.data; }
353  else
354  {
355  REPORT // not accessed
356  ColCopy = new Real [nrows_val]; MatrixErrorNoSpace(ColCopy);
357  MONITOR_REAL_NEW("Make (LT GetCol)",nrows_val,ColCopy)
358  mrc.cw += HaveStore; mrc.data = ColCopy;
359  }
360 
361  if (+(mrc.cw*LoadOnEntry))
362  {
363  REPORT
364  Real* Mstore = store+(col*(col+3))/2;
365  // while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
366  if (i) for (;;)
367  { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore += ++col; }
368  }
369 }
370 
372 {
373  REPORT
374  int col=mrc.rowcol; mrc.skip=col; mrc.length=nrows_val;
375  int i=nrows_val-col; mrc.storage=i; mrc.data = mrc.store + col;
376 
377  if (+(mrc.cw*LoadOnEntry))
378  {
379  REPORT Real* ColCopy = mrc.data;
380  Real* Mstore = store+(col*(col+3))/2;
381  // while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
382  if (i) for (;;)
383  { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore += ++col; }
384  }
385 }
386 
388 {
389  REPORT
390  int col=mrc.rowcol; Real* Cstore = mrc.data;
391  Real* Mstore = store+(col*(col+3))/2; int i=nrows_val-col;
392  //while (i--) { *Mstore = *Cstore++; Mstore += ++col; }
393  if (i) for (;;)
394  { *Mstore = *Cstore++; if (!(--i)) break; Mstore += ++col; }
395 }
396 
398  //712
399 
400 // routines for symmetric matrix
401 
403 {
404  REPORT //571
405  mrc.skip=0; int row=mrc.rowcol; mrc.length=ncols_val;
406  if (+(mrc.cw*DirectPart))
407  { REPORT mrc.storage=row+1; mrc.data=store+(row*(row+1))/2; }
408  else
409  {
410  // do not allow StoreOnExit and !DirectPart
411  if (+(mrc.cw*StoreOnExit))
412  Throw(InternalException("SymmetricMatrix::GetRow(MatrixRowCol&)"));
413  mrc.storage=ncols_val; Real* RowCopy;
414  if (!(mrc.cw*HaveStore))
415  {
416  REPORT
417  RowCopy = new Real [ncols_val]; MatrixErrorNoSpace(RowCopy);
418  MONITOR_REAL_NEW("Make (SymGetRow)",ncols_val,RowCopy)
419  mrc.cw += HaveStore; mrc.data = RowCopy;
420  }
421  else { REPORT RowCopy = mrc.data; }
422  if (+(mrc.cw*LoadOnEntry))
423  {
424  REPORT // 544
425  Real* Mstore = store+(row*(row+1))/2; int i = row;
426  while (i--) *RowCopy++ = *Mstore++;
427  i = ncols_val-row;
428  // while (i--) { *RowCopy++ = *Mstore; Mstore += ++row; }
429  if (i) for (;;)
430  { *RowCopy++ = *Mstore; if (!(--i)) break; Mstore += ++row; }
431  }
432  }
433 }
434 
436 {
437  // do not allow StoreHere
438  if (+(mrc.cw*StoreHere))
439  Throw(InternalException("SymmetricMatrix::GetCol(MatrixRowCol&)"));
440 
441  int col=mrc.rowcol; mrc.length=nrows_val;
442  REPORT
443  mrc.skip=0;
444  if (+(mrc.cw*DirectPart)) // actually get row ??
445  { REPORT mrc.storage=col+1; mrc.data=store+(col*(col+1))/2; }
446  else
447  {
448  // do not allow StoreOnExit and !DirectPart
449  if (+(mrc.cw*StoreOnExit))
450  Throw(InternalException("SymmetricMatrix::GetCol(MatrixRowCol&)"));
451 
452  mrc.storage=ncols_val; Real* ColCopy;
453  if ( +(mrc.cw*HaveStore)) { REPORT ColCopy = mrc.data; }
454  else
455  {
456  REPORT // not accessed
457  ColCopy = new Real [ncols_val]; MatrixErrorNoSpace(ColCopy);
458  MONITOR_REAL_NEW("Make (SymGetCol)",ncols_val,ColCopy)
459  mrc.cw += HaveStore; mrc.data = ColCopy;
460  }
461  if (+(mrc.cw*LoadOnEntry))
462  {
463  REPORT
464  Real* Mstore = store+(col*(col+1))/2; int i = col;
465  while (i--) *ColCopy++ = *Mstore++;
466  i = ncols_val-col;
467  // while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
468  if (i) for (;;)
469  { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore += ++col; }
470  }
471  }
472 }
473 
475 {
476  int col=mrc.rowcol; mrc.length=nrows_val;
477  if (+(mrc.cw*DirectPart))
478  {
479  REPORT
480  mrc.skip=col; int i=nrows_val-col; mrc.storage=i;
481  mrc.data = mrc.store+col;
482  if (+(mrc.cw*LoadOnEntry))
483  {
484  REPORT // not accessed
485  Real* ColCopy = mrc.data;
486  Real* Mstore = store+(col*(col+3))/2;
487  // while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
488  if (i) for (;;)
489  { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore += ++col; }
490  }
491  }
492  else
493  {
494  REPORT
495  // do not allow StoreOnExit and !DirectPart
496  if (+(mrc.cw*StoreOnExit))
497  Throw(InternalException("SymmetricMatrix::GetCol(MatrixColX&)"));
498 
499  mrc.skip=0; mrc.storage=ncols_val;
500  if (+(mrc.cw*LoadOnEntry))
501  {
502  REPORT
503  Real* ColCopy = mrc.data;
504  Real* Mstore = store+(col*(col+1))/2; int i = col;
505  while (i--) *ColCopy++ = *Mstore++;
506  i = ncols_val-col;
507  // while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
508  if (i) for (;;)
509  { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore += ++col; }
510  }
511  }
512 }
513 
514 // Do not need RestoreRow because we do not allow !DirectPart && StoreOnExit
515 
517 {
518  REPORT
519  // Really do restore column
520  int col=mrc.rowcol; Real* Cstore = mrc.data;
521  Real* Mstore = store+(col*(col+3))/2; int i = nrows_val-col;
522  // while (i--) { *Mstore = *Cstore++; Mstore+= ++col; }
523  if (i) for (;;)
524  { *Mstore = *Cstore++; if (!(--i)) break; Mstore+= ++col; }
525 
526 }
527 
528 // routines for row vector
529 
531 {
532  REPORT
533  // do not allow StoreHere
534  if (+(mrc.cw*StoreHere))
535  Throw(InternalException("RowVector::GetCol(MatrixRowCol&)"));
536 
537  mrc.skip=0; mrc.storage=1; mrc.length=nrows_val;
538  mrc.data = store+mrc.rowcol;
539 
540 }
541 
543 {
544  REPORT
545  mrc.skip=0; mrc.storage=1; mrc.length=nrows_val;
546  if (mrc.cw >= LoadOnEntry)
547  { REPORT *(mrc.data) = *(store+mrc.rowcol); }
548 
549 }
550 
552 { REPORT mrc.rowcol++; mrc.data++; }
553 
555 {
556  if (+(mrc.cw*StoreOnExit)) { REPORT *(store+mrc.rowcol)=*(mrc.data); }
557 
558  mrc.rowcol++;
559  if (mrc.rowcol < ncols_val)
560  {
561  if (+(mrc.cw*LoadOnEntry)) { REPORT *(mrc.data)=*(store+mrc.rowcol); }
562  }
563  else { REPORT mrc.cw -= StoreOnExit; }
564 }
565 
567  { REPORT *(store+mrc.rowcol)=*(mrc.data); } // not accessed
568 
569 
570 // routines for band matrices
571 
573 {
574  REPORT
575  int r = mrc.rowcol; int w = lower_val+1+upper_val; mrc.length=ncols_val;
576  int s = r-lower_val;
577  if (s<0) { mrc.data = store+(r*w-s); w += s; s = 0; }
578  else mrc.data = store+r*w;
579  mrc.skip = s; s += w-ncols_val; if (s>0) w -= s; mrc.storage = w;
580 }
581 
582 // should make special versions of this for upper and lower band matrices
583 
585 {
586  REPORT
587  int r = ++mrc.rowcol;
588  if (r<=lower_val) { mrc.storage++; mrc.data += lower_val+upper_val; }
589  else { mrc.skip++; mrc.data += lower_val+upper_val+1; }
590  if (r>=ncols_val-upper_val) mrc.storage--;
591 }
592 
594 {
595  REPORT
596  int c = mrc.rowcol; int n = lower_val+upper_val; int w = n+1;
597  mrc.length=nrows_val; Real* ColCopy;
598  int b; int s = c-upper_val;
599  if (s<=0) { w += s; s = 0; b = c+lower_val; } else b = s*w+n;
600  mrc.skip = s; s += w-nrows_val; if (s>0) w -= s; mrc.storage = w;
601  if ( +(mrc.cw*(StoreHere+HaveStore)) )
602  { REPORT ColCopy = mrc.data; }
603  else
604  {
605  REPORT
606  ColCopy = new Real [n+1]; MatrixErrorNoSpace(ColCopy);
607  MONITOR_REAL_NEW("Make (BMGetCol)",n+1,ColCopy)
608  mrc.cw += HaveStore; mrc.data = ColCopy;
609  }
610 
611  if (+(mrc.cw*LoadOnEntry))
612  {
613  REPORT
614  Real* Mstore = store+b;
615  // while (w--) { *ColCopy++ = *Mstore; Mstore+=n; }
616  if (w) for (;;)
617  { *ColCopy++ = *Mstore; if (!(--w)) break; Mstore+=n; }
618  }
619 }
620 
622 {
623  REPORT
624  int c = mrc.rowcol; int n = lower_val+upper_val; int w = n+1;
625  mrc.length=nrows_val; int b; int s = c-upper_val;
626  if (s<=0) { w += s; s = 0; b = c+lower_val; } else b = s*w+n;
627  mrc.skip = s; s += w-nrows_val; if (s>0) w -= s; mrc.storage = w;
628  mrc.data = mrc.store+mrc.skip;
629 
630  if (+(mrc.cw*LoadOnEntry))
631  {
632  REPORT
633  Real* ColCopy = mrc.data; Real* Mstore = store+b;
634  // while (w--) { *ColCopy++ = *Mstore; Mstore+=n; }
635  if (w) for (;;)
636  { *ColCopy++ = *Mstore; if (!(--w)) break; Mstore+=n; }
637  }
638 }
639 
641 {
642  REPORT
643  int c = mrc.rowcol; int n = lower_val+upper_val; int s = c-upper_val;
644  Real* Mstore = store + ((s<=0) ? c+lower_val : s*n+s+n);
645  Real* Cstore = mrc.data;
646  int w = mrc.storage;
647  // while (w--) { *Mstore = *Cstore++; Mstore += n; }
648  if (w) for (;;)
649  { *Mstore = *Cstore++; if (!(--w)) break; Mstore += n; }
650 }
651 
652 // routines for symmetric band matrix
653 
655 {
656  REPORT
657  int r=mrc.rowcol; int s = r-lower_val; int w1 = lower_val+1; int o = r*w1;
658  mrc.length = ncols_val;
659  if (s<0) { w1 += s; o -= s; s = 0; }
660  mrc.skip = s;
661 
662  if (+(mrc.cw*DirectPart))
663  { REPORT mrc.data = store+o; mrc.storage = w1; }
664  else
665  {
666  // do not allow StoreOnExit and !DirectPart
667  if (+(mrc.cw*StoreOnExit))
668  Throw(InternalException("SymmetricBandMatrix::GetRow(MatrixRowCol&)"));
669  int w = w1+lower_val; s += w-ncols_val; Real* RowCopy;
670  if (s>0) w -= s; mrc.storage = w; int w2 = w-w1;
671  if (!(mrc.cw*HaveStore))
672  {
673  REPORT
674  RowCopy = new Real [2*lower_val+1]; MatrixErrorNoSpace(RowCopy);
675  MONITOR_REAL_NEW("Make (SmBGetRow)",2*lower_val+1,RowCopy)
676  mrc.cw += HaveStore; mrc.data = RowCopy;
677  }
678  else { REPORT RowCopy = mrc.data; }
679 
680  if (+(mrc.cw*LoadOnEntry) && ncols_val > 0)
681  {
682  REPORT
683  Real* Mstore = store+o;
684  while (w1--) *RowCopy++ = *Mstore++;
685  Mstore--;
686  while (w2--) { Mstore += lower_val; *RowCopy++ = *Mstore; }
687  }
688  }
689 }
690 
692 {
693  // do not allow StoreHere
694  if (+(mrc.cw*StoreHere))
695  Throw(InternalException("SymmetricBandMatrix::GetCol(MatrixRowCol&)"));
696 
697  int c=mrc.rowcol; int w1 = lower_val+1; mrc.length=nrows_val;
698  REPORT
699  int s = c-lower_val; int o = c*w1;
700  if (s<0) { w1 += s; o -= s; s = 0; }
701  mrc.skip = s;
702 
703  if (+(mrc.cw*DirectPart))
704  { REPORT mrc.data = store+o; mrc.storage = w1; }
705  else
706  {
707  // do not allow StoreOnExit and !DirectPart
708  if (+(mrc.cw*StoreOnExit))
709  Throw(InternalException("SymmetricBandMatrix::GetCol(MatrixRowCol&)"));
710  int w = w1+lower_val; s += w-ncols_val; Real* ColCopy;
711  if (s>0) w -= s; mrc.storage = w; int w2 = w-w1;
712 
713  if ( +(mrc.cw*HaveStore) ) { REPORT ColCopy = mrc.data; }
714  else
715  {
716  REPORT ColCopy = new Real [2*lower_val+1]; MatrixErrorNoSpace(ColCopy);
717  MONITOR_REAL_NEW("Make (SmBGetCol)",2*lower_val+1,ColCopy)
718  mrc.cw += HaveStore; mrc.data = ColCopy;
719  }
720 
721  if (+(mrc.cw*LoadOnEntry))
722  {
723  REPORT
724  Real* Mstore = store+o;
725  while (w1--) *ColCopy++ = *Mstore++;
726  Mstore--;
727  while (w2--) { Mstore += lower_val; *ColCopy++ = *Mstore; }
728  }
729  }
730 }
731 
733 {
734  int c=mrc.rowcol; int w1 = lower_val+1; mrc.length=nrows_val;
735  if (+(mrc.cw*DirectPart))
736  {
737  REPORT
738  int b = c*w1+lower_val;
739  mrc.skip = c; c += w1-nrows_val; w1 -= c; mrc.storage = w1;
740  Real* ColCopy = mrc.data = mrc.store+mrc.skip;
741  if (+(mrc.cw*LoadOnEntry))
742  {
743  REPORT
744  Real* Mstore = store+b;
745  // while (w1--) { *ColCopy++ = *Mstore; Mstore += lower; }
746  if (w1) for (;;)
747  { *ColCopy++ = *Mstore; if (!(--w1)) break; Mstore += lower_val; }
748  }
749  }
750  else
751  {
752  REPORT
753  // do not allow StoreOnExit and !DirectPart
754  if (+(mrc.cw*StoreOnExit))
755  Throw(InternalException("SymmetricBandMatrix::GetCol(MatrixColX&)"));
756  int s = c-lower_val; int o = c*w1;
757  if (s<0) { w1 += s; o -= s; s = 0; }
758  mrc.skip = s;
759 
760  int w = w1+lower_val; s += w-ncols_val;
761  if (s>0) w -= s; mrc.storage = w; int w2 = w-w1;
762 
763  Real* ColCopy = mrc.data = mrc.store+mrc.skip;
764 
765  if (+(mrc.cw*LoadOnEntry))
766  {
767  REPORT
768  Real* Mstore = store+o;
769  while (w1--) *ColCopy++ = *Mstore++;
770  Mstore--;
771  while (w2--) { Mstore += lower_val; *ColCopy++ = *Mstore; }
772  }
773 
774  }
775 }
776 
778 {
779  REPORT
780  int c = mrc.rowcol;
781  Real* Mstore = store + c*lower_val+c+lower_val;
782  Real* Cstore = mrc.data; int w = mrc.storage;
783  // while (w--) { *Mstore = *Cstore++; Mstore += lower_val; }
784  if (w) for (;;)
785  { *Mstore = *Cstore++; if (!(--w)) break; Mstore += lower_val; }
786 }
787 
788 // routines for identity matrix
789 
791 {
792  REPORT
793  mrc.skip=mrc.rowcol; mrc.storage=1; mrc.data=store; mrc.length=ncols_val;
794 }
795 
797 {
798  REPORT
799  mrc.skip=mrc.rowcol; mrc.storage=1; mrc.length=nrows_val;
800  if (+(mrc.cw*StoreHere)) // should not happen
801  Throw(InternalException("IdentityMatrix::GetCol(MatrixRowCol&)"));
802  else { REPORT mrc.data=store; }
803 }
804 
806 {
807  REPORT
808  mrc.skip=mrc.rowcol; mrc.storage=1; mrc.length=nrows_val;
809  mrc.data = mrc.store+mrc.skip; *(mrc.data)=*store;
810 }
811 
813 
815 
817 {
818  REPORT
819  if (+(mrc.cw*StoreOnExit)) { REPORT *store=*(mrc.data); }
820  mrc.IncrDiag(); // must increase mrc.data so need IncrDiag
821  int t1 = +(mrc.cw*LoadOnEntry);
822  if (t1 && mrc.rowcol < ncols_val) { REPORT *(mrc.data)=*store; }
823 }
824 
825 
826 
827 
828 // *************************** destructors *******************************
829 
831 {
832  if (+(cw*HaveStore))
833  {
834  MONITOR_REAL_DELETE("Free (RowCol)",-1,data) // do not know length
835  delete [] data;
836  }
837 }
838 
839 MatrixRow::~MatrixRow() { if (+(cw*StoreOnExit)) gm->RestoreRow(*this); }
840 
841 MatrixCol::~MatrixCol() { if (+(cw*StoreOnExit)) gm->RestoreCol(*this); }
842 
843 MatrixColX::~MatrixColX() { if (+(cw*StoreOnExit)) gm->RestoreCol(*this); }
844 
845 #ifdef use_namespace
846 }
847 #endif
848 
void RestoreCol(MatrixRowCol &)
Definition: newmat.h:1285
void IncrMat()
Definition: newmatrc.h:53
void IncrId()
Definition: newmatrc.h:55
void NextRow(MatrixRowCol &)
Definition: newmat3.cpp:334
void GetCol(MatrixRowCol &)
Definition: newmat3.cpp:530
Internal newmat exception - shouldn&#39;t happen.
Definition: newmat.h:2025
void NextRow(MatrixRowCol &)
Definition: newmat3.cpp:189
~MatrixCol()
Definition: newmat3.cpp:841
void NextRow(MatrixRowCol &)
Definition: newmat3.cpp:258
double Real
Definition: include.h:307
void GetCol(MatrixRowCol &)
Definition: newmat3.cpp:285
void RestoreCol(MatrixRowCol &)
Definition: newmat3.cpp:640
void GetCol(MatrixRowCol &)
Definition: newmat3.cpp:123
void GetRow(MatrixRowCol &)
Definition: newmat3.cpp:572
void GetCol(MatrixRowCol &)
Definition: newmat3.cpp:691
void RestoreCol(MatrixRowCol &)
Definition: newmat3.cpp:165
void NextRow(MatrixRowCol &)
Definition: newmat3.cpp:812
void GetCol(MatrixRowCol &)
Definition: newmat3.cpp:346
void IncrDiag()
Definition: newmatrc.h:54
void GetRow(MatrixRowCol &)
Definition: newmat3.cpp:115
void GetCol(MatrixRowCol &)
Definition: newmat3.cpp:240
#define REPORT
Definition: newmat3.cpp:24
virtual void NextRow(MatrixRowCol &)
Definition: newmat3.cpp:85
virtual void NextCol(MatrixRowCol &)
Definition: newmat3.cpp:94
void IncrUT()
Definition: newmatrc.h:56
void GetCol(MatrixRowCol &)
Definition: newmat3.cpp:593
Real * data
Definition: newmatrc.h:51
int storage
Definition: newmatrc.h:48
#define MONITOR_REAL_DELETE(Operation, Size, Pointer)
Definition: myexcept.h:331
Real * store
Definition: newmatrc.h:151
void NextCol(MatrixRowCol &)
Definition: newmat3.cpp:551
void GetRow(MatrixRowCol &)
Definition: newmat3.cpp:277
void NextCol(MatrixRowCol &)
Definition: newmat3.cpp:191
void RestoreCol(MatrixRowCol &)
Definition: newmat.h:782
LoadAndStoreFlag cw
Definition: newmatrc.h:52
#define Throw(E)
Definition: myexcept.h:191
#define MONITOR_REAL_NEW(Operation, Size, Pointer)
Definition: myexcept.h:329
void NextCol(MatrixRowCol &)
Definition: newmat3.cpp:261
void RestoreCol(MatrixRowCol &)
Definition: newmat3.cpp:324
void RestoreCol(MatrixRowCol &)
Definition: newmat3.cpp:387
void GetCol(MatrixRowCol &)
Definition: newmat3.cpp:796
void MatrixErrorNoSpace(const void *)
test for allocation fails
Definition: newmatex.cpp:301
void NextCol(MatrixRowCol &)
Definition: newmat3.cpp:814
void RestoreCol(MatrixRowCol &)
Definition: newmat.h:981
~MatrixRow()
Definition: newmat3.cpp:839
void GetCol(MatrixRowCol &)
Definition: newmat3.cpp:435
void GetRow(MatrixRowCol &)
Definition: newmat3.cpp:654
void GetRow(MatrixRowCol &)
Definition: newmat3.cpp:402
int rowcol
Definition: newmatrc.h:49
void IncrLT()
Definition: newmatrc.h:57
void NextRow(MatrixRowCol &)
Definition: newmat3.cpp:584
void GetRow(MatrixRowCol &)
Definition: newmat3.cpp:233
void GetRow(MatrixRowCol &)
Definition: newmat3.cpp:339
void GetRow(MatrixRowCol &)
Definition: newmat3.cpp:790
void NextRow(MatrixRowCol &)
Definition: newmat3.cpp:397
int length
Definition: newmatrc.h:46


kni
Author(s): Martin Günther
autogenerated on Fri Jan 3 2020 04:01:16