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