source: nanovis/branches/1.1/newmat11/newmat7.cpp @ 4906

Last change on this file since 4906 was 2096, checked in by ldelgass, 13 years ago

Normalize line endings, set eol-style to native on *.cpp, *.h files

  • Property svn:eol-style set to native
File size: 30.4 KB
Line 
1/// \ingroup newmat
2///@{
3
4/// \file newmat7.cpp
5/// Invert, solve, binary operations.
6
7// Copyright (C) 1991,2,3,4: R B Davies
8
9#include "include.h"
10
11#include "newmat.h"
12#include "newmatrc.h"
13
14#ifdef use_namespace
15namespace NEWMAT {
16#endif
17
18
19#ifdef DO_REPORT
20#define REPORT { static ExeCounter ExeCount(__LINE__,7); ++ExeCount; }
21#else
22#define REPORT {}
23#endif
24
25
26//***************************** solve routines ******************************/
27
28GeneralMatrix* GeneralMatrix::MakeSolver()
29{
30   REPORT
31   GeneralMatrix* gm = new CroutMatrix(*this);
32   MatrixErrorNoSpace(gm); gm->ReleaseAndDelete(); return gm;
33}
34
35GeneralMatrix* Matrix::MakeSolver()
36{
37   REPORT
38   GeneralMatrix* gm = new CroutMatrix(*this);
39   MatrixErrorNoSpace(gm); gm->ReleaseAndDelete(); return gm;
40}
41
42void CroutMatrix::Solver(MatrixColX& mcout, const MatrixColX& mcin)
43{
44   REPORT
45   int i = mcin.skip; Real* el = mcin.data-i; Real* el1 = el;
46   while (i--) *el++ = 0.0;
47   el += mcin.storage; i = nrows_val - mcin.skip - mcin.storage;
48   while (i--) *el++ = 0.0;
49   lubksb(el1, mcout.skip);
50}
51
52
53// Do we need check for entirely zero output?
54
55void UpperTriangularMatrix::Solver(MatrixColX& mcout,
56   const MatrixColX& mcin)
57{
58   REPORT
59   int i = mcin.skip-mcout.skip; Real* elx = mcin.data-i;
60   while (i-- > 0) *elx++ = 0.0;
61   int nr = mcin.skip+mcin.storage;
62   elx = mcin.data+mcin.storage; Real* el = elx;
63   int j = mcout.skip+mcout.storage-nr;
64   int nc = ncols_val-nr; i = nr-mcout.skip;
65   while (j-- > 0) *elx++ = 0.0;
66   Real* Ael = store + (nr*(2*ncols_val-nr+1))/2; j = 0;
67   while (i-- > 0)
68   {
69      elx = el; Real sum = 0.0; int jx = j++; Ael -= nc;
70      while (jx--) sum += *(--Ael) * *(--elx);
71      elx--; *elx = (*elx - sum) / *(--Ael);
72   }
73}
74
75void LowerTriangularMatrix::Solver(MatrixColX& mcout,
76   const MatrixColX& mcin)
77{
78   REPORT
79   int i = mcin.skip-mcout.skip; Real* elx = mcin.data-i;
80   while (i-- > 0) *elx++ = 0.0;
81   int nc = mcin.skip; i = nc+mcin.storage; elx = mcin.data+mcin.storage;
82   int nr = mcout.skip+mcout.storage; int j = nr-i; i = nr-nc;
83   while (j-- > 0) *elx++ = 0.0;
84   Real* el = mcin.data; Real* Ael = store + (nc*(nc+1))/2; j = 0;
85   while (i-- > 0)
86   {
87      elx = el; Real sum = 0.0; int jx = j++; Ael += nc;
88      while (jx--) sum += *Ael++ * *elx++;
89      *elx = (*elx - sum) / *Ael++;
90   }
91}
92
93//******************* carry out binary operations *************************/
94
95static GeneralMatrix*
96   GeneralMult(GeneralMatrix*,GeneralMatrix*,MultipliedMatrix*,MatrixType);
97static GeneralMatrix*
98   GeneralSolv(GeneralMatrix*,GeneralMatrix*,BaseMatrix*,MatrixType);
99static GeneralMatrix*
100   GeneralSolvI(GeneralMatrix*,BaseMatrix*,MatrixType);
101static GeneralMatrix*
102   GeneralKP(GeneralMatrix*,GeneralMatrix*,KPMatrix*,MatrixType);
103
104GeneralMatrix* MultipliedMatrix::Evaluate(MatrixType mt)
105{
106   REPORT
107   gm2 = ((BaseMatrix*&)bm2)->Evaluate();
108   gm2 = gm2->Evaluate(gm2->type().MultRHS());     // no symmetric on RHS
109   gm1 = ((BaseMatrix*&)bm1)->Evaluate();
110   return GeneralMult(gm1, gm2, this, mt);
111}
112
113GeneralMatrix* SolvedMatrix::Evaluate(MatrixType mt)
114{
115   REPORT
116   gm1 = ((BaseMatrix*&)bm1)->Evaluate();
117   gm2 = ((BaseMatrix*&)bm2)->Evaluate();
118   return GeneralSolv(gm1,gm2,this,mt);
119}
120
121GeneralMatrix* KPMatrix::Evaluate(MatrixType mt)
122{
123   REPORT
124   gm1 = ((BaseMatrix*&)bm1)->Evaluate();
125   gm2 = ((BaseMatrix*&)bm2)->Evaluate();
126   return GeneralKP(gm1,gm2,this,mt);
127}
128
129// routines for adding or subtracting matrices of identical storage structure
130
131static void Add(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
132{
133   REPORT
134   Real* s1=gm1->Store(); Real* s2=gm2->Store();
135   Real* s=gm->Store(); int i=gm->Storage() >> 2;
136   while (i--)
137   {
138       *s++ = *s1++ + *s2++; *s++ = *s1++ + *s2++;
139       *s++ = *s1++ + *s2++; *s++ = *s1++ + *s2++;
140   }
141   i=gm->Storage() & 3; while (i--) *s++ = *s1++ + *s2++;
142}
143
144static void AddTo(GeneralMatrix* gm, const GeneralMatrix* gm2)
145{
146   REPORT
147   const Real* s2=gm2->Store(); Real* s=gm->Store(); int i=gm->Storage() >> 2;
148   while (i--)
149   { *s++ += *s2++; *s++ += *s2++; *s++ += *s2++; *s++ += *s2++; }
150   i=gm->Storage() & 3; while (i--) *s++ += *s2++;
151}
152
153void GeneralMatrix::PlusEqual(const GeneralMatrix& gm)
154{
155   REPORT
156   if (nrows_val != gm.nrows_val || ncols_val != gm.ncols_val)
157      Throw(IncompatibleDimensionsException(*this, gm));
158   AddTo(this, &gm);
159}
160
161static void Subtract(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
162{
163   REPORT
164   Real* s1=gm1->Store(); Real* s2=gm2->Store();
165   Real* s=gm->Store(); int i=gm->Storage() >> 2;
166   while (i--)
167   {
168       *s++ = *s1++ - *s2++; *s++ = *s1++ - *s2++;
169       *s++ = *s1++ - *s2++; *s++ = *s1++ - *s2++;
170   }
171   i=gm->Storage() & 3; while (i--) *s++ = *s1++ - *s2++;
172}
173
174static void SubtractFrom(GeneralMatrix* gm, const GeneralMatrix* gm2)
175{
176   REPORT
177   const Real* s2=gm2->Store(); Real* s=gm->Store(); int i=gm->Storage() >> 2;
178   while (i--)
179   { *s++ -= *s2++; *s++ -= *s2++; *s++ -= *s2++; *s++ -= *s2++; }
180   i=gm->Storage() & 3; while (i--) *s++ -= *s2++;
181}
182
183void GeneralMatrix::MinusEqual(const GeneralMatrix& gm)
184{
185   REPORT
186   if (nrows_val != gm.nrows_val || ncols_val != gm.ncols_val)
187      Throw(IncompatibleDimensionsException(*this, gm));
188   SubtractFrom(this, &gm);
189}
190
191static void ReverseSubtract(GeneralMatrix* gm, const GeneralMatrix* gm2)
192{
193   REPORT
194   const Real* s2=gm2->Store(); Real* s=gm->Store(); int i=gm->Storage() >> 2;
195   while (i--)
196   {
197      *s = *s2++ - *s; s++; *s = *s2++ - *s; s++;
198      *s = *s2++ - *s; s++; *s = *s2++ - *s; s++;
199   }
200   i=gm->Storage() & 3; while (i--) { *s = *s2++ - *s; s++; }
201}
202
203static void SP(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
204{
205   REPORT
206   Real* s1=gm1->Store(); Real* s2=gm2->Store();
207   Real* s=gm->Store(); int i=gm->Storage() >> 2;
208   while (i--)
209   {
210       *s++ = *s1++ * *s2++; *s++ = *s1++ * *s2++;
211       *s++ = *s1++ * *s2++; *s++ = *s1++ * *s2++;
212   }
213   i=gm->Storage() & 3; while (i--) *s++ = *s1++ * *s2++;
214}
215
216static void SP(GeneralMatrix* gm, GeneralMatrix* gm2)
217{
218   REPORT
219   Real* s2=gm2->Store(); Real* s=gm->Store(); int i=gm->Storage() >> 2;
220   while (i--)
221   { *s++ *= *s2++; *s++ *= *s2++; *s++ *= *s2++; *s++ *= *s2++; }
222   i=gm->Storage() & 3; while (i--) *s++ *= *s2++;
223}
224
225// routines for adding or subtracting matrices of different storage structure
226
227static void AddDS(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
228{
229   REPORT
230   int nr = gm->Nrows();
231   MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
232   MatrixRow mr(gm, StoreOnExit+DirectPart);
233   while (nr--) { mr.Add(mr1,mr2); mr1.Next(); mr2.Next(); mr.Next(); }
234}
235
236static void AddDS(GeneralMatrix* gm, GeneralMatrix* gm2)
237// Add into first argument
238{
239   REPORT
240   int nr = gm->Nrows();
241   MatrixRow mr(gm, StoreOnExit+LoadOnEntry+DirectPart);
242   MatrixRow mr2(gm2, LoadOnEntry);
243   while (nr--) { mr.Add(mr2); mr.Next(); mr2.Next(); }
244}
245
246static void SubtractDS
247   (GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
248{
249   REPORT
250   int nr = gm->Nrows();
251   MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
252   MatrixRow mr(gm, StoreOnExit+DirectPart);
253   while (nr--) { mr.Sub(mr1,mr2); mr1.Next(); mr2.Next(); mr.Next(); }
254}
255
256static void SubtractDS(GeneralMatrix* gm, GeneralMatrix* gm2)
257{
258   REPORT
259   int nr = gm->Nrows();
260   MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart);
261   MatrixRow mr2(gm2, LoadOnEntry);
262   while (nr--) { mr.Sub(mr2); mr.Next(); mr2.Next(); }
263}
264
265static void ReverseSubtractDS(GeneralMatrix* gm, GeneralMatrix* gm2)
266{
267   REPORT
268   int nr = gm->Nrows();
269   MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart);
270   MatrixRow mr2(gm2, LoadOnEntry);
271   while (nr--) { mr.RevSub(mr2); mr2.Next(); mr.Next(); }
272}
273
274static void SPDS(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
275{
276   REPORT
277   int nr = gm->Nrows();
278   MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
279   MatrixRow mr(gm, StoreOnExit+DirectPart);
280   while (nr--) { mr.Multiply(mr1,mr2); mr1.Next(); mr2.Next(); mr.Next(); }
281}
282
283static void SPDS(GeneralMatrix* gm, GeneralMatrix* gm2)
284// SP into first argument
285{
286   REPORT
287   int nr = gm->Nrows();
288   MatrixRow mr(gm, StoreOnExit+LoadOnEntry+DirectPart);
289   MatrixRow mr2(gm2, LoadOnEntry);
290   while (nr--) { mr.Multiply(mr2); mr.Next(); mr2.Next(); }
291}
292
293static GeneralMatrix* GeneralMult1(GeneralMatrix* gm1, GeneralMatrix* gm2,
294   MultipliedMatrix* mm, MatrixType mtx)
295{
296   REPORT
297   Tracer tr("GeneralMult1");
298   int nr=gm1->Nrows(); int nc=gm2->Ncols();
299   if (gm1->Ncols() !=gm2->Nrows())
300      Throw(IncompatibleDimensionsException(*gm1, *gm2));
301   GeneralMatrix* gmx = mtx.New(nr,nc,mm);
302
303   MatrixCol mcx(gmx, StoreOnExit+DirectPart);
304   MatrixCol mc2(gm2, LoadOnEntry);
305   while (nc--)
306   {
307      MatrixRow mr1(gm1, LoadOnEntry, mcx.Skip());
308      Real* el = mcx.Data();                         // pointer to an element
309      int n = mcx.Storage();
310      while (n--) { *(el++) = DotProd(mr1,mc2); mr1.Next(); }
311      mc2.Next(); mcx.Next();
312   }
313   gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gmx;
314}
315
316static GeneralMatrix* GeneralMult2(GeneralMatrix* gm1, GeneralMatrix* gm2,
317   MultipliedMatrix* mm, MatrixType mtx)
318{
319   // version that accesses by row only - not good for thin matrices
320   // or column vectors in right hand term.
321   REPORT
322   Tracer tr("GeneralMult2");
323   int nr=gm1->Nrows(); int nc=gm2->Ncols();
324   if (gm1->Ncols() !=gm2->Nrows())
325      Throw(IncompatibleDimensionsException(*gm1, *gm2));
326   GeneralMatrix* gmx = mtx.New(nr,nc,mm);
327
328   MatrixRow mrx(gmx, LoadOnEntry+StoreOnExit+DirectPart);
329   MatrixRow mr1(gm1, LoadOnEntry);
330   while (nr--)
331   {
332      MatrixRow mr2(gm2, LoadOnEntry, mr1.Skip());
333      Real* el = mr1.Data();                         // pointer to an element
334      int n = mr1.Storage();
335      mrx.Zero();
336      while (n--) { mrx.AddScaled(mr2, *el++); mr2.Next(); }
337      mr1.Next(); mrx.Next();
338   }
339   gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gmx;
340}
341
342static GeneralMatrix* mmMult(GeneralMatrix* gm1, GeneralMatrix* gm2)
343{
344   // matrix multiplication for type Matrix only
345   REPORT
346   Tracer tr("MatrixMult");
347
348   int nr=gm1->Nrows(); int ncr=gm1->Ncols(); int nc=gm2->Ncols();
349   if (ncr != gm2->Nrows()) Throw(IncompatibleDimensionsException(*gm1,*gm2));
350
351   Matrix* gm = new Matrix(nr,nc); MatrixErrorNoSpace(gm);
352
353   Real* s1=gm1->Store(); Real* s2=gm2->Store(); Real* s=gm->Store();
354
355   if (ncr)
356   {
357      while (nr--)
358      {
359         Real* s2x = s2; int j = ncr;
360         Real* sx = s; Real f = *s1++; int k = nc;
361         while (k--) *sx++ = f * *s2x++;
362         while (--j)
363            { sx = s; f = *s1++; k = nc; while (k--) *sx++ += f * *s2x++; }
364         s = sx;
365      }
366   }
367   else *gm = 0.0;
368
369   gm->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gm;
370}
371
372static GeneralMatrix* GeneralMult(GeneralMatrix* gm1, GeneralMatrix* gm2,
373   MultipliedMatrix* mm, MatrixType mtx)
374{
375   if ( Rectangular(gm1->type(), gm2->type(), mtx))
376      { REPORT return mmMult(gm1, gm2); }
377   Compare(gm1->type() * gm2->type(),mtx);
378   int nr = gm2->Nrows(); int nc = gm2->Ncols();
379   if (nc <= 5 && nr > nc) { REPORT return GeneralMult1(gm1, gm2, mm, mtx); }
380   REPORT return GeneralMult2(gm1, gm2, mm, mtx);
381}
382
383static GeneralMatrix* GeneralKP(GeneralMatrix* gm1, GeneralMatrix* gm2,
384   KPMatrix* kp, MatrixType mtx)
385{
386   REPORT
387   Tracer tr("GeneralKP");
388   int nr1 = gm1->Nrows(); int nc1 = gm1->Ncols();
389   int nr2 = gm2->Nrows(); int nc2 = gm2->Ncols();
390   Compare((gm1->type()).KP(gm2->type()),mtx);
391   GeneralMatrix* gmx = mtx.New(nr1*nr2, nc1*nc2, kp);
392   MatrixRow mrx(gmx, LoadOnEntry+StoreOnExit+DirectPart);
393   MatrixRow mr1(gm1, LoadOnEntry);
394   for (int i = 1; i <= nr1; ++i)
395   {
396      MatrixRow mr2(gm2, LoadOnEntry);
397      for (int j = 1; j <= nr2; ++j)
398         { mrx.KP(mr1,mr2); mr2.Next(); mrx.Next(); }
399      mr1.Next();
400   }
401   gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gmx;
402}
403
404static GeneralMatrix* GeneralSolv(GeneralMatrix* gm1, GeneralMatrix* gm2,
405   BaseMatrix* sm, MatrixType mtx)
406{
407   REPORT
408   Tracer tr("GeneralSolv");
409   Compare(gm1->type().i() * gm2->type(),mtx);
410   int nr = gm1->Nrows();
411   if (nr != gm1->Ncols()) Throw(NotSquareException(*gm1));
412   int nc = gm2->Ncols();
413   if (gm1->Ncols() != gm2->Nrows())
414      Throw(IncompatibleDimensionsException(*gm1, *gm2));
415   GeneralMatrix* gmx = mtx.New(nr,nc,sm); MatrixErrorNoSpace(gmx);
416   Real* r = new Real [nr]; MatrixErrorNoSpace(r);
417   MONITOR_REAL_NEW("Make   (GenSolv)",nr,r)
418   GeneralMatrix* gms = gm1->MakeSolver();
419   Try
420   {
421
422      MatrixColX mcx(gmx, r, StoreOnExit+DirectPart);   // copy to and from r
423         // this must be inside Try so mcx is destroyed before gmx
424      MatrixColX mc2(gm2, r, LoadOnEntry);
425      while (nc--) { gms->Solver(mcx, mc2); mcx.Next(); mc2.Next(); }
426   }
427   CatchAll
428   {
429      if (gms) gms->tDelete();
430      delete gmx;                   // <--------------------
431      gm2->tDelete();
432      MONITOR_REAL_DELETE("Delete (GenSolv)",nr,r)
433                          // AT&T version 2.1 gives an internal error
434      delete [] r;
435      ReThrow;
436   }
437   gms->tDelete(); gmx->ReleaseAndDelete(); gm2->tDelete();
438   MONITOR_REAL_DELETE("Delete (GenSolv)",nr,r)
439                          // AT&T version 2.1 gives an internal error
440   delete [] r;
441   return gmx;
442}
443
444// version for inverses - gm2 is identity
445static GeneralMatrix* GeneralSolvI(GeneralMatrix* gm1, BaseMatrix* sm,
446   MatrixType mtx)
447{
448   REPORT
449   Tracer tr("GeneralSolvI");
450   Compare(gm1->type().i(),mtx);
451   int nr = gm1->Nrows();
452   if (nr != gm1->Ncols()) Throw(NotSquareException(*gm1));
453   int nc = nr;
454   // DiagonalMatrix I(nr); I = 1;
455   IdentityMatrix I(nr);
456   GeneralMatrix* gmx = mtx.New(nr,nc,sm); MatrixErrorNoSpace(gmx);
457   Real* r = new Real [nr]; MatrixErrorNoSpace(r);
458   MONITOR_REAL_NEW("Make   (GenSolvI)",nr,r)
459   GeneralMatrix* gms = gm1->MakeSolver();
460   Try
461   {
462
463      MatrixColX mcx(gmx, r, StoreOnExit+DirectPart);   // copy to and from r
464         // this must be inside Try so mcx is destroyed before gmx
465      MatrixColX mc2(&I, r, LoadOnEntry);
466      while (nc--) { gms->Solver(mcx, mc2); mcx.Next(); mc2.Next(); }
467   }
468   CatchAll
469   {
470      if (gms) gms->tDelete();
471      delete gmx;
472      MONITOR_REAL_DELETE("Delete (GenSolvI)",nr,r)
473                          // AT&T version 2.1 gives an internal error
474      delete [] r;
475      ReThrow;
476   }
477   gms->tDelete(); gmx->ReleaseAndDelete();
478   MONITOR_REAL_DELETE("Delete (GenSolvI)",nr,r)
479                          // AT&T version 2.1 gives an internal error
480   delete [] r;
481   return gmx;
482}
483
484GeneralMatrix* InvertedMatrix::Evaluate(MatrixType mtx)
485{
486   // Matrix Inversion - use solve routines
487   Tracer tr("InvertedMatrix::Evaluate");
488   REPORT
489   gm=((BaseMatrix*&)bm)->Evaluate();
490   return GeneralSolvI(gm,this,mtx);
491}
492
493//*************************** New versions ************************
494
495GeneralMatrix* AddedMatrix::Evaluate(MatrixType mtd)
496{
497   REPORT
498   Tracer tr("AddedMatrix::Evaluate");
499   gm1=((BaseMatrix*&)bm1)->Evaluate(); gm2=((BaseMatrix*&)bm2)->Evaluate();
500   int nr=gm1->Nrows(); int nc=gm1->Ncols();
501   if (nr!=gm2->Nrows() || nc!=gm2->Ncols())
502   {
503      Try { Throw(IncompatibleDimensionsException(*gm1, *gm2)); }
504      CatchAll
505      {
506         gm1->tDelete(); gm2->tDelete();
507         ReThrow;
508      }
509   }
510   MatrixType mt1 = gm1->type(), mt2 = gm2->type(); MatrixType mts = mt1 + mt2;
511   if (!mtd) { REPORT mtd = mts; }
512   else if (!(mtd.DataLossOK || mtd >= mts))
513   {
514      REPORT
515      gm1->tDelete(); gm2->tDelete();
516      Throw(ProgramException("Illegal Conversion", mts, mtd));
517   }
518   GeneralMatrix* gmx;
519   bool c1 = (mtd == mt1), c2 = (mtd == mt2);
520   if ( c1 && c2 && (gm1->SimpleAddOK(gm2) == 0) )
521   {
522      if (gm1->reuse()) { REPORT AddTo(gm1,gm2); gm2->tDelete(); gmx = gm1; }
523      else if (gm2->reuse()) { REPORT AddTo(gm2,gm1); gmx = gm2; }
524      else
525      {
526         REPORT
527         // what if new throws an exception
528         Try { gmx = mt1.New(nr,nc,this); }
529         CatchAll
530         {
531            ReThrow;
532         }
533         gmx->ReleaseAndDelete(); Add(gmx,gm1,gm2);
534      }
535   }
536   else
537   {
538      if (c1 && c2)
539      {
540         short SAO = gm1->SimpleAddOK(gm2);
541         if (SAO & 1) { REPORT c1 = false; }
542         if (SAO & 2) { REPORT c2 = false; }
543      }
544      if (c1 && gm1->reuse() )               // must have type test first
545         { REPORT AddDS(gm1,gm2); gm2->tDelete(); gmx = gm1; }
546      else if (c2 && gm2->reuse() )
547         { REPORT AddDS(gm2,gm1); if (!c1) gm1->tDelete(); gmx = gm2; }
548      else
549      {
550         REPORT
551         Try { gmx = mtd.New(nr,nc,this); }
552         CatchAll
553         {
554            if (!c1) gm1->tDelete(); if (!c2) gm2->tDelete();
555            ReThrow;
556         }
557         AddDS(gmx,gm1,gm2);
558         if (!c1) gm1->tDelete(); if (!c2) gm2->tDelete();
559         gmx->ReleaseAndDelete();
560      }
561   }
562   return gmx;
563}
564
565GeneralMatrix* SubtractedMatrix::Evaluate(MatrixType mtd)
566{
567   REPORT
568   Tracer tr("SubtractedMatrix::Evaluate");
569   gm1=((BaseMatrix*&)bm1)->Evaluate(); gm2=((BaseMatrix*&)bm2)->Evaluate();
570   int nr=gm1->Nrows(); int nc=gm1->Ncols();
571   if (nr!=gm2->Nrows() || nc!=gm2->Ncols())
572   {
573      Try { Throw(IncompatibleDimensionsException(*gm1, *gm2)); }
574      CatchAll
575      {
576         gm1->tDelete(); gm2->tDelete();
577         ReThrow;
578      }
579   }
580   MatrixType mt1 = gm1->type(), mt2 = gm2->type(); MatrixType mts = mt1 + mt2;
581   if (!mtd) { REPORT mtd = mts; }
582   else if (!(mtd.DataLossOK || mtd >= mts))
583   {
584      gm1->tDelete(); gm2->tDelete();
585      Throw(ProgramException("Illegal Conversion", mts, mtd));
586   }
587   GeneralMatrix* gmx;
588   bool c1 = (mtd == mt1), c2 = (mtd == mt2);
589   if ( c1 && c2 && (gm1->SimpleAddOK(gm2) == 0) )
590   {
591      if (gm1->reuse())
592         { REPORT SubtractFrom(gm1,gm2); gm2->tDelete(); gmx = gm1; }
593      else if (gm2->reuse()) { REPORT ReverseSubtract(gm2,gm1); gmx = gm2; }
594      else
595      {
596         REPORT
597         Try { gmx = mt1.New(nr,nc,this); }
598         CatchAll
599         {
600            ReThrow;
601         }
602         gmx->ReleaseAndDelete(); Subtract(gmx,gm1,gm2);
603      }
604   }
605   else
606   {
607      if (c1 && c2)
608      {
609         short SAO = gm1->SimpleAddOK(gm2);
610         if (SAO & 1) { REPORT c1 = false; }
611         if (SAO & 2) { REPORT c2 = false; }
612      }
613      if (c1 && gm1->reuse() )               // must have type test first
614         { REPORT SubtractDS(gm1,gm2); gm2->tDelete(); gmx = gm1; }
615      else if (c2 && gm2->reuse() )
616      {
617         REPORT ReverseSubtractDS(gm2,gm1);
618         if (!c1) gm1->tDelete(); gmx = gm2;
619      }
620      else
621      {
622         REPORT
623         // what if New throws and exception
624         Try { gmx = mtd.New(nr,nc,this); }
625         CatchAll
626         {
627            if (!c1) gm1->tDelete(); if (!c2) gm2->tDelete();
628            ReThrow;
629         }
630         SubtractDS(gmx,gm1,gm2);
631         if (!c1) gm1->tDelete(); if (!c2) gm2->tDelete();
632         gmx->ReleaseAndDelete();
633      }
634   }
635   return gmx;
636}
637
638GeneralMatrix* SPMatrix::Evaluate(MatrixType mtd)
639{
640   REPORT
641   Tracer tr("SPMatrix::Evaluate");
642   gm1=((BaseMatrix*&)bm1)->Evaluate(); gm2=((BaseMatrix*&)bm2)->Evaluate();
643   int nr=gm1->Nrows(); int nc=gm1->Ncols();
644   if (nr!=gm2->Nrows() || nc!=gm2->Ncols())
645   {
646      Try { Throw(IncompatibleDimensionsException(*gm1, *gm2)); }
647      CatchAll
648      {
649         gm1->tDelete(); gm2->tDelete();
650         ReThrow;
651      }
652   }
653   MatrixType mt1 = gm1->type(), mt2 = gm2->type();
654   MatrixType mts = mt1.SP(mt2);
655   if (!mtd) { REPORT mtd = mts; }
656   else if (!(mtd.DataLossOK || mtd >= mts))
657   {
658      gm1->tDelete(); gm2->tDelete();
659      Throw(ProgramException("Illegal Conversion", mts, mtd));
660   }
661   GeneralMatrix* gmx;
662   bool c1 = (mtd == mt1), c2 = (mtd == mt2);
663   if ( c1 && c2 && (gm1->SimpleAddOK(gm2) == 0) )
664   {
665      if (gm1->reuse()) { REPORT SP(gm1,gm2); gm2->tDelete(); gmx = gm1; }
666      else if (gm2->reuse()) { REPORT SP(gm2,gm1); gmx = gm2; }
667      else
668      {
669         REPORT
670         Try { gmx = mt1.New(nr,nc,this); }
671         CatchAll
672         {
673            ReThrow;
674         }
675         gmx->ReleaseAndDelete(); SP(gmx,gm1,gm2);
676      }
677   }
678   else
679   {
680      if (c1 && c2)
681      {
682         short SAO = gm1->SimpleAddOK(gm2);
683         if (SAO & 1) { REPORT c2 = false; }    // c1 and c2 swapped
684         if (SAO & 2) { REPORT c1 = false; }
685      }
686      if (c1 && gm1->reuse() )               // must have type test first
687         { REPORT SPDS(gm1,gm2); gm2->tDelete(); gmx = gm1; }
688      else if (c2 && gm2->reuse() )
689         { REPORT SPDS(gm2,gm1); if (!c1) gm1->tDelete(); gmx = gm2; }
690      else
691      {
692         REPORT
693         // what if New throws and exception
694         Try { gmx = mtd.New(nr,nc,this); }
695         CatchAll
696         {
697            if (!c1) gm1->tDelete(); if (!c2) gm2->tDelete();
698            ReThrow;
699         }
700         SPDS(gmx,gm1,gm2);
701         if (!c1) gm1->tDelete(); if (!c2) gm2->tDelete();
702         gmx->ReleaseAndDelete();
703      }
704   }
705   return gmx;
706}
707
708
709
710//*************************** norm functions ****************************/
711
712Real BaseMatrix::norm1() const
713{
714   // maximum of sum of absolute values of a column
715   REPORT
716   GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
717   int nc = gm->Ncols(); Real value = 0.0;
718   MatrixCol mc(gm, LoadOnEntry);
719   while (nc--)
720      { Real v = mc.SumAbsoluteValue(); if (value < v) value = v; mc.Next(); }
721   gm->tDelete(); return value;
722}
723
724Real BaseMatrix::norm_infinity() const
725{
726   // maximum of sum of absolute values of a row
727   REPORT
728   GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
729   int nr = gm->Nrows(); Real value = 0.0;
730   MatrixRow mr(gm, LoadOnEntry);
731   while (nr--)
732      { Real v = mr.SumAbsoluteValue(); if (value < v) value = v; mr.Next(); }
733   gm->tDelete(); return value;
734}
735
736//********************** Concatenation and stacking *************************/
737
738GeneralMatrix* ConcatenatedMatrix::Evaluate(MatrixType mtx)
739{
740   REPORT
741   Tracer tr("Concatenate");
742      gm2 = ((BaseMatrix*&)bm2)->Evaluate();
743      gm1 = ((BaseMatrix*&)bm1)->Evaluate();
744      Compare(gm1->type() | gm2->type(),mtx);
745      int nr=gm1->Nrows(); int nc = gm1->Ncols() + gm2->Ncols();
746      if (nr != gm2->Nrows())
747         Throw(IncompatibleDimensionsException(*gm1, *gm2));
748      GeneralMatrix* gmx = mtx.New(nr,nc,this);
749      MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
750      MatrixRow mr(gmx, StoreOnExit+DirectPart);
751      while (nr--) { mr.ConCat(mr1,mr2); mr1.Next(); mr2.Next(); mr.Next(); }
752      gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gmx;
753}
754
755GeneralMatrix* StackedMatrix::Evaluate(MatrixType mtx)
756{
757   REPORT
758   Tracer tr("Stack");
759      gm2 = ((BaseMatrix*&)bm2)->Evaluate();
760      gm1 = ((BaseMatrix*&)bm1)->Evaluate();
761      Compare(gm1->type() & gm2->type(),mtx);
762      int nc=gm1->Ncols();
763      int nr1 = gm1->Nrows(); int nr2 = gm2->Nrows();
764      if (nc != gm2->Ncols())
765         Throw(IncompatibleDimensionsException(*gm1, *gm2));
766      GeneralMatrix* gmx = mtx.New(nr1+nr2,nc,this);
767      MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
768      MatrixRow mr(gmx, StoreOnExit+DirectPart);
769      while (nr1--) { mr.Copy(mr1); mr1.Next(); mr.Next(); }
770      while (nr2--) { mr.Copy(mr2); mr2.Next(); mr.Next(); }
771      gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gmx;
772}
773
774// ************************* equality of matrices ******************** //
775
776static bool RealEqual(Real* s1, Real* s2, int n)
777{
778   int i = n >> 2;
779   while (i--)
780   {
781      if (*s1++ != *s2++) return false; if (*s1++ != *s2++) return false;
782      if (*s1++ != *s2++) return false; if (*s1++ != *s2++) return false;
783   }
784   i = n & 3; while (i--) if (*s1++ != *s2++) return false;
785   return true;
786}
787
788static bool intEqual(int* s1, int* s2, int n)
789{
790   int i = n >> 2;
791   while (i--)
792   {
793      if (*s1++ != *s2++) return false; if (*s1++ != *s2++) return false;
794      if (*s1++ != *s2++) return false; if (*s1++ != *s2++) return false;
795   }
796   i = n & 3; while (i--) if (*s1++ != *s2++) return false;
797   return true;
798}
799
800
801bool operator==(const BaseMatrix& A, const BaseMatrix& B)
802{
803   Tracer tr("BaseMatrix ==");
804   REPORT
805   GeneralMatrix* gmA = ((BaseMatrix&)A).Evaluate();
806   GeneralMatrix* gmB = ((BaseMatrix&)B).Evaluate();
807
808   if (gmA == gmB)                            // same matrix
809      { REPORT gmA->tDelete(); return true; }
810
811   if ( gmA->Nrows() != gmB->Nrows() || gmA->Ncols() != gmB->Ncols() )
812                                              // different dimensions
813      { REPORT gmA->tDelete(); gmB->tDelete(); return false; }
814
815   // check for CroutMatrix or BandLUMatrix
816   MatrixType AType = gmA->type(); MatrixType BType = gmB->type();
817   if (AType.CannotConvert() || BType.CannotConvert() )
818   {
819      REPORT
820      bool bx = gmA->IsEqual(*gmB);
821      gmA->tDelete(); gmB->tDelete();
822      return bx;
823   }
824
825   // is matrix storage the same
826   // will need to modify if further matrix structures are introduced
827   if (AType == BType && gmA->bandwidth() == gmB->bandwidth())
828   {                                          // compare store
829      REPORT
830      bool bx = RealEqual(gmA->Store(),gmB->Store(),gmA->Storage());
831      gmA->tDelete(); gmB->tDelete();
832      return bx;
833   }
834
835   // matrix storage different - just subtract
836   REPORT  return is_zero(*gmA-*gmB);
837}
838
839bool operator==(const GeneralMatrix& A, const GeneralMatrix& B)
840{
841   Tracer tr("GeneralMatrix ==");
842   // May or may not call tDeletes
843   REPORT
844
845   if (&A == &B)                              // same matrix
846      { REPORT return true; }
847
848   if ( A.Nrows() != B.Nrows() || A.Ncols() != B.Ncols() )
849      { REPORT return false; }                // different dimensions
850
851   // check for CroutMatrix or BandLUMatrix
852   MatrixType AType = A.Type(); MatrixType BType = B.Type();
853   if (AType.CannotConvert() || BType.CannotConvert() )
854      { REPORT  return A.IsEqual(B); }
855
856   // is matrix storage the same
857   // will need to modify if further matrix structures are introduced
858   if (AType == BType && A.bandwidth() == B.bandwidth())
859      { REPORT return RealEqual(A.Store(),B.Store(),A.Storage()); }
860
861   // matrix storage different - just subtract
862   REPORT  return is_zero(A-B);
863}
864
865bool GeneralMatrix::is_zero() const
866{
867   REPORT
868   Real* s=store; int i = storage >> 2;
869   while (i--)
870   {
871      if (*s++) return false; if (*s++) return false;
872      if (*s++) return false; if (*s++) return false;
873   }
874   i = storage & 3; while (i--) if (*s++) return false;
875   return true;
876}
877
878bool is_zero(const BaseMatrix& A)
879{
880   Tracer tr("BaseMatrix::is_zero");
881   REPORT
882   GeneralMatrix* gm1 = 0; bool bx;
883   Try { gm1=((BaseMatrix&)A).Evaluate(); bx = gm1->is_zero(); }
884   CatchAll { if (gm1) gm1->tDelete(); ReThrow; }
885   gm1->tDelete();
886   return bx;
887}
888
889// IsEqual functions - insist matrices are of same type
890// as well as equal values to be equal
891
892bool GeneralMatrix::IsEqual(const GeneralMatrix& A) const
893{
894   Tracer tr("GeneralMatrix IsEqual");
895   if (A.type() != type())                       // not same types
896      { REPORT return false; }
897   if (&A == this)                               // same matrix
898      { REPORT  return true; }
899   if (A.nrows_val != nrows_val || A.ncols_val != ncols_val)
900                                                 // different dimensions
901   { REPORT return false; }
902   // is matrix storage the same - compare store
903   REPORT
904   return RealEqual(A.store,store,storage);
905}
906
907bool CroutMatrix::IsEqual(const GeneralMatrix& A) const
908{
909   Tracer tr("CroutMatrix IsEqual");
910   if (A.type() != type())                       // not same types
911      { REPORT return false; }
912   if (&A == this)                               // same matrix
913      { REPORT  return true; }
914   if (A.nrows_val != nrows_val || A.ncols_val != ncols_val)
915                                                 // different dimensions
916   { REPORT return false; }
917   // is matrix storage the same - compare store
918   REPORT
919   return RealEqual(A.store,store,storage)
920      && intEqual(((CroutMatrix&)A).indx, indx, nrows_val);
921}
922
923
924bool BandLUMatrix::IsEqual(const GeneralMatrix& A) const
925{
926   Tracer tr("BandLUMatrix IsEqual");
927   if (A.type() != type())                       // not same types
928      { REPORT  return false; }
929   if (&A == this)                               // same matrix
930      { REPORT  return true; }
931   if ( A.Nrows() != nrows_val || A.Ncols() != ncols_val
932      || ((BandLUMatrix&)A).m1 != m1 || ((BandLUMatrix&)A).m2 != m2 )
933                                                 // different dimensions
934   { REPORT  return false; }
935
936   // matrix storage the same - compare store
937   REPORT
938   return RealEqual(A.Store(),store,storage)
939      && RealEqual(((BandLUMatrix&)A).store2,store2,storage2)
940      && intEqual(((BandLUMatrix&)A).indx, indx, nrows_val);
941}
942
943
944// ************************* cross products ******************** //
945
946inline void crossproduct_body(Real* a, Real* b, Real* c)
947{
948   c[0] = a[1] * b[2] - a[2] * b[1];
949   c[1] = a[2] * b[0] - a[0] * b[2];
950   c[2] = a[0] * b[1] - a[1] * b[0];
951}
952
953Matrix crossproduct(const Matrix& A, const Matrix& B)
954{
955   REPORT
956   int ac = A.Ncols(); int ar = A.Nrows();
957   int bc = B.Ncols(); int br = B.Nrows();
958   Real* a = A.Store(); Real* b = B.Store();
959   if (ac == 3)
960   {
961      if (bc != 3 || ar != 1 || br != 1)
962         { Tracer et("crossproduct"); IncompatibleDimensionsException(A, B); }
963      REPORT
964      RowVector C(3);  Real* c = C.Store(); crossproduct_body(a, b, c);
965      return (Matrix&)C;
966   }
967   else
968   {
969      if (ac != 1 || bc != 1 || ar != 3 || br != 3)
970         { Tracer et("crossproduct"); IncompatibleDimensionsException(A, B); }
971      REPORT
972      ColumnVector C(3);  Real* c = C.Store(); crossproduct_body(a, b, c);
973      return (Matrix&)C;
974   }
975}
976
977ReturnMatrix crossproduct_rows(const Matrix& A, const Matrix& B)
978{
979   REPORT
980   int n = A.Nrows();
981   if (A.Ncols() != 3 || B.Ncols() != 3 || n != B.Nrows())
982   {
983      Tracer et("crossproduct_rows"); IncompatibleDimensionsException(A, B);
984   }
985   Matrix C(n, 3);
986   Real* a = A.Store(); Real* b = B.Store(); Real* c = C.Store();
987   if (n--)
988   {
989      for (;;)
990      {
991         crossproduct_body(a, b, c);
992         if (!(n--)) break;
993         a += 3; b += 3; c += 3;
994      }
995   }
996
997   return C.ForReturn();
998}
999
1000ReturnMatrix crossproduct_columns(const Matrix& A, const Matrix& B)
1001{
1002   REPORT
1003   int n = A.Ncols();
1004   if (A.Nrows() != 3 || B.Nrows() != 3 || n != B.Ncols())
1005   {
1006      Tracer et("crossproduct_columns");
1007      IncompatibleDimensionsException(A, B);
1008   }
1009   Matrix C(3, n);
1010   Real* a = A.Store(); Real* b = B.Store(); Real* c = C.Store();
1011   Real* an = a+n; Real* an2 = an+n;
1012   Real* bn = b+n; Real* bn2 = bn+n;
1013   Real* cn = c+n; Real* cn2 = cn+n;
1014
1015   int i = n;
1016   while (i--)
1017   {
1018      *c++   = *an    * *bn2   - *an2   * *bn;
1019      *cn++  = *an2++ * *b     - *a     * *bn2++;
1020      *cn2++ = *a++   * *bn++  - *an++  * *b++;
1021   }
1022
1023   return C.ForReturn();
1024}
1025
1026
1027#ifdef use_namespace
1028}
1029#endif
1030
1031///@}
1032
Note: See TracBrowser for help on using the repository browser.