source: ntrip/trunk/BNC/newmat/newmat4.cpp@ 8209

Last change on this file since 8209 was 2013, checked in by mervart, 15 years ago

* empty log message *
File size: 34.2 KB
Line 
1/// \ingroup newmat
2///@{
3
4/// \file newmat4.cpp
5/// Constructors, resize, basic utilities, SimpleIntArray.
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
18namespace 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
34static int tristore(int n) // elements in triangular matrix
35{ return (n*(n+1))/2; }
36
37
38// **************************** constructors ***************************/
39
40GeneralMatrix::GeneralMatrix()
41{ store=0; storage=0; nrows_val=0; ncols_val=0; tag_val=-1; }
42
43GeneralMatrix::GeneralMatrix(ArrayLengthSpecifier s)
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
55Matrix::Matrix(int m, int n) : GeneralMatrix(m*n)
56{ REPORT nrows_val=m; ncols_val=n; }
57
58SquareMatrix::SquareMatrix(ArrayLengthSpecifier n)
59 : Matrix(n.Value(),n.Value())
60{ REPORT }
61
62SymmetricMatrix::SymmetricMatrix(ArrayLengthSpecifier n)
63 : GeneralMatrix(tristore(n.Value()))
64{ REPORT nrows_val=n.Value(); ncols_val=n.Value(); }
65
66UpperTriangularMatrix::UpperTriangularMatrix(ArrayLengthSpecifier n)
67 : GeneralMatrix(tristore(n.Value()))
68{ REPORT nrows_val=n.Value(); ncols_val=n.Value(); }
69
70LowerTriangularMatrix::LowerTriangularMatrix(ArrayLengthSpecifier n)
71 : GeneralMatrix(tristore(n.Value()))
72{ REPORT nrows_val=n.Value(); ncols_val=n.Value(); }
73
74DiagonalMatrix::DiagonalMatrix(ArrayLengthSpecifier m) : GeneralMatrix(m)
75{ REPORT nrows_val=m.Value(); ncols_val=m.Value(); }
76
77Matrix::Matrix(const BaseMatrix& M)
78{
79 REPORT // CheckConversion(M);
80 // MatrixConversionCheck mcc;
81 GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::Rt);
82 GetMatrix(gmx);
83}
84
85SquareMatrix::SquareMatrix(const BaseMatrix& M) : Matrix(M)
86{
87 REPORT
88 if (ncols_val != nrows_val)
89 {
90 Tracer tr("SquareMatrix");
91 Throw(NotSquareException(*this));
92 }
93}
94
95
96SquareMatrix::SquareMatrix(const Matrix& gm)
97{
98 REPORT
99 if (gm.ncols_val != gm.nrows_val)
100 {
101 Tracer tr("SquareMatrix(Matrix)");
102 Throw(NotSquareException(gm));
103 }
104 GetMatrix(&gm);
105}
106
107
108RowVector::RowVector(const BaseMatrix& M) : Matrix(M)
109{
110 REPORT
111 if (nrows_val!=1)
112 {
113 Tracer tr("RowVector");
114 Throw(VectorException(*this));
115 }
116}
117
118ColumnVector::ColumnVector(const BaseMatrix& M) : Matrix(M)
119{
120 REPORT
121 if (ncols_val!=1)
122 {
123 Tracer tr("ColumnVector");
124 Throw(VectorException(*this));
125 }
126}
127
128SymmetricMatrix::SymmetricMatrix(const BaseMatrix& M)
129{
130 REPORT // CheckConversion(M);
131 // MatrixConversionCheck mcc;
132 GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::Sm);
133 GetMatrix(gmx);
134}
135
136UpperTriangularMatrix::UpperTriangularMatrix(const BaseMatrix& M)
137{
138 REPORT // CheckConversion(M);
139 // MatrixConversionCheck mcc;
140 GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::UT);
141 GetMatrix(gmx);
142}
143
144LowerTriangularMatrix::LowerTriangularMatrix(const BaseMatrix& M)
145{
146 REPORT // CheckConversion(M);
147 // MatrixConversionCheck mcc;
148 GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::LT);
149 GetMatrix(gmx);
150}
151
152DiagonalMatrix::DiagonalMatrix(const BaseMatrix& M)
153{
154 REPORT //CheckConversion(M);
155 // MatrixConversionCheck mcc;
156 GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::Dg);
157 GetMatrix(gmx);
158}
159
160IdentityMatrix::IdentityMatrix(const BaseMatrix& M)
161{
162 REPORT //CheckConversion(M);
163 // MatrixConversionCheck mcc;
164 GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::Id);
165 GetMatrix(gmx);
166}
167
168GeneralMatrix::~GeneralMatrix()
169{
170 if (store)
171 {
172 MONITOR_REAL_DELETE("Free (GenMatrix)",storage,store)
173 delete [] store;
174 }
175}
176
177CroutMatrix::CroutMatrix(const BaseMatrix& m)
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
190 GeneralMatrix* gm1 = gm->Evaluate(MatrixType::Rt);
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
200void CroutMatrix::get_aux(CroutMatrix& X)
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
219CroutMatrix::CroutMatrix(const CroutMatrix& gm) : GeneralMatrix()
220{
221 REPORT
222 Tracer tr("CroutMatrix(const CroutMatrix&)");
223 ((CroutMatrix&)gm).get_aux(*this);
224 GetMatrix(&gm);
225}
226
227CroutMatrix::~CroutMatrix()
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
240GeneralMatrix::operator ReturnMatrix() const
241{
242 REPORT
243 GeneralMatrix* gm = Image(); gm->ReleaseAndDelete();
244 return ReturnMatrix(gm);
245}
246
247
248
249ReturnMatrix GeneralMatrix::for_return() const
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
260Matrix::Matrix(Real a, int m, int n) : GeneralMatrix(m * n)
261 { REPORT nrows_val=m; ncols_val=n; operator=(a); }
262
263Matrix::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
272void 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 {
283 store = new Real [storage]; MatrixErrorNoSpace(store);
284 MONITOR_REAL_NEW("Make (ReDimensi)",storage,store)
285 }
286 else store = 0;
287}
288
289void Matrix::resize(int nr, int nc)
290{ REPORT GeneralMatrix::resize(nr,nc,nr*nc); }
291
292void SquareMatrix::resize(int n)
293{ REPORT GeneralMatrix::resize(n,n,n*n); }
294
295void 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
303void SymmetricMatrix::resize(int nr)
304{ REPORT GeneralMatrix::resize(nr,nr,tristore(nr)); }
305
306void UpperTriangularMatrix::resize(int nr)
307{ REPORT GeneralMatrix::resize(nr,nr,tristore(nr)); }
308
309void LowerTriangularMatrix::resize(int nr)
310{ REPORT GeneralMatrix::resize(nr,nr,tristore(nr)); }
311
312void DiagonalMatrix::resize(int nr)
313{ REPORT GeneralMatrix::resize(nr,nr,nr); }
314
315void RowVector::resize(int nc)
316{ REPORT GeneralMatrix::resize(1,nc,nc); }
317
318void ColumnVector::resize(int nr)
319{ REPORT GeneralMatrix::resize(nr,1,nr); }
320
321void 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
328void 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
335void IdentityMatrix::resize(int nr)
336{ REPORT GeneralMatrix::resize(nr,nr,1); *store = 1; }
337
338
339void Matrix::resize(const GeneralMatrix& A)
340{ REPORT resize(A.Nrows(), A.Ncols()); }
341
342void SquareMatrix::resize(const GeneralMatrix& A)
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
354void nricMatrix::resize(const GeneralMatrix& A)
355{ REPORT resize(A.Nrows(), A.Ncols()); }
356
357void ColumnVector::resize(const GeneralMatrix& A)
358{ REPORT resize(A.Nrows(), A.Ncols()); }
359
360void RowVector::resize(const GeneralMatrix& A)
361{ REPORT resize(A.Nrows(), A.Ncols()); }
362
363void SymmetricMatrix::resize(const GeneralMatrix& A)
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
375void DiagonalMatrix::resize(const GeneralMatrix& A)
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
387void UpperTriangularMatrix::resize(const GeneralMatrix& A)
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
399void LowerTriangularMatrix::resize(const GeneralMatrix& A)
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
411void IdentityMatrix::resize(const GeneralMatrix& A)
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
423void GeneralMatrix::resize(const GeneralMatrix&)
424{
425 REPORT
426 Tracer tr("GeneralMatrix::resize(GM)");
427 Throw(NotDefinedException("resize", "this type of matrix"));
428}
429
430//*********************** resize_keep *******************************
431
432void 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
461void SquareMatrix::resize_keep(int nr)
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
479void 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
488void SymmetricMatrix::resize_keep(int nr)
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
506void UpperTriangularMatrix::resize_keep(int nr)
507{
508 Tracer tr("UpperTriangularMatrix::resize_keep");
509 if (nr < nrows_val)
510 {
511 REPORT
512 UpperTriangularMatrix X = sym_submatrix(1,nr);
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
524void LowerTriangularMatrix::resize_keep(int nr)
525{
526 Tracer tr("LowerTriangularMatrix::resize_keep");
527 if (nr < nrows_val)
528 {
529 REPORT
530 LowerTriangularMatrix X = sym_submatrix(1,nr);
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
542void DiagonalMatrix::resize_keep(int nr)
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
560void RowVector::resize_keep(int nc)
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
578void 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
586void ColumnVector::resize_keep(int nr)
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
604void 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/*
614void GeneralMatrix::resizeForAdd(const GeneralMatrix& A, const GeneralMatrix&)
615{ REPORT resize(A); }
616
617void 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
626bool GeneralMatrix::SameStorageType(const GeneralMatrix& A) const
627{
628 REPORT
629 return type() == A.type();
630}
631*/
632
633// ******************* manipulate types, storage **************************/
634
635int GeneralMatrix::search(const BaseMatrix* s) const
636{ REPORT return (s==this) ? 1 : 0; }
637
638int GenericMatrix::search(const BaseMatrix* s) const
639{ REPORT return gm->search(s); }
640
641int MultipliedMatrix::search(const BaseMatrix* s) const
642{ REPORT return bm1->search(s) + bm2->search(s); }
643
644int ShiftedMatrix::search(const BaseMatrix* s) const
645{ REPORT return bm->search(s); }
646
647int NegatedMatrix::search(const BaseMatrix* s) const
648{ REPORT return bm->search(s); }
649
650int ReturnMatrix::search(const BaseMatrix* s) const
651{ REPORT return (s==gm) ? 1 : 0; }
652
653MatrixType Matrix::type() const { return MatrixType::Rt; }
654MatrixType SquareMatrix::type() const { return MatrixType::Sq; }
655MatrixType SymmetricMatrix::type() const { return MatrixType::Sm; }
656MatrixType UpperTriangularMatrix::type() const { return MatrixType::UT; }
657MatrixType LowerTriangularMatrix::type() const { return MatrixType::LT; }
658MatrixType DiagonalMatrix::type() const { return MatrixType::Dg; }
659MatrixType RowVector::type() const { return MatrixType::RV; }
660MatrixType ColumnVector::type() const { return MatrixType::CV; }
661MatrixType CroutMatrix::type() const { return MatrixType::Ct; }
662MatrixType BandMatrix::type() const { return MatrixType::BM; }
663MatrixType UpperBandMatrix::type() const { return MatrixType::UB; }
664MatrixType LowerBandMatrix::type() const { return MatrixType::LB; }
665MatrixType SymmetricBandMatrix::type() const { return MatrixType::SB; }
666
667MatrixType IdentityMatrix::type() const { return MatrixType::Id; }
668
669
670
671MatrixBandWidth BaseMatrix::bandwidth() const { REPORT return -1; }
672MatrixBandWidth DiagonalMatrix::bandwidth() const { REPORT return 0; }
673MatrixBandWidth IdentityMatrix::bandwidth() const { REPORT return 0; }
674
675MatrixBandWidth UpperTriangularMatrix::bandwidth() const
676 { REPORT return MatrixBandWidth(0,-1); }
677
678MatrixBandWidth LowerTriangularMatrix::bandwidth() const
679 { REPORT return MatrixBandWidth(-1,0); }
680
681MatrixBandWidth BandMatrix::bandwidth() const
682 { REPORT return MatrixBandWidth(lower_val,upper_val); }
683
684MatrixBandWidth BandLUMatrix::bandwidth() const
685 { REPORT return MatrixBandWidth(m1,m2); }
686
687MatrixBandWidth GenericMatrix::bandwidth()const
688 { REPORT return gm->bandwidth(); }
689
690MatrixBandWidth AddedMatrix::bandwidth() const
691 { REPORT return gm1->bandwidth() + gm2->bandwidth(); }
692
693MatrixBandWidth SPMatrix::bandwidth() const
694 { REPORT return gm1->bandwidth().minimum(gm2->bandwidth()); }
695
696MatrixBandWidth KPMatrix::bandwidth() const
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
725MatrixBandWidth MultipliedMatrix::bandwidth() const
726{ REPORT return gm1->bandwidth() * gm2->bandwidth(); }
727
728MatrixBandWidth ConcatenatedMatrix::bandwidth() const { REPORT return -1; }
729
730MatrixBandWidth SolvedMatrix::bandwidth() const
731{
732 if (+gm1->type() & MatrixType::Diagonal)
733 { REPORT return gm2->bandwidth(); }
734 else { REPORT return -1; }
735}
736
737MatrixBandWidth ScaledMatrix::bandwidth() const
738 { REPORT return gm->bandwidth(); }
739
740MatrixBandWidth NegatedMatrix::bandwidth() const
741 { REPORT return gm->bandwidth(); }
742
743MatrixBandWidth TransposedMatrix::bandwidth() const
744 { REPORT return gm->bandwidth().t(); }
745
746MatrixBandWidth InvertedMatrix::bandwidth() const
747{
748 if (+gm->type() & MatrixType::Diagonal)
749 { REPORT return MatrixBandWidth(0,0); }
750 else { REPORT return -1; }
751}
752
753MatrixBandWidth RowedMatrix::bandwidth() const { REPORT return -1; }
754MatrixBandWidth ColedMatrix::bandwidth() const { REPORT return -1; }
755MatrixBandWidth DiagedMatrix::bandwidth() const { REPORT return 0; }
756MatrixBandWidth MatedMatrix::bandwidth() const { REPORT return -1; }
757MatrixBandWidth ReturnMatrix::bandwidth() const
758 { REPORT return gm->bandwidth(); }
759
760MatrixBandWidth GetSubMatrix::bandwidth() const
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
786void GeneralMatrix::tDelete()
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
807void 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
819bool GeneralMatrix::reuse()
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)
828 newmat_block_copy(storage, store, s); store = 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
838Real* GeneralMatrix::GetStore()
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)
847 newmat_block_copy(storage, store, 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
861void GeneralMatrix::GetMatrix(const GeneralMatrix* gmx)
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
868GeneralMatrix* GeneralMatrix::BorrowStore(GeneralMatrix* gmx, MatrixType mt)
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
888void GeneralMatrix::Eq(const BaseMatrix& X, MatrixType mt)
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
928void GeneralMatrix::Eq(const GeneralMatrix& X)
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<<
946void GeneralMatrix::Eq(const BaseMatrix& X, MatrixType mt, bool ldok)
947{
948 REPORT
949 if (ldok) mt.SetDataLossOK();
950 Eq(X, mt);
951}
952
953void GeneralMatrix::Eq2(const BaseMatrix& X, MatrixType mt)
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
965void GeneralMatrix::inject(const GeneralMatrix& X)
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)
971 Throw(IncompatibleDimensionsException());
972 MatrixRow mr((GeneralMatrix*)&X, LoadOnEntry);
973 MatrixRow mrx(this, LoadOnEntry+StoreOnExit+DirectPart);
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
980bool 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
991GeneralMatrix* Matrix::Image() const
992{
993 REPORT
994 GeneralMatrix* gm = new Matrix(*this); MatrixErrorNoSpace(gm);
995 return gm;
996}
997
998GeneralMatrix* SquareMatrix::Image() const
999{
1000 REPORT
1001 GeneralMatrix* gm = new SquareMatrix(*this); MatrixErrorNoSpace(gm);
1002 return gm;
1003}
1004
1005GeneralMatrix* SymmetricMatrix::Image() const
1006{
1007 REPORT
1008 GeneralMatrix* gm = new SymmetricMatrix(*this); MatrixErrorNoSpace(gm);
1009 return gm;
1010}
1011
1012GeneralMatrix* UpperTriangularMatrix::Image() const
1013{
1014 REPORT
1015 GeneralMatrix* gm = new UpperTriangularMatrix(*this);
1016 MatrixErrorNoSpace(gm); return gm;
1017}
1018
1019GeneralMatrix* LowerTriangularMatrix::Image() const
1020{
1021 REPORT
1022 GeneralMatrix* gm = new LowerTriangularMatrix(*this);
1023 MatrixErrorNoSpace(gm); return gm;
1024}
1025
1026GeneralMatrix* DiagonalMatrix::Image() const
1027{
1028 REPORT
1029 GeneralMatrix* gm = new DiagonalMatrix(*this); MatrixErrorNoSpace(gm);
1030 return gm;
1031}
1032
1033GeneralMatrix* RowVector::Image() const
1034{
1035 REPORT
1036 GeneralMatrix* gm = new RowVector(*this); MatrixErrorNoSpace(gm);
1037 return gm;
1038}
1039
1040GeneralMatrix* ColumnVector::Image() const
1041{
1042 REPORT
1043 GeneralMatrix* gm = new ColumnVector(*this); MatrixErrorNoSpace(gm);
1044 return gm;
1045}
1046
1047GeneralMatrix* nricMatrix::Image() const
1048{
1049 REPORT
1050 GeneralMatrix* gm = new nricMatrix(*this); MatrixErrorNoSpace(gm);
1051 return gm;
1052}
1053
1054GeneralMatrix* IdentityMatrix::Image() const
1055{
1056 REPORT
1057 GeneralMatrix* gm = new IdentityMatrix(*this); MatrixErrorNoSpace(gm);
1058 return gm;
1059}
1060
1061GeneralMatrix* CroutMatrix::Image() const
1062{
1063 REPORT
1064 GeneralMatrix* gm = new CroutMatrix(*this); MatrixErrorNoSpace(gm);
1065 return gm;
1066}
1067
1068GeneralMatrix* GeneralMatrix::Image() const
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
1079void nricMatrix::MakeRowPointer()
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
1091void nricMatrix::DeleteRowPointer()
1092 { REPORT if (nrows_val) delete [] row_pointer; }
1093
1094void GeneralMatrix::CheckStore() const
1095{
1096 REPORT
1097 if (!store)
1098 Throw(ProgramException("NRIC accessing matrix with unset dimensions"));
1099}
1100
1101
1102// *************************** CleanUp routines *****************************/
1103
1104void GeneralMatrix::cleanup()
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
1116void nricMatrix::cleanup()
1117 { REPORT DeleteRowPointer(); GeneralMatrix::cleanup(); }
1118
1119void nricMatrix::MiniCleanUp()
1120 { REPORT DeleteRowPointer(); GeneralMatrix::MiniCleanUp(); }
1121
1122void RowVector::cleanup()
1123 { REPORT GeneralMatrix::cleanup(); nrows_val=1; }
1124
1125void ColumnVector::cleanup()
1126 { REPORT GeneralMatrix::cleanup(); ncols_val=1; }
1127
1128void CroutMatrix::cleanup()
1129{
1130 REPORT
1131 if (nrows_val) delete [] indx;
1132 GeneralMatrix::cleanup();
1133}
1134
1135void CroutMatrix::MiniCleanUp()
1136{
1137 REPORT
1138 if (nrows_val) delete [] indx;
1139 GeneralMatrix::MiniCleanUp();
1140}
1141
1142void BandLUMatrix::cleanup()
1143{
1144 REPORT
1145 if (nrows_val) delete [] indx;
1146 if (storage2) delete [] store2;
1147 GeneralMatrix::cleanup();
1148}
1149
1150void BandLUMatrix::MiniCleanUp()
1151{
1152 REPORT
1153 if (nrows_val) delete [] indx;
1154 if (storage2) delete [] store2;
1155 GeneralMatrix::MiniCleanUp();
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
1163SimpleIntArray::SimpleIntArray(int xn) : n(xn)
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
1172SimpleIntArray::~SimpleIntArray() { REPORT if (a) delete [] a; }
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
1179int& SimpleIntArray::operator[](int i)
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
1189int SimpleIntArray::operator[](int i) const
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
1198void SimpleIntArray::operator=(int ai)
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
1204void SimpleIntArray::operator=(const SimpleIntArray& b)
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
1214SimpleIntArray::SimpleIntArray(const SimpleIntArray& b) : Janitor(), n(b.n)
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
1228void 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
1264void GeneralMatrix::swap(GeneralMatrix& gm)
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
1275void nricMatrix::swap(nricMatrix& gm)
1276{
1277 REPORT
1278 GeneralMatrix::swap((GeneralMatrix&)gm);
1279 Real** rp = row_pointer; row_pointer = gm.row_pointer; gm.row_pointer = rp;
1280}
1281
1282void CroutMatrix::swap(CroutMatrix& gm)
1283{
1284 REPORT
1285 GeneralMatrix::swap((GeneralMatrix&)gm);
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
1292void BandMatrix::swap(BandMatrix& gm)
1293{
1294 REPORT
1295 GeneralMatrix::swap((GeneralMatrix&)gm);
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
1301void SymmetricBandMatrix::swap(SymmetricBandMatrix& gm)
1302{
1303 REPORT
1304 GeneralMatrix::swap((GeneralMatrix&)gm);
1305 int i;
1306 i = lower_val; lower_val = gm.lower_val; gm.lower_val = i;
1307}
1308
1309void BandLUMatrix::swap(BandLUMatrix& gm)
1310{
1311 REPORT
1312 GeneralMatrix::swap((GeneralMatrix&)gm);
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
1324void GenericMatrix::swap(GenericMatrix& x)
1325{
1326 REPORT
1327 GeneralMatrix* tgm = gm; gm = x.gm; x.gm = tgm;
1328}
1329
1330// ********************** C subscript classes ****************************
1331
1332RealStarStar::RealStarStar(Matrix& A)
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
1344ConstRealStarStar::ConstRealStarStar(const Matrix& A)
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
1363/// \fn GeneralMatrix::SimpleAddOK(const GeneralMatrix* gm)
1364/// Can we add two matrices with simple vector add.
1365/// SimpleAddOK shows when we can add two matrices by a simple vector add
1366/// and when we can add one matrix into another
1367///
1368/// *gm must be the same type as *this
1369/// - return 0 if simple add is OK
1370/// - return 1 if we can add into *gm only
1371/// - return 2 if we can add into *this only
1372/// - return 3 if we can't add either way
1373///
1374/// Also applies to subtract;
1375/// for SP this will still be valid if we swap 1 and 2
1376///
1377/// For types Matrix, DiagonalMatrix, UpperTriangularMatrix,
1378/// LowerTriangularMatrix, SymmetricMatrix etc return 0.
1379/// For band matrices we will need to check bandwidths.
1380
1381
1382
1383
1384
1385
1386
1387///@}
Note: See TracBrowser for help on using the repository browser.