newmat4.cpp
Go to the documentation of this file.
1 
6 
7 
8 // Copyright (C) 1991,2,3,4,8,9: R B Davies
9 
10 //#define WANT_STREAM
11 
12 #include "include.h"
13 
14 #include "newmat.h"
15 #include "newmatrc.h"
16 
17 #ifdef use_namespace
18 namespace NEWMAT {
19 #endif
20 
21 
22 
23 #ifdef DO_REPORT
24 #define REPORT { static ExeCounter ExeCount(__LINE__,4); ++ExeCount; }
25 #else
26 #define REPORT {}
27 #endif
28 
29 
30 #define DO_SEARCH // search for LHS of = in RHS
31 
32 // ************************* general utilities *************************/
33 
34 static int tristore(int n) // elements in triangular matrix
35 { return (n*(n+1))/2; }
36 
37 
38 // **************************** constructors ***************************/
39 
41 { store=0; storage=0; nrows_val=0; ncols_val=0; tag_val=-1; }
42 
44 {
45  REPORT
46  storage=s.Value(); tag_val=-1;
47  if (storage)
48  {
49  store = new Real [storage]; MatrixErrorNoSpace(store);
50  MONITOR_REAL_NEW("Make (GenMatrix)",storage,store)
51  }
52  else store = 0;
53 }
54 
55 Matrix::Matrix(int m, int n) : GeneralMatrix(m*n)
57 
59  : Matrix(n.Value(),n.Value())
60 { REPORT }
61 
63  : GeneralMatrix(tristore(n.Value()))
64 { REPORT nrows_val=n.Value(); ncols_val=n.Value(); }
65 
67  : GeneralMatrix(tristore(n.Value()))
68 { REPORT nrows_val=n.Value(); ncols_val=n.Value(); }
69 
71  : GeneralMatrix(tristore(n.Value()))
72 { REPORT nrows_val=n.Value(); ncols_val=n.Value(); }
73 
75 { REPORT nrows_val=m.Value(); ncols_val=m.Value(); }
76 
78 {
79  REPORT // CheckConversion(M);
80  // MatrixConversionCheck mcc;
82  GetMatrix(gmx);
83 }
84 
86 {
87  REPORT
88  if (ncols_val != nrows_val)
89  {
90  Tracer tr("SquareMatrix");
91  Throw(NotSquareException(*this));
92  }
93 }
94 
95 
97 {
98  REPORT
99  if (gm.ncols_val != gm.nrows_val)
100  {
101  Tracer tr("SquareMatrix(Matrix)");
103  }
104  GetMatrix(&gm);
105 }
106 
107 
109 {
110  REPORT
111  if (nrows_val!=1)
112  {
113  Tracer tr("RowVector");
114  Throw(VectorException(*this));
115  }
116 }
117 
119 {
120  REPORT
121  if (ncols_val!=1)
122  {
123  Tracer tr("ColumnVector");
124  Throw(VectorException(*this));
125  }
126 }
127 
129 {
130  REPORT // CheckConversion(M);
131  // MatrixConversionCheck mcc;
133  GetMatrix(gmx);
134 }
135 
137 {
138  REPORT // CheckConversion(M);
139  // MatrixConversionCheck mcc;
141  GetMatrix(gmx);
142 }
143 
145 {
146  REPORT // CheckConversion(M);
147  // MatrixConversionCheck mcc;
149  GetMatrix(gmx);
150 }
151 
153 {
154  REPORT //CheckConversion(M);
155  // MatrixConversionCheck mcc;
157  GetMatrix(gmx);
158 }
159 
161 {
162  REPORT //CheckConversion(M);
163  // MatrixConversionCheck mcc;
165  GetMatrix(gmx);
166 }
167 
169 {
170  if (store)
171  {
172  MONITOR_REAL_DELETE("Free (GenMatrix)",storage,store)
173  delete [] store;
174  }
175 }
176 
178 {
179  REPORT
180  Tracer tr("CroutMatrix");
181  indx = 0; // in case of exception at next line
182  GeneralMatrix* gm = ((BaseMatrix&)m).Evaluate();
183  if (gm->nrows_val!=gm->ncols_val)
184  { gm->tDelete(); Throw(NotSquareException(*gm)); }
185  if (gm->type() == MatrixType::Ct)
186  { REPORT ((CroutMatrix*)gm)->get_aux(*this); GetMatrix(gm); }
187  else
188  {
189  REPORT
191  GetMatrix(gm1);
192  d=true; sing=false;
193  indx=new int [nrows_val]; MatrixErrorNoSpace(indx);
194  MONITOR_INT_NEW("Index (CroutMat)",nrows_val,indx)
195  ludcmp();
196  }
197 }
198 
199 // could we use SetParameters instead of this
201 {
202  X.d = d; X.sing = sing;
203  if (tag_val == 0 || tag_val == 1) // reuse the array
204  { REPORT X.indx = indx; indx = 0; d = true; sing = true; return; }
205  else if (nrows_val == 0)
206  { REPORT indx = 0; d = true; sing = true; return; }
207  else // copy the array
208  {
209  REPORT
210  Tracer tr("CroutMatrix::get_aux");
211  int *ix = new int [nrows_val]; MatrixErrorNoSpace(ix);
212  MONITOR_INT_NEW("Index (CM::get_aux)", nrows_val, ix)
213  int n = nrows_val; int* i = ix; int* j = indx;
214  while(n--) *i++ = *j++;
215  X.indx = ix;
216  }
217 }
218 
220 {
221  REPORT
222  Tracer tr("CroutMatrix(const CroutMatrix&)");
223  ((CroutMatrix&)gm).get_aux(*this);
224  GetMatrix(&gm);
225 }
226 
228 {
229  MONITOR_INT_DELETE("Index (CroutMat)",nrows_val,indx)
230  delete [] indx;
231 }
232 
233 //ReturnMatrix::ReturnMatrix(GeneralMatrix& gmx)
234 //{
235 // REPORT
236 // gm = gmx.Image(); gm->ReleaseAndDelete();
237 //}
238 
239 
240 GeneralMatrix::operator ReturnMatrix() const
241 {
242  REPORT
243  GeneralMatrix* gm = Image(); gm->ReleaseAndDelete();
244  return ReturnMatrix(gm);
245 }
246 
247 
248 
250 {
251  REPORT
252  GeneralMatrix* gm = Image(); gm->ReleaseAndDelete();
253  return ReturnMatrix(gm);
254 }
255 
256 // ************ Constructors for use with NR in C++ interface ***********
257 
258 #ifdef SETUP_C_SUBSCRIPTS
259 
260 Matrix::Matrix(Real a, int m, int n) : GeneralMatrix(m * n)
261  { REPORT nrows_val=m; ncols_val=n; operator=(a); }
262 
263 Matrix::Matrix(const Real* a, int m, int n) : GeneralMatrix(m * n)
264  { REPORT nrows_val=m; ncols_val=n; *this << a; }
265 
266 #endif
267 
268 
269 
270 // ************************** resize matrices ***************************/
271 
272 void GeneralMatrix::resize(int nr, int nc, int s)
273 {
274  REPORT
275  if (store)
276  {
277  MONITOR_REAL_DELETE("Free (ReDimensi)",storage,store)
278  delete [] store;
279  }
280  storage=s; nrows_val=nr; ncols_val=nc; tag_val=-1;
281  if (s)
282  {
284  MONITOR_REAL_NEW("Make (ReDimensi)",storage,store)
285  }
286  else store = 0;
287 }
288 
289 void Matrix::resize(int nr, int nc)
290 { REPORT GeneralMatrix::resize(nr,nc,nr*nc); }
291 
293 { REPORT GeneralMatrix::resize(n,n,n*n); }
294 
295 void SquareMatrix::resize(int nr, int nc)
296 {
297  REPORT
298  Tracer tr("SquareMatrix::resize");
299  if (nc != nr) Throw(NotSquareException(*this));
300  GeneralMatrix::resize(nr,nc,nr*nc);
301 }
302 
304 { REPORT GeneralMatrix::resize(nr,nr,tristore(nr)); }
305 
307 { REPORT GeneralMatrix::resize(nr,nr,tristore(nr)); }
308 
310 { REPORT GeneralMatrix::resize(nr,nr,tristore(nr)); }
311 
313 { REPORT GeneralMatrix::resize(nr,nr,nr); }
314 
315 void RowVector::resize(int nc)
316 { REPORT GeneralMatrix::resize(1,nc,nc); }
317 
319 { REPORT GeneralMatrix::resize(nr,1,nr); }
320 
321 void RowVector::resize(int nr, int nc)
322 {
323  Tracer tr("RowVector::resize");
324  if (nr != 1) Throw(VectorException(*this));
325  REPORT GeneralMatrix::resize(1,nc,nc);
326 }
327 
328 void ColumnVector::resize(int nr, int nc)
329 {
330  Tracer tr("ColumnVector::resize");
331  if (nc != 1) Throw(VectorException(*this));
332  REPORT GeneralMatrix::resize(nr,1,nr);
333 }
334 
336 { REPORT GeneralMatrix::resize(nr,nr,1); *store = 1; }
337 
338 
340 { REPORT resize(A.Nrows(), A.Ncols()); }
341 
343 {
344  REPORT
345  int n = A.Nrows();
346  if (n != A.Ncols())
347  {
348  Tracer tr("SquareMatrix::resize(GM)");
349  Throw(NotSquareException(*this));
350  }
351  resize(n);
352 }
353 
355 { REPORT resize(A.Nrows(), A.Ncols()); }
356 
358 { REPORT resize(A.Nrows(), A.Ncols()); }
359 
361 { REPORT resize(A.Nrows(), A.Ncols()); }
362 
364 {
365  REPORT
366  int n = A.Nrows();
367  if (n != A.Ncols())
368  {
369  Tracer tr("SymmetricMatrix::resize(GM)");
370  Throw(NotSquareException(*this));
371  }
372  resize(n);
373 }
374 
376 {
377  REPORT
378  int n = A.Nrows();
379  if (n != A.Ncols())
380  {
381  Tracer tr("DiagonalMatrix::resize(GM)");
382  Throw(NotSquareException(*this));
383  }
384  resize(n);
385 }
386 
388 {
389  REPORT
390  int n = A.Nrows();
391  if (n != A.Ncols())
392  {
393  Tracer tr("UpperTriangularMatrix::resize(GM)");
394  Throw(NotSquareException(*this));
395  }
396  resize(n);
397 }
398 
400 {
401  REPORT
402  int n = A.Nrows();
403  if (n != A.Ncols())
404  {
405  Tracer tr("LowerTriangularMatrix::resize(GM)");
406  Throw(NotSquareException(*this));
407  }
408  resize(n);
409 }
410 
412 {
413  REPORT
414  int n = A.Nrows();
415  if (n != A.Ncols())
416  {
417  Tracer tr("IdentityMatrix::resize(GM)");
418  Throw(NotSquareException(*this));
419  }
420  resize(n);
421 }
422 
424 {
425  REPORT
426  Tracer tr("GeneralMatrix::resize(GM)");
427  Throw(NotDefinedException("resize", "this type of matrix"));
428 }
429 
430 //*********************** resize_keep *******************************
431 
432 void Matrix::resize_keep(int nr, int nc)
433 {
434  Tracer tr("Matrix::resize_keep");
435  if (nr == nrows_val && nc == ncols_val) { REPORT return; }
436 
437  if (nr <= nrows_val && nc <= ncols_val)
438  {
439  REPORT
440  Matrix X = submatrix(1,nr,1,nc);
441  swap(X);
442  }
443  else if (nr >= nrows_val && nc >= ncols_val)
444  {
445  REPORT
446  Matrix X(nr, nc); X = 0;
447  X.submatrix(1,nrows_val,1,ncols_val) = *this;
448  swap(X);
449  }
450  else
451  {
452  REPORT
453  Matrix X(nr, nc); X = 0;
454  if (nr > nrows_val) nr = nrows_val;
455  if (nc > ncols_val) nc = ncols_val;
456  X.submatrix(1,nr,1,nc) = submatrix(1,nr,1,nc);
457  swap(X);
458  }
459 }
460 
462 {
463  Tracer tr("SquareMatrix::resize_keep");
464  if (nr < nrows_val)
465  {
466  REPORT
467  SquareMatrix X = sym_submatrix(1,nr);
468  swap(X);
469  }
470  else if (nr > nrows_val)
471  {
472  REPORT
473  SquareMatrix X(nr); X = 0;
474  X.sym_submatrix(1,nrows_val) = *this;
475  swap(X);
476  }
477 }
478 
479 void SquareMatrix::resize_keep(int nr, int nc)
480 {
481  Tracer tr("SquareMatrix::resize_keep 2");
482  REPORT
483  if (nr != nc) Throw(NotSquareException(*this));
484  resize_keep(nr);
485 }
486 
487 
489 {
490  Tracer tr("SymmetricMatrix::resize_keep");
491  if (nr < nrows_val)
492  {
493  REPORT
494  SymmetricMatrix X = sym_submatrix(1,nr);
495  swap(X);
496  }
497  else if (nr > nrows_val)
498  {
499  REPORT
500  SymmetricMatrix X(nr); X = 0;
501  X.sym_submatrix(1,nrows_val) = *this;
502  swap(X);
503  }
504 }
505 
507 {
508  Tracer tr("UpperTriangularMatrix::resize_keep");
509  if (nr < nrows_val)
510  {
511  REPORT
513  swap(X);
514  }
515  else if (nr > nrows_val)
516  {
517  REPORT
518  UpperTriangularMatrix X(nr); X = 0;
519  X.sym_submatrix(1,nrows_val) = *this;
520  swap(X);
521  }
522 }
523 
525 {
526  Tracer tr("LowerTriangularMatrix::resize_keep");
527  if (nr < nrows_val)
528  {
529  REPORT
531  swap(X);
532  }
533  else if (nr > nrows_val)
534  {
535  REPORT
536  LowerTriangularMatrix X(nr); X = 0;
537  X.sym_submatrix(1,nrows_val) = *this;
538  swap(X);
539  }
540 }
541 
543 {
544  Tracer tr("DiagonalMatrix::resize_keep");
545  if (nr < nrows_val)
546  {
547  REPORT
548  DiagonalMatrix X = sym_submatrix(1,nr);
549  swap(X);
550  }
551  else if (nr > nrows_val)
552  {
553  REPORT
554  DiagonalMatrix X(nr); X = 0;
555  X.sym_submatrix(1,nrows_val) = *this;
556  swap(X);
557  }
558 }
559 
561 {
562  Tracer tr("RowVector::resize_keep");
563  if (nc < ncols_val)
564  {
565  REPORT
566  RowVector X = columns(1,nc);
567  swap(X);
568  }
569  else if (nc > ncols_val)
570  {
571  REPORT
572  RowVector X(nc); X = 0;
573  X.columns(1,ncols_val) = *this;
574  swap(X);
575  }
576 }
577 
578 void RowVector::resize_keep(int nr, int nc)
579 {
580  Tracer tr("RowVector::resize_keep 2");
581  REPORT
582  if (nr != 1) Throw(VectorException(*this));
583  resize_keep(nc);
584 }
585 
587 {
588  Tracer tr("ColumnVector::resize_keep");
589  if (nr < nrows_val)
590  {
591  REPORT
592  ColumnVector X = rows(1,nr);
593  swap(X);
594  }
595  else if (nr > nrows_val)
596  {
597  REPORT
598  ColumnVector X(nr); X = 0;
599  X.rows(1,nrows_val) = *this;
600  swap(X);
601  }
602 }
603 
604 void ColumnVector::resize_keep(int nr, int nc)
605 {
606  Tracer tr("ColumnVector::resize_keep 2");
607  REPORT
608  if (nc != 1) Throw(VectorException(*this));
609  resize_keep(nr);
610 }
611 
612 
613 /*
614 void GeneralMatrix::resizeForAdd(const GeneralMatrix& A, const GeneralMatrix&)
615 { REPORT resize(A); }
616 
617 void GeneralMatrix::resizeForSP(const GeneralMatrix& A, const GeneralMatrix&)
618 { REPORT resize(A); }
619 
620 
621 // ************************* SameStorageType ******************************
622 
623 // SameStorageType checks A and B have same storage type including bandwidth
624 // It does not check same dimensions since we assume this is already done
625 
626 bool GeneralMatrix::SameStorageType(const GeneralMatrix& A) const
627 {
628  REPORT
629  return type() == A.type();
630 }
631 */
632 
633 // ******************* manipulate types, storage **************************/
634 
636 { REPORT return (s==this) ? 1 : 0; }
637 
639 { REPORT return gm->search(s); }
640 
642 { REPORT return bm1->search(s) + bm2->search(s); }
643 
645 { REPORT return bm->search(s); }
646 
648 { REPORT return bm->search(s); }
649 
650 int ReturnMatrix::search(const BaseMatrix* s) const
651 { REPORT return (s==gm) ? 1 : 0; }
652 
666 
668 
669 
670 
674 
676  { REPORT return MatrixBandWidth(0,-1); }
677 
679  { REPORT return MatrixBandWidth(-1,0); }
680 
682  { REPORT return MatrixBandWidth(lower_val,upper_val); }
683 
685  { REPORT return MatrixBandWidth(m1,m2); }
686 
688  { REPORT return gm->bandwidth(); }
689 
691  { REPORT return gm1->bandwidth() + gm2->bandwidth(); }
692 
694  { REPORT return gm1->bandwidth().minimum(gm2->bandwidth()); }
695 
697 {
698  int lower, upper;
699  MatrixBandWidth bw1 = gm1->bandwidth(), bw2 = gm2->bandwidth();
700  if (bw1.Lower() < 0)
701  {
702  if (bw2.Lower() < 0) { REPORT lower = -1; }
703  else { REPORT lower = bw2.Lower() + (gm1->Nrows() - 1) * gm2->Nrows(); }
704  }
705  else
706  {
707  if (bw2.Lower() < 0)
708  { REPORT lower = (1 + bw1.Lower()) * gm2->Nrows() - 1; }
709  else { REPORT lower = bw2.Lower() + bw1.Lower() * gm2->Nrows(); }
710  }
711  if (bw1.Upper() < 0)
712  {
713  if (bw2.Upper() < 0) { REPORT upper = -1; }
714  else { REPORT upper = bw2.Upper() + (gm1->Nrows() - 1) * gm2->Nrows(); }
715  }
716  else
717  {
718  if (bw2.Upper() < 0)
719  { REPORT upper = (1 + bw1.Upper()) * gm2->Nrows() - 1; }
720  else { REPORT upper = bw2.Upper() + bw1.Upper() * gm2->Nrows(); }
721  }
722  return MatrixBandWidth(lower, upper);
723 }
724 
726 { REPORT return gm1->bandwidth() * gm2->bandwidth(); }
727 
729 
731 {
732  if (+gm1->type() & MatrixType::Diagonal)
733  { REPORT return gm2->bandwidth(); }
734  else { REPORT return -1; }
735 }
736 
738  { REPORT return gm->bandwidth(); }
739 
741  { REPORT return gm->bandwidth(); }
742 
744  { REPORT return gm->bandwidth().t(); }
745 
747 {
748  if (+gm->type() & MatrixType::Diagonal)
749  { REPORT return MatrixBandWidth(0,0); }
750  else { REPORT return -1; }
751 }
752 
758  { REPORT return gm->bandwidth(); }
759 
761 {
762 
763  if (row_skip==col_skip && row_number==col_number)
764  { REPORT return gm->bandwidth(); }
765  else { REPORT return MatrixBandWidth(-1); }
766 }
767 
768 // ********************** the memory managment tools **********************/
769 
770 // Rules regarding tDelete, reuse, GetStore, BorrowStore
771 // All matrices processed during expression evaluation must be subject
772 // to exactly one of reuse(), tDelete(), GetStore() or BorrowStore().
773 // If reuse returns true the matrix must be reused.
774 // GetMatrix(gm) always calls gm->GetStore()
775 // gm->Evaluate obeys rules
776 // bm->Evaluate obeys rules for matrices in bm structure
777 
778 // Meaning of tag_val
779 // tag_val = -1 memory cannot be reused (default situation)
780 // tag_val = -2 memory has been borrowed from another matrix
781 // (don't change values)
782 // tag_val = i > 0 delete or reuse memory after i operations
783 // tag_val = 0 like value 1 but matrix was created by new
784 // so delete it
785 
787 {
788  if (tag_val<0)
789  {
790  if (tag_val<-1) { REPORT store = 0; delete this; return; } // borrowed
791  else { REPORT return; } // not a temporary matrix - leave alone
792  }
793  if (tag_val==1)
794  {
795  if (store)
796  {
797  REPORT MONITOR_REAL_DELETE("Free (tDelete)",storage,store)
798  delete [] store;
799  }
800  MiniCleanUp(); return; // CleanUp
801  }
802  if (tag_val==0) { REPORT delete this; return; }
803 
804  REPORT tag_val--; return;
805 }
806 
807 void newmat_block_copy(int n, Real* from, Real* to)
808 {
809  REPORT
810  int i = (n >> 3);
811  while (i--)
812  {
813  *to++ = *from++; *to++ = *from++; *to++ = *from++; *to++ = *from++;
814  *to++ = *from++; *to++ = *from++; *to++ = *from++; *to++ = *from++;
815  }
816  i = n & 7; while (i--) *to++ = *from++;
817 }
818 
820 {
821  if (tag_val < -1) // borrowed storage
822  {
823  if (storage)
824  {
825  REPORT
826  Real* s = new Real [storage]; MatrixErrorNoSpace(s);
827  MONITOR_REAL_NEW("Make (reuse)",storage,s)
829  }
830  else { REPORT MiniCleanUp(); } // CleanUp
831  tag_val = 0; return true;
832  }
833  if (tag_val < 0 ) { REPORT return false; }
834  if (tag_val <= 1 ) { REPORT return true; }
835  REPORT tag_val--; return false;
836 }
837 
839 {
840  if (tag_val<0 || tag_val>1)
841  {
842  Real* s;
843  if (storage)
844  {
845  s = new Real [storage]; MatrixErrorNoSpace(s);
846  MONITOR_REAL_NEW("Make (GetStore)",storage,s)
848  }
849  else s = 0;
850  if (tag_val > 1) { REPORT tag_val--; }
851  else if (tag_val < -1) { REPORT store = 0; delete this; } // borrowed store
852  else { REPORT }
853  return s;
854  }
855  Real* s = store; // cleanup - done later
856  if (tag_val==0) { REPORT store = 0; delete this; }
857  else { REPORT MiniCleanUp(); }
858  return s;
859 }
860 
862 {
863  REPORT tag_val=-1; nrows_val=gmx->Nrows(); ncols_val=gmx->Ncols();
864  storage=gmx->storage; SetParameters(gmx);
865  store=((GeneralMatrix*)gmx)->GetStore();
866 }
867 
869 // Copy storage of *this to storage of *gmx. Then convert to type mt.
870 // If mt == 0 just let *gmx point to storage of *this if tag_val==-1.
871 {
872  if (!mt)
873  {
874  if (tag_val == -1) { REPORT gmx->tag_val = -2; gmx->store = store; }
875  else { REPORT gmx->tag_val = 0; gmx->store = GetStore(); }
876  }
877  else if (Compare(gmx->type(),mt))
878  { REPORT gmx->tag_val = 0; gmx->store = GetStore(); }
879  else
880  {
881  REPORT gmx->tag_val = -2; gmx->store = store;
882  gmx = gmx->Evaluate(mt); gmx->tag_val = 0; tDelete();
883  }
884 
885  return gmx;
886 }
887 
889 // Count number of references to this in X.
890 // If zero delete storage in this;
891 // otherwise tag this to show when storage can be deleted
892 // evaluate X and copy to this
893 {
894 #ifdef DO_SEARCH
895  int counter=X.search(this);
896  if (counter==0)
897  {
898  REPORT
899  if (store)
900  {
901  MONITOR_REAL_DELETE("Free (operator=)",storage,store)
902  REPORT delete [] store; storage = 0; store = 0;
903  }
904  }
905  else { REPORT Release(counter); }
906  GeneralMatrix* gmx = ((BaseMatrix&)X).Evaluate(mt);
907  if (gmx!=this) { REPORT GetMatrix(gmx); }
908  else { REPORT }
909  Protect();
910 #else
911  GeneralMatrix* gmx = ((BaseMatrix&)X).Evaluate(mt);
912  if (gmx!=this)
913  {
914  REPORT
915  if (store)
916  {
917  MONITOR_REAL_DELETE("Free (operator=)",storage,store)
918  REPORT delete [] store; storage = 0; store = 0;
919  }
920  GetMatrix(gmx);
921  }
922  else { REPORT }
923  Protect();
924 #endif
925 }
926 
927 // version with no conversion
929 {
930  GeneralMatrix* gmx = (GeneralMatrix*)&X;
931  if (gmx!=this)
932  {
933  REPORT
934  if (store)
935  {
936  MONITOR_REAL_DELETE("Free (operator=)",storage,store)
937  REPORT delete [] store; storage = 0; store = 0;
938  }
939  GetMatrix(gmx);
940  }
941  else { REPORT }
942  Protect();
943 }
944 
945 // version to work with operator<<
946 void GeneralMatrix::Eq(const BaseMatrix& X, MatrixType mt, bool ldok)
947 {
948  REPORT
949  if (ldok) mt.SetDataLossOK();
950  Eq(X, mt);
951 }
952 
954 // a cut down version of Eq for use with += etc.
955 // we know BaseMatrix points to two GeneralMatrix objects,
956 // the first being this (may be the same).
957 // we know tag_val has been set correctly in each.
958 {
959  GeneralMatrix* gmx = ((BaseMatrix&)X).Evaluate(mt);
960  if (gmx!=this) { REPORT GetMatrix(gmx); } // simplify GetMatrix ?
961  else { REPORT }
962  Protect();
963 }
964 
966 // copy stored values of X; otherwise leave els of *this unchanged
967 {
968  REPORT
969  Tracer tr("inject");
970  if (nrows_val != X.nrows_val || ncols_val != X.ncols_val)
974  int i=nrows_val;
975  while (i--) { mrx.Inject(mr); mrx.Next(); mr.Next(); }
976 }
977 
978 // ************* checking for data loss during conversion *******************/
979 
980 bool Compare(const MatrixType& source, MatrixType& destination)
981 {
982  if (!destination) { destination=source; return true; }
983  if (destination==source) return true;
984  if (!destination.DataLossOK && !(destination>=source))
985  Throw(ProgramException("Illegal Conversion", source, destination));
986  return false;
987 }
988 
989 // ************* Make a copy of a matrix on the heap *********************/
990 
992 {
993  REPORT
994  GeneralMatrix* gm = new Matrix(*this); MatrixErrorNoSpace(gm);
995  return gm;
996 }
997 
999 {
1000  REPORT
1001  GeneralMatrix* gm = new SquareMatrix(*this); MatrixErrorNoSpace(gm);
1002  return gm;
1003 }
1004 
1006 {
1007  REPORT
1008  GeneralMatrix* gm = new SymmetricMatrix(*this); MatrixErrorNoSpace(gm);
1009  return gm;
1010 }
1011 
1013 {
1014  REPORT
1015  GeneralMatrix* gm = new UpperTriangularMatrix(*this);
1016  MatrixErrorNoSpace(gm); return gm;
1017 }
1018 
1020 {
1021  REPORT
1022  GeneralMatrix* gm = new LowerTriangularMatrix(*this);
1023  MatrixErrorNoSpace(gm); return gm;
1024 }
1025 
1027 {
1028  REPORT
1029  GeneralMatrix* gm = new DiagonalMatrix(*this); MatrixErrorNoSpace(gm);
1030  return gm;
1031 }
1032 
1034 {
1035  REPORT
1036  GeneralMatrix* gm = new RowVector(*this); MatrixErrorNoSpace(gm);
1037  return gm;
1038 }
1039 
1041 {
1042  REPORT
1043  GeneralMatrix* gm = new ColumnVector(*this); MatrixErrorNoSpace(gm);
1044  return gm;
1045 }
1046 
1048 {
1049  REPORT
1050  GeneralMatrix* gm = new nricMatrix(*this); MatrixErrorNoSpace(gm);
1051  return gm;
1052 }
1053 
1055 {
1056  REPORT
1057  GeneralMatrix* gm = new IdentityMatrix(*this); MatrixErrorNoSpace(gm);
1058  return gm;
1059 }
1060 
1062 {
1063  REPORT
1064  GeneralMatrix* gm = new CroutMatrix(*this); MatrixErrorNoSpace(gm);
1065  return gm;
1066 }
1067 
1069 {
1070  bool dummy = true;
1071  if (dummy) // get rid of warning message
1072  Throw(InternalException("Cannot apply Image to this matrix type"));
1073  return 0;
1074 }
1075 
1076 
1077 // *********************** nricMatrix routines *****************************/
1078 
1080 {
1081  REPORT
1082  if (nrows_val > 0)
1083  {
1084  row_pointer = new Real* [nrows_val]; MatrixErrorNoSpace(row_pointer);
1085  Real* s = Store() - 1; int i = nrows_val; Real** rp = row_pointer;
1086  if (i) for (;;) { *rp++ = s; if (!(--i)) break; s+=ncols_val; }
1087  }
1088  else row_pointer = 0;
1089 }
1090 
1092  { REPORT if (nrows_val) delete [] row_pointer; }
1093 
1095 {
1096  REPORT
1097  if (!store)
1098  Throw(ProgramException("NRIC accessing matrix with unset dimensions"));
1099 }
1100 
1101 
1102 // *************************** CleanUp routines *****************************/
1103 
1105 {
1106  // set matrix dimensions to zero, delete storage
1107  REPORT
1108  if (store && storage)
1109  {
1110  MONITOR_REAL_DELETE("Free (cleanup) ",storage,store)
1111  REPORT delete [] store;
1112  }
1113  store=0; storage=0; nrows_val=0; ncols_val=0; tag_val = -1;
1114 }
1115 
1117  { REPORT DeleteRowPointer(); GeneralMatrix::cleanup(); }
1118 
1120  { REPORT DeleteRowPointer(); GeneralMatrix::MiniCleanUp(); }
1121 
1124 
1127 
1129 {
1130  REPORT
1131  if (nrows_val) delete [] indx;
1133 }
1134 
1136 {
1137  REPORT
1138  if (nrows_val) delete [] indx;
1140 }
1141 
1143 {
1144  REPORT
1145  if (nrows_val) delete [] indx;
1146  if (storage2) delete [] store2;
1148 }
1149 
1151 {
1152  REPORT
1153  if (nrows_val) delete [] indx;
1154  if (storage2) delete [] store2;
1156 }
1157 
1158 // ************************ simple integer array class ***********************
1159 
1160 // construct a new array of length xn. Check that xn is non-negative and
1161 // that space is available
1162 
1164 {
1165  if (n < 0) Throw(Logic_error("invalid array length"));
1166  else if (n == 0) { REPORT a = 0; }
1167  else { REPORT a = new int [n]; if (!a) Throw(Bad_alloc()); }
1168 }
1169 
1170 // destroy an array - return its space to memory
1171 
1173 
1174 // access an element of an array; return a "reference" so elements
1175 // can be modified.
1176 // check index is within range
1177 // in this array class the index runs from 0 to n-1
1178 
1180 {
1181  REPORT
1182  if (i < 0 || i >= n) Throw(Logic_error("array index out of range"));
1183  return a[i];
1184 }
1185 
1186 // same thing again but for arrays declared constant so we can't
1187 // modify its elements
1188 
1190 {
1191  REPORT
1192  if (i < 0 || i >= n) Throw(Logic_error("array index out of range"));
1193  return a[i];
1194 }
1195 
1196 // set all the elements equal to a given value
1197 
1199  { REPORT for (int i = 0; i < n; i++) a[i] = ai; }
1200 
1201 // set the elements equal to those of another array.
1202 // now allow length to be changed
1203 
1205 {
1206  REPORT
1207  if (b.n != n) resize(b.n);
1208  for (int i = 0; i < n; i++) a[i] = b.a[i];
1209 }
1210 
1211 // construct a new array equal to an existing array
1212 // check that space is available
1213 
1215 {
1216  if (n == 0) { REPORT a = 0; }
1217  else
1218  {
1219  REPORT
1220  a = new int [n]; if (!a) Throw(Bad_alloc());
1221  for (int i = 0; i < n; i++) a[i] = b.a[i];
1222  }
1223 }
1224 
1225 // change the size of an array; optionally copy data from old array to
1226 // new array
1227 
1228 void SimpleIntArray::resize(int n1, bool keep)
1229 {
1230  if (n1 == n) { REPORT return; }
1231  else if (n1 == 0) { REPORT n = 0; delete [] a; a = 0; }
1232  else if (n == 0)
1233  {
1234  REPORT
1235  a = new int [n1]; if (!a) Throw(Bad_alloc());
1236  n = n1;
1237  if (keep) operator=(0);
1238  }
1239  else
1240  {
1241  int* a1 = a;
1242  if (keep)
1243  {
1244  REPORT
1245  int i;
1246  a = new int [n1]; if (!a) Throw(Bad_alloc());
1247  if (n > n1) n = n1;
1248  else for (i = n; i < n1; i++) a[i] = 0;
1249  for (i = 0; i < n; i++) a[i] = a1[i];
1250  n = n1; delete [] a1;
1251  }
1252  else
1253  {
1254  REPORT n = n1; delete [] a1;
1255  a = new int [n]; if (!a) Throw(Bad_alloc());
1256  }
1257  }
1258 }
1259 
1260 //************************** swap values ********************************
1261 
1262 // swap values
1263 
1265 {
1266  REPORT
1267  int t;
1268  t = tag_val; tag_val = gm.tag_val; gm.tag_val = t;
1269  t = nrows_val; nrows_val = gm.nrows_val; gm.nrows_val = t;
1270  t = ncols_val; ncols_val = gm.ncols_val; gm.ncols_val = t;
1271  t = storage; storage = gm.storage; gm.storage = t;
1272  Real* s = store; store = gm.store; gm.store = s;
1273 }
1274 
1276 {
1277  REPORT
1279  Real** rp = row_pointer; row_pointer = gm.row_pointer; gm.row_pointer = rp;
1280 }
1281 
1283 {
1284  REPORT
1286  int* i = indx; indx = gm.indx; gm.indx = i;
1287  bool b;
1288  b = d; d = gm.d; gm.d = b;
1289  b = sing; sing = gm.sing; gm.sing = b;
1290 }
1291 
1293 {
1294  REPORT
1296  int i;
1297  i = lower_val; lower_val = gm.lower_val; gm.lower_val = i;
1298  i = upper_val; upper_val = gm.upper_val; gm.upper_val = i;
1299 }
1300 
1302 {
1303  REPORT
1305  int i;
1306  i = lower_val; lower_val = gm.lower_val; gm.lower_val = i;
1307 }
1308 
1310 {
1311  REPORT
1313  int* i = indx; indx = gm.indx; gm.indx = i;
1314  bool b;
1315  b = d; d = gm.d; gm.d = b;
1316  b = sing; sing = gm.sing; gm.sing = b;
1317  int m;
1318  m = storage2; storage2 = gm.storage2; gm.storage2 = m;
1319  m = m1; m1 = gm.m1; gm.m1 = m;
1320  m = m2; m2 = gm.m2; gm.m2 = m;
1321  Real* s = store2; store2 = gm.store2; gm.store2 = s;
1322 }
1323 
1325 {
1326  REPORT
1327  GeneralMatrix* tgm = gm; gm = x.gm; x.gm = tgm;
1328 }
1329 
1330 // ********************** C subscript classes ****************************
1331 
1333 {
1334  REPORT
1335  Tracer tr("RealStarStar");
1336  int n = A.ncols();
1337  int m = A.nrows();
1338  a = new Real*[m];
1339  MatrixErrorNoSpace(a);
1340  Real* d = A.data();
1341  for (int i = 0; i < m; ++i) a[i] = d + i * n;
1342 }
1343 
1345 {
1346  REPORT
1347  Tracer tr("ConstRealStarStar");
1348  int n = A.ncols();
1349  int m = A.nrows();
1350  a = new const Real*[m];
1351  MatrixErrorNoSpace(a);
1352  const Real* d = A.data();
1353  for (int i = 0; i < m; ++i) a[i] = d + i * n;
1354 }
1355 
1356 
1357 
1358 #ifdef use_namespace
1359 }
1360 #endif
1361 
1362 
1380 
1381 
1382 
1383 
1384 
1385 
1386 
MatrixBandWidth bandwidth() const
Definition: newmat4.cpp:740
friend class LowerTriangularMatrix
Definition: newmat.h:587
GeneralMatrix * Image() const
Definition: newmat4.cpp:1054
int search(const BaseMatrix *) const
Definition: newmat4.cpp:635
MatrixType type() const
Definition: newmat4.cpp:664
void resize(int)
Definition: newmat4.cpp:312
void MakeRowPointer()
Definition: newmat4.cpp:1079
Real * GetStore()
Definition: newmat4.cpp:838
GeneralMatrix * Image() const
Definition: newmat4.cpp:1040
void swap(SymmetricBandMatrix &gm)
Definition: newmat4.cpp:1301
MatrixBandWidth bandwidth() const
Definition: newmat4.cpp:746
MatrixType type() const
Definition: newmat4.cpp:662
int search(const BaseMatrix *) const
Definition: newmat4.cpp:644
void Release()
Definition: newmat.h:514
friend class SquareMatrix
Definition: newmat.h:583
void tDelete()
Definition: newmat4.cpp:786
void MiniCleanUp()
Definition: newmat4.cpp:1150
void operator=(const CroutMatrix &)
Definition: newmat6.cpp:427
void cleanup()
Definition: newmat4.cpp:1125
friend class UpperTriangularMatrix
Definition: newmat.h:586
MatrixType type() const
Definition: newmat4.cpp:667
void operator=(int ai)
set the array equal to a constant
Definition: newmat4.cpp:1198
Miscellaneous exception (details in character string).
Definition: newmat.h:1947
MatrixType type() const
Definition: newmat4.cpp:659
MatrixBandWidth bandwidth() const
Definition: newmat4.cpp:737
MatrixType type() const
Definition: newmat4.cpp:657
void inject(const GeneralMatrix &)
Definition: newmat4.cpp:965
#define MONITOR_INT_DELETE(Operation, Size, Pointer)
Definition: myexcept.h:332
Matrix()
Definition: newmat.h:629
ReturnMatrix for_return() const
Definition: newmat4.cpp:249
virtual int search(const BaseMatrix *) const =0
MatrixBandWidth bandwidth() const
Definition: newmat4.cpp:754
MatrixType type() const
Definition: newmat4.cpp:665
friend class nricMatrix
Definition: newmat.h:584
GetSubMatrix rows(int, int) const
Definition: submat.cpp:58
MatrixType type() const
Definition: newmat4.cpp:663
int * a
pointer to the array
Definition: newmat.h:1855
void resize_keep(int)
Definition: newmat4.cpp:524
Internal newmat exception - shouldn&#39;t happen.
Definition: newmat.h:2025
ColumnVector()
Definition: newmat.h:1012
void resize(int n)
Definition: newmat4.cpp:335
void DeleteRowPointer()
Definition: newmat4.cpp:1091
GeneralMatrix * gm
Definition: newmat.h:1399
friend class ColumnVector
Definition: newmat.h:591
GetSubMatrix columns(int, int) const
Definition: submat.cpp:77
virtual void MiniCleanUp()
Definition: newmat.h:482
MatrixType type() const
Definition: newmat4.cpp:654
void swap(GenericMatrix &x)
Definition: newmat4.cpp:1324
void resize_keep(int)
Definition: newmat4.cpp:560
double Real
Definition: include.h:307
GeneralMatrix * Image() const
Definition: newmat4.cpp:1047
void MiniCleanUp()
Definition: newmat4.cpp:1135
friend class DiagonalMatrix
Definition: newmat.h:588
int n
length of the array
Definition: newmat.h:1856
virtual MatrixType type() const =0
virtual MatrixBandWidth bandwidth() const
Definition: newmat4.cpp:671
GeneralMatrix * Evaluate(MatrixType mt=MatrixTypeUnSp)
Definition: newmat5.cpp:85
void Eq2(const BaseMatrix &, MatrixType)
Definition: newmat4.cpp:953
void swap(nricMatrix &gm)
Definition: newmat4.cpp:1275
void get_aux(CroutMatrix &)
Definition: newmat4.cpp:200
int Nrows() const
Definition: newmat.h:494
GeneralMatrix * Image() const
Definition: newmat4.cpp:998
~SimpleIntArray()
return the space to memory
Definition: newmat4.cpp:1172
MatrixType type() const
Definition: newmat4.cpp:661
bool sing
Definition: newmat.h:1063
CroutMatrix()
Definition: newmat.h:1069
MatrixBandWidth bandwidth() const
Definition: newmat4.cpp:684
MatrixBandWidth bandwidth() const
Definition: newmat4.cpp:690
#define MONITOR_INT_NEW(Operation, Size, Pointer)
Definition: myexcept.h:330
Real * data()
Definition: newmat.h:502
SquareMatrix()
Definition: newmat.h:683
MatrixBandWidth bandwidth() const
Definition: newmat4.cpp:728
int nrows_val
Definition: newmat.h:452
A matrix is not square exception.
Definition: newmat.h:1981
int tag_val
Definition: newmat.h:451
Upper triangular matrix.
Definition: newmat.h:799
void resize_keep(int)
Definition: newmat4.cpp:506
int Value() const
Definition: newmat.h:234
Square matrix.
Definition: newmat.h:679
MatrixType type() const
Definition: newmat4.cpp:653
MatrixType type() const
Definition: newmat4.cpp:660
MatrixBandWidth bandwidth() const
Definition: newmat4.cpp:730
int Upper() const
Definition: newmat.h:216
void resize(int)
Definition: newmat4.cpp:318
void resize(int, int, int)
Definition: newmat4.cpp:272
int search(const BaseMatrix *) const
Definition: newmat4.cpp:647
FloatVector * d
#define REPORT
Definition: newmat4.cpp:26
bool reuse()
Definition: newmat4.cpp:819
void resize(int)
Definition: newmat4.cpp:315
void swap(GeneralMatrix &gm)
Definition: newmat4.cpp:1264
void resize_keep(int)
Definition: newmat4.cpp:461
Band matrix.
Definition: newmat.h:1096
void cleanup()
Definition: newmat4.cpp:1116
GeneralMatrix * Image() const
Definition: newmat4.cpp:1026
Real * Store() const
Definition: newmat.h:497
void resize(int)
Definition: newmat4.cpp:292
MatrixBandWidth bandwidth() const
Definition: newmat4.cpp:757
MatrixBandWidth bandwidth() const
Definition: newmat4.cpp:756
GeneralMatrix * Image() const
Definition: newmat4.cpp:1012
bool Compare(const MatrixType &source, MatrixType &destination)
Definition: newmat4.cpp:980
GetSubMatrix sym_submatrix(int, int) const
Definition: submat.cpp:39
#define MONITOR_REAL_DELETE(Operation, Size, Pointer)
Definition: myexcept.h:331
MatrixBandWidth bandwidth() const
Definition: newmat4.cpp:678
int search(const BaseMatrix *) const
Definition: newmat4.cpp:650
void MiniCleanUp()
Definition: newmat4.cpp:1119
MatrixBandWidth bandwidth() const
Definition: newmat4.cpp:725
ConstRealStarStar(const Matrix &A)
Definition: newmat4.cpp:1344
void resize_keep(int)
Definition: newmat4.cpp:542
MatrixBandWidth bandwidth() const
Definition: newmat4.cpp:687
void resize_keep(int)
Definition: newmat4.cpp:488
MatrixBandWidth bandwidth() const
Definition: newmat4.cpp:673
Base of the matrix classes.
Definition: newmat.h:292
bool DataLossOK
Definition: newmat.h:140
#define Throw(E)
Definition: myexcept.h:191
MatrixType type() const
Definition: newmat4.cpp:655
MatrixBandWidth bandwidth() const
Definition: newmat4.cpp:693
void Eq(const BaseMatrix &, MatrixType)
Definition: newmat4.cpp:888
MatrixBandWidth bandwidth() const
Definition: newmat4.cpp:675
MatrixBandWidth bandwidth() const
Definition: newmat4.cpp:755
The usual rectangular matrix.
Definition: newmat.h:625
Real * store
Definition: newmat.h:454
#define MONITOR_REAL_NEW(Operation, Size, Pointer)
Definition: myexcept.h:329
InvertedMatrix i() const
Definition: newmat6.cpp:329
friend class Matrix
Definition: newmat.h:582
GeneralMatrix * Image() const
Definition: newmat4.cpp:1005
int ncols_val
Definition: newmat.h:452
int search(const BaseMatrix *) const
Definition: newmat4.cpp:641
MatrixBandWidth bandwidth() const
Definition: newmat4.cpp:753
Incompatible dimensions exception.
Definition: newmat.h:1998
Rectangular matrix for use with Numerical Recipes in C.
Definition: newmat.h:711
MatrixBandWidth bandwidth() const
Definition: newmat4.cpp:681
GeneralMatrix * Image() const
Definition: newmat4.cpp:1019
void resize(int m, int n)
Definition: newmat.h:732
FloatVector FloatVector * a
MatrixType type() const
Definition: newmat4.cpp:656
int Ncols() const
Definition: newmat.h:495
LU decomposition of a band matrix.
Definition: newmat.h:1306
void resize(int i, bool keep=false)
change length, keep data if keep = true
Definition: newmat4.cpp:1228
GeneralMatrix * Image() const
Definition: newmat4.cpp:991
void MatrixErrorNoSpace(const void *)
test for allocation fails
Definition: newmatex.cpp:301
void swap(CroutMatrix &gm)
Definition: newmat4.cpp:1282
Diagonal matrix.
Definition: newmat.h:896
MatrixType type() const
Definition: newmat4.cpp:658
void resize(int)
Definition: newmat4.cpp:303
bool d
Definition: newmat.h:1062
MatrixBandWidth bandwidth() const
Definition: newmat4.cpp:696
Lower triangular matrix.
Definition: newmat.h:848
void cleanup()
Definition: newmat4.cpp:1128
void cleanup()
Definition: newmat4.cpp:1104
int search(const BaseMatrix *bm) const
Definition: newmat4.cpp:638
void swap(BandLUMatrix &gm)
Definition: newmat4.cpp:1309
int ncols() const
Definition: newmat.h:500
Symmetric band matrix.
Definition: newmat.h:1245
static int tristore(int n)
Definition: newmat4.cpp:34
virtual void resize(int, int)
Definition: newmat4.cpp:289
void cleanup()
Definition: newmat4.cpp:1142
SimpleIntArray()
build an array length 0
Definition: newmat.h:1859
void CheckStore() const
Definition: newmat4.cpp:1094
void newmat_block_copy(int n, Real *from, Real *to)
Definition: newmat4.cpp:807
Not defined exception.
Definition: newmat.h:2008
int & operator[](int i)
access element of the array - start at 0
Definition: newmat4.cpp:1179
void SetDataLossOK()
Definition: newmat.h:150
GeneralMatrix * BorrowStore(GeneralMatrix *, MatrixType)
Definition: newmat4.cpp:868
int nrows() const
Definition: newmat.h:499
DiagonalMatrix()
Definition: newmat.h:900
Row vector.
Definition: newmat.h:953
void resize_keep(int)
Definition: newmat4.cpp:586
void Protect()
Definition: newmat.h:509
int Lower() const
Definition: newmat.h:218
RealStarStar(Matrix &A)
Definition: newmat4.cpp:1332
int storage
Definition: newmat.h:453
int * indx
Definition: newmat.h:1061
Column vector.
Definition: newmat.h:1008
void swap(BandMatrix &gm)
Definition: newmat4.cpp:1292
RowVector()
Definition: newmat.h:957
virtual void resize_keep(int, int)
Definition: newmat4.cpp:432
The classes for matrices that can contain data are derived from this.
Definition: newmat.h:447
MatrixBandWidth bandwidth() const
Definition: newmat4.cpp:743
friend class RowVector
Definition: newmat.h:590
virtual GeneralMatrix * Image() const
Definition: newmat4.cpp:1068
void Inject(const MatrixRowCol &)
Definition: newmat2.cpp:69
friend class SymmetricMatrix
Definition: newmat.h:585
GeneralMatrix * Image() const
Definition: newmat4.cpp:1033
Cannot convert to vector exception.
Definition: newmat.h:1972
GetSubMatrix submatrix(int, int, int, int) const
Definition: submat.cpp:27
A matrix which can be of any GeneralMatrix type.
Definition: newmat.h:1397
virtual void SetParameters(const GeneralMatrix *)
Definition: newmat.h:566
void Next()
Definition: newmatrc.h:161
GeneralMatrix * Image() const
Definition: newmat4.cpp:1061
void GetMatrix(const GeneralMatrix *)
Definition: newmat4.cpp:861
MatrixBandWidth bandwidth() const
Definition: newmat4.cpp:672
virtual ~GeneralMatrix()
Definition: newmat4.cpp:168
friend class ReturnMatrix
Definition: newmat.h:617
Identity matrix.
Definition: newmat.h:1350
void cleanup()
Definition: newmat4.cpp:1122
MatrixBandWidth bandwidth() const
Definition: newmat4.cpp:760
Symmetric matrix.
Definition: newmat.h:753


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