newmat7.cpp
Go to the documentation of this file.
1 
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
15 namespace 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 
29 {
30  REPORT
31  GeneralMatrix* gm = new CroutMatrix(*this);
32  MatrixErrorNoSpace(gm); gm->ReleaseAndDelete(); return gm;
33 }
34 
36 {
37  REPORT
38  GeneralMatrix* gm = new CroutMatrix(*this);
39  MatrixErrorNoSpace(gm); gm->ReleaseAndDelete(); return gm;
40 }
41 
42 void 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 
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 
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 
95 static GeneralMatrix*
97 static GeneralMatrix*
99 static GeneralMatrix*
101 static GeneralMatrix*
103 
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 
114 {
115  REPORT
116  gm1 = ((BaseMatrix*&)bm1)->Evaluate();
117  gm2 = ((BaseMatrix*&)bm2)->Evaluate();
118  return GeneralSolv(gm1,gm2,this,mt);
119 }
120 
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 
131 static 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 
144 static 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 
154 {
155  REPORT
156  if (nrows_val != gm.nrows_val || ncols_val != gm.ncols_val)
158  AddTo(this, &gm);
159 }
160 
161 static 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 
174 static 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 
184 {
185  REPORT
186  if (nrows_val != gm.nrows_val || ncols_val != gm.ncols_val)
188  SubtractFrom(this, &gm);
189 }
190 
191 static 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 
203 static 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 
216 static 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 
227 static 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);
233  while (nr--) { mr.Add(mr1,mr2); mr1.Next(); mr2.Next(); mr.Next(); }
234 }
235 
236 static void AddDS(GeneralMatrix* gm, GeneralMatrix* gm2)
237 // Add into first argument
238 {
239  REPORT
240  int nr = gm->Nrows();
242  MatrixRow mr2(gm2, LoadOnEntry);
243  while (nr--) { mr.Add(mr2); mr.Next(); mr2.Next(); }
244 }
245 
246 static void SubtractDS
248 {
249  REPORT
250  int nr = gm->Nrows();
251  MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
253  while (nr--) { mr.Sub(mr1,mr2); mr1.Next(); mr2.Next(); mr.Next(); }
254 }
255 
256 static void SubtractDS(GeneralMatrix* gm, GeneralMatrix* gm2)
257 {
258  REPORT
259  int nr = gm->Nrows();
261  MatrixRow mr2(gm2, LoadOnEntry);
262  while (nr--) { mr.Sub(mr2); mr.Next(); mr2.Next(); }
263 }
264 
266 {
267  REPORT
268  int nr = gm->Nrows();
270  MatrixRow mr2(gm2, LoadOnEntry);
271  while (nr--) { mr.RevSub(mr2); mr2.Next(); mr.Next(); }
272 }
273 
274 static 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);
280  while (nr--) { mr.Multiply(mr1,mr2); mr1.Next(); mr2.Next(); mr.Next(); }
281 }
282 
283 static void SPDS(GeneralMatrix* gm, GeneralMatrix* gm2)
284 // SP into first argument
285 {
286  REPORT
287  int nr = gm->Nrows();
289  MatrixRow mr2(gm2, LoadOnEntry);
290  while (nr--) { mr.Multiply(mr2); mr.Next(); mr2.Next(); }
291 }
292 
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())
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 
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())
326  GeneralMatrix* gmx = mtx.New(nr,nc,mm);
327 
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 
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 
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 
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);
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 
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())
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
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 
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 
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 
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 
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 
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 
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 
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())
748  GeneralMatrix* gmx = mtx.New(nr,nc,this);
749  MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
751  while (nr--) { mr.ConCat(mr1,mr2); mr1.Next(); mr2.Next(); mr.Next(); }
752  gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gmx;
753 }
754 
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())
766  GeneralMatrix* gmx = mtx.New(nr1+nr2,nc,this);
767  MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
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 
776 static 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 
788 static 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 
801 bool 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 
839 bool 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 
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 
878 bool 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 
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 
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 
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 
946 inline 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 
953 Matrix 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 
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 
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");
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 
1032 
#define ReThrow
Definition: myexcept.h:192
void Solver(MatrixColX &, const MatrixColX &)
Definition: newmat7.cpp:42
MatrixType Type() const
Definition: newmat.h:493
virtual short SimpleAddOK(const GeneralMatrix *)
Definition: newmat.h:480
Real DotProd(const MatrixRowCol &mrc1, const MatrixRowCol &mrc2)
Definition: newmat2.cpp:80
GeneralMatrix * Evaluate(MatrixType mt=MatrixTypeUnSp)
Definition: newmat7.cpp:104
void Copy(const MatrixRowCol &)
Definition: newmat2.cpp:419
int Storage()
Definition: newmatrc.h:101
void tDelete()
Definition: newmat4.cpp:786
#define Try
Definition: myexcept.h:190
GeneralMatrix * New() const
new matrix of given type
ReturnMatrix crossproduct_rows(const Matrix &A, const Matrix &B)
Definition: newmat7.cpp:977
Miscellaneous exception (details in character string).
Definition: newmat.h:1947
MatrixType i() const
type of inverse
Definition: newmat1.cpp:82
static GeneralMatrix * GeneralMult2(GeneralMatrix *gm1, GeneralMatrix *gm2, MultipliedMatrix *mm, MatrixType mtx)
Definition: newmat7.cpp:316
void KP(const MatrixRowCol &, const MatrixRowCol &)
Definition: newmat2.cpp:352
bool IsEqual(const GeneralMatrix &) const
Definition: newmat7.cpp:924
static void SPDS(GeneralMatrix *gm, GeneralMatrix *gm1, GeneralMatrix *gm2)
Definition: newmat7.cpp:274
static bool intEqual(int *s1, int *s2, int n)
Definition: newmat7.cpp:788
Matrix crossproduct(const Matrix &A, const Matrix &B)
Definition: newmat7.cpp:953
GeneralMatrix * Evaluate(MatrixType mt=MatrixTypeUnSp)
Definition: newmat7.cpp:638
void RevSub(const MatrixRowCol &)
Definition: newmat2.cpp:271
void Add(const MatrixRowCol &)
Definition: newmat2.cpp:36
GeneralMatrix * MakeSolver()
Definition: newmat7.cpp:35
void ConCat(const MatrixRowCol &, const MatrixRowCol &)
Definition: newmat2.cpp:287
static GeneralMatrix * mmMult(GeneralMatrix *gm1, GeneralMatrix *gm2)
Definition: newmat7.cpp:342
static void ReverseSubtractDS(GeneralMatrix *gm, GeneralMatrix *gm2)
Definition: newmat7.cpp:265
void AddScaled(const MatrixRowCol &, Real)
Definition: newmat2.cpp:47
static GeneralMatrix * GeneralSolvI(GeneralMatrix *, BaseMatrix *, MatrixType)
Definition: newmat7.cpp:445
double Real
Definition: include.h:307
virtual MatrixType type() const =0
virtual MatrixBandWidth bandwidth() const
Definition: newmat4.cpp:671
GeneralMatrix * Evaluate(MatrixType mt=MatrixTypeUnSp)
Definition: newmat5.cpp:85
static void Add(GeneralMatrix *gm, GeneralMatrix *gm1, GeneralMatrix *gm2)
Definition: newmat7.cpp:131
Real norm1() const
Definition: newmat7.cpp:712
int Skip()
Definition: newmatrc.h:100
bool is_zero() const
Definition: newmat7.cpp:865
static void ReverseSubtract(GeneralMatrix *gm, const GeneralMatrix *gm2)
Definition: newmat7.cpp:191
#define REPORT
Definition: newmat7.cpp:22
int Nrows() const
Definition: newmat.h:494
KPMatrix KP(const BaseMatrix &, const BaseMatrix &)
Definition: newmat6.cpp:281
void PlusEqual(const GeneralMatrix &gm)
Definition: newmat7.cpp:153
bool CannotConvert() const
Definition: newmat.h:190
#define CatchAll
Definition: myexcept.h:194
int Storage() const
Definition: newmat.h:496
int nrows_val
Definition: newmat.h:452
A matrix is not square exception.
Definition: newmat.h:1981
GeneralMatrix * Evaluate(MatrixType mt=MatrixTypeUnSp)
Definition: newmat7.cpp:495
bool is_zero(const BaseMatrix &A)
Definition: newmat7.cpp:878
void crossproduct_body(Real *a, Real *b, Real *c)
Definition: newmat7.cpp:946
void ReleaseAndDelete()
Definition: newmat.h:516
GeneralMatrix * Evaluate(MatrixType mt=MatrixTypeUnSp)
Definition: newmat7.cpp:484
static void SubtractFrom(GeneralMatrix *gm, const GeneralMatrix *gm2)
Definition: newmat7.cpp:174
void Next()
Definition: newmatrc.h:177
bool reuse()
Definition: newmat4.cpp:819
Real * data
Definition: newmatrc.h:51
int storage
Definition: newmatrc.h:48
MatrixType SP(const MatrixType &) const
Definition: newmat1.cpp:45
Real * Data()
Definition: newmatrc.h:99
void Multiply(const MatrixRowCol &)
Definition: newmat2.cpp:314
Real * Store() const
Definition: newmat.h:497
static GeneralMatrix * GeneralMult(GeneralMatrix *, GeneralMatrix *, MultipliedMatrix *, MatrixType)
Definition: newmat7.cpp:372
bool Compare(const MatrixType &, MatrixType &)
Definition: newmat4.cpp:980
void Sub(const MatrixRowCol &)
Definition: newmat2.cpp:58
#define MONITOR_REAL_DELETE(Operation, Size, Pointer)
Definition: myexcept.h:331
bool IsEqual(const GeneralMatrix &) const
Definition: newmat7.cpp:907
static bool RealEqual(Real *s1, Real *s2, int n)
Definition: newmat7.cpp:776
ReturnMatrix crossproduct_columns(const Matrix &A, const Matrix &B)
Definition: newmat7.cpp:1000
static void SP(GeneralMatrix *gm, GeneralMatrix *gm1, GeneralMatrix *gm2)
Definition: newmat7.cpp:203
Base of the matrix classes.
Definition: newmat.h:292
bool DataLossOK
Definition: newmat.h:140
virtual bool IsEqual(const GeneralMatrix &) const
Definition: newmat7.cpp:892
#define Throw(E)
Definition: myexcept.h:191
Real norm_infinity() const
Definition: newmat7.cpp:724
The usual rectangular matrix.
Definition: newmat.h:625
static void SubtractDS(GeneralMatrix *gm, GeneralMatrix *gm1, GeneralMatrix *gm2)
Definition: newmat7.cpp:247
Real * store
Definition: newmat.h:454
#define MONITOR_REAL_NEW(Operation, Size, Pointer)
Definition: myexcept.h:329
GeneralMatrix * Evaluate(MatrixType mt=MatrixTypeUnSp)
Definition: newmat7.cpp:113
static GeneralMatrix * GeneralSolv(GeneralMatrix *, GeneralMatrix *, BaseMatrix *, MatrixType)
Definition: newmat7.cpp:404
int ncols_val
Definition: newmat.h:452
void MinusEqual(const GeneralMatrix &gm)
Definition: newmat7.cpp:183
void Solver(MatrixColX &, const MatrixColX &)
Definition: newmat7.cpp:75
static void Subtract(GeneralMatrix *gm, GeneralMatrix *gm1, GeneralMatrix *gm2)
Definition: newmat7.cpp:161
Incompatible dimensions exception.
Definition: newmat.h:1998
FloatVector FloatVector * a
bool Rectangular(MatrixType a, MatrixType b, MatrixType c)
Definition: newmat1.cpp:108
void Next()
Definition: newmatrc.h:175
virtual GeneralMatrix * MakeSolver()
Definition: newmat7.cpp:28
int Ncols() const
Definition: newmat.h:495
LU decomposition of a band matrix.
Definition: newmat.h:1306
void MatrixErrorNoSpace(const void *)
test for allocation fails
Definition: newmatex.cpp:301
Real SumAbsoluteValue()
Definition: newmat2.cpp:585
static GeneralMatrix * GeneralMult1(GeneralMatrix *gm1, GeneralMatrix *gm2, MultipliedMatrix *mm, MatrixType mtx)
Definition: newmat7.cpp:293
ReturnMatrix ForReturn() const
Definition: newmat.h:2157
GeneralMatrix * Evaluate(MatrixType mt=MatrixTypeUnSp)
Definition: newmat7.cpp:738
Real sum(const BaseMatrix &B)
Definition: newmat.h:2105
static GeneralMatrix * GeneralKP(GeneralMatrix *, GeneralMatrix *, KPMatrix *, MatrixType)
Definition: newmat7.cpp:383
bool operator==(const BaseMatrix &A, const BaseMatrix &B)
Definition: newmat7.cpp:801
GeneralMatrix * Evaluate(MatrixType mt=MatrixTypeUnSp)
Definition: newmat7.cpp:565
Row vector.
Definition: newmat.h:953
virtual void Solver(MatrixColX &, const MatrixColX &)
Definition: newmat.h:537
Column vector.
Definition: newmat.h:1008
The classes for matrices that can contain data are derived from this.
Definition: newmat.h:447
static void AddDS(GeneralMatrix *gm, GeneralMatrix *gm1, GeneralMatrix *gm2)
Definition: newmat7.cpp:227
void Solver(MatrixColX &, const MatrixColX &)
Definition: newmat7.cpp:55
void Zero()
Definition: newmat2.cpp:566
void Next()
Definition: newmatrc.h:161
GeneralMatrix * Evaluate(MatrixType mt=MatrixTypeUnSp)
Definition: newmat7.cpp:755
static void AddTo(GeneralMatrix *gm, const GeneralMatrix *gm2)
Definition: newmat7.cpp:144
GeneralMatrix * Evaluate(MatrixType mt=MatrixTypeUnSp)
Definition: newmat7.cpp:121
Identity matrix.
Definition: newmat.h:1350


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