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 |
---|
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 | |
---|
28 | GeneralMatrix* GeneralMatrix::MakeSolver() |
---|
29 | { |
---|
30 | REPORT |
---|
31 | GeneralMatrix* gm = new CroutMatrix(*this); |
---|
32 | MatrixErrorNoSpace(gm); gm->ReleaseAndDelete(); return gm; |
---|
33 | } |
---|
34 | |
---|
35 | GeneralMatrix* Matrix::MakeSolver() |
---|
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 | |
---|
55 | void 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 | |
---|
75 | void 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 | |
---|
95 | static GeneralMatrix* |
---|
96 | GeneralMult(GeneralMatrix*,GeneralMatrix*,MultipliedMatrix*,MatrixType); |
---|
97 | static GeneralMatrix* |
---|
98 | GeneralSolv(GeneralMatrix*,GeneralMatrix*,BaseMatrix*,MatrixType); |
---|
99 | static GeneralMatrix* |
---|
100 | GeneralSolvI(GeneralMatrix*,BaseMatrix*,MatrixType); |
---|
101 | static GeneralMatrix* |
---|
102 | GeneralKP(GeneralMatrix*,GeneralMatrix*,KPMatrix*,MatrixType); |
---|
103 | |
---|
104 | GeneralMatrix* 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 | |
---|
113 | GeneralMatrix* 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 | |
---|
121 | GeneralMatrix* 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 | |
---|
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 | |
---|
153 | void 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 | |
---|
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 | |
---|
183 | void 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 | |
---|
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); |
---|
232 | MatrixRow mr(gm, StoreOnExit+DirectPart); |
---|
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(); |
---|
241 | MatrixRow mr(gm, StoreOnExit+LoadOnEntry+DirectPart); |
---|
242 | MatrixRow mr2(gm2, LoadOnEntry); |
---|
243 | while (nr--) { mr.Add(mr2); mr.Next(); mr2.Next(); } |
---|
244 | } |
---|
245 | |
---|
246 | static 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 | |
---|
256 | static 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 | |
---|
265 | static 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 | |
---|
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); |
---|
279 | MatrixRow mr(gm, StoreOnExit+DirectPart); |
---|
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(); |
---|
288 | MatrixRow mr(gm, StoreOnExit+LoadOnEntry+DirectPart); |
---|
289 | MatrixRow mr2(gm2, LoadOnEntry); |
---|
290 | while (nr--) { mr.Multiply(mr2); mr.Next(); mr2.Next(); } |
---|
291 | } |
---|
292 | |
---|
293 | static 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 | |
---|
316 | static 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 | |
---|
342 | static 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 | |
---|
372 | static 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 | |
---|
383 | static 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 | |
---|
404 | static 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 |
---|
445 | static 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 | |
---|
484 | GeneralMatrix* 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 | |
---|
495 | GeneralMatrix* 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 | |
---|
565 | GeneralMatrix* 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 | |
---|
638 | GeneralMatrix* 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 | |
---|
712 | Real 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 | |
---|
724 | Real 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 | |
---|
738 | GeneralMatrix* 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 | |
---|
755 | GeneralMatrix* 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 | |
---|
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 | |
---|
865 | bool 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 | |
---|
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 | |
---|
892 | bool 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 | |
---|
907 | bool 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 | |
---|
924 | bool 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 | |
---|
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 | |
---|
977 | ReturnMatrix 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 | |
---|
1000 | ReturnMatrix 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 | |
---|