1 | /// \ingroup newmat |
---|
2 | ///@{ |
---|
3 | |
---|
4 | /// \file newmat6.cpp |
---|
5 | /// Operators, element access. |
---|
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 | |
---|
20 | #ifdef DO_REPORT |
---|
21 | #define REPORT { static ExeCounter ExeCount(__LINE__,6); ++ExeCount; } |
---|
22 | #else |
---|
23 | #define REPORT {} |
---|
24 | #endif |
---|
25 | |
---|
26 | /*************************** general utilities *************************/ |
---|
27 | |
---|
28 | static int tristore(int n) // els in triangular matrix |
---|
29 | { return (n*(n+1))/2; } |
---|
30 | |
---|
31 | |
---|
32 | /****************************** operators *******************************/ |
---|
33 | |
---|
34 | Real& Matrix::operator()(int m, int n) |
---|
35 | { |
---|
36 | REPORT |
---|
37 | if (m<=0 || m>nrows_val || n<=0 || n>ncols_val) |
---|
38 | Throw(IndexException(m,n,*this)); |
---|
39 | return store[(m-1)*ncols_val+n-1]; |
---|
40 | } |
---|
41 | |
---|
42 | Real& SymmetricMatrix::operator()(int m, int n) |
---|
43 | { |
---|
44 | REPORT |
---|
45 | if (m<=0 || n<=0 || m>nrows_val || n>ncols_val) |
---|
46 | Throw(IndexException(m,n,*this)); |
---|
47 | if (m>=n) return store[tristore(m-1)+n-1]; |
---|
48 | else return store[tristore(n-1)+m-1]; |
---|
49 | } |
---|
50 | |
---|
51 | Real& UpperTriangularMatrix::operator()(int m, int n) |
---|
52 | { |
---|
53 | REPORT |
---|
54 | if (m<=0 || n<m || n>ncols_val) |
---|
55 | Throw(IndexException(m,n,*this)); |
---|
56 | return store[(m-1)*ncols_val+n-1-tristore(m-1)]; |
---|
57 | } |
---|
58 | |
---|
59 | Real& LowerTriangularMatrix::operator()(int m, int n) |
---|
60 | { |
---|
61 | REPORT |
---|
62 | if (n<=0 || m<n || m>nrows_val) |
---|
63 | Throw(IndexException(m,n,*this)); |
---|
64 | return store[tristore(m-1)+n-1]; |
---|
65 | } |
---|
66 | |
---|
67 | Real& DiagonalMatrix::operator()(int m, int n) |
---|
68 | { |
---|
69 | REPORT |
---|
70 | if (n<=0 || m!=n || m>nrows_val || n>ncols_val) |
---|
71 | Throw(IndexException(m,n,*this)); |
---|
72 | return store[n-1]; |
---|
73 | } |
---|
74 | |
---|
75 | Real& DiagonalMatrix::operator()(int m) |
---|
76 | { |
---|
77 | REPORT |
---|
78 | if (m<=0 || m>nrows_val) Throw(IndexException(m,*this)); |
---|
79 | return store[m-1]; |
---|
80 | } |
---|
81 | |
---|
82 | Real& ColumnVector::operator()(int m) |
---|
83 | { |
---|
84 | REPORT |
---|
85 | if (m<=0 || m> nrows_val) Throw(IndexException(m,*this)); |
---|
86 | return store[m-1]; |
---|
87 | } |
---|
88 | |
---|
89 | Real& RowVector::operator()(int n) |
---|
90 | { |
---|
91 | REPORT |
---|
92 | if (n<=0 || n> ncols_val) Throw(IndexException(n,*this)); |
---|
93 | return store[n-1]; |
---|
94 | } |
---|
95 | |
---|
96 | Real& BandMatrix::operator()(int m, int n) |
---|
97 | { |
---|
98 | REPORT |
---|
99 | int w = upper_val+lower_val+1; int i = lower_val+n-m; |
---|
100 | if (m<=0 || m>nrows_val || n<=0 || n>ncols_val || i<0 || i>=w) |
---|
101 | Throw(IndexException(m,n,*this)); |
---|
102 | return store[w*(m-1)+i]; |
---|
103 | } |
---|
104 | |
---|
105 | Real& UpperBandMatrix::operator()(int m, int n) |
---|
106 | { |
---|
107 | REPORT |
---|
108 | int w = upper_val+1; int i = n-m; |
---|
109 | if (m<=0 || m>nrows_val || n<=0 || n>ncols_val || i<0 || i>=w) |
---|
110 | Throw(IndexException(m,n,*this)); |
---|
111 | return store[w*(m-1)+i]; |
---|
112 | } |
---|
113 | |
---|
114 | Real& LowerBandMatrix::operator()(int m, int n) |
---|
115 | { |
---|
116 | REPORT |
---|
117 | int w = lower_val+1; int i = lower_val+n-m; |
---|
118 | if (m<=0 || m>nrows_val || n<=0 || n>ncols_val || i<0 || i>=w) |
---|
119 | Throw(IndexException(m,n,*this)); |
---|
120 | return store[w*(m-1)+i]; |
---|
121 | } |
---|
122 | |
---|
123 | Real& SymmetricBandMatrix::operator()(int m, int n) |
---|
124 | { |
---|
125 | REPORT |
---|
126 | int w = lower_val+1; |
---|
127 | if (m>=n) |
---|
128 | { |
---|
129 | REPORT |
---|
130 | int i = lower_val+n-m; |
---|
131 | if ( m>nrows_val || n<=0 || i<0 ) |
---|
132 | Throw(IndexException(m,n,*this)); |
---|
133 | return store[w*(m-1)+i]; |
---|
134 | } |
---|
135 | else |
---|
136 | { |
---|
137 | REPORT |
---|
138 | int i = lower_val+m-n; |
---|
139 | if ( n>nrows_val || m<=0 || i<0 ) |
---|
140 | Throw(IndexException(m,n,*this)); |
---|
141 | return store[w*(n-1)+i]; |
---|
142 | } |
---|
143 | } |
---|
144 | |
---|
145 | |
---|
146 | Real Matrix::operator()(int m, int n) const |
---|
147 | { |
---|
148 | REPORT |
---|
149 | if (m<=0 || m>nrows_val || n<=0 || n>ncols_val) |
---|
150 | Throw(IndexException(m,n,*this)); |
---|
151 | return store[(m-1)*ncols_val+n-1]; |
---|
152 | } |
---|
153 | |
---|
154 | Real SymmetricMatrix::operator()(int m, int n) const |
---|
155 | { |
---|
156 | REPORT |
---|
157 | if (m<=0 || n<=0 || m>nrows_val || n>ncols_val) |
---|
158 | Throw(IndexException(m,n,*this)); |
---|
159 | if (m>=n) return store[tristore(m-1)+n-1]; |
---|
160 | else return store[tristore(n-1)+m-1]; |
---|
161 | } |
---|
162 | |
---|
163 | Real UpperTriangularMatrix::operator()(int m, int n) const |
---|
164 | { |
---|
165 | REPORT |
---|
166 | if (m<=0 || n<m || n>ncols_val) |
---|
167 | Throw(IndexException(m,n,*this)); |
---|
168 | return store[(m-1)*ncols_val+n-1-tristore(m-1)]; |
---|
169 | } |
---|
170 | |
---|
171 | Real LowerTriangularMatrix::operator()(int m, int n) const |
---|
172 | { |
---|
173 | REPORT |
---|
174 | if (n<=0 || m<n || m>nrows_val) |
---|
175 | Throw(IndexException(m,n,*this)); |
---|
176 | return store[tristore(m-1)+n-1]; |
---|
177 | } |
---|
178 | |
---|
179 | Real DiagonalMatrix::operator()(int m, int n) const |
---|
180 | { |
---|
181 | REPORT |
---|
182 | if (n<=0 || m!=n || m>nrows_val || n>ncols_val) |
---|
183 | Throw(IndexException(m,n,*this)); |
---|
184 | return store[n-1]; |
---|
185 | } |
---|
186 | |
---|
187 | Real DiagonalMatrix::operator()(int m) const |
---|
188 | { |
---|
189 | REPORT |
---|
190 | if (m<=0 || m>nrows_val) Throw(IndexException(m,*this)); |
---|
191 | return store[m-1]; |
---|
192 | } |
---|
193 | |
---|
194 | Real ColumnVector::operator()(int m) const |
---|
195 | { |
---|
196 | REPORT |
---|
197 | if (m<=0 || m> nrows_val) Throw(IndexException(m,*this)); |
---|
198 | return store[m-1]; |
---|
199 | } |
---|
200 | |
---|
201 | Real RowVector::operator()(int n) const |
---|
202 | { |
---|
203 | REPORT |
---|
204 | if (n<=0 || n> ncols_val) Throw(IndexException(n,*this)); |
---|
205 | return store[n-1]; |
---|
206 | } |
---|
207 | |
---|
208 | Real BandMatrix::operator()(int m, int n) const |
---|
209 | { |
---|
210 | REPORT |
---|
211 | int w = upper_val+lower_val+1; int i = lower_val+n-m; |
---|
212 | if (m<=0 || m>nrows_val || n<=0 || n>ncols_val || i<0 || i>=w) |
---|
213 | Throw(IndexException(m,n,*this)); |
---|
214 | return store[w*(m-1)+i]; |
---|
215 | } |
---|
216 | |
---|
217 | Real UpperBandMatrix::operator()(int m, int n) const |
---|
218 | { |
---|
219 | REPORT |
---|
220 | int w = upper_val+1; int i = n-m; |
---|
221 | if (m<=0 || m>nrows_val || n<=0 || n>ncols_val || i<0 || i>=w) |
---|
222 | Throw(IndexException(m,n,*this)); |
---|
223 | return store[w*(m-1)+i]; |
---|
224 | } |
---|
225 | |
---|
226 | Real LowerBandMatrix::operator()(int m, int n) const |
---|
227 | { |
---|
228 | REPORT |
---|
229 | int w = lower_val+1; int i = lower_val+n-m; |
---|
230 | if (m<=0 || m>nrows_val || n<=0 || n>ncols_val || i<0 || i>=w) |
---|
231 | Throw(IndexException(m,n,*this)); |
---|
232 | return store[w*(m-1)+i]; |
---|
233 | } |
---|
234 | |
---|
235 | Real SymmetricBandMatrix::operator()(int m, int n) const |
---|
236 | { |
---|
237 | REPORT |
---|
238 | int w = lower_val+1; |
---|
239 | if (m>=n) |
---|
240 | { |
---|
241 | REPORT |
---|
242 | int i = lower_val+n-m; |
---|
243 | if ( m>nrows_val || n<=0 || i<0 ) |
---|
244 | Throw(IndexException(m,n,*this)); |
---|
245 | return store[w*(m-1)+i]; |
---|
246 | } |
---|
247 | else |
---|
248 | { |
---|
249 | REPORT |
---|
250 | int i = lower_val+m-n; |
---|
251 | if ( n>nrows_val || m<=0 || i<0 ) |
---|
252 | Throw(IndexException(m,n,*this)); |
---|
253 | return store[w*(n-1)+i]; |
---|
254 | } |
---|
255 | } |
---|
256 | |
---|
257 | |
---|
258 | Real BaseMatrix::as_scalar() const |
---|
259 | { |
---|
260 | REPORT |
---|
261 | GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate(); |
---|
262 | |
---|
263 | if (gm->nrows_val!=1 || gm->ncols_val!=1) |
---|
264 | { |
---|
265 | Tracer tr("as_scalar"); |
---|
266 | Try |
---|
267 | { Throw(ProgramException("Cannot convert to scalar", *gm)); } |
---|
268 | CatchAll { gm->tDelete(); ReThrow; } |
---|
269 | } |
---|
270 | |
---|
271 | Real x = *(gm->store); gm->tDelete(); return x; |
---|
272 | } |
---|
273 | |
---|
274 | |
---|
275 | AddedMatrix BaseMatrix::operator+(const BaseMatrix& bm) const |
---|
276 | { REPORT return AddedMatrix(this, &bm); } |
---|
277 | |
---|
278 | SPMatrix SP(const BaseMatrix& bm1,const BaseMatrix& bm2) |
---|
279 | { REPORT return SPMatrix(&bm1, &bm2); } |
---|
280 | |
---|
281 | KPMatrix KP(const BaseMatrix& bm1,const BaseMatrix& bm2) |
---|
282 | { REPORT return KPMatrix(&bm1, &bm2); } |
---|
283 | |
---|
284 | MultipliedMatrix BaseMatrix::operator*(const BaseMatrix& bm) const |
---|
285 | { REPORT return MultipliedMatrix(this, &bm); } |
---|
286 | |
---|
287 | ConcatenatedMatrix BaseMatrix::operator|(const BaseMatrix& bm) const |
---|
288 | { REPORT return ConcatenatedMatrix(this, &bm); } |
---|
289 | |
---|
290 | StackedMatrix BaseMatrix::operator&(const BaseMatrix& bm) const |
---|
291 | { REPORT return StackedMatrix(this, &bm); } |
---|
292 | |
---|
293 | SolvedMatrix InvertedMatrix::operator*(const BaseMatrix& bmx) const |
---|
294 | { REPORT return SolvedMatrix(bm, &bmx); } |
---|
295 | |
---|
296 | SubtractedMatrix BaseMatrix::operator-(const BaseMatrix& bm) const |
---|
297 | { REPORT return SubtractedMatrix(this, &bm); } |
---|
298 | |
---|
299 | ShiftedMatrix BaseMatrix::operator+(Real f) const |
---|
300 | { REPORT return ShiftedMatrix(this, f); } |
---|
301 | |
---|
302 | ShiftedMatrix operator+(Real f, const BaseMatrix& BM) |
---|
303 | { REPORT return ShiftedMatrix(&BM, f); } |
---|
304 | |
---|
305 | NegShiftedMatrix operator-(Real f, const BaseMatrix& bm) |
---|
306 | { REPORT return NegShiftedMatrix(f, &bm); } |
---|
307 | |
---|
308 | ScaledMatrix BaseMatrix::operator*(Real f) const |
---|
309 | { REPORT return ScaledMatrix(this, f); } |
---|
310 | |
---|
311 | ScaledMatrix BaseMatrix::operator/(Real f) const |
---|
312 | { REPORT return ScaledMatrix(this, 1.0/f); } |
---|
313 | |
---|
314 | ScaledMatrix operator*(Real f, const BaseMatrix& BM) |
---|
315 | { REPORT return ScaledMatrix(&BM, f); } |
---|
316 | |
---|
317 | ShiftedMatrix BaseMatrix::operator-(Real f) const |
---|
318 | { REPORT return ShiftedMatrix(this, -f); } |
---|
319 | |
---|
320 | TransposedMatrix BaseMatrix::t() const |
---|
321 | { REPORT return TransposedMatrix(this); } |
---|
322 | |
---|
323 | NegatedMatrix BaseMatrix::operator-() const |
---|
324 | { REPORT return NegatedMatrix(this); } |
---|
325 | |
---|
326 | ReversedMatrix BaseMatrix::reverse() const |
---|
327 | { REPORT return ReversedMatrix(this); } |
---|
328 | |
---|
329 | InvertedMatrix BaseMatrix::i() const |
---|
330 | { REPORT return InvertedMatrix(this); } |
---|
331 | |
---|
332 | |
---|
333 | RowedMatrix BaseMatrix::as_row() const |
---|
334 | { REPORT return RowedMatrix(this); } |
---|
335 | |
---|
336 | ColedMatrix BaseMatrix::as_column() const |
---|
337 | { REPORT return ColedMatrix(this); } |
---|
338 | |
---|
339 | DiagedMatrix BaseMatrix::as_diagonal() const |
---|
340 | { REPORT return DiagedMatrix(this); } |
---|
341 | |
---|
342 | MatedMatrix BaseMatrix::as_matrix(int nrx, int ncx) const |
---|
343 | { REPORT return MatedMatrix(this,nrx,ncx); } |
---|
344 | |
---|
345 | |
---|
346 | void GeneralMatrix::operator=(Real f) |
---|
347 | { REPORT int i=storage; Real* s=store; while (i--) { *s++ = f; } } |
---|
348 | |
---|
349 | void Matrix::operator=(const BaseMatrix& X) |
---|
350 | { |
---|
351 | REPORT //CheckConversion(X); |
---|
352 | // MatrixConversionCheck mcc; |
---|
353 | Eq(X,MatrixType::Rt); |
---|
354 | } |
---|
355 | |
---|
356 | void SquareMatrix::operator=(const BaseMatrix& X) |
---|
357 | { |
---|
358 | REPORT //CheckConversion(X); |
---|
359 | // MatrixConversionCheck mcc; |
---|
360 | Eq(X,MatrixType::Rt); |
---|
361 | if (nrows_val != ncols_val) |
---|
362 | { Tracer tr("SquareMatrix(=)"); Throw(NotSquareException(*this)); } |
---|
363 | } |
---|
364 | |
---|
365 | void SquareMatrix::operator=(const Matrix& m) |
---|
366 | { |
---|
367 | REPORT |
---|
368 | if (m.nrows_val != m.ncols_val) |
---|
369 | { Tracer tr("SquareMatrix(=Matrix)"); Throw(NotSquareException(*this)); } |
---|
370 | Eq(m); |
---|
371 | } |
---|
372 | |
---|
373 | void RowVector::operator=(const BaseMatrix& X) |
---|
374 | { |
---|
375 | REPORT // CheckConversion(X); |
---|
376 | // MatrixConversionCheck mcc; |
---|
377 | Eq(X,MatrixType::RV); |
---|
378 | if (nrows_val!=1) |
---|
379 | { Tracer tr("RowVector(=)"); Throw(VectorException(*this)); } |
---|
380 | } |
---|
381 | |
---|
382 | void ColumnVector::operator=(const BaseMatrix& X) |
---|
383 | { |
---|
384 | REPORT //CheckConversion(X); |
---|
385 | // MatrixConversionCheck mcc; |
---|
386 | Eq(X,MatrixType::CV); |
---|
387 | if (ncols_val!=1) |
---|
388 | { Tracer tr("ColumnVector(=)"); Throw(VectorException(*this)); } |
---|
389 | } |
---|
390 | |
---|
391 | void SymmetricMatrix::operator=(const BaseMatrix& X) |
---|
392 | { |
---|
393 | REPORT // CheckConversion(X); |
---|
394 | // MatrixConversionCheck mcc; |
---|
395 | Eq(X,MatrixType::Sm); |
---|
396 | } |
---|
397 | |
---|
398 | void UpperTriangularMatrix::operator=(const BaseMatrix& X) |
---|
399 | { |
---|
400 | REPORT //CheckConversion(X); |
---|
401 | // MatrixConversionCheck mcc; |
---|
402 | Eq(X,MatrixType::UT); |
---|
403 | } |
---|
404 | |
---|
405 | void LowerTriangularMatrix::operator=(const BaseMatrix& X) |
---|
406 | { |
---|
407 | REPORT //CheckConversion(X); |
---|
408 | // MatrixConversionCheck mcc; |
---|
409 | Eq(X,MatrixType::LT); |
---|
410 | } |
---|
411 | |
---|
412 | void DiagonalMatrix::operator=(const BaseMatrix& X) |
---|
413 | { |
---|
414 | REPORT // CheckConversion(X); |
---|
415 | // MatrixConversionCheck mcc; |
---|
416 | Eq(X,MatrixType::Dg); |
---|
417 | } |
---|
418 | |
---|
419 | void IdentityMatrix::operator=(const BaseMatrix& X) |
---|
420 | { |
---|
421 | REPORT // CheckConversion(X); |
---|
422 | // MatrixConversionCheck mcc; |
---|
423 | Eq(X,MatrixType::Id); |
---|
424 | } |
---|
425 | |
---|
426 | |
---|
427 | void CroutMatrix::operator=(const CroutMatrix& gm) |
---|
428 | { |
---|
429 | if (&gm == this) { REPORT tag_val = -1; return; } |
---|
430 | REPORT |
---|
431 | if (indx > 0) { delete [] indx; indx = 0; } |
---|
432 | ((CroutMatrix&)gm).get_aux(*this); |
---|
433 | Eq(gm); |
---|
434 | } |
---|
435 | |
---|
436 | |
---|
437 | |
---|
438 | |
---|
439 | |
---|
440 | void GeneralMatrix::operator<<(const double* r) |
---|
441 | { |
---|
442 | REPORT |
---|
443 | int i = storage; Real* s=store; |
---|
444 | while(i--) *s++ = (Real)*r++; |
---|
445 | } |
---|
446 | |
---|
447 | |
---|
448 | void GeneralMatrix::operator<<(const float* r) |
---|
449 | { |
---|
450 | REPORT |
---|
451 | int i = storage; Real* s=store; |
---|
452 | while(i--) *s++ = (Real)*r++; |
---|
453 | } |
---|
454 | |
---|
455 | |
---|
456 | void GeneralMatrix::operator<<(const int* r) |
---|
457 | { |
---|
458 | REPORT |
---|
459 | int i = storage; Real* s=store; |
---|
460 | while(i--) *s++ = (Real)*r++; |
---|
461 | } |
---|
462 | |
---|
463 | |
---|
464 | void GenericMatrix::operator=(const GenericMatrix& bmx) |
---|
465 | { |
---|
466 | if (&bmx != this) { REPORT if (gm) delete gm; gm = bmx.gm->Image();} |
---|
467 | else { REPORT } |
---|
468 | gm->Protect(); |
---|
469 | } |
---|
470 | |
---|
471 | void GenericMatrix::operator=(const BaseMatrix& bmx) |
---|
472 | { |
---|
473 | if (gm) |
---|
474 | { |
---|
475 | int counter=bmx.search(gm); |
---|
476 | if (counter==0) { REPORT delete gm; gm=0; } |
---|
477 | else { REPORT gm->Release(counter); } |
---|
478 | } |
---|
479 | else { REPORT } |
---|
480 | GeneralMatrix* gmx = ((BaseMatrix&)bmx).Evaluate(); |
---|
481 | if (gmx != gm) { REPORT if (gm) delete gm; gm = gmx->Image(); } |
---|
482 | else { REPORT } |
---|
483 | gm->Protect(); |
---|
484 | } |
---|
485 | |
---|
486 | |
---|
487 | /*************************** += etc ***************************************/ |
---|
488 | |
---|
489 | |
---|
490 | // GeneralMatrix operators |
---|
491 | |
---|
492 | void GeneralMatrix::operator+=(const BaseMatrix& X) |
---|
493 | { |
---|
494 | REPORT |
---|
495 | Tracer tr("GeneralMatrix::operator+="); |
---|
496 | // MatrixConversionCheck mcc; |
---|
497 | Protect(); // so it cannot get deleted |
---|
498 | // during Evaluate |
---|
499 | GeneralMatrix* gm = ((BaseMatrix&)X).Evaluate(); |
---|
500 | AddedMatrix am(this,gm); |
---|
501 | if (gm==this) Release(2); else Release(); |
---|
502 | Eq2(am,type()); |
---|
503 | } |
---|
504 | |
---|
505 | void GeneralMatrix::operator-=(const BaseMatrix& X) |
---|
506 | { |
---|
507 | REPORT |
---|
508 | Tracer tr("GeneralMatrix::operator-="); |
---|
509 | // MatrixConversionCheck mcc; |
---|
510 | Protect(); // so it cannot get deleted |
---|
511 | // during Evaluate |
---|
512 | GeneralMatrix* gm = ((BaseMatrix&)X).Evaluate(); |
---|
513 | SubtractedMatrix am(this,gm); |
---|
514 | if (gm==this) Release(2); else Release(); |
---|
515 | Eq2(am,type()); |
---|
516 | } |
---|
517 | |
---|
518 | void GeneralMatrix::operator*=(const BaseMatrix& X) |
---|
519 | { |
---|
520 | REPORT |
---|
521 | Tracer tr("GeneralMatrix::operator*="); |
---|
522 | // MatrixConversionCheck mcc; |
---|
523 | Protect(); // so it cannot get deleted |
---|
524 | // during Evaluate |
---|
525 | GeneralMatrix* gm = ((BaseMatrix&)X).Evaluate(); |
---|
526 | MultipliedMatrix am(this,gm); |
---|
527 | if (gm==this) Release(2); else Release(); |
---|
528 | Eq2(am,type()); |
---|
529 | } |
---|
530 | |
---|
531 | void GeneralMatrix::operator|=(const BaseMatrix& X) |
---|
532 | { |
---|
533 | REPORT |
---|
534 | Tracer tr("GeneralMatrix::operator|="); |
---|
535 | // MatrixConversionCheck mcc; |
---|
536 | Protect(); // so it cannot get deleted |
---|
537 | // during Evaluate |
---|
538 | GeneralMatrix* gm = ((BaseMatrix&)X).Evaluate(); |
---|
539 | ConcatenatedMatrix am(this,gm); |
---|
540 | if (gm==this) Release(2); else Release(); |
---|
541 | Eq2(am,type()); |
---|
542 | } |
---|
543 | |
---|
544 | void GeneralMatrix::operator&=(const BaseMatrix& X) |
---|
545 | { |
---|
546 | REPORT |
---|
547 | Tracer tr("GeneralMatrix::operator&="); |
---|
548 | // MatrixConversionCheck mcc; |
---|
549 | Protect(); // so it cannot get deleted |
---|
550 | // during Evaluate |
---|
551 | GeneralMatrix* gm = ((BaseMatrix&)X).Evaluate(); |
---|
552 | StackedMatrix am(this,gm); |
---|
553 | if (gm==this) Release(2); else Release(); |
---|
554 | Eq2(am,type()); |
---|
555 | } |
---|
556 | |
---|
557 | void GeneralMatrix::operator+=(Real r) |
---|
558 | { |
---|
559 | REPORT |
---|
560 | Tracer tr("GeneralMatrix::operator+=(Real)"); |
---|
561 | // MatrixConversionCheck mcc; |
---|
562 | ShiftedMatrix am(this,r); |
---|
563 | Release(); Eq2(am,type()); |
---|
564 | } |
---|
565 | |
---|
566 | void GeneralMatrix::operator*=(Real r) |
---|
567 | { |
---|
568 | REPORT |
---|
569 | Tracer tr("GeneralMatrix::operator*=(Real)"); |
---|
570 | // MatrixConversionCheck mcc; |
---|
571 | ScaledMatrix am(this,r); |
---|
572 | Release(); Eq2(am,type()); |
---|
573 | } |
---|
574 | |
---|
575 | |
---|
576 | // Generic matrix operators |
---|
577 | |
---|
578 | void GenericMatrix::operator+=(const BaseMatrix& X) |
---|
579 | { |
---|
580 | REPORT |
---|
581 | Tracer tr("GenericMatrix::operator+="); |
---|
582 | if (!gm) Throw(ProgramException("GenericMatrix is null")); |
---|
583 | gm->Protect(); // so it cannot get deleted during Evaluate |
---|
584 | GeneralMatrix* gmx = ((BaseMatrix&)X).Evaluate(); |
---|
585 | AddedMatrix am(gm,gmx); |
---|
586 | if (gmx==gm) gm->Release(2); else gm->Release(); |
---|
587 | GeneralMatrix* gmy = am.Evaluate(); |
---|
588 | if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); } |
---|
589 | else { REPORT } |
---|
590 | gm->Protect(); |
---|
591 | } |
---|
592 | |
---|
593 | void GenericMatrix::operator-=(const BaseMatrix& X) |
---|
594 | { |
---|
595 | REPORT |
---|
596 | Tracer tr("GenericMatrix::operator-="); |
---|
597 | if (!gm) Throw(ProgramException("GenericMatrix is null")); |
---|
598 | gm->Protect(); // so it cannot get deleted during Evaluate |
---|
599 | GeneralMatrix* gmx = ((BaseMatrix&)X).Evaluate(); |
---|
600 | SubtractedMatrix am(gm,gmx); |
---|
601 | if (gmx==gm) gm->Release(2); else gm->Release(); |
---|
602 | GeneralMatrix* gmy = am.Evaluate(); |
---|
603 | if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); } |
---|
604 | else { REPORT } |
---|
605 | gm->Protect(); |
---|
606 | } |
---|
607 | |
---|
608 | void GenericMatrix::operator*=(const BaseMatrix& X) |
---|
609 | { |
---|
610 | REPORT |
---|
611 | Tracer tr("GenericMatrix::operator*="); |
---|
612 | if (!gm) Throw(ProgramException("GenericMatrix is null")); |
---|
613 | gm->Protect(); // so it cannot get deleted during Evaluate |
---|
614 | GeneralMatrix* gmx = ((BaseMatrix&)X).Evaluate(); |
---|
615 | MultipliedMatrix am(gm,gmx); |
---|
616 | if (gmx==gm) gm->Release(2); else gm->Release(); |
---|
617 | GeneralMatrix* gmy = am.Evaluate(); |
---|
618 | if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); } |
---|
619 | else { REPORT } |
---|
620 | gm->Protect(); |
---|
621 | } |
---|
622 | |
---|
623 | void GenericMatrix::operator|=(const BaseMatrix& X) |
---|
624 | { |
---|
625 | REPORT |
---|
626 | Tracer tr("GenericMatrix::operator|="); |
---|
627 | if (!gm) Throw(ProgramException("GenericMatrix is null")); |
---|
628 | gm->Protect(); // so it cannot get deleted during Evaluate |
---|
629 | GeneralMatrix* gmx = ((BaseMatrix&)X).Evaluate(); |
---|
630 | ConcatenatedMatrix am(gm,gmx); |
---|
631 | if (gmx==gm) gm->Release(2); else gm->Release(); |
---|
632 | GeneralMatrix* gmy = am.Evaluate(); |
---|
633 | if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); } |
---|
634 | else { REPORT } |
---|
635 | gm->Protect(); |
---|
636 | } |
---|
637 | |
---|
638 | void GenericMatrix::operator&=(const BaseMatrix& X) |
---|
639 | { |
---|
640 | REPORT |
---|
641 | Tracer tr("GenericMatrix::operator&="); |
---|
642 | if (!gm) Throw(ProgramException("GenericMatrix is null")); |
---|
643 | gm->Protect(); // so it cannot get deleted during Evaluate |
---|
644 | GeneralMatrix* gmx = ((BaseMatrix&)X).Evaluate(); |
---|
645 | StackedMatrix am(gm,gmx); |
---|
646 | if (gmx==gm) gm->Release(2); else gm->Release(); |
---|
647 | GeneralMatrix* gmy = am.Evaluate(); |
---|
648 | if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); } |
---|
649 | else { REPORT } |
---|
650 | gm->Protect(); |
---|
651 | } |
---|
652 | |
---|
653 | void GenericMatrix::operator+=(Real r) |
---|
654 | { |
---|
655 | REPORT |
---|
656 | Tracer tr("GenericMatrix::operator+= (Real)"); |
---|
657 | if (!gm) Throw(ProgramException("GenericMatrix is null")); |
---|
658 | ShiftedMatrix am(gm,r); |
---|
659 | gm->Release(); |
---|
660 | GeneralMatrix* gmy = am.Evaluate(); |
---|
661 | if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); } |
---|
662 | else { REPORT } |
---|
663 | gm->Protect(); |
---|
664 | } |
---|
665 | |
---|
666 | void GenericMatrix::operator*=(Real r) |
---|
667 | { |
---|
668 | REPORT |
---|
669 | Tracer tr("GenericMatrix::operator*= (Real)"); |
---|
670 | if (!gm) Throw(ProgramException("GenericMatrix is null")); |
---|
671 | ScaledMatrix am(gm,r); |
---|
672 | gm->Release(); |
---|
673 | GeneralMatrix* gmy = am.Evaluate(); |
---|
674 | if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); } |
---|
675 | else { REPORT } |
---|
676 | gm->Protect(); |
---|
677 | } |
---|
678 | |
---|
679 | |
---|
680 | /************************* element access *********************************/ |
---|
681 | |
---|
682 | Real& Matrix::element(int m, int n) |
---|
683 | { |
---|
684 | REPORT |
---|
685 | if (m<0 || m>= nrows_val || n<0 || n>= ncols_val) |
---|
686 | Throw(IndexException(m,n,*this,true)); |
---|
687 | return store[m*ncols_val+n]; |
---|
688 | } |
---|
689 | |
---|
690 | Real Matrix::element(int m, int n) const |
---|
691 | { |
---|
692 | REPORT |
---|
693 | if (m<0 || m>= nrows_val || n<0 || n>= ncols_val) |
---|
694 | Throw(IndexException(m,n,*this,true)); |
---|
695 | return store[m*ncols_val+n]; |
---|
696 | } |
---|
697 | |
---|
698 | Real& SymmetricMatrix::element(int m, int n) |
---|
699 | { |
---|
700 | REPORT |
---|
701 | if (m<0 || n<0 || m >= nrows_val || n>=ncols_val) |
---|
702 | Throw(IndexException(m,n,*this,true)); |
---|
703 | if (m>=n) return store[tristore(m)+n]; |
---|
704 | else return store[tristore(n)+m]; |
---|
705 | } |
---|
706 | |
---|
707 | Real SymmetricMatrix::element(int m, int n) const |
---|
708 | { |
---|
709 | REPORT |
---|
710 | if (m<0 || n<0 || m >= nrows_val || n>=ncols_val) |
---|
711 | Throw(IndexException(m,n,*this,true)); |
---|
712 | if (m>=n) return store[tristore(m)+n]; |
---|
713 | else return store[tristore(n)+m]; |
---|
714 | } |
---|
715 | |
---|
716 | Real& UpperTriangularMatrix::element(int m, int n) |
---|
717 | { |
---|
718 | REPORT |
---|
719 | if (m<0 || n<m || n>=ncols_val) |
---|
720 | Throw(IndexException(m,n,*this,true)); |
---|
721 | return store[m*ncols_val+n-tristore(m)]; |
---|
722 | } |
---|
723 | |
---|
724 | Real UpperTriangularMatrix::element(int m, int n) const |
---|
725 | { |
---|
726 | REPORT |
---|
727 | if (m<0 || n<m || n>=ncols_val) |
---|
728 | Throw(IndexException(m,n,*this,true)); |
---|
729 | return store[m*ncols_val+n-tristore(m)]; |
---|
730 | } |
---|
731 | |
---|
732 | Real& LowerTriangularMatrix::element(int m, int n) |
---|
733 | { |
---|
734 | REPORT |
---|
735 | if (n<0 || m<n || m>=nrows_val) |
---|
736 | Throw(IndexException(m,n,*this,true)); |
---|
737 | return store[tristore(m)+n]; |
---|
738 | } |
---|
739 | |
---|
740 | Real LowerTriangularMatrix::element(int m, int n) const |
---|
741 | { |
---|
742 | REPORT |
---|
743 | if (n<0 || m<n || m>=nrows_val) |
---|
744 | Throw(IndexException(m,n,*this,true)); |
---|
745 | return store[tristore(m)+n]; |
---|
746 | } |
---|
747 | |
---|
748 | Real& DiagonalMatrix::element(int m, int n) |
---|
749 | { |
---|
750 | REPORT |
---|
751 | if (n<0 || m!=n || m>=nrows_val || n>=ncols_val) |
---|
752 | Throw(IndexException(m,n,*this,true)); |
---|
753 | return store[n]; |
---|
754 | } |
---|
755 | |
---|
756 | Real DiagonalMatrix::element(int m, int n) const |
---|
757 | { |
---|
758 | REPORT |
---|
759 | if (n<0 || m!=n || m>=nrows_val || n>=ncols_val) |
---|
760 | Throw(IndexException(m,n,*this,true)); |
---|
761 | return store[n]; |
---|
762 | } |
---|
763 | |
---|
764 | Real& DiagonalMatrix::element(int m) |
---|
765 | { |
---|
766 | REPORT |
---|
767 | if (m<0 || m>=nrows_val) Throw(IndexException(m,*this,true)); |
---|
768 | return store[m]; |
---|
769 | } |
---|
770 | |
---|
771 | Real DiagonalMatrix::element(int m) const |
---|
772 | { |
---|
773 | REPORT |
---|
774 | if (m<0 || m>=nrows_val) Throw(IndexException(m,*this,true)); |
---|
775 | return store[m]; |
---|
776 | } |
---|
777 | |
---|
778 | Real& ColumnVector::element(int m) |
---|
779 | { |
---|
780 | REPORT |
---|
781 | if (m<0 || m>= nrows_val) Throw(IndexException(m,*this,true)); |
---|
782 | return store[m]; |
---|
783 | } |
---|
784 | |
---|
785 | Real ColumnVector::element(int m) const |
---|
786 | { |
---|
787 | REPORT |
---|
788 | if (m<0 || m>= nrows_val) Throw(IndexException(m,*this,true)); |
---|
789 | return store[m]; |
---|
790 | } |
---|
791 | |
---|
792 | Real& RowVector::element(int n) |
---|
793 | { |
---|
794 | REPORT |
---|
795 | if (n<0 || n>= ncols_val) Throw(IndexException(n,*this,true)); |
---|
796 | return store[n]; |
---|
797 | } |
---|
798 | |
---|
799 | Real RowVector::element(int n) const |
---|
800 | { |
---|
801 | REPORT |
---|
802 | if (n<0 || n>= ncols_val) Throw(IndexException(n,*this,true)); |
---|
803 | return store[n]; |
---|
804 | } |
---|
805 | |
---|
806 | Real& BandMatrix::element(int m, int n) |
---|
807 | { |
---|
808 | REPORT |
---|
809 | int w = upper_val+lower_val+1; int i = lower_val+n-m; |
---|
810 | if (m<0 || m>= nrows_val || n<0 || n>= ncols_val || i<0 || i>=w) |
---|
811 | Throw(IndexException(m,n,*this,true)); |
---|
812 | return store[w*m+i]; |
---|
813 | } |
---|
814 | |
---|
815 | Real BandMatrix::element(int m, int n) const |
---|
816 | { |
---|
817 | REPORT |
---|
818 | int w = upper_val+lower_val+1; int i = lower_val+n-m; |
---|
819 | if (m<0 || m>= nrows_val || n<0 || n>= ncols_val || i<0 || i>=w) |
---|
820 | Throw(IndexException(m,n,*this,true)); |
---|
821 | return store[w*m+i]; |
---|
822 | } |
---|
823 | |
---|
824 | Real& UpperBandMatrix::element(int m, int n) |
---|
825 | { |
---|
826 | REPORT |
---|
827 | int w = upper_val+1; int i = n-m; |
---|
828 | if (m<0 || m>= nrows_val || n<0 || n>= ncols_val || i<0 || i>=w) |
---|
829 | Throw(IndexException(m,n,*this,true)); |
---|
830 | return store[w*m+i]; |
---|
831 | } |
---|
832 | |
---|
833 | Real UpperBandMatrix::element(int m, int n) const |
---|
834 | { |
---|
835 | REPORT |
---|
836 | int w = upper_val+1; int i = n-m; |
---|
837 | if (m<0 || m>= nrows_val || n<0 || n>= ncols_val || i<0 || i>=w) |
---|
838 | Throw(IndexException(m,n,*this,true)); |
---|
839 | return store[w*m+i]; |
---|
840 | } |
---|
841 | |
---|
842 | Real& LowerBandMatrix::element(int m, int n) |
---|
843 | { |
---|
844 | REPORT |
---|
845 | int w = lower_val+1; int i = lower_val+n-m; |
---|
846 | if (m<0 || m>= nrows_val || n<0 || n>= ncols_val || i<0 || i>=w) |
---|
847 | Throw(IndexException(m,n,*this,true)); |
---|
848 | return store[w*m+i]; |
---|
849 | } |
---|
850 | |
---|
851 | Real LowerBandMatrix::element(int m, int n) const |
---|
852 | { |
---|
853 | REPORT |
---|
854 | int w = lower_val+1; int i = lower_val+n-m; |
---|
855 | if (m<0 || m>= nrows_val || n<0 || n>= ncols_val || i<0 || i>=w) |
---|
856 | Throw(IndexException(m,n,*this,true)); |
---|
857 | return store[w*m+i]; |
---|
858 | } |
---|
859 | |
---|
860 | Real& SymmetricBandMatrix::element(int m, int n) |
---|
861 | { |
---|
862 | REPORT |
---|
863 | int w = lower_val+1; |
---|
864 | if (m>=n) |
---|
865 | { |
---|
866 | REPORT |
---|
867 | int i = lower_val+n-m; |
---|
868 | if ( m>=nrows_val || n<0 || i<0 ) |
---|
869 | Throw(IndexException(m,n,*this,true)); |
---|
870 | return store[w*m+i]; |
---|
871 | } |
---|
872 | else |
---|
873 | { |
---|
874 | REPORT |
---|
875 | int i = lower_val+m-n; |
---|
876 | if ( n>=nrows_val || m<0 || i<0 ) |
---|
877 | Throw(IndexException(m,n,*this,true)); |
---|
878 | return store[w*n+i]; |
---|
879 | } |
---|
880 | } |
---|
881 | |
---|
882 | Real SymmetricBandMatrix::element(int m, int n) const |
---|
883 | { |
---|
884 | REPORT |
---|
885 | int w = lower_val+1; |
---|
886 | if (m>=n) |
---|
887 | { |
---|
888 | REPORT |
---|
889 | int i = lower_val+n-m; |
---|
890 | if ( m>=nrows_val || n<0 || i<0 ) |
---|
891 | Throw(IndexException(m,n,*this,true)); |
---|
892 | return store[w*m+i]; |
---|
893 | } |
---|
894 | else |
---|
895 | { |
---|
896 | REPORT |
---|
897 | int i = lower_val+m-n; |
---|
898 | if ( n>=nrows_val || m<0 || i<0 ) |
---|
899 | Throw(IndexException(m,n,*this,true)); |
---|
900 | return store[w*n+i]; |
---|
901 | } |
---|
902 | } |
---|
903 | |
---|
904 | #ifdef use_namespace |
---|
905 | } |
---|
906 | #endif |
---|
907 | |
---|
908 | |
---|
909 | ///} |
---|