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

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

upgrade to newmat11 library
File size: 23.6 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 during Evaluate
498 GeneralMatrix* gm = ((BaseMatrix&)X).Evaluate();
499 AddedMatrix am(this,gm);
500 if (gm==this) Release(2); else Release();
501 Eq2(am,type());
502}
503
504// GeneralMatrix operators
505
506void GeneralMatrix::SP_eq(const BaseMatrix& X)
507{
508 REPORT
509 Tracer tr("GeneralMatrix::SP_eq");
510 // MatrixConversionCheck mcc;
511 Protect(); // so it cannot get deleted during Evaluate
512 GeneralMatrix* gm = ((BaseMatrix&)X).Evaluate();
513 SPMatrix spm(this,gm);
514 if (gm==this) Release(2); else Release();
515 Eq2(spm,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 SubtractedMatrix 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 MultipliedMatrix 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 ConcatenatedMatrix am(this,gm);
553 if (gm==this) Release(2); else Release();
554 Eq2(am,type());
555}
556
557void GeneralMatrix::operator&=(const BaseMatrix& X)
558{
559 REPORT
560 Tracer tr("GeneralMatrix::operator&=");
561 // MatrixConversionCheck mcc;
562 Protect(); // so it cannot get deleted
563 // during Evaluate
564 GeneralMatrix* gm = ((BaseMatrix&)X).Evaluate();
565 StackedMatrix am(this,gm);
566 if (gm==this) Release(2); else Release();
567 Eq2(am,type());
568}
569
570void GeneralMatrix::operator+=(Real r)
571{
572 REPORT
573 Tracer tr("GeneralMatrix::operator+=(Real)");
574 // MatrixConversionCheck mcc;
575 ShiftedMatrix am(this,r);
576 Release(); Eq2(am,type());
577}
578
579void GeneralMatrix::operator*=(Real r)
580{
581 REPORT
582 Tracer tr("GeneralMatrix::operator*=(Real)");
583 // MatrixConversionCheck mcc;
584 ScaledMatrix am(this,r);
585 Release(); Eq2(am,type());
586}
587
588
589// Generic matrix operators
590
591void GenericMatrix::operator+=(const BaseMatrix& X)
592{
593 REPORT
594 Tracer tr("GenericMatrix::operator+=");
595 if (!gm) Throw(ProgramException("GenericMatrix is null"));
596 gm->Protect(); // so it cannot get deleted during Evaluate
597 GeneralMatrix* gmx = ((BaseMatrix&)X).Evaluate();
598 AddedMatrix am(gm,gmx);
599 if (gmx==gm) gm->Release(2); else gm->Release();
600 GeneralMatrix* gmy = am.Evaluate();
601 if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); }
602 else { REPORT }
603 gm->Protect();
604}
605
606void GenericMatrix::SP_eq(const BaseMatrix& X)
607{
608 REPORT
609 Tracer tr("GenericMatrix::SP_eq");
610 if (!gm) Throw(ProgramException("GenericMatrix is null"));
611 gm->Protect(); // so it cannot get deleted during Evaluate
612 GeneralMatrix* gmx = ((BaseMatrix&)X).Evaluate();
613 SPMatrix spm(gm,gmx);
614 if (gmx==gm) gm->Release(2); else gm->Release();
615 GeneralMatrix* gmy = spm.Evaluate();
616 if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); }
617 else { REPORT }
618 gm->Protect();
619}
620
621void GenericMatrix::operator-=(const BaseMatrix& X)
622{
623 REPORT
624 Tracer tr("GenericMatrix::operator-=");
625 if (!gm) Throw(ProgramException("GenericMatrix is null"));
626 gm->Protect(); // so it cannot get deleted during Evaluate
627 GeneralMatrix* gmx = ((BaseMatrix&)X).Evaluate();
628 SubtractedMatrix am(gm,gmx);
629 if (gmx==gm) gm->Release(2); else gm->Release();
630 GeneralMatrix* gmy = am.Evaluate();
631 if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); }
632 else { REPORT }
633 gm->Protect();
634}
635
636void GenericMatrix::operator*=(const BaseMatrix& X)
637{
638 REPORT
639 Tracer tr("GenericMatrix::operator*=");
640 if (!gm) Throw(ProgramException("GenericMatrix is null"));
641 gm->Protect(); // so it cannot get deleted during Evaluate
642 GeneralMatrix* gmx = ((BaseMatrix&)X).Evaluate();
643 MultipliedMatrix am(gm,gmx);
644 if (gmx==gm) gm->Release(2); else gm->Release();
645 GeneralMatrix* gmy = am.Evaluate();
646 if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); }
647 else { REPORT }
648 gm->Protect();
649}
650
651void GenericMatrix::operator|=(const BaseMatrix& X)
652{
653 REPORT
654 Tracer tr("GenericMatrix::operator|=");
655 if (!gm) Throw(ProgramException("GenericMatrix is null"));
656 gm->Protect(); // so it cannot get deleted during Evaluate
657 GeneralMatrix* gmx = ((BaseMatrix&)X).Evaluate();
658 ConcatenatedMatrix am(gm,gmx);
659 if (gmx==gm) gm->Release(2); else 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&=(const BaseMatrix& X)
667{
668 REPORT
669 Tracer tr("GenericMatrix::operator&=");
670 if (!gm) Throw(ProgramException("GenericMatrix is null"));
671 gm->Protect(); // so it cannot get deleted during Evaluate
672 GeneralMatrix* gmx = ((BaseMatrix&)X).Evaluate();
673 StackedMatrix am(gm,gmx);
674 if (gmx==gm) gm->Release(2); else gm->Release();
675 GeneralMatrix* gmy = am.Evaluate();
676 if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); }
677 else { REPORT }
678 gm->Protect();
679}
680
681void GenericMatrix::operator+=(Real r)
682{
683 REPORT
684 Tracer tr("GenericMatrix::operator+= (Real)");
685 if (!gm) Throw(ProgramException("GenericMatrix is null"));
686 ShiftedMatrix am(gm,r);
687 gm->Release();
688 GeneralMatrix* gmy = am.Evaluate();
689 if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); }
690 else { REPORT }
691 gm->Protect();
692}
693
694void GenericMatrix::operator*=(Real r)
695{
696 REPORT
697 Tracer tr("GenericMatrix::operator*= (Real)");
698 if (!gm) Throw(ProgramException("GenericMatrix is null"));
699 ScaledMatrix am(gm,r);
700 gm->Release();
701 GeneralMatrix* gmy = am.Evaluate();
702 if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); }
703 else { REPORT }
704 gm->Protect();
705}
706
707
708/************************* element access *********************************/
709
710Real& Matrix::element(int m, int n)
711{
712 REPORT
713 if (m<0 || m>= nrows_val || n<0 || n>= ncols_val)
714 Throw(IndexException(m,n,*this,true));
715 return store[m*ncols_val+n];
716}
717
718Real Matrix::element(int m, int n) const
719{
720 REPORT
721 if (m<0 || m>= nrows_val || n<0 || n>= ncols_val)
722 Throw(IndexException(m,n,*this,true));
723 return store[m*ncols_val+n];
724}
725
726Real& SymmetricMatrix::element(int m, int n)
727{
728 REPORT
729 if (m<0 || n<0 || m >= nrows_val || n>=ncols_val)
730 Throw(IndexException(m,n,*this,true));
731 if (m>=n) return store[tristore(m)+n];
732 else return store[tristore(n)+m];
733}
734
735Real SymmetricMatrix::element(int m, int n) const
736{
737 REPORT
738 if (m<0 || n<0 || m >= nrows_val || n>=ncols_val)
739 Throw(IndexException(m,n,*this,true));
740 if (m>=n) return store[tristore(m)+n];
741 else return store[tristore(n)+m];
742}
743
744Real& UpperTriangularMatrix::element(int m, int n)
745{
746 REPORT
747 if (m<0 || n<m || n>=ncols_val)
748 Throw(IndexException(m,n,*this,true));
749 return store[m*ncols_val+n-tristore(m)];
750}
751
752Real UpperTriangularMatrix::element(int m, int n) const
753{
754 REPORT
755 if (m<0 || n<m || n>=ncols_val)
756 Throw(IndexException(m,n,*this,true));
757 return store[m*ncols_val+n-tristore(m)];
758}
759
760Real& LowerTriangularMatrix::element(int m, int n)
761{
762 REPORT
763 if (n<0 || m<n || m>=nrows_val)
764 Throw(IndexException(m,n,*this,true));
765 return store[tristore(m)+n];
766}
767
768Real LowerTriangularMatrix::element(int m, int n) const
769{
770 REPORT
771 if (n<0 || m<n || m>=nrows_val)
772 Throw(IndexException(m,n,*this,true));
773 return store[tristore(m)+n];
774}
775
776Real& DiagonalMatrix::element(int m, int n)
777{
778 REPORT
779 if (n<0 || m!=n || m>=nrows_val || n>=ncols_val)
780 Throw(IndexException(m,n,*this,true));
781 return store[n];
782}
783
784Real DiagonalMatrix::element(int m, int n) const
785{
786 REPORT
787 if (n<0 || m!=n || m>=nrows_val || n>=ncols_val)
788 Throw(IndexException(m,n,*this,true));
789 return store[n];
790}
791
792Real& DiagonalMatrix::element(int m)
793{
794 REPORT
795 if (m<0 || m>=nrows_val) Throw(IndexException(m,*this,true));
796 return store[m];
797}
798
799Real DiagonalMatrix::element(int m) const
800{
801 REPORT
802 if (m<0 || m>=nrows_val) Throw(IndexException(m,*this,true));
803 return store[m];
804}
805
806Real& ColumnVector::element(int m)
807{
808 REPORT
809 if (m<0 || m>= nrows_val) Throw(IndexException(m,*this,true));
810 return store[m];
811}
812
813Real ColumnVector::element(int m) const
814{
815 REPORT
816 if (m<0 || m>= nrows_val) Throw(IndexException(m,*this,true));
817 return store[m];
818}
819
820Real& RowVector::element(int n)
821{
822 REPORT
823 if (n<0 || n>= ncols_val) Throw(IndexException(n,*this,true));
824 return store[n];
825}
826
827Real RowVector::element(int n) const
828{
829 REPORT
830 if (n<0 || n>= ncols_val) Throw(IndexException(n,*this,true));
831 return store[n];
832}
833
834Real& BandMatrix::element(int m, int n)
835{
836 REPORT
837 int w = upper_val+lower_val+1; int i = lower_val+n-m;
838 if (m<0 || m>= nrows_val || n<0 || n>= ncols_val || i<0 || i>=w)
839 Throw(IndexException(m,n,*this,true));
840 return store[w*m+i];
841}
842
843Real BandMatrix::element(int m, int n) const
844{
845 REPORT
846 int w = upper_val+lower_val+1; int i = lower_val+n-m;
847 if (m<0 || m>= nrows_val || n<0 || n>= ncols_val || i<0 || i>=w)
848 Throw(IndexException(m,n,*this,true));
849 return store[w*m+i];
850}
851
852Real& UpperBandMatrix::element(int m, int n)
853{
854 REPORT
855 int w = upper_val+1; int i = n-m;
856 if (m<0 || m>= nrows_val || n<0 || n>= ncols_val || i<0 || i>=w)
857 Throw(IndexException(m,n,*this,true));
858 return store[w*m+i];
859}
860
861Real UpperBandMatrix::element(int m, int n) const
862{
863 REPORT
864 int w = upper_val+1; int i = n-m;
865 if (m<0 || m>= nrows_val || n<0 || n>= ncols_val || i<0 || i>=w)
866 Throw(IndexException(m,n,*this,true));
867 return store[w*m+i];
868}
869
870Real& LowerBandMatrix::element(int m, int n)
871{
872 REPORT
873 int w = lower_val+1; int i = lower_val+n-m;
874 if (m<0 || m>= nrows_val || n<0 || n>= ncols_val || i<0 || i>=w)
875 Throw(IndexException(m,n,*this,true));
876 return store[w*m+i];
877}
878
879Real LowerBandMatrix::element(int m, int n) const
880{
881 REPORT
882 int w = lower_val+1; int i = lower_val+n-m;
883 if (m<0 || m>= nrows_val || n<0 || n>= ncols_val || i<0 || i>=w)
884 Throw(IndexException(m,n,*this,true));
885 return store[w*m+i];
886}
887
888Real& SymmetricBandMatrix::element(int m, int n)
889{
890 REPORT
891 int w = lower_val+1;
892 if (m>=n)
893 {
894 REPORT
895 int i = lower_val+n-m;
896 if ( m>=nrows_val || n<0 || i<0 )
897 Throw(IndexException(m,n,*this,true));
898 return store[w*m+i];
899 }
900 else
901 {
902 REPORT
903 int i = lower_val+m-n;
904 if ( n>=nrows_val || m<0 || i<0 )
905 Throw(IndexException(m,n,*this,true));
906 return store[w*n+i];
907 }
908}
909
910Real SymmetricBandMatrix::element(int m, int n) const
911{
912 REPORT
913 int w = lower_val+1;
914 if (m>=n)
915 {
916 REPORT
917 int i = lower_val+n-m;
918 if ( m>=nrows_val || n<0 || i<0 )
919 Throw(IndexException(m,n,*this,true));
920 return store[w*m+i];
921 }
922 else
923 {
924 REPORT
925 int i = lower_val+m-n;
926 if ( n>=nrows_val || m<0 || i<0 )
927 Throw(IndexException(m,n,*this,true));
928 return store[w*n+i];
929 }
930}
931
932#ifdef use_namespace
933}
934#endif
935
936
937///}
Note: See TracBrowser for help on using the repository browser.