Planet
navi homePPSaboutscreenshotsdownloaddevelopmentforum

source: orxonox.OLD/orxonox/trunk/src/util/newmat/newmat6.cpp @ 4565

Last change on this file since 4565 was 4565, checked in by patrick, 19 years ago

orxonox/trunk: added the newmat library to the project. needs some translation in directory, temp under util/newmat. is needed by the collision detection engine to perform lin alg operations such as eigenvector decomposition. perhaps we will make our own library to do that later.
File size: 26.5 KB
Line 
1//$$ newmat6.cpp            Operators, element access, submatrices
2
3// Copyright (C) 1991,2,3,4: R B Davies
4
5#include "include.h"
6
7#include "newmat.h"
8#include "newmatrc.h"
9
10#ifdef use_namespace
11namespace NEWMAT {
12#endif
13
14
15
16#ifdef DO_REPORT
17#define REPORT { static ExeCounter ExeCount(__LINE__,6); ++ExeCount; }
18#else
19#define REPORT {}
20#endif
21
22/*************************** general utilities *************************/
23
24static int tristore(int n)                      // els in triangular matrix
25{ return (n*(n+1))/2; }
26
27
28/****************************** operators *******************************/
29
30Real& Matrix::operator()(int m, int n)
31{
32   REPORT
33   if (m<=0 || m>nrows || n<=0 || n>ncols)
34      Throw(IndexException(m,n,*this));
35   return store[(m-1)*ncols+n-1];
36}
37
38Real& SymmetricMatrix::operator()(int m, int n)
39{
40   REPORT
41   if (m<=0 || n<=0 || m>nrows || n>ncols)
42      Throw(IndexException(m,n,*this));
43   if (m>=n) return store[tristore(m-1)+n-1];
44   else return store[tristore(n-1)+m-1];
45}
46
47Real& UpperTriangularMatrix::operator()(int m, int n)
48{
49   REPORT
50   if (m<=0 || n<m || n>ncols)
51      Throw(IndexException(m,n,*this));
52   return store[(m-1)*ncols+n-1-tristore(m-1)];
53}
54
55Real& LowerTriangularMatrix::operator()(int m, int n)
56{
57   REPORT
58   if (n<=0 || m<n || m>nrows)
59      Throw(IndexException(m,n,*this));
60   return store[tristore(m-1)+n-1];
61}
62
63Real& DiagonalMatrix::operator()(int m, int n)
64{
65   REPORT
66   if (n<=0 || m!=n || m>nrows || n>ncols)
67      Throw(IndexException(m,n,*this));
68   return store[n-1];
69}
70
71Real& DiagonalMatrix::operator()(int m)
72{
73   REPORT
74   if (m<=0 || m>nrows) Throw(IndexException(m,*this));
75   return store[m-1];
76}
77
78Real& ColumnVector::operator()(int m)
79{
80   REPORT
81   if (m<=0 || m> nrows) Throw(IndexException(m,*this));
82   return store[m-1];
83}
84
85Real& RowVector::operator()(int n)
86{
87   REPORT
88   if (n<=0 || n> ncols) Throw(IndexException(n,*this));
89   return store[n-1];
90}
91
92Real& BandMatrix::operator()(int m, int n)
93{
94   REPORT
95   int w = upper+lower+1; int i = lower+n-m;
96   if (m<=0 || m>nrows || n<=0 || n>ncols || i<0 || i>=w)
97      Throw(IndexException(m,n,*this));
98   return store[w*(m-1)+i];
99}
100
101Real& UpperBandMatrix::operator()(int m, int n)
102{
103   REPORT
104   int w = upper+1; int i = n-m;
105   if (m<=0 || m>nrows || n<=0 || n>ncols || i<0 || i>=w)
106      Throw(IndexException(m,n,*this));
107   return store[w*(m-1)+i];
108}
109
110Real& LowerBandMatrix::operator()(int m, int n)
111{
112   REPORT
113   int w = lower+1; int i = lower+n-m;
114   if (m<=0 || m>nrows || n<=0 || n>ncols || i<0 || i>=w)
115      Throw(IndexException(m,n,*this));
116   return store[w*(m-1)+i];
117}
118
119Real& SymmetricBandMatrix::operator()(int m, int n)
120{
121   REPORT
122   int w = lower+1;
123   if (m>=n)
124   {
125      REPORT
126      int i = lower+n-m;
127      if ( m>nrows || n<=0 || i<0 )
128         Throw(IndexException(m,n,*this));
129      return store[w*(m-1)+i];
130   }
131   else
132   {
133      REPORT
134      int i = lower+m-n;
135      if ( n>nrows || m<=0 || i<0 )
136         Throw(IndexException(m,n,*this));
137      return store[w*(n-1)+i];
138   }
139}
140
141
142Real Matrix::operator()(int m, int n) const
143{
144   REPORT
145   if (m<=0 || m>nrows || n<=0 || n>ncols)
146      Throw(IndexException(m,n,*this));
147   return store[(m-1)*ncols+n-1];
148}
149
150Real SymmetricMatrix::operator()(int m, int n) const
151{
152   REPORT
153   if (m<=0 || n<=0 || m>nrows || n>ncols)
154      Throw(IndexException(m,n,*this));
155   if (m>=n) return store[tristore(m-1)+n-1];
156   else return store[tristore(n-1)+m-1];
157}
158
159Real UpperTriangularMatrix::operator()(int m, int n) const
160{
161   REPORT
162   if (m<=0 || n<m || n>ncols)
163      Throw(IndexException(m,n,*this));
164   return store[(m-1)*ncols+n-1-tristore(m-1)];
165}
166
167Real LowerTriangularMatrix::operator()(int m, int n) const
168{
169   REPORT
170   if (n<=0 || m<n || m>nrows)
171      Throw(IndexException(m,n,*this));
172   return store[tristore(m-1)+n-1];
173}
174
175Real DiagonalMatrix::operator()(int m, int n) const
176{
177   REPORT
178   if (n<=0 || m!=n || m>nrows || n>ncols)
179      Throw(IndexException(m,n,*this));
180   return store[n-1];
181}
182
183Real DiagonalMatrix::operator()(int m) const
184{
185   REPORT
186   if (m<=0 || m>nrows) Throw(IndexException(m,*this));
187   return store[m-1];
188}
189
190Real ColumnVector::operator()(int m) const
191{
192   REPORT
193   if (m<=0 || m> nrows) Throw(IndexException(m,*this));
194   return store[m-1];
195}
196
197Real RowVector::operator()(int n) const
198{
199   REPORT
200   if (n<=0 || n> ncols) Throw(IndexException(n,*this));
201   return store[n-1];
202}
203
204Real BandMatrix::operator()(int m, int n) const
205{
206   REPORT
207   int w = upper+lower+1; int i = lower+n-m;
208   if (m<=0 || m>nrows || n<=0 || n>ncols || i<0 || i>=w)
209      Throw(IndexException(m,n,*this));
210   return store[w*(m-1)+i];
211}
212
213Real UpperBandMatrix::operator()(int m, int n) const
214{
215   REPORT
216   int w = upper+1; int i = n-m;
217   if (m<=0 || m>nrows || n<=0 || n>ncols || i<0 || i>=w)
218      Throw(IndexException(m,n,*this));
219   return store[w*(m-1)+i];
220}
221
222Real LowerBandMatrix::operator()(int m, int n) const
223{
224   REPORT
225   int w = lower+1; int i = lower+n-m;
226   if (m<=0 || m>nrows || n<=0 || n>ncols || i<0 || i>=w)
227      Throw(IndexException(m,n,*this));
228   return store[w*(m-1)+i];
229}
230
231Real SymmetricBandMatrix::operator()(int m, int n) const
232{
233   REPORT
234   int w = lower+1;
235   if (m>=n)
236   {
237      REPORT
238      int i = lower+n-m;
239      if ( m>nrows || n<=0 || i<0 )
240         Throw(IndexException(m,n,*this));
241      return store[w*(m-1)+i];
242   }
243   else
244   {
245      REPORT
246      int i = lower+m-n;
247      if ( n>nrows || m<=0 || i<0 )
248         Throw(IndexException(m,n,*this));
249      return store[w*(n-1)+i];
250   }
251}
252
253
254Real BaseMatrix::AsScalar() const
255{
256   REPORT
257   GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
258
259   if (gm->nrows!=1 || gm->ncols!=1)
260   {
261      Tracer tr("AsScalar");
262      Try
263         { Throw(ProgramException("Cannot convert to scalar", *gm)); }
264      CatchAll { gm->tDelete(); ReThrow; }
265   }
266
267   Real x = *(gm->store); gm->tDelete(); return x;
268}
269
270#ifdef TEMPS_DESTROYED_QUICKLY
271
272AddedMatrix& BaseMatrix::operator+(const BaseMatrix& bm) const
273{
274   REPORT
275   AddedMatrix* x = new AddedMatrix(this, &bm);
276   MatrixErrorNoSpace(x);
277   return *x;
278}
279
280SPMatrix& SP(const BaseMatrix& bm1,const BaseMatrix& bm2)
281{
282   REPORT
283   SPMatrix* x = new SPMatrix(&bm1, &bm2);
284   MatrixErrorNoSpace(x);
285   return *x;
286}
287
288KPMatrix& KP(const BaseMatrix& bm1,const BaseMatrix& bm2)
289{
290   REPORT
291   KPMatrix* x = new KPMatrix(&bm1, &bm2);
292   MatrixErrorNoSpace(x);
293   return *x;
294}
295
296MultipliedMatrix& BaseMatrix::operator*(const BaseMatrix& bm) const
297{
298   REPORT
299   MultipliedMatrix* x = new MultipliedMatrix(this, &bm);
300   MatrixErrorNoSpace(x);
301   return *x;
302}
303
304ConcatenatedMatrix& BaseMatrix::operator|(const BaseMatrix& bm) const
305{
306   REPORT
307   ConcatenatedMatrix* x = new ConcatenatedMatrix(this, &bm);
308   MatrixErrorNoSpace(x);
309   return *x;
310}
311
312StackedMatrix& BaseMatrix::operator&(const BaseMatrix& bm) const
313{
314   REPORT
315   StackedMatrix* x = new StackedMatrix(this, &bm);
316   MatrixErrorNoSpace(x);
317   return *x;
318}
319
320//SolvedMatrix& InvertedMatrix::operator*(const BaseMatrix& bmx) const
321SolvedMatrix& InvertedMatrix::operator*(const BaseMatrix& bmx)
322{
323   REPORT
324   SolvedMatrix* x;
325   Try { x = new SolvedMatrix(bm, &bmx); MatrixErrorNoSpace(x); }
326   CatchAll { delete this; ReThrow; }
327   delete this;                // since we are using bm rather than this
328   return *x;
329}
330
331SubtractedMatrix& BaseMatrix::operator-(const BaseMatrix& bm) const
332{
333   REPORT
334   SubtractedMatrix* x = new SubtractedMatrix(this, &bm);
335   MatrixErrorNoSpace(x);
336   return *x;
337}
338
339ShiftedMatrix& BaseMatrix::operator+(Real f) const
340{
341   REPORT
342   ShiftedMatrix* x = new ShiftedMatrix(this, f);
343   MatrixErrorNoSpace(x);
344   return *x;
345}
346
347NegShiftedMatrix& operator-(Real f,const BaseMatrix& bm1)
348{
349   REPORT
350   NegShiftedMatrix* x = new NegShiftedMatrix(f, &bm1);
351   MatrixErrorNoSpace(x);
352   return *x;
353}
354
355ScaledMatrix& BaseMatrix::operator*(Real f) const
356{
357   REPORT
358   ScaledMatrix* x = new ScaledMatrix(this, f);
359   MatrixErrorNoSpace(x);
360   return *x;
361}
362
363ScaledMatrix& BaseMatrix::operator/(Real f) const
364{
365   REPORT
366   ScaledMatrix* x = new ScaledMatrix(this, 1.0/f);
367   MatrixErrorNoSpace(x);
368   return *x;
369}
370
371ShiftedMatrix& BaseMatrix::operator-(Real f) const
372{
373   REPORT
374   ShiftedMatrix* x = new ShiftedMatrix(this, -f);
375   MatrixErrorNoSpace(x);
376   return *x;
377}
378
379TransposedMatrix& BaseMatrix::t() const
380{
381   REPORT
382   TransposedMatrix* x = new TransposedMatrix(this);
383   MatrixErrorNoSpace(x);
384   return *x;
385}
386
387NegatedMatrix& BaseMatrix::operator-() const
388{
389   REPORT
390   NegatedMatrix* x = new NegatedMatrix(this);
391   MatrixErrorNoSpace(x);
392   return *x;
393}
394
395ReversedMatrix& BaseMatrix::Reverse() const
396{
397   REPORT
398   ReversedMatrix* x = new ReversedMatrix(this);
399   MatrixErrorNoSpace(x);
400   return *x;
401}
402
403InvertedMatrix& BaseMatrix::i() const
404{
405   REPORT
406   InvertedMatrix* x = new InvertedMatrix(this);
407   MatrixErrorNoSpace(x);
408   return *x;
409}
410
411RowedMatrix& BaseMatrix::AsRow() const
412{
413   REPORT
414   RowedMatrix* x = new RowedMatrix(this);
415   MatrixErrorNoSpace(x);
416   return *x;
417}
418
419ColedMatrix& BaseMatrix::AsColumn() const
420{
421   REPORT
422   ColedMatrix* x = new ColedMatrix(this);
423   MatrixErrorNoSpace(x);
424   return *x;
425}
426
427DiagedMatrix& BaseMatrix::AsDiagonal() const
428{
429   REPORT
430   DiagedMatrix* x = new DiagedMatrix(this);
431   MatrixErrorNoSpace(x);
432   return *x;
433}
434
435MatedMatrix& BaseMatrix::AsMatrix(int nrx, int ncx) const
436{
437   REPORT
438   MatedMatrix* x = new MatedMatrix(this,nrx,ncx);
439   MatrixErrorNoSpace(x);
440   return *x;
441}
442
443#else
444
445AddedMatrix BaseMatrix::operator+(const BaseMatrix& bm) const
446{ REPORT return AddedMatrix(this, &bm); }
447
448SPMatrix SP(const BaseMatrix& bm1,const BaseMatrix& bm2)
449{ REPORT return SPMatrix(&bm1, &bm2); }
450
451KPMatrix KP(const BaseMatrix& bm1,const BaseMatrix& bm2)
452{ REPORT return KPMatrix(&bm1, &bm2); }
453
454MultipliedMatrix BaseMatrix::operator*(const BaseMatrix& bm) const
455{ REPORT return MultipliedMatrix(this, &bm); }
456
457ConcatenatedMatrix BaseMatrix::operator|(const BaseMatrix& bm) const
458{ REPORT return ConcatenatedMatrix(this, &bm); }
459
460StackedMatrix BaseMatrix::operator&(const BaseMatrix& bm) const
461{ REPORT return StackedMatrix(this, &bm); }
462
463SolvedMatrix InvertedMatrix::operator*(const BaseMatrix& bmx) const
464{ REPORT return SolvedMatrix(bm, &bmx); }
465
466SubtractedMatrix BaseMatrix::operator-(const BaseMatrix& bm) const
467{ REPORT return SubtractedMatrix(this, &bm); }
468
469ShiftedMatrix BaseMatrix::operator+(Real f) const
470{ REPORT return ShiftedMatrix(this, f); }
471
472NegShiftedMatrix operator-(Real f, const BaseMatrix& bm)
473{ REPORT return NegShiftedMatrix(f, &bm); }
474
475ScaledMatrix BaseMatrix::operator*(Real f) const
476{ REPORT return ScaledMatrix(this, f); }
477
478ScaledMatrix BaseMatrix::operator/(Real f) const
479{ REPORT return ScaledMatrix(this, 1.0/f); }
480
481ShiftedMatrix BaseMatrix::operator-(Real f) const
482{ REPORT return ShiftedMatrix(this, -f); }
483
484TransposedMatrix BaseMatrix::t() const
485{ REPORT return TransposedMatrix(this); }
486
487NegatedMatrix BaseMatrix::operator-() const
488{ REPORT return NegatedMatrix(this); }
489
490ReversedMatrix BaseMatrix::Reverse() const
491{ REPORT return ReversedMatrix(this); }
492
493InvertedMatrix BaseMatrix::i() const
494{ REPORT return InvertedMatrix(this); }
495
496
497RowedMatrix BaseMatrix::AsRow() const
498{ REPORT return RowedMatrix(this); }
499
500ColedMatrix BaseMatrix::AsColumn() const
501{ REPORT return ColedMatrix(this); }
502
503DiagedMatrix BaseMatrix::AsDiagonal() const
504{ REPORT return DiagedMatrix(this); }
505
506MatedMatrix BaseMatrix::AsMatrix(int nrx, int ncx) const
507{ REPORT return MatedMatrix(this,nrx,ncx); }
508
509#endif
510
511void GeneralMatrix::operator=(Real f)
512{ REPORT int i=storage; Real* s=store; while (i--) { *s++ = f; } }
513
514void Matrix::operator=(const BaseMatrix& X)
515{
516   REPORT //CheckConversion(X);
517   // MatrixConversionCheck mcc;
518   Eq(X,MatrixType::Rt);
519} 
520
521void RowVector::operator=(const BaseMatrix& X)
522{
523   REPORT  // CheckConversion(X);
524   // MatrixConversionCheck mcc;
525   Eq(X,MatrixType::RV);
526   if (nrows!=1)
527      { Tracer tr("RowVector(=)"); Throw(VectorException(*this)); }
528}
529
530void ColumnVector::operator=(const BaseMatrix& X)
531{
532   REPORT //CheckConversion(X);
533   // MatrixConversionCheck mcc;
534   Eq(X,MatrixType::CV);
535   if (ncols!=1)
536      { Tracer tr("ColumnVector(=)"); Throw(VectorException(*this)); }
537}
538
539void SymmetricMatrix::operator=(const BaseMatrix& X)
540{
541   REPORT // CheckConversion(X);
542   // MatrixConversionCheck mcc;
543   Eq(X,MatrixType::Sm);
544}
545
546void UpperTriangularMatrix::operator=(const BaseMatrix& X)
547{
548   REPORT //CheckConversion(X);
549   // MatrixConversionCheck mcc;
550   Eq(X,MatrixType::UT);
551}
552
553void LowerTriangularMatrix::operator=(const BaseMatrix& X)
554{
555   REPORT //CheckConversion(X);
556   // MatrixConversionCheck mcc;
557   Eq(X,MatrixType::LT);
558}
559
560void DiagonalMatrix::operator=(const BaseMatrix& X)
561{
562   REPORT // CheckConversion(X);
563   // MatrixConversionCheck mcc;
564   Eq(X,MatrixType::Dg);
565}
566
567void IdentityMatrix::operator=(const BaseMatrix& X)
568{
569   REPORT // CheckConversion(X);
570   // MatrixConversionCheck mcc;
571   Eq(X,MatrixType::Id);
572}
573
574void GeneralMatrix::operator<<(const Real* r)
575{
576   REPORT
577   int i = storage; Real* s=store;
578   while(i--) *s++ = *r++;
579}
580
581
582void GenericMatrix::operator=(const GenericMatrix& bmx)
583{
584   if (&bmx != this) { REPORT if (gm) delete gm; gm = bmx.gm->Image();}
585   else { REPORT }
586   gm->Protect();
587}
588
589void GenericMatrix::operator=(const BaseMatrix& bmx)
590{
591   if (gm)
592   {
593      int counter=bmx.search(gm);
594      if (counter==0) { REPORT delete gm; gm=0; }
595      else { REPORT gm->Release(counter); }
596   }
597   else { REPORT }
598   GeneralMatrix* gmx = ((BaseMatrix&)bmx).Evaluate();
599   if (gmx != gm) { REPORT if (gm) delete gm; gm = gmx->Image(); }
600   else { REPORT }
601   gm->Protect();
602}
603
604
605/*************************** += etc ***************************************/
606
607// will also need versions for SubMatrix
608
609
610
611// GeneralMatrix operators
612
613void GeneralMatrix::operator+=(const BaseMatrix& X)
614{
615   REPORT
616   Tracer tr("GeneralMatrix::operator+=");
617   // MatrixConversionCheck mcc;
618   Protect();                                   // so it cannot get deleted
619                                                // during Evaluate
620   GeneralMatrix* gm = ((BaseMatrix&)X).Evaluate();
621#ifdef TEMPS_DESTROYED_QUICKLY
622   AddedMatrix* am = new AddedMatrix(this,gm);
623   MatrixErrorNoSpace(am);
624   if (gm==this) Release(2); else Release();
625   Eq2(*am,Type());
626#else
627   AddedMatrix am(this,gm);
628   if (gm==this) Release(2); else Release();
629   Eq2(am,Type());
630#endif
631}
632
633void GeneralMatrix::operator-=(const BaseMatrix& X)
634{
635   REPORT
636   Tracer tr("GeneralMatrix::operator-=");
637   // MatrixConversionCheck mcc;
638   Protect();                                   // so it cannot get deleted
639                                                // during Evaluate
640   GeneralMatrix* gm = ((BaseMatrix&)X).Evaluate();
641#ifdef TEMPS_DESTROYED_QUICKLY
642   SubtractedMatrix* am = new SubtractedMatrix(this,gm);
643   MatrixErrorNoSpace(am);
644   if (gm==this) Release(2); else Release();
645   Eq2(*am,Type());
646#else
647   SubtractedMatrix am(this,gm);
648   if (gm==this) Release(2); else Release();
649   Eq2(am,Type());
650#endif
651}
652
653void GeneralMatrix::operator*=(const BaseMatrix& X)
654{
655   REPORT
656   Tracer tr("GeneralMatrix::operator*=");
657   // MatrixConversionCheck mcc;
658   Protect();                                   // so it cannot get deleted
659                                                // during Evaluate
660   GeneralMatrix* gm = ((BaseMatrix&)X).Evaluate();
661#ifdef TEMPS_DESTROYED_QUICKLY
662   MultipliedMatrix* am = new MultipliedMatrix(this,gm);
663   MatrixErrorNoSpace(am);
664   if (gm==this) Release(2); else Release();
665   Eq2(*am,Type());
666#else
667   MultipliedMatrix am(this,gm);
668   if (gm==this) Release(2); else Release();
669   Eq2(am,Type());
670#endif
671}
672
673void GeneralMatrix::operator|=(const BaseMatrix& X)
674{
675   REPORT
676   Tracer tr("GeneralMatrix::operator|=");
677   // MatrixConversionCheck mcc;
678   Protect();                                   // so it cannot get deleted
679                                                // during Evaluate
680   GeneralMatrix* gm = ((BaseMatrix&)X).Evaluate();
681#ifdef TEMPS_DESTROYED_QUICKLY
682   ConcatenatedMatrix* am = new ConcatenatedMatrix(this,gm);
683   MatrixErrorNoSpace(am);
684   if (gm==this) Release(2); else Release();
685   Eq2(*am,Type());
686#else
687   ConcatenatedMatrix am(this,gm);
688   if (gm==this) Release(2); else Release();
689   Eq2(am,Type());
690#endif
691}
692
693void GeneralMatrix::operator&=(const BaseMatrix& X)
694{
695   REPORT
696   Tracer tr("GeneralMatrix::operator&=");
697   // MatrixConversionCheck mcc;
698   Protect();                                   // so it cannot get deleted
699                                                // during Evaluate
700   GeneralMatrix* gm = ((BaseMatrix&)X).Evaluate();
701#ifdef TEMPS_DESTROYED_QUICKLY
702   StackedMatrix* am = new StackedMatrix(this,gm);
703   MatrixErrorNoSpace(am);
704   if (gm==this) Release(2); else Release();
705   Eq2(*am,Type());
706#else
707   StackedMatrix am(this,gm);
708   if (gm==this) Release(2); else Release();
709   Eq2(am,Type());
710#endif
711}
712
713void GeneralMatrix::operator+=(Real r)
714{
715   REPORT
716   Tracer tr("GeneralMatrix::operator+=(Real)");
717   // MatrixConversionCheck mcc;
718#ifdef TEMPS_DESTROYED_QUICKLY
719   ShiftedMatrix* am = new ShiftedMatrix(this,r);
720   MatrixErrorNoSpace(am);
721   Release(); Eq2(*am,Type());
722#else
723   ShiftedMatrix am(this,r);
724   Release(); Eq2(am,Type());
725#endif
726}
727
728void GeneralMatrix::operator*=(Real r)
729{
730   REPORT
731   Tracer tr("GeneralMatrix::operator*=(Real)");
732   // MatrixConversionCheck mcc;
733#ifdef TEMPS_DESTROYED_QUICKLY
734   ScaledMatrix* am = new ScaledMatrix(this,r);
735   MatrixErrorNoSpace(am);
736   Release(); Eq2(*am,Type());
737#else
738   ScaledMatrix am(this,r);
739   Release(); Eq2(am,Type());
740#endif
741}
742
743
744// Generic matrix operators
745
746void GenericMatrix::operator+=(const BaseMatrix& X)
747{
748   REPORT
749   Tracer tr("GenericMatrix::operator+=");
750   if (!gm) Throw(ProgramException("GenericMatrix is null"));
751   gm->Protect();            // so it cannot get deleted during Evaluate
752   GeneralMatrix* gmx = ((BaseMatrix&)X).Evaluate();
753#ifdef TEMPS_DESTROYED_QUICKLY
754   AddedMatrix* am = new AddedMatrix(gm,gmx);
755   MatrixErrorNoSpace(am);
756   if (gmx==gm) gm->Release(2); else gm->Release();
757   GeneralMatrix* gmy = am->Evaluate();
758#else
759   AddedMatrix am(gm,gmx);
760   if (gmx==gm) gm->Release(2); else gm->Release();
761   GeneralMatrix* gmy = am.Evaluate();
762#endif
763   if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); }
764   else { REPORT }
765   gm->Protect();
766}
767
768void GenericMatrix::operator-=(const BaseMatrix& X)
769{
770   REPORT
771   Tracer tr("GenericMatrix::operator-=");
772   if (!gm) Throw(ProgramException("GenericMatrix is null"));
773   gm->Protect();            // so it cannot get deleted during Evaluate
774   GeneralMatrix* gmx = ((BaseMatrix&)X).Evaluate();
775#ifdef TEMPS_DESTROYED_QUICKLY
776   SubtractedMatrix* am = new SubtractedMatrix(gm,gmx);
777   MatrixErrorNoSpace(am);
778   if (gmx==gm) gm->Release(2); else gm->Release();
779   GeneralMatrix* gmy = am->Evaluate();
780#else
781   SubtractedMatrix am(gm,gmx);
782   if (gmx==gm) gm->Release(2); else gm->Release();
783   GeneralMatrix* gmy = am.Evaluate();
784#endif
785   if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); }
786   else { REPORT }
787   gm->Protect();
788}
789
790void GenericMatrix::operator*=(const BaseMatrix& X)
791{
792   REPORT
793   Tracer tr("GenericMatrix::operator*=");
794   if (!gm) Throw(ProgramException("GenericMatrix is null"));
795   gm->Protect();            // so it cannot get deleted during Evaluate
796   GeneralMatrix* gmx = ((BaseMatrix&)X).Evaluate();
797#ifdef TEMPS_DESTROYED_QUICKLY
798   MultipliedMatrix* am = new MultipliedMatrix(gm,gmx);
799   MatrixErrorNoSpace(am);
800   if (gmx==gm) gm->Release(2); else gm->Release();
801   GeneralMatrix* gmy = am->Evaluate();
802#else
803   MultipliedMatrix am(gm,gmx);
804   if (gmx==gm) gm->Release(2); else gm->Release();
805   GeneralMatrix* gmy = am.Evaluate();
806#endif
807   if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); }
808   else { REPORT }
809   gm->Protect();
810}
811
812void GenericMatrix::operator|=(const BaseMatrix& X)
813{
814   REPORT
815   Tracer tr("GenericMatrix::operator|=");
816   if (!gm) Throw(ProgramException("GenericMatrix is null"));
817   gm->Protect();            // so it cannot get deleted during Evaluate
818   GeneralMatrix* gmx = ((BaseMatrix&)X).Evaluate();
819#ifdef TEMPS_DESTROYED_QUICKLY
820   ConcatenatedMatrix* am = new ConcatenatedMatrix(gm,gmx);
821   MatrixErrorNoSpace(am);
822   if (gmx==gm) gm->Release(2); else gm->Release();
823   GeneralMatrix* gmy = am->Evaluate();
824#else
825   ConcatenatedMatrix am(gm,gmx);
826   if (gmx==gm) gm->Release(2); else gm->Release();
827   GeneralMatrix* gmy = am.Evaluate();
828#endif
829   if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); }
830   else { REPORT }
831   gm->Protect();
832}
833
834void GenericMatrix::operator&=(const BaseMatrix& X)
835{
836   REPORT
837   Tracer tr("GenericMatrix::operator&=");
838   if (!gm) Throw(ProgramException("GenericMatrix is null"));
839   gm->Protect();            // so it cannot get deleted during Evaluate
840   GeneralMatrix* gmx = ((BaseMatrix&)X).Evaluate();
841#ifdef TEMPS_DESTROYED_QUICKLY
842   StackedMatrix* am = new StackedMatrix(gm,gmx);
843   MatrixErrorNoSpace(am);
844   if (gmx==gm) gm->Release(2); else gm->Release();
845   GeneralMatrix* gmy = am->Evaluate();
846#else
847   StackedMatrix am(gm,gmx);
848   if (gmx==gm) gm->Release(2); else gm->Release();
849   GeneralMatrix* gmy = am.Evaluate();
850#endif
851   if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); }
852   else { REPORT }
853   gm->Protect();
854}
855
856void GenericMatrix::operator+=(Real r)
857{
858   REPORT
859   Tracer tr("GenericMatrix::operator+= (Real)");
860   if (!gm) Throw(ProgramException("GenericMatrix is null"));
861#ifdef TEMPS_DESTROYED_QUICKLY
862   ShiftedMatrix* am = new ShiftedMatrix(gm,r);
863   MatrixErrorNoSpace(am);
864   gm->Release();
865   GeneralMatrix* gmy = am->Evaluate();
866#else
867   ShiftedMatrix am(gm,r);
868   gm->Release();
869   GeneralMatrix* gmy = am.Evaluate();
870#endif
871   if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); }
872   else { REPORT }
873   gm->Protect();
874}
875
876void GenericMatrix::operator*=(Real r)
877{
878   REPORT
879   Tracer tr("GenericMatrix::operator*= (Real)");
880   if (!gm) Throw(ProgramException("GenericMatrix is null"));
881#ifdef TEMPS_DESTROYED_QUICKLY
882   ScaledMatrix* am = new ScaledMatrix(gm,r);
883   MatrixErrorNoSpace(am);
884   gm->Release();
885   GeneralMatrix* gmy = am->Evaluate();
886#else
887   ScaledMatrix am(gm,r);
888   gm->Release();
889   GeneralMatrix* gmy = am.Evaluate();
890#endif
891   if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); }
892   else { REPORT }
893   gm->Protect();
894}
895
896
897/************************* element access *********************************/
898
899Real& Matrix::element(int m, int n)
900{
901   REPORT
902   if (m<0 || m>= nrows || n<0 || n>= ncols)
903      Throw(IndexException(m,n,*this,true));
904   return store[m*ncols+n];
905}
906
907Real Matrix::element(int m, int n) const
908{
909   REPORT
910   if (m<0 || m>= nrows || n<0 || n>= ncols)
911      Throw(IndexException(m,n,*this,true));
912   return store[m*ncols+n];
913}
914
915Real& SymmetricMatrix::element(int m, int n)
916{
917   REPORT
918   if (m<0 || n<0 || m >= nrows || n>=ncols)
919      Throw(IndexException(m,n,*this,true));
920   if (m>=n) return store[tristore(m)+n];
921   else return store[tristore(n)+m];
922}
923
924Real SymmetricMatrix::element(int m, int n) const
925{
926   REPORT
927   if (m<0 || n<0 || m >= nrows || n>=ncols)
928      Throw(IndexException(m,n,*this,true));
929   if (m>=n) return store[tristore(m)+n];
930   else return store[tristore(n)+m];
931}
932
933Real& UpperTriangularMatrix::element(int m, int n)
934{
935   REPORT
936   if (m<0 || n<m || n>=ncols)
937      Throw(IndexException(m,n,*this,true));
938   return store[m*ncols+n-tristore(m)];
939}
940
941Real UpperTriangularMatrix::element(int m, int n) const
942{
943   REPORT
944   if (m<0 || n<m || n>=ncols)
945      Throw(IndexException(m,n,*this,true));
946   return store[m*ncols+n-tristore(m)];
947}
948
949Real& LowerTriangularMatrix::element(int m, int n)
950{
951   REPORT
952   if (n<0 || m<n || m>=nrows)
953      Throw(IndexException(m,n,*this,true));
954   return store[tristore(m)+n];
955}
956
957Real LowerTriangularMatrix::element(int m, int n) const
958{
959   REPORT
960   if (n<0 || m<n || m>=nrows)
961      Throw(IndexException(m,n,*this,true));
962   return store[tristore(m)+n];
963}
964
965Real& DiagonalMatrix::element(int m, int n)
966{
967   REPORT
968   if (n<0 || m!=n || m>=nrows || n>=ncols)
969      Throw(IndexException(m,n,*this,true));
970   return store[n];
971}
972
973Real DiagonalMatrix::element(int m, int n) const
974{
975   REPORT
976   if (n<0 || m!=n || m>=nrows || n>=ncols)
977      Throw(IndexException(m,n,*this,true));
978   return store[n];
979}
980
981Real& DiagonalMatrix::element(int m)
982{
983   REPORT
984   if (m<0 || m>=nrows) Throw(IndexException(m,*this,true));
985   return store[m];
986}
987
988Real DiagonalMatrix::element(int m) const
989{
990   REPORT
991   if (m<0 || m>=nrows) Throw(IndexException(m,*this,true));
992   return store[m];
993}
994
995Real& ColumnVector::element(int m)
996{
997   REPORT
998   if (m<0 || m>= nrows) Throw(IndexException(m,*this,true));
999   return store[m];
1000}
1001
1002Real ColumnVector::element(int m) const
1003{
1004   REPORT
1005   if (m<0 || m>= nrows) Throw(IndexException(m,*this,true));
1006   return store[m];
1007}
1008
1009Real& RowVector::element(int n)
1010{
1011   REPORT
1012   if (n<0 || n>= ncols)  Throw(IndexException(n,*this,true));
1013   return store[n];
1014}
1015
1016Real RowVector::element(int n) const
1017{
1018   REPORT
1019   if (n<0 || n>= ncols)  Throw(IndexException(n,*this,true));
1020   return store[n];
1021}
1022
1023Real& BandMatrix::element(int m, int n)
1024{
1025   REPORT
1026   int w = upper+lower+1; int i = lower+n-m;
1027   if (m<0 || m>= nrows || n<0 || n>= ncols || i<0 || i>=w)
1028      Throw(IndexException(m,n,*this,true));
1029   return store[w*m+i];
1030}
1031
1032Real BandMatrix::element(int m, int n) const
1033{
1034   REPORT
1035   int w = upper+lower+1; int i = lower+n-m;
1036   if (m<0 || m>= nrows || n<0 || n>= ncols || i<0 || i>=w)
1037      Throw(IndexException(m,n,*this,true));
1038   return store[w*m+i];
1039}
1040
1041Real& UpperBandMatrix::element(int m, int n)
1042{
1043   REPORT
1044   int w = upper+1; int i = n-m;
1045   if (m<0 || m>= nrows || n<0 || n>= ncols || i<0 || i>=w)
1046      Throw(IndexException(m,n,*this,true));
1047   return store[w*m+i];
1048}
1049
1050Real UpperBandMatrix::element(int m, int n) const
1051{
1052   REPORT
1053   int w = upper+1; int i = n-m;
1054   if (m<0 || m>= nrows || n<0 || n>= ncols || i<0 || i>=w)
1055      Throw(IndexException(m,n,*this,true));
1056   return store[w*m+i];
1057}
1058
1059Real& LowerBandMatrix::element(int m, int n)
1060{
1061   REPORT
1062   int w = lower+1; int i = lower+n-m;
1063   if (m<0 || m>= nrows || n<0 || n>= ncols || i<0 || i>=w)
1064      Throw(IndexException(m,n,*this,true));
1065   return store[w*m+i];
1066}
1067
1068Real LowerBandMatrix::element(int m, int n) const
1069{
1070   REPORT
1071   int w = lower+1; int i = lower+n-m;
1072   if (m<0 || m>= nrows || n<0 || n>= ncols || i<0 || i>=w)
1073      Throw(IndexException(m,n,*this,true));
1074   return store[w*m+i];
1075}
1076
1077Real& SymmetricBandMatrix::element(int m, int n)
1078{
1079   REPORT
1080   int w = lower+1;
1081   if (m>=n)
1082   {
1083      REPORT
1084      int i = lower+n-m;
1085      if ( m>=nrows || n<0 || i<0 )
1086         Throw(IndexException(m,n,*this,true));
1087      return store[w*m+i];
1088   }
1089   else
1090   {
1091      REPORT
1092      int i = lower+m-n;
1093      if ( n>=nrows || m<0 || i<0 )
1094         Throw(IndexException(m,n,*this,true));
1095      return store[w*n+i];
1096   }
1097}
1098
1099Real SymmetricBandMatrix::element(int m, int n) const
1100{
1101   REPORT
1102   int w = lower+1;
1103   if (m>=n)
1104   {
1105      REPORT
1106      int i = lower+n-m;
1107      if ( m>=nrows || n<0 || i<0 )
1108         Throw(IndexException(m,n,*this,true));
1109      return store[w*m+i];
1110   }
1111   else
1112   {
1113      REPORT
1114      int i = lower+m-n;
1115      if ( n>=nrows || m<0 || i<0 )
1116         Throw(IndexException(m,n,*this,true));
1117      return store[w*n+i];
1118   }
1119}
1120
1121#ifdef use_namespace
1122}
1123#endif
1124
Note: See TracBrowser for help on using the repository browser.