source: ntrip/trunk/BNC/newmat/newmat7.cpp@ 10381

Last change on this file since 10381 was 9434, checked in by stuerze, 3 years ago

newmat is reverted to its stable version 10 because version 11beta was not working with mac compilers
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.