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

Last change on this file since 9012 was 8901, checked in by stuerze, 5 years ago

upgrade to newmat11 library
File size: 31.6 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, const GeneralMatrix* gm2)
217{
218 REPORT
219 const 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
225void GeneralMatrix::SP_Equal(const GeneralMatrix& gm)
226{
227 REPORT
228 if (nrows_val != gm.nrows_val || ncols_val != gm.ncols_val)
229 Throw(IncompatibleDimensionsException(*this, gm));
230 SP(this, &gm);
231}
232
233// routines for adding or subtracting matrices of different storage structure
234
235static void AddDS(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
236{
237 REPORT
238 int nr = gm->Nrows();
239 MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
240 MatrixRow mr(gm, StoreOnExit+DirectPart);
241 while (nr--) { mr.Add(mr1,mr2); mr1.Next(); mr2.Next(); mr.Next(); }
242}
243
244static void AddDS(GeneralMatrix* gm, GeneralMatrix* gm2)
245// Add into first argument
246{
247 REPORT
248 int nr = gm->Nrows();
249 MatrixRow mr(gm, StoreOnExit+LoadOnEntry+DirectPart);
250 MatrixRow mr2(gm2, LoadOnEntry);
251 while (nr--) { mr.Add(mr2); mr.Next(); mr2.Next(); }
252}
253
254static void SubtractDS
255 (GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
256{
257 REPORT
258 int nr = gm->Nrows();
259 MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
260 MatrixRow mr(gm, StoreOnExit+DirectPart);
261 while (nr--) { mr.Sub(mr1,mr2); mr1.Next(); mr2.Next(); mr.Next(); }
262}
263
264static void SubtractDS(GeneralMatrix* gm, GeneralMatrix* gm2)
265{
266 REPORT
267 int nr = gm->Nrows();
268 MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart);
269 MatrixRow mr2(gm2, LoadOnEntry);
270 while (nr--) { mr.Sub(mr2); mr.Next(); mr2.Next(); }
271}
272
273static void ReverseSubtractDS(GeneralMatrix* gm, GeneralMatrix* gm2)
274{
275 REPORT
276 int nr = gm->Nrows();
277 MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart);
278 MatrixRow mr2(gm2, LoadOnEntry);
279 while (nr--) { mr.RevSub(mr2); mr2.Next(); mr.Next(); }
280}
281
282static void SPDS(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
283{
284 REPORT
285 int nr = gm->Nrows();
286 MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
287 MatrixRow mr(gm, StoreOnExit+DirectPart);
288 while (nr--) { mr.Multiply(mr1,mr2); mr1.Next(); mr2.Next(); mr.Next(); }
289}
290
291static void SPDS(GeneralMatrix* gm, GeneralMatrix* gm2)
292// SP into first argument
293{
294 REPORT
295 int nr = gm->Nrows();
296 MatrixRow mr(gm, StoreOnExit+LoadOnEntry+DirectPart);
297 MatrixRow mr2(gm2, LoadOnEntry);
298 while (nr--) { mr.Multiply(mr2); mr.Next(); mr2.Next(); }
299}
300
301static GeneralMatrix* GeneralMult1(GeneralMatrix* gm1, GeneralMatrix* gm2,
302 MultipliedMatrix* mm, MatrixType mtx)
303{
304 REPORT
305 Tracer tr("GeneralMult1");
306 int nr=gm1->Nrows(); int nc=gm2->Ncols();
307 if (gm1->Ncols() !=gm2->Nrows())
308 Throw(IncompatibleDimensionsException(*gm1, *gm2));
309 GeneralMatrix* gmx = mtx.New(nr,nc,mm);
310
311 MatrixCol mcx(gmx, StoreOnExit+DirectPart);
312 MatrixCol mc2(gm2, LoadOnEntry);
313 while (nc--)
314 {
315 MatrixRow mr1(gm1, LoadOnEntry, mcx.Skip());
316 Real* el = mcx.Data(); // pointer to an element
317 int n = mcx.Storage();
318 while (n--) { *(el++) = DotProd(mr1,mc2); mr1.Next(); }
319 mc2.Next(); mcx.Next();
320 }
321 gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gmx;
322}
323
324static GeneralMatrix* GeneralMult2(GeneralMatrix* gm1, GeneralMatrix* gm2,
325 MultipliedMatrix* mm, MatrixType mtx)
326{
327 // version that accesses by row only - not good for thin matrices
328 // or column vectors in right hand term.
329 REPORT
330 Tracer tr("GeneralMult2");
331 int nr=gm1->Nrows(); int nc=gm2->Ncols();
332 if (gm1->Ncols() !=gm2->Nrows())
333 Throw(IncompatibleDimensionsException(*gm1, *gm2));
334 GeneralMatrix* gmx = mtx.New(nr,nc,mm);
335
336 MatrixRow mrx(gmx, LoadOnEntry+StoreOnExit+DirectPart);
337 MatrixRow mr1(gm1, LoadOnEntry);
338 while (nr--)
339 {
340 MatrixRow mr2(gm2, LoadOnEntry, mr1.Skip());
341 Real* el = mr1.Data(); // pointer to an element
342 int n = mr1.Storage();
343 mrx.Zero();
344 while (n--) { mrx.AddScaled(mr2, *el++); mr2.Next(); }
345 mr1.Next(); mrx.Next();
346 }
347 gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gmx;
348}
349
350static GeneralMatrix* mmMult(GeneralMatrix* gm1, GeneralMatrix* gm2)
351{
352 // matrix multiplication for type Matrix only
353 REPORT
354 Tracer tr("MatrixMult");
355
356 int nr=gm1->Nrows(); int ncr=gm1->Ncols(); int nc=gm2->Ncols();
357 if (ncr != gm2->Nrows()) Throw(IncompatibleDimensionsException(*gm1,*gm2));
358
359 Matrix* gm = new Matrix(nr,nc); MatrixErrorNoSpace(gm);
360
361 Real* s1=gm1->Store(); Real* s2=gm2->Store(); Real* s=gm->Store();
362
363 if (ncr)
364 {
365 while (nr--)
366 {
367 Real* s2x = s2; int j = ncr;
368 Real* sx = s; Real f = *s1++; int k = nc;
369 while (k--) *sx++ = f * *s2x++;
370 while (--j)
371 { sx = s; f = *s1++; k = nc; while (k--) *sx++ += f * *s2x++; }
372 s = sx;
373 }
374 }
375 else *gm = 0.0;
376
377 gm->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gm;
378}
379
380static GeneralMatrix* GeneralMult(GeneralMatrix* gm1, GeneralMatrix* gm2,
381 MultipliedMatrix* mm, MatrixType mtx)
382{
383 if ( Rectangular(gm1->type(), gm2->type(), mtx))
384 { REPORT return mmMult(gm1, gm2); }
385 Compare(gm1->type() * gm2->type(),mtx);
386 int nr = gm2->Nrows(); int nc = gm2->Ncols();
387 if (nc <= 5 && nr > nc) { REPORT return GeneralMult1(gm1, gm2, mm, mtx); }
388 REPORT return GeneralMult2(gm1, gm2, mm, mtx);
389}
390
391static GeneralMatrix* GeneralKP(GeneralMatrix* gm1, GeneralMatrix* gm2,
392 KPMatrix* kp, MatrixType mtx)
393{
394 REPORT
395 Tracer tr("GeneralKP");
396 int nr1 = gm1->Nrows(); int nc1 = gm1->Ncols();
397 int nr2 = gm2->Nrows(); int nc2 = gm2->Ncols();
398 Compare((gm1->type()).KP(gm2->type()),mtx);
399 GeneralMatrix* gmx = mtx.New(nr1*nr2, nc1*nc2, kp);
400 MatrixRow mrx(gmx, LoadOnEntry+StoreOnExit+DirectPart);
401 MatrixRow mr1(gm1, LoadOnEntry);
402 for (int i = 1; i <= nr1; ++i)
403 {
404 MatrixRow mr2(gm2, LoadOnEntry);
405 for (int j = 1; j <= nr2; ++j)
406 { mrx.KP(mr1,mr2); mr2.Next(); mrx.Next(); }
407 mr1.Next();
408 }
409 gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gmx;
410}
411
412static GeneralMatrix* GeneralSolv(GeneralMatrix* gm1, GeneralMatrix* gm2,
413 BaseMatrix* sm, MatrixType mtx)
414{
415 REPORT
416 Tracer tr("GeneralSolv");
417 Compare(gm1->type().i() * gm2->type(),mtx);
418 int nr = gm1->Nrows();
419 if (nr != gm1->Ncols()) Throw(NotSquareException(*gm1));
420 int nc = gm2->Ncols();
421 if (gm1->Ncols() != gm2->Nrows())
422 Throw(IncompatibleDimensionsException(*gm1, *gm2));
423 GeneralMatrix* gmx = mtx.New(nr,nc,sm); MatrixErrorNoSpace(gmx);
424 Real* r = new Real [nr]; MatrixErrorNoSpace(r);
425 MONITOR_REAL_NEW("Make (GenSolv)",nr,r)
426 GeneralMatrix* gms = gm1->MakeSolver();
427 Try
428 {
429
430 MatrixColX mcx(gmx, r, StoreOnExit+DirectPart); // copy to and from r
431 // this must be inside Try so mcx is destroyed before gmx
432 MatrixColX mc2(gm2, r, LoadOnEntry);
433 while (nc--) { gms->Solver(mcx, mc2); mcx.Next(); mc2.Next(); }
434 }
435 CatchAll
436 {
437 if (gms) gms->tDelete();
438 delete gmx; // <--------------------
439 gm2->tDelete();
440 MONITOR_REAL_DELETE("Delete (GenSolv)",nr,r)
441 // AT&T version 2.1 gives an internal error
442 delete [] r;
443 ReThrow;
444 }
445 gms->tDelete(); gmx->ReleaseAndDelete(); gm2->tDelete();
446 MONITOR_REAL_DELETE("Delete (GenSolv)",nr,r)
447 // AT&T version 2.1 gives an internal error
448 delete [] r;
449 return gmx;
450}
451
452// version for inverses - gm2 is identity
453static GeneralMatrix* GeneralSolvI(GeneralMatrix* gm1, BaseMatrix* sm,
454 MatrixType mtx)
455{
456 REPORT
457 Tracer tr("GeneralSolvI");
458 Compare(gm1->type().i(),mtx);
459 int nr = gm1->Nrows();
460 if (nr != gm1->Ncols()) Throw(NotSquareException(*gm1));
461 int nc = nr;
462 // DiagonalMatrix I(nr); I = 1;
463 IdentityMatrix I(nr);
464 GeneralMatrix* gmx = mtx.New(nr,nc,sm); MatrixErrorNoSpace(gmx);
465 Real* r = new Real [nr]; MatrixErrorNoSpace(r);
466 MONITOR_REAL_NEW("Make (GenSolvI)",nr,r)
467 GeneralMatrix* gms = gm1->MakeSolver();
468 Try
469 {
470
471 MatrixColX mcx(gmx, r, StoreOnExit+DirectPart); // copy to and from r
472 // this must be inside Try so mcx is destroyed before gmx
473 MatrixColX mc2(&I, r, LoadOnEntry);
474 while (nc--) { gms->Solver(mcx, mc2); mcx.Next(); mc2.Next(); }
475 }
476 CatchAll
477 {
478 if (gms) gms->tDelete();
479 delete gmx;
480 MONITOR_REAL_DELETE("Delete (GenSolvI)",nr,r)
481 // AT&T version 2.1 gives an internal error
482 delete [] r;
483 ReThrow;
484 }
485 gms->tDelete(); gmx->ReleaseAndDelete();
486 MONITOR_REAL_DELETE("Delete (GenSolvI)",nr,r)
487 // AT&T version 2.1 gives an internal error
488 delete [] r;
489 return gmx;
490}
491
492GeneralMatrix* InvertedMatrix::Evaluate(MatrixType mtx)
493{
494 // Matrix Inversion - use solve routines
495 Tracer tr("InvertedMatrix::Evaluate");
496 REPORT
497 gm=((BaseMatrix*&)bm)->Evaluate();
498 return GeneralSolvI(gm,this,mtx);
499}
500
501//*************************** New versions ************************
502
503GeneralMatrix* AddedMatrix::Evaluate(MatrixType mtd)
504{
505 REPORT
506 Tracer tr("AddedMatrix::Evaluate");
507 gm1=((BaseMatrix*&)bm1)->Evaluate(); gm2=((BaseMatrix*&)bm2)->Evaluate();
508 int nr=gm1->Nrows(); int nc=gm1->Ncols();
509 if (nr!=gm2->Nrows() || nc!=gm2->Ncols())
510 {
511 Try { Throw(IncompatibleDimensionsException(*gm1, *gm2)); }
512 CatchAll
513 {
514 gm1->tDelete(); gm2->tDelete();
515 ReThrow;
516 }
517 }
518 MatrixType mt1 = gm1->type(), mt2 = gm2->type(); MatrixType mts = mt1 + mt2;
519 if (!mtd) { REPORT mtd = mts; }
520 else if (!(mtd.DataLossOK || mtd >= mts))
521 {
522 REPORT
523 gm1->tDelete(); gm2->tDelete();
524 Throw(ProgramException("Illegal Conversion", mts, mtd));
525 }
526 GeneralMatrix* gmx;
527 bool c1 = (mtd == mt1), c2 = (mtd == mt2);
528 if ( c1 && c2 && (gm1->SimpleAddOK(gm2) == 0) )
529 {
530 if (gm1->reuse()) { REPORT AddTo(gm1,gm2); gm2->tDelete(); gmx = gm1; }
531 else if (gm2->reuse()) { REPORT AddTo(gm2,gm1); gmx = gm2; }
532 else
533 {
534 REPORT
535 // what if new throws an exception
536 Try { gmx = mt1.New(nr,nc,this); }
537 CatchAll
538 {
539 ReThrow;
540 }
541 gmx->ReleaseAndDelete(); Add(gmx,gm1,gm2);
542 }
543 }
544 else
545 {
546 if (c1 && c2)
547 {
548 short SAO = gm1->SimpleAddOK(gm2);
549 if (SAO & 1) { REPORT c1 = false; }
550 if (SAO & 2) { REPORT c2 = false; }
551 }
552 if (c1 && gm1->reuse() ) // must have type test first
553 { REPORT AddDS(gm1,gm2); gm2->tDelete(); gmx = gm1; }
554 else if (c2 && gm2->reuse() )
555 { REPORT AddDS(gm2,gm1); if (!c1) gm1->tDelete(); gmx = gm2; }
556 else
557 {
558 REPORT
559 Try { gmx = mtd.New(nr,nc,this); }
560 CatchAll
561 {
562 if (!c1) gm1->tDelete(); if (!c2) gm2->tDelete();
563 ReThrow;
564 }
565 AddDS(gmx,gm1,gm2);
566 if (!c1) gm1->tDelete(); if (!c2) gm2->tDelete();
567 gmx->ReleaseAndDelete();
568 }
569 }
570 return gmx;
571}
572
573GeneralMatrix* SubtractedMatrix::Evaluate(MatrixType mtd)
574{
575 REPORT
576 Tracer tr("SubtractedMatrix::Evaluate");
577 gm1=((BaseMatrix*&)bm1)->Evaluate(); gm2=((BaseMatrix*&)bm2)->Evaluate();
578 int nr=gm1->Nrows(); int nc=gm1->Ncols();
579 if (nr!=gm2->Nrows() || nc!=gm2->Ncols())
580 {
581 Try { Throw(IncompatibleDimensionsException(*gm1, *gm2)); }
582 CatchAll
583 {
584 gm1->tDelete(); gm2->tDelete();
585 ReThrow;
586 }
587 }
588 MatrixType mt1 = gm1->type(), mt2 = gm2->type(); MatrixType mts = mt1 + mt2;
589 if (!mtd) { REPORT mtd = mts; }
590 else if (!(mtd.DataLossOK || mtd >= mts))
591 {
592 gm1->tDelete(); gm2->tDelete();
593 Throw(ProgramException("Illegal Conversion", mts, mtd));
594 }
595 GeneralMatrix* gmx;
596 bool c1 = (mtd == mt1), c2 = (mtd == mt2);
597 if ( c1 && c2 && (gm1->SimpleAddOK(gm2) == 0) )
598 {
599 if (gm1->reuse())
600 { REPORT SubtractFrom(gm1,gm2); gm2->tDelete(); gmx = gm1; }
601 else if (gm2->reuse()) { REPORT ReverseSubtract(gm2,gm1); gmx = gm2; }
602 else
603 {
604 REPORT
605 Try { gmx = mt1.New(nr,nc,this); }
606 CatchAll
607 {
608 ReThrow;
609 }
610 gmx->ReleaseAndDelete(); Subtract(gmx,gm1,gm2);
611 }
612 }
613 else
614 {
615 if (c1 && c2)
616 {
617 short SAO = gm1->SimpleAddOK(gm2);
618 if (SAO & 1) { REPORT c1 = false; }
619 if (SAO & 2) { REPORT c2 = false; }
620 }
621 if (c1 && gm1->reuse() ) // must have type test first
622 { REPORT SubtractDS(gm1,gm2); gm2->tDelete(); gmx = gm1; }
623 else if (c2 && gm2->reuse() )
624 {
625 REPORT ReverseSubtractDS(gm2,gm1);
626 if (!c1) gm1->tDelete(); gmx = gm2;
627 }
628 else
629 {
630 REPORT
631 // what if New throws and exception
632 Try { gmx = mtd.New(nr,nc,this); }
633 CatchAll
634 {
635 if (!c1) gm1->tDelete(); if (!c2) gm2->tDelete();
636 ReThrow;
637 }
638 SubtractDS(gmx,gm1,gm2);
639 if (!c1) gm1->tDelete(); if (!c2) gm2->tDelete();
640 gmx->ReleaseAndDelete();
641 }
642 }
643 return gmx;
644}
645
646GeneralMatrix* SPMatrix::Evaluate(MatrixType mtd)
647{
648 REPORT
649 Tracer tr("SPMatrix::Evaluate");
650 gm1=((BaseMatrix*&)bm1)->Evaluate(); gm2=((BaseMatrix*&)bm2)->Evaluate();
651 int nr=gm1->Nrows(); int nc=gm1->Ncols();
652 if (nr!=gm2->Nrows() || nc!=gm2->Ncols())
653 {
654 Try { Throw(IncompatibleDimensionsException(*gm1, *gm2)); }
655 CatchAll
656 {
657 gm1->tDelete(); gm2->tDelete();
658 ReThrow;
659 }
660 }
661 MatrixType mt1 = gm1->type(), mt2 = gm2->type();
662 MatrixType mts = mt1.SP(mt2);
663 if (!mtd) { REPORT mtd = mts; }
664 else if (!(mtd.DataLossOK || mtd >= mts))
665 {
666 gm1->tDelete(); gm2->tDelete();
667 Throw(ProgramException("Illegal Conversion", mts, mtd));
668 }
669 GeneralMatrix* gmx;
670 bool c1 = (mtd == mt1), c2 = (mtd == mt2);
671 if ( c1 && c2 && (gm1->SimpleAddOK(gm2) == 0) )
672 {
673 if (gm1->reuse()) { REPORT SP(gm1,gm2); gm2->tDelete(); gmx = gm1; }
674 else if (gm2->reuse()) { REPORT SP(gm2,gm1); gmx = gm2; }
675 else
676 {
677 REPORT
678 Try { gmx = mt1.New(nr,nc,this); }
679 CatchAll
680 {
681 ReThrow;
682 }
683 gmx->ReleaseAndDelete(); SP(gmx,gm1,gm2);
684 }
685 }
686 else
687 {
688 if (c1 && c2)
689 {
690 short SAO = gm1->SimpleAddOK(gm2);
691 if (SAO & 1) { REPORT c2 = false; } // c1 and c2 swapped
692 if (SAO & 2) { REPORT c1 = false; }
693 }
694 if (c1 && gm1->reuse() ) // must have type test first
695 { REPORT SPDS(gm1,gm2); gm2->tDelete(); gmx = gm1; }
696 else if (c2 && gm2->reuse() )
697 { REPORT SPDS(gm2,gm1); if (!c1) gm1->tDelete(); gmx = gm2; }
698 else
699 {
700 REPORT
701 // what if New throws and exception
702 Try { gmx = mtd.New(nr,nc,this); }
703 CatchAll
704 {
705 if (!c1) gm1->tDelete(); if (!c2) gm2->tDelete();
706 ReThrow;
707 }
708 SPDS(gmx,gm1,gm2);
709 if (!c1) gm1->tDelete(); if (!c2) gm2->tDelete();
710 gmx->ReleaseAndDelete();
711 }
712 }
713 return gmx;
714}
715
716
717
718//*************************** norm functions ****************************/
719
720Real BaseMatrix::norm1() const
721{
722 // maximum of sum of absolute values of a column
723 REPORT
724 GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
725 int nc = gm->Ncols(); Real value = 0.0;
726 MatrixCol mc(gm, LoadOnEntry);
727 while (nc--)
728 { Real v = mc.SumAbsoluteValue(); if (value < v) value = v; mc.Next(); }
729 gm->tDelete(); return value;
730}
731
732Real BaseMatrix::norm_infinity() const
733{
734 // maximum of sum of absolute values of a row
735 REPORT
736 GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
737 int nr = gm->Nrows(); Real value = 0.0;
738 MatrixRow mr(gm, LoadOnEntry);
739 while (nr--)
740 { Real v = mr.SumAbsoluteValue(); if (value < v) value = v; mr.Next(); }
741 gm->tDelete(); return value;
742}
743
744//********************** Concatenation and stacking *************************/
745
746GeneralMatrix* ConcatenatedMatrix::Evaluate(MatrixType mtx)
747{
748 REPORT
749 Tracer tr("Concatenate");
750 gm2 = ((BaseMatrix*&)bm2)->Evaluate();
751 gm1 = ((BaseMatrix*&)bm1)->Evaluate();
752 Compare(gm1->type() | gm2->type(),mtx);
753 int nr=gm1->Nrows(); int nc = gm1->Ncols() + gm2->Ncols();
754 if (nr != gm2->Nrows())
755 Throw(IncompatibleDimensionsException(*gm1, *gm2));
756 GeneralMatrix* gmx = mtx.New(nr,nc,this);
757 MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
758 MatrixRow mr(gmx, StoreOnExit+DirectPart);
759 while (nr--) { mr.ConCat(mr1,mr2); mr1.Next(); mr2.Next(); mr.Next(); }
760 gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gmx;
761}
762
763GeneralMatrix* StackedMatrix::Evaluate(MatrixType mtx)
764{
765 REPORT
766 Tracer tr("Stack");
767 gm2 = ((BaseMatrix*&)bm2)->Evaluate();
768 gm1 = ((BaseMatrix*&)bm1)->Evaluate();
769 Compare(gm1->type() & gm2->type(),mtx);
770 int nc=gm1->Ncols();
771 int nr1 = gm1->Nrows(); int nr2 = gm2->Nrows();
772 if (nc != gm2->Ncols())
773 Throw(IncompatibleDimensionsException(*gm1, *gm2));
774 GeneralMatrix* gmx = mtx.New(nr1+nr2,nc,this);
775 MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
776 MatrixRow mr(gmx, StoreOnExit+DirectPart);
777 while (nr1--) { mr.Copy(mr1); mr1.Next(); mr.Next(); }
778 while (nr2--) { mr.Copy(mr2); mr2.Next(); mr.Next(); }
779 gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gmx;
780}
781
782// ************************* equality of matrices ******************** //
783
784static bool RealEqual(Real* s1, Real* s2, int n)
785{
786 int i = n >> 2;
787 while (i--)
788 {
789 if (*s1++ != *s2++) return false; if (*s1++ != *s2++) return false;
790 if (*s1++ != *s2++) return false; if (*s1++ != *s2++) return false;
791 }
792 i = n & 3; while (i--) if (*s1++ != *s2++) return false;
793 return true;
794}
795
796static bool intEqual(int* s1, int* s2, int n)
797{
798 int i = n >> 2;
799 while (i--)
800 {
801 if (*s1++ != *s2++) return false; if (*s1++ != *s2++) return false;
802 if (*s1++ != *s2++) return false; if (*s1++ != *s2++) return false;
803 }
804 i = n & 3; while (i--) if (*s1++ != *s2++) return false;
805 return true;
806}
807
808
809bool operator==(const BaseMatrix& A, const BaseMatrix& B)
810{
811 Tracer tr("BaseMatrix ==");
812 REPORT
813 GeneralMatrix* gmA = ((BaseMatrix&)A).Evaluate();
814 GeneralMatrix* gmB = ((BaseMatrix&)B).Evaluate();
815
816 if (gmA == gmB) // same matrix
817 { REPORT gmA->tDelete(); return true; }
818
819 if ( gmA->Nrows() != gmB->Nrows() || gmA->Ncols() != gmB->Ncols() )
820 // different dimensions
821 { REPORT gmA->tDelete(); gmB->tDelete(); return false; }
822
823 // check for CroutMatrix or BandLUMatrix
824 MatrixType AType = gmA->type(); MatrixType BType = gmB->type();
825 if (AType.CannotConvert() || BType.CannotConvert() )
826 {
827 REPORT
828 bool bx = gmA->IsEqual(*gmB);
829 gmA->tDelete(); gmB->tDelete();
830 return bx;
831 }
832
833 // is matrix storage the same
834 // will need to modify if further matrix structures are introduced
835 if (AType == BType && gmA->bandwidth() == gmB->bandwidth())
836 { // compare store
837 REPORT
838 bool bx = RealEqual(gmA->Store(),gmB->Store(),gmA->Storage());
839 gmA->tDelete(); gmB->tDelete();
840 return bx;
841 }
842
843 // matrix storage different - just subtract
844 REPORT return is_zero(*gmA-*gmB);
845}
846
847bool operator==(const GeneralMatrix& A, const GeneralMatrix& B)
848{
849 Tracer tr("GeneralMatrix ==");
850 // May or may not call tDeletes
851 REPORT
852
853 if (&A == &B) // same matrix
854 { REPORT return true; }
855
856 if ( A.Nrows() != B.Nrows() || A.Ncols() != B.Ncols() )
857 { REPORT return false; } // different dimensions
858
859 // check for CroutMatrix or BandLUMatrix
860 MatrixType AType = A.Type(); MatrixType BType = B.Type();
861 if (AType.CannotConvert() || BType.CannotConvert() )
862 { REPORT return A.IsEqual(B); }
863
864 // is matrix storage the same
865 // will need to modify if further matrix structures are introduced
866 if (AType == BType && A.bandwidth() == B.bandwidth())
867 { REPORT return RealEqual(A.Store(),B.Store(),A.Storage()); }
868
869 // matrix storage different - just subtract
870 REPORT return is_zero(A-B);
871}
872
873bool GeneralMatrix::is_zero() const
874{
875 REPORT
876 Real* s=store; int i = storage >> 2;
877 while (i--)
878 {
879 if (*s++) return false; if (*s++) return false;
880 if (*s++) return false; if (*s++) return false;
881 }
882 i = storage & 3; while (i--) if (*s++) return false;
883 return true;
884}
885
886bool is_zero(const BaseMatrix& A)
887{
888 Tracer tr("BaseMatrix::is_zero");
889 REPORT
890 GeneralMatrix* gm1 = 0; bool bx;
891 Try { gm1=((BaseMatrix&)A).Evaluate(); bx = gm1->is_zero(); }
892 CatchAll { if (gm1) gm1->tDelete(); ReThrow; }
893 gm1->tDelete();
894 return bx;
895}
896
897// IsEqual functions - insist matrices are of same type
898// as well as equal values to be equal
899
900bool GeneralMatrix::IsEqual(const GeneralMatrix& A) const
901{
902 Tracer tr("GeneralMatrix IsEqual");
903 if (A.type() != type()) // not same types
904 { REPORT return false; }
905 if (&A == this) // same matrix
906 { REPORT return true; }
907 if (A.nrows_val != nrows_val || A.ncols_val != ncols_val)
908 // different dimensions
909 { REPORT return false; }
910 // is matrix storage the same - compare store
911 REPORT
912 return RealEqual(A.store,store,storage);
913}
914
915bool CroutMatrix::IsEqual(const GeneralMatrix& A) const
916{
917 Tracer tr("CroutMatrix IsEqual");
918 if (A.type() != type()) // not same types
919 { REPORT return false; }
920 if (&A == this) // same matrix
921 { REPORT return true; }
922 if (A.nrows_val != nrows_val || A.ncols_val != ncols_val)
923 // different dimensions
924 { REPORT return false; }
925 // is matrix storage the same - compare store
926 REPORT
927 return RealEqual(A.store,store,storage)
928 && intEqual(((CroutMatrix&)A).indx, indx, nrows_val);
929}
930
931
932bool BandLUMatrix::IsEqual(const GeneralMatrix& A) const
933{
934 Tracer tr("BandLUMatrix IsEqual");
935 if (A.type() != type()) // not same types
936 { REPORT return false; }
937 if (&A == this) // same matrix
938 { REPORT return true; }
939 if ( A.Nrows() != nrows_val || A.Ncols() != ncols_val
940 || ((BandLUMatrix&)A).m1 != m1 || ((BandLUMatrix&)A).m2 != m2 )
941 // different dimensions
942 { REPORT return false; }
943
944 // matrix storage the same - compare store
945 REPORT
946 return RealEqual(A.Store(),store,storage)
947 && RealEqual(((BandLUMatrix&)A).store2,store2,storage2)
948 && intEqual(((BandLUMatrix&)A).indx, indx, nrows_val);
949}
950
951
952// ************************* cross products ******************** //
953
954inline void crossproduct_body(Real* a, Real* b, Real* c)
955{
956 c[0] = a[1] * b[2] - a[2] * b[1];
957 c[1] = a[2] * b[0] - a[0] * b[2];
958 c[2] = a[0] * b[1] - a[1] * b[0];
959}
960
961Matrix crossproduct(const Matrix& A, const Matrix& B)
962{
963 REPORT
964 int ac = A.Ncols(); int ar = A.Nrows();
965 int bc = B.Ncols(); int br = B.Nrows();
966 Real* a = A.Store(); Real* b = B.Store();
967 if (ac == 3)
968 {
969 if (bc != 3 || ar != 1 || br != 1)
970 { Tracer et("crossproduct"); IncompatibleDimensionsException(A, B); }
971 REPORT
972 RowVector C(3); Real* c = C.Store(); crossproduct_body(a, b, c);
973 return (Matrix&)C;
974 }
975 else
976 {
977 if (ac != 1 || bc != 1 || ar != 3 || br != 3)
978 { Tracer et("crossproduct"); IncompatibleDimensionsException(A, B); }
979 REPORT
980 ColumnVector C(3); Real* c = C.Store(); crossproduct_body(a, b, c);
981 return (Matrix&)C;
982 }
983}
984
985ReturnMatrix crossproduct_rows(const Matrix& A, const Matrix& B)
986{
987 REPORT
988 int n = A.Nrows();
989 if (A.Ncols() != 3 || B.Ncols() != 3 || n != B.Nrows())
990 {
991 Tracer et("crossproduct_rows"); IncompatibleDimensionsException(A, B);
992 }
993 Matrix C(n, 3);
994 Real* a = A.Store(); Real* b = B.Store(); Real* c = C.Store();
995 if (n--)
996 {
997 for (;;)
998 {
999 crossproduct_body(a, b, c);
1000 if (!(n--)) break;
1001 a += 3; b += 3; c += 3;
1002 }
1003 }
1004
1005 return C.ForReturn();
1006}
1007
1008ReturnMatrix crossproduct_columns(const Matrix& A, const Matrix& B)
1009{
1010 REPORT
1011 int n = A.Ncols();
1012 if (A.Nrows() != 3 || B.Nrows() != 3 || n != B.Ncols())
1013 {
1014 Tracer et("crossproduct_columns");
1015 IncompatibleDimensionsException(A, B);
1016 }
1017 Matrix C(3, n);
1018 Real* a = A.Store(); Real* b = B.Store(); Real* c = C.Store();
1019 Real* an = a+n; Real* an2 = an+n;
1020 Real* bn = b+n; Real* bn2 = bn+n;
1021 Real* cn = c+n; Real* cn2 = cn+n;
1022
1023 int i = n;
1024 while (i--)
1025 {
1026 *c++ = *an * *bn2 - *an2 * *bn;
1027 *cn++ = *an2++ * *b - *a * *bn2++;
1028 *cn2++ = *a++ * *bn++ - *an++ * *b++;
1029 }
1030
1031 return C.ForReturn();
1032}
1033
1034
1035#ifdef use_namespace
1036}
1037#endif
1038
1039///@}
1040
Note: See TracBrowser for help on using the repository browser.