source: ntrip/trunk/BNC/newmat/newmat6.cpp@ 8443

Last change on this file since 8443 was 2013, checked in by mervart, 15 years ago

* empty log message *
File size: 21.9 KB
Line 
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
15namespace 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
28static int tristore(int n) // els in triangular matrix
29{ return (n*(n+1))/2; }
30
31
32/****************************** operators *******************************/
33
34Real& 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
42Real& 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
51Real& 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
59Real& 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
67Real& 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
75Real& 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
82Real& 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
89Real& 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
96Real& 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
105Real& 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
114Real& 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
123Real& 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
146Real 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
154Real 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
163Real 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
171Real 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
179Real 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
187Real 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
194Real 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
201Real 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
208Real 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
217Real 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
226Real 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
235Real 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
258Real 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
275AddedMatrix BaseMatrix::operator+(const BaseMatrix& bm) const
276{ REPORT return AddedMatrix(this, &bm); }
277
278SPMatrix SP(const BaseMatrix& bm1,const BaseMatrix& bm2)
279{ REPORT return SPMatrix(&bm1, &bm2); }
280
281KPMatrix KP(const BaseMatrix& bm1,const BaseMatrix& bm2)
282{ REPORT return KPMatrix(&bm1, &bm2); }
283
284MultipliedMatrix BaseMatrix::operator*(const BaseMatrix& bm) const
285{ REPORT return MultipliedMatrix(this, &bm); }
286
287ConcatenatedMatrix BaseMatrix::operator|(const BaseMatrix& bm) const
288{ REPORT return ConcatenatedMatrix(this, &bm); }
289
290StackedMatrix BaseMatrix::operator&(const BaseMatrix& bm) const
291{ REPORT return StackedMatrix(this, &bm); }
292
293SolvedMatrix InvertedMatrix::operator*(const BaseMatrix& bmx) const
294{ REPORT return SolvedMatrix(bm, &bmx); }
295
296SubtractedMatrix BaseMatrix::operator-(const BaseMatrix& bm) const
297{ REPORT return SubtractedMatrix(this, &bm); }
298
299ShiftedMatrix BaseMatrix::operator+(Real f) const
300{ REPORT return ShiftedMatrix(this, f); }
301
302ShiftedMatrix operator+(Real f, const BaseMatrix& BM)
303{ REPORT return ShiftedMatrix(&BM, f); }
304
305NegShiftedMatrix operator-(Real f, const BaseMatrix& bm)
306{ REPORT return NegShiftedMatrix(f, &bm); }
307
308ScaledMatrix BaseMatrix::operator*(Real f) const
309{ REPORT return ScaledMatrix(this, f); }
310
311ScaledMatrix BaseMatrix::operator/(Real f) const
312{ REPORT return ScaledMatrix(this, 1.0/f); }
313
314ScaledMatrix operator*(Real f, const BaseMatrix& BM)
315{ REPORT return ScaledMatrix(&BM, f); }
316
317ShiftedMatrix BaseMatrix::operator-(Real f) const
318{ REPORT return ShiftedMatrix(this, -f); }
319
320TransposedMatrix BaseMatrix::t() const
321{ REPORT return TransposedMatrix(this); }
322
323NegatedMatrix BaseMatrix::operator-() const
324{ REPORT return NegatedMatrix(this); }
325
326ReversedMatrix BaseMatrix::reverse() const
327{ REPORT return ReversedMatrix(this); }
328
329InvertedMatrix BaseMatrix::i() const
330{ REPORT return InvertedMatrix(this); }
331
332
333RowedMatrix BaseMatrix::as_row() const
334{ REPORT return RowedMatrix(this); }
335
336ColedMatrix BaseMatrix::as_column() const
337{ REPORT return ColedMatrix(this); }
338
339DiagedMatrix BaseMatrix::as_diagonal() const
340{ REPORT return DiagedMatrix(this); }
341
342MatedMatrix BaseMatrix::as_matrix(int nrx, int ncx) const
343{ REPORT return MatedMatrix(this,nrx,ncx); }
344
345
346void GeneralMatrix::operator=(Real f)
347{ REPORT int i=storage; Real* s=store; while (i--) { *s++ = f; } }
348
349void Matrix::operator=(const BaseMatrix& X)
350{
351 REPORT //CheckConversion(X);
352 // MatrixConversionCheck mcc;
353 Eq(X,MatrixType::Rt);
354}
355
356void 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
365void 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
373void 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
382void 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
391void SymmetricMatrix::operator=(const BaseMatrix& X)
392{
393 REPORT // CheckConversion(X);
394 // MatrixConversionCheck mcc;
395 Eq(X,MatrixType::Sm);
396}
397
398void UpperTriangularMatrix::operator=(const BaseMatrix& X)
399{
400 REPORT //CheckConversion(X);
401 // MatrixConversionCheck mcc;
402 Eq(X,MatrixType::UT);
403}
404
405void LowerTriangularMatrix::operator=(const BaseMatrix& X)
406{
407 REPORT //CheckConversion(X);
408 // MatrixConversionCheck mcc;
409 Eq(X,MatrixType::LT);
410}
411
412void DiagonalMatrix::operator=(const BaseMatrix& X)
413{
414 REPORT // CheckConversion(X);
415 // MatrixConversionCheck mcc;
416 Eq(X,MatrixType::Dg);
417}
418
419void IdentityMatrix::operator=(const BaseMatrix& X)
420{
421 REPORT // CheckConversion(X);
422 // MatrixConversionCheck mcc;
423 Eq(X,MatrixType::Id);
424}
425
426
427void 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
440void GeneralMatrix::operator<<(const double* r)
441{
442 REPORT
443 int i = storage; Real* s=store;
444 while(i--) *s++ = (Real)*r++;
445}
446
447
448void GeneralMatrix::operator<<(const float* r)
449{
450 REPORT
451 int i = storage; Real* s=store;
452 while(i--) *s++ = (Real)*r++;
453}
454
455
456void GeneralMatrix::operator<<(const int* r)
457{
458 REPORT
459 int i = storage; Real* s=store;
460 while(i--) *s++ = (Real)*r++;
461}
462
463
464void 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
471void 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
492void 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
505void 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
518void 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
531void 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
544void 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
557void 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
566void 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
578void 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
593void 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
608void 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
623void 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
638void 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
653void 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
666void 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
682Real& 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
690Real 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
698Real& 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
707Real 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
716Real& 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
724Real 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
732Real& 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
740Real 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
748Real& 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
756Real 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
764Real& 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
771Real 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
778Real& 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
785Real 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
792Real& 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
799Real 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
806Real& 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
815Real 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
824Real& 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
833Real 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
842Real& 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
851Real 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
860Real& 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
882Real 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///}
Note: See TracBrowser for help on using the repository browser.