source: ntrip/trunk/BNC/newmat/newmat.h@ 6698

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

* empty log message *
File size: 84.1 KB
Line 
1/// \defgroup newmat Newmat matrix manipulation library
2///@{
3
4/// \file newmat.h
5/// Definition file for matrix library.
6
7// Copyright (C) 2004: R B Davies
8
9#ifndef NEWMAT_LIB
10#define NEWMAT_LIB 0
11
12#include "include.h"
13
14#include "myexcept.h"
15
16
17#ifdef use_namespace
18namespace NEWMAT { using namespace RBD_COMMON; }
19namespace RBD_LIBRARIES { using namespace NEWMAT; }
20namespace NEWMAT {
21#endif
22
23//#define DO_REPORT // to activate REPORT
24
25#ifdef NO_LONG_NAMES
26#define UpperTriangularMatrix UTMatrix
27#define LowerTriangularMatrix LTMatrix
28#define SymmetricMatrix SMatrix
29#define DiagonalMatrix DMatrix
30#define BandMatrix BMatrix
31#define UpperBandMatrix UBMatrix
32#define LowerBandMatrix LBMatrix
33#define SymmetricBandMatrix SBMatrix
34#define BandLUMatrix BLUMatrix
35#endif
36
37// ************************** general utilities ****************************/
38
39class GeneralMatrix; // defined later
40class BaseMatrix; // defined later
41class MatrixInput; // defined later
42
43void MatrixErrorNoSpace(const void*); ///< test for allocation fails
44
45/// Return from LogDeterminant function.
46/// Members are the log of the absolute value and the sign (+1, -1 or 0)
47class LogAndSign
48{
49 Real log_val;
50 int sign_val;
51public:
52 LogAndSign() { log_val=0.0; sign_val=1; }
53 LogAndSign(Real);
54 void operator*=(Real); ///< multiply by a real
55 void pow_eq(int k); ///< raise to power of k
56 void PowEq(int k) { pow_eq(k); }
57 void ChangeSign() { sign_val = -sign_val; }
58 void change_sign() { sign_val = -sign_val; } ///< change sign
59 Real LogValue() const { return log_val; }
60 Real log_value() const { return log_val; } ///< log of the absolute value
61 int Sign() const { return sign_val; }
62 int sign() const { return sign_val; } ///< sign of the value
63 Real value() const; ///< the value
64 Real Value() const { return value(); }
65 FREE_CHECK(LogAndSign)
66};
67
68// the following class is for counting the number of times a piece of code
69// is executed. It is used for locating any code not executed by test
70// routines. Use turbo GREP locate all places this code is called and
71// check which ones are not accessed.
72// Somewhat implementation dependent as it relies on "cout" still being
73// present when ExeCounter objects are destructed.
74
75#ifdef DO_REPORT
76
77class ExeCounter
78{
79 int line; // code line number
80 int fileid; // file identifier
81 long nexe; // number of executions
82 static int nreports; // number of reports
83public:
84 ExeCounter(int,int);
85 void operator++() { nexe++; }
86 ~ExeCounter(); // prints out reports
87};
88
89#endif
90
91
92// ************************** class MatrixType *****************************
93
94/// Find the type of a matrix resulting from matrix operations.
95/// Also identify what conversions are permissible.
96/// This class must be updated when new matrix types are added.
97
98class MatrixType
99{
100public:
101 enum Attribute { Valid = 1,
102 Diagonal = 2, // order of these is important
103 Symmetric = 4,
104 Band = 8,
105 Lower = 16,
106 Upper = 32,
107 Square = 64,
108 Skew = 128,
109 LUDeco = 256,
110 Ones = 512 };
111
112 enum { US = 0,
113 UT = Valid + Upper + Square,
114 LT = Valid + Lower + Square,
115 Rt = Valid,
116 Sq = Valid + Square,
117 Sm = Valid + Symmetric + Square,
118 Sk = Valid + Skew + Square,
119 Dg = Valid + Diagonal + Band + Lower + Upper + Symmetric
120 + Square,
121 Id = Valid + Diagonal + Band + Lower + Upper + Symmetric
122 + Square + Ones,
123 RV = Valid, // do not separate out
124 CV = Valid, // vectors
125 BM = Valid + Band + Square,
126 UB = Valid + Band + Upper + Square,
127 LB = Valid + Band + Lower + Square,
128 SB = Valid + Band + Symmetric + Square,
129 KB = Valid + Band + Skew + Square,
130 Ct = Valid + LUDeco + Square,
131 BC = Valid + Band + LUDeco + Square,
132 Mask = ~Square
133 };
134
135
136 static int nTypes() { return 13; } // number of different types
137 // exclude Ct, US, BC
138public:
139 int attribute;
140 bool DataLossOK; // true if data loss is OK when
141 // this represents a destination
142public:
143 MatrixType () : DataLossOK(false) {}
144 MatrixType (int i) : attribute(i), DataLossOK(false) {}
145 MatrixType (int i, bool dlok) : attribute(i), DataLossOK(dlok) {}
146 MatrixType (const MatrixType& mt)
147 : attribute(mt.attribute), DataLossOK(mt.DataLossOK) {}
148 void operator=(const MatrixType& mt)
149 { attribute = mt.attribute; DataLossOK = mt.DataLossOK; }
150 void SetDataLossOK() { DataLossOK = true; }
151 int operator+() const { return attribute; }
152 MatrixType operator+(MatrixType mt) const
153 { return MatrixType(attribute & mt.attribute); }
154 MatrixType operator*(const MatrixType&) const;
155 MatrixType SP(const MatrixType&) const;
156 MatrixType KP(const MatrixType&) const;
157 MatrixType operator|(const MatrixType& mt) const
158 { return MatrixType(attribute & mt.attribute & Valid); }
159 MatrixType operator&(const MatrixType& mt) const
160 { return MatrixType(attribute & mt.attribute & Valid); }
161 bool operator>=(MatrixType mt) const
162 { return ( attribute & ~mt.attribute & Mask ) == 0; }
163 bool operator<(MatrixType mt) const // for MS Visual C++ 4
164 { return ( attribute & ~mt.attribute & Mask ) != 0; }
165 bool operator==(MatrixType t) const
166 { return (attribute == t.attribute); }
167 bool operator!=(MatrixType t) const
168 { return (attribute != t.attribute); }
169 bool operator!() const { return (attribute & Valid) == 0; }
170 MatrixType i() const; ///< type of inverse
171 MatrixType t() const; ///< type of transpose
172 MatrixType AddEqualEl() const ///< add constant to matrix
173 { return MatrixType(attribute & (Valid + Symmetric + Square)); }
174 MatrixType MultRHS() const; ///< type for rhs of multiply
175 MatrixType sub() const ///< type of submatrix
176 { return MatrixType(attribute & Valid); }
177 MatrixType ssub() const ///< type of sym submatrix
178 { return MatrixType(attribute); } // not for selection matrix
179 GeneralMatrix* New() const; ///< new matrix of given type
180 GeneralMatrix* New(int,int,BaseMatrix*) const;
181 ///< new matrix of given type
182 const char* value() const; ///< type as char string
183 const char* Value() const { return value(); }
184 friend bool Rectangular(MatrixType a, MatrixType b, MatrixType c);
185 friend bool Compare(const MatrixType&, MatrixType&);
186 ///< compare and check conversion
187 bool is_band() const { return (attribute & Band) != 0; }
188 bool is_diagonal() const { return (attribute & Diagonal) != 0; }
189 bool is_symmetric() const { return (attribute & Symmetric) != 0; }
190 bool CannotConvert() const { return (attribute & LUDeco) != 0; }
191 // used by operator==
192 FREE_CHECK(MatrixType)
193};
194
195
196// *********************** class MatrixBandWidth ***********************/
197
198///Upper and lower bandwidths of a matrix.
199///That is number of diagonals strictly above or below main diagonal,
200///e.g. diagonal matrix has 0 upper and lower bandwiths.
201///-1 means the matrix may have the maximum bandwidth.
202class MatrixBandWidth
203{
204public:
205 int lower_val;
206 int upper_val;
207 MatrixBandWidth(const int l, const int u) : lower_val(l), upper_val(u) {}
208 MatrixBandWidth(const int i) : lower_val(i), upper_val(i) {}
209 MatrixBandWidth operator+(const MatrixBandWidth&) const;
210 MatrixBandWidth operator*(const MatrixBandWidth&) const;
211 MatrixBandWidth minimum(const MatrixBandWidth&) const;
212 MatrixBandWidth t() const { return MatrixBandWidth(upper_val,lower_val); }
213 bool operator==(const MatrixBandWidth& bw) const
214 { return (lower_val == bw.lower_val) && (upper_val == bw.upper_val); }
215 bool operator!=(const MatrixBandWidth& bw) const { return !operator==(bw); }
216 int Upper() const { return upper_val; }
217 int upper() const { return upper_val; }
218 int Lower() const { return lower_val; }
219 int lower() const { return lower_val; }
220 FREE_CHECK(MatrixBandWidth)
221};
222
223
224// ********************* Array length specifier ************************/
225
226/// This class is used to avoid constructors such as
227/// ColumnVector(int) being used for conversions.
228/// Eventually this should be replaced by the use of the keyword "explicit".
229
230class ArrayLengthSpecifier
231{
232 int v;
233public:
234 int Value() const { return v; }
235 int value() const { return v; }
236 ArrayLengthSpecifier(int l) : v(l) {}
237};
238
239// ************************* Matrix routines ***************************/
240
241
242class MatrixRowCol; // defined later
243class MatrixRow;
244class MatrixCol;
245class MatrixColX;
246
247class GeneralMatrix; // defined later
248class AddedMatrix;
249class MultipliedMatrix;
250class SubtractedMatrix;
251class SPMatrix;
252class KPMatrix;
253class ConcatenatedMatrix;
254class StackedMatrix;
255class SolvedMatrix;
256class ShiftedMatrix;
257class NegShiftedMatrix;
258class ScaledMatrix;
259class TransposedMatrix;
260class ReversedMatrix;
261class NegatedMatrix;
262class InvertedMatrix;
263class RowedMatrix;
264class ColedMatrix;
265class DiagedMatrix;
266class MatedMatrix;
267class GetSubMatrix;
268class ReturnMatrix;
269class Matrix;
270class SquareMatrix;
271class nricMatrix;
272class RowVector;
273class ColumnVector;
274class SymmetricMatrix;
275class UpperTriangularMatrix;
276class LowerTriangularMatrix;
277class DiagonalMatrix;
278class CroutMatrix;
279class BandMatrix;
280class LowerBandMatrix;
281class UpperBandMatrix;
282class SymmetricBandMatrix;
283class LinearEquationSolver;
284class GenericMatrix;
285
286
287#define MatrixTypeUnSp 0
288//static MatrixType MatrixTypeUnSp(MatrixType::US);
289// // AT&T needs this
290
291/// Base of the matrix classes.
292class BaseMatrix : public Janitor
293{
294protected:
295 virtual int search(const BaseMatrix*) const = 0;
296 // count number of times matrix is referred to
297public:
298 virtual GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp) = 0;
299 // evaluate temporary
300 // for old version of G++
301 // virtual GeneralMatrix* Evaluate(MatrixType mt) = 0;
302 // GeneralMatrix* Evaluate() { return Evaluate(MatrixTypeUnSp); }
303 AddedMatrix operator+(const BaseMatrix&) const; // results of operations
304 MultipliedMatrix operator*(const BaseMatrix&) const;
305 SubtractedMatrix operator-(const BaseMatrix&) const;
306 ConcatenatedMatrix operator|(const BaseMatrix&) const;
307 StackedMatrix operator&(const BaseMatrix&) const;
308 ShiftedMatrix operator+(Real) const;
309 ScaledMatrix operator*(Real) const;
310 ScaledMatrix operator/(Real) const;
311 ShiftedMatrix operator-(Real) const;
312 TransposedMatrix t() const;
313// TransposedMatrix t;
314 NegatedMatrix operator-() const; // change sign of elements
315 ReversedMatrix reverse() const;
316 ReversedMatrix Reverse() const;
317 InvertedMatrix i() const;
318// InvertedMatrix i;
319 RowedMatrix as_row() const;
320 RowedMatrix AsRow() const;
321 ColedMatrix as_column() const;
322 ColedMatrix AsColumn() const;
323 DiagedMatrix as_diagonal() const;
324 DiagedMatrix AsDiagonal() const;
325 MatedMatrix as_matrix(int,int) const;
326 MatedMatrix AsMatrix(int m, int n) const;
327 GetSubMatrix submatrix(int,int,int,int) const;
328 GetSubMatrix SubMatrix(int fr, int lr, int fc, int lc) const;
329 GetSubMatrix sym_submatrix(int,int) const;
330 GetSubMatrix SymSubMatrix(int f, int l) const;
331 GetSubMatrix row(int) const;
332 GetSubMatrix rows(int,int) const;
333 GetSubMatrix column(int) const;
334 GetSubMatrix columns(int,int) const;
335 GetSubMatrix Row(int f) const;
336 GetSubMatrix Rows(int f, int l) const;
337 GetSubMatrix Column(int f) const;
338 GetSubMatrix Columns(int f, int l) const;
339 Real as_scalar() const; // conversion of 1 x 1 matrix
340 Real AsScalar() const;
341 virtual LogAndSign log_determinant() const;
342 LogAndSign LogDeterminant() const { return log_determinant(); }
343 Real determinant() const;
344 Real Determinant() const { return determinant(); }
345 virtual Real sum_square() const;
346 Real SumSquare() const { return sum_square(); }
347 Real norm_Frobenius() const;
348 Real norm_frobenius() const { return norm_Frobenius(); }
349 Real NormFrobenius() const { return norm_Frobenius(); }
350 virtual Real sum_absolute_value() const;
351 Real SumAbsoluteValue() const { return sum_absolute_value(); }
352 virtual Real sum() const;
353 virtual Real Sum() const { return sum(); }
354 virtual Real maximum_absolute_value() const;
355 Real MaximumAbsoluteValue() const { return maximum_absolute_value(); }
356 virtual Real maximum_absolute_value1(int& i) const;
357 Real MaximumAbsoluteValue1(int& i) const
358 { return maximum_absolute_value1(i); }
359 virtual Real maximum_absolute_value2(int& i, int& j) const;
360 Real MaximumAbsoluteValue2(int& i, int& j) const
361 { return maximum_absolute_value2(i,j); }
362 virtual Real minimum_absolute_value() const;
363 Real MinimumAbsoluteValue() const { return minimum_absolute_value(); }
364 virtual Real minimum_absolute_value1(int& i) const;
365 Real MinimumAbsoluteValue1(int& i) const
366 { return minimum_absolute_value1(i); }
367 virtual Real minimum_absolute_value2(int& i, int& j) const;
368 Real MinimumAbsoluteValue2(int& i, int& j) const
369 { return minimum_absolute_value2(i,j); }
370 virtual Real maximum() const;
371 Real Maximum() const { return maximum(); }
372 virtual Real maximum1(int& i) const;
373 Real Maximum1(int& i) const { return maximum1(i); }
374 virtual Real maximum2(int& i, int& j) const;
375 Real Maximum2(int& i, int& j) const { return maximum2(i,j); }
376 virtual Real minimum() const;
377 Real Minimum() const { return minimum(); }
378 virtual Real minimum1(int& i) const;
379 Real Minimum1(int& i) const { return minimum1(i); }
380 virtual Real minimum2(int& i, int& j) const;
381 Real Minimum2(int& i, int& j) const { return minimum2(i,j); }
382 virtual Real trace() const;
383 Real Trace() const { return trace(); }
384 Real norm1() const;
385 Real Norm1() const { return norm1(); }
386 Real norm_infinity() const;
387 Real NormInfinity() const { return norm_infinity(); }
388 virtual MatrixBandWidth bandwidth() const; // bandwidths of band matrix
389 virtual MatrixBandWidth BandWidth() const { return bandwidth(); }
390 void IEQND() const; // called by ineq. ops
391 ReturnMatrix sum_square_columns() const;
392 ReturnMatrix sum_square_rows() const;
393 ReturnMatrix sum_columns() const;
394 ReturnMatrix sum_rows() const;
395 virtual void cleanup() {}
396 void CleanUp() { cleanup(); }
397
398// virtual ReturnMatrix Reverse() const; // reverse order of elements
399//protected:
400// BaseMatrix() : t(this), i(this) {}
401
402 friend class GeneralMatrix;
403 friend class Matrix;
404 friend class SquareMatrix;
405 friend class nricMatrix;
406 friend class RowVector;
407 friend class ColumnVector;
408 friend class SymmetricMatrix;
409 friend class UpperTriangularMatrix;
410 friend class LowerTriangularMatrix;
411 friend class DiagonalMatrix;
412 friend class CroutMatrix;
413 friend class BandMatrix;
414 friend class LowerBandMatrix;
415 friend class UpperBandMatrix;
416 friend class SymmetricBandMatrix;
417 friend class AddedMatrix;
418 friend class MultipliedMatrix;
419 friend class SubtractedMatrix;
420 friend class SPMatrix;
421 friend class KPMatrix;
422 friend class ConcatenatedMatrix;
423 friend class StackedMatrix;
424 friend class SolvedMatrix;
425 friend class ShiftedMatrix;
426 friend class NegShiftedMatrix;
427 friend class ScaledMatrix;
428 friend class TransposedMatrix;
429 friend class ReversedMatrix;
430 friend class NegatedMatrix;
431 friend class InvertedMatrix;
432 friend class RowedMatrix;
433 friend class ColedMatrix;
434 friend class DiagedMatrix;
435 friend class MatedMatrix;
436 friend class GetSubMatrix;
437 friend class ReturnMatrix;
438 friend class LinearEquationSolver;
439 friend class GenericMatrix;
440 NEW_DELETE(BaseMatrix)
441};
442
443
444// ***************************** working classes **************************/
445
446/// The classes for matrices that can contain data are derived from this.
447class GeneralMatrix : public BaseMatrix // declarable matrix types
448{
449 virtual GeneralMatrix* Image() const; // copy of matrix
450protected:
451 int tag_val; // shows whether can reuse
452 int nrows_val, ncols_val; // dimensions
453 int storage; // total store required
454 Real* store; // point to store (0=not set)
455 GeneralMatrix(); // initialise with no store
456 GeneralMatrix(ArrayLengthSpecifier); // constructor getting store
457 void Add(GeneralMatrix*, Real); // sum of GM and Real
458 void Add(Real); // add Real to this
459 void NegAdd(GeneralMatrix*, Real); // Real - GM
460 void NegAdd(Real); // this = this - Real
461 void Multiply(GeneralMatrix*, Real); // product of GM and Real
462 void Multiply(Real); // multiply this by Real
463 void Negate(GeneralMatrix*); // change sign
464 void Negate(); // change sign
465 void ReverseElements(); // internal reverse of elements
466 void ReverseElements(GeneralMatrix*); // reverse order of elements
467 void operator=(Real); // set matrix to constant
468 Real* GetStore(); // get store or copy
469 GeneralMatrix* BorrowStore(GeneralMatrix*, MatrixType);
470 // temporarily access store
471 void GetMatrix(const GeneralMatrix*); // used by = and initialise
472 void Eq(const BaseMatrix&, MatrixType); // used by =
473 void Eq(const GeneralMatrix&); // version with no conversion
474 void Eq(const BaseMatrix&, MatrixType, bool);// used by <<
475 void Eq2(const BaseMatrix&, MatrixType); // cut down version of Eq
476 int search(const BaseMatrix*) const;
477 virtual GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
478 void CheckConversion(const BaseMatrix&); // check conversion OK
479 void resize(int, int, int); // change dimensions
480 virtual short SimpleAddOK(const GeneralMatrix*) { return 0; }
481 // see bandmat.cpp for explanation
482 virtual void MiniCleanUp()
483 { store = 0; storage = 0; nrows_val = 0; ncols_val = 0; tag_val = -1;}
484 // CleanUp when the data array has already been deleted
485 void PlusEqual(const GeneralMatrix& gm);
486 void MinusEqual(const GeneralMatrix& gm);
487 void PlusEqual(Real f);
488 void MinusEqual(Real f);
489 void swap(GeneralMatrix& gm); // swap values
490public:
491 GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
492 virtual MatrixType type() const = 0; // type of a matrix
493 MatrixType Type() const { return type(); }
494 int Nrows() const { return nrows_val; } // get dimensions
495 int Ncols() const { return ncols_val; }
496 int Storage() const { return storage; }
497 Real* Store() const { return store; }
498 // updated names
499 int nrows() const { return nrows_val; } // get dimensions
500 int ncols() const { return ncols_val; }
501 int size() const { return storage; }
502 Real* data() { return store; }
503 const Real* data() const { return store; }
504 const Real* const_data() const { return store; }
505 virtual ~GeneralMatrix(); // delete store if set
506 void tDelete(); // delete if tag_val permits
507 bool reuse(); // true if tag_val allows reuse
508 void protect() { tag_val=-1; } // cannot delete or reuse
509 void Protect() { tag_val=-1; } // cannot delete or reuse
510 int tag() const { return tag_val; }
511 int Tag() const { return tag_val; }
512 bool is_zero() const; // test matrix has all zeros
513 bool IsZero() const { return is_zero(); } // test matrix has all zeros
514 void Release() { tag_val=1; } // del store after next use
515 void Release(int t) { tag_val=t; } // del store after t accesses
516 void ReleaseAndDelete() { tag_val=0; } // delete matrix after use
517 void release() { tag_val=1; } // del store after next use
518 void release(int t) { tag_val=t; } // del store after t accesses
519 void release_and_delete() { tag_val=0; } // delete matrix after use
520 void operator<<(const double*); // assignment from an array
521 void operator<<(const float*); // assignment from an array
522 void operator<<(const int*); // assignment from an array
523 void operator<<(const BaseMatrix& X)
524 { Eq(X,this->type(),true); } // = without checking type
525 void inject(const GeneralMatrix&); // copy stored els only
526 void Inject(const GeneralMatrix& GM) { inject(GM); }
527 void operator+=(const BaseMatrix&);
528 void operator-=(const BaseMatrix&);
529 void operator*=(const BaseMatrix&);
530 void operator|=(const BaseMatrix&);
531 void operator&=(const BaseMatrix&);
532 void operator+=(Real);
533 void operator-=(Real r) { operator+=(-r); }
534 void operator*=(Real);
535 void operator/=(Real r) { operator*=(1.0/r); }
536 virtual GeneralMatrix* MakeSolver(); // for solving
537 virtual void Solver(MatrixColX&, const MatrixColX&) {}
538 virtual void GetRow(MatrixRowCol&) = 0; // Get matrix row
539 virtual void RestoreRow(MatrixRowCol&) {} // Restore matrix row
540 virtual void NextRow(MatrixRowCol&); // Go to next row
541 virtual void GetCol(MatrixRowCol&) = 0; // Get matrix col
542 virtual void GetCol(MatrixColX&) = 0; // Get matrix col
543 virtual void RestoreCol(MatrixRowCol&) {} // Restore matrix col
544 virtual void RestoreCol(MatrixColX&) {} // Restore matrix col
545 virtual void NextCol(MatrixRowCol&); // Go to next col
546 virtual void NextCol(MatrixColX&); // Go to next col
547 Real sum_square() const;
548 Real sum_absolute_value() const;
549 Real sum() const;
550 Real maximum_absolute_value1(int& i) const;
551 Real minimum_absolute_value1(int& i) const;
552 Real maximum1(int& i) const;
553 Real minimum1(int& i) const;
554 Real maximum_absolute_value() const;
555 Real maximum_absolute_value2(int& i, int& j) const;
556 Real minimum_absolute_value() const;
557 Real minimum_absolute_value2(int& i, int& j) const;
558 Real maximum() const;
559 Real maximum2(int& i, int& j) const;
560 Real minimum() const;
561 Real minimum2(int& i, int& j) const;
562 LogAndSign log_determinant() const;
563 virtual bool IsEqual(const GeneralMatrix&) const;
564 // same type, same values
565 void CheckStore() const; // check store is non-zero
566 virtual void SetParameters(const GeneralMatrix*) {}
567 // set parameters in GetMatrix
568 operator ReturnMatrix() const; // for building a ReturnMatrix
569 ReturnMatrix for_return() const;
570 ReturnMatrix ForReturn() const;
571 //virtual bool SameStorageType(const GeneralMatrix& A) const;
572 //virtual void ReSizeForAdd(const GeneralMatrix& A, const GeneralMatrix& B);
573 //virtual void ReSizeForSP(const GeneralMatrix& A, const GeneralMatrix& B);
574 virtual void resize(const GeneralMatrix& A);
575 virtual void ReSize(const GeneralMatrix& A) { resize(A); }
576 MatrixInput operator<<(double); // for loading a list
577 MatrixInput operator<<(float); // for loading a list
578 MatrixInput operator<<(int f);
579// ReturnMatrix Reverse() const; // reverse order of elements
580 void cleanup(); // to clear store
581
582 friend class Matrix;
583 friend class SquareMatrix;
584 friend class nricMatrix;
585 friend class SymmetricMatrix;
586 friend class UpperTriangularMatrix;
587 friend class LowerTriangularMatrix;
588 friend class DiagonalMatrix;
589 friend class CroutMatrix;
590 friend class RowVector;
591 friend class ColumnVector;
592 friend class BandMatrix;
593 friend class LowerBandMatrix;
594 friend class UpperBandMatrix;
595 friend class SymmetricBandMatrix;
596 friend class BaseMatrix;
597 friend class AddedMatrix;
598 friend class MultipliedMatrix;
599 friend class SubtractedMatrix;
600 friend class SPMatrix;
601 friend class KPMatrix;
602 friend class ConcatenatedMatrix;
603 friend class StackedMatrix;
604 friend class SolvedMatrix;
605 friend class ShiftedMatrix;
606 friend class NegShiftedMatrix;
607 friend class ScaledMatrix;
608 friend class TransposedMatrix;
609 friend class ReversedMatrix;
610 friend class NegatedMatrix;
611 friend class InvertedMatrix;
612 friend class RowedMatrix;
613 friend class ColedMatrix;
614 friend class DiagedMatrix;
615 friend class MatedMatrix;
616 friend class GetSubMatrix;
617 friend class ReturnMatrix;
618 friend class LinearEquationSolver;
619 friend class GenericMatrix;
620 NEW_DELETE(GeneralMatrix)
621};
622
623
624/// The usual rectangular matrix.
625class Matrix : public GeneralMatrix
626{
627 GeneralMatrix* Image() const; // copy of matrix
628public:
629 Matrix() {}
630 ~Matrix() {}
631 Matrix(int, int); // standard declaration
632 Matrix(const BaseMatrix&); // evaluate BaseMatrix
633 void operator=(const BaseMatrix&);
634 void operator=(Real f) { GeneralMatrix::operator=(f); }
635 void operator=(const Matrix& m) { Eq(m); }
636 MatrixType type() const;
637 Real& operator()(int, int); // access element
638 Real& element(int, int); // access element
639 Real operator()(int, int) const; // access element
640 Real element(int, int) const; // access element
641#ifdef SETUP_C_SUBSCRIPTS
642 Real* operator[](int m) { return store+m*ncols_val; }
643 const Real* operator[](int m) const { return store+m*ncols_val; }
644 // following for Numerical Recipes in C++
645 Matrix(Real, int, int);
646 Matrix(const Real*, int, int);
647#endif
648 Matrix(const Matrix& gm) : GeneralMatrix() { GetMatrix(&gm); }
649 GeneralMatrix* MakeSolver();
650 Real trace() const;
651 void GetRow(MatrixRowCol&);
652 void GetCol(MatrixRowCol&);
653 void GetCol(MatrixColX&);
654 void RestoreCol(MatrixRowCol&);
655 void RestoreCol(MatrixColX&);
656 void NextRow(MatrixRowCol&);
657 void NextCol(MatrixRowCol&);
658 void NextCol(MatrixColX&);
659 virtual void resize(int,int); // change dimensions
660 // virtual so we will catch it being used in a vector called as a matrix
661 virtual void resize_keep(int,int);
662 virtual void ReSize(int m, int n) { resize(m, n); }
663 void resize(const GeneralMatrix& A);
664 void ReSize(const GeneralMatrix& A) { resize(A); }
665 Real maximum_absolute_value2(int& i, int& j) const;
666 Real minimum_absolute_value2(int& i, int& j) const;
667 Real maximum2(int& i, int& j) const;
668 Real minimum2(int& i, int& j) const;
669 void operator+=(const Matrix& M) { PlusEqual(M); }
670 void operator-=(const Matrix& M) { MinusEqual(M); }
671 void operator+=(Real f) { GeneralMatrix::Add(f); }
672 void operator-=(Real f) { GeneralMatrix::Add(-f); }
673 void swap(Matrix& gm) { GeneralMatrix::swap((GeneralMatrix&)gm); }
674 friend Real dotproduct(const Matrix& A, const Matrix& B);
675 NEW_DELETE(Matrix)
676};
677
678/// Square matrix.
679class SquareMatrix : public Matrix
680{
681 GeneralMatrix* Image() const; // copy of matrix
682public:
683 SquareMatrix() {}
684 ~SquareMatrix() {}
685 SquareMatrix(ArrayLengthSpecifier); // standard declaration
686 SquareMatrix(const BaseMatrix&); // evaluate BaseMatrix
687 void operator=(const BaseMatrix&);
688 void operator=(Real f) { GeneralMatrix::operator=(f); }
689 void operator=(const SquareMatrix& m) { Eq(m); }
690 void operator=(const Matrix& m);
691 MatrixType type() const;
692 SquareMatrix(const SquareMatrix& gm) : Matrix() { GetMatrix(&gm); }
693 SquareMatrix(const Matrix& gm);
694 void resize(int); // change dimensions
695 void ReSize(int m) { resize(m); }
696 void resize_keep(int);
697 void resize_keep(int,int);
698 void resize(int,int); // change dimensions
699 void ReSize(int m, int n) { resize(m, n); }
700 void resize(const GeneralMatrix& A);
701 void ReSize(const GeneralMatrix& A) { resize(A); }
702 void operator+=(const Matrix& M) { PlusEqual(M); }
703 void operator-=(const Matrix& M) { MinusEqual(M); }
704 void operator+=(Real f) { GeneralMatrix::Add(f); }
705 void operator-=(Real f) { GeneralMatrix::Add(-f); }
706 void swap(SquareMatrix& gm) { GeneralMatrix::swap((GeneralMatrix&)gm); }
707 NEW_DELETE(SquareMatrix)
708};
709
710/// Rectangular matrix for use with Numerical Recipes in C.
711class nricMatrix : public Matrix
712{
713 GeneralMatrix* Image() const; // copy of matrix
714 Real** row_pointer; // points to rows
715 void MakeRowPointer(); // build rowpointer
716 void DeleteRowPointer();
717public:
718 nricMatrix() {}
719 nricMatrix(int m, int n) // standard declaration
720 : Matrix(m,n) { MakeRowPointer(); }
721 nricMatrix(const BaseMatrix& bm) // evaluate BaseMatrix
722 : Matrix(bm) { MakeRowPointer(); }
723 void operator=(const BaseMatrix& bm)
724 { DeleteRowPointer(); Matrix::operator=(bm); MakeRowPointer(); }
725 void operator=(Real f) { GeneralMatrix::operator=(f); }
726 void operator=(const nricMatrix& m)
727 { DeleteRowPointer(); Eq(m); MakeRowPointer(); }
728 void operator<<(const BaseMatrix& X)
729 { DeleteRowPointer(); Eq(X,this->type(),true); MakeRowPointer(); }
730 nricMatrix(const nricMatrix& gm) : Matrix()
731 { GetMatrix(&gm); MakeRowPointer(); }
732 void resize(int m, int n) // change dimensions
733 { DeleteRowPointer(); Matrix::resize(m,n); MakeRowPointer(); }
734 void resize_keep(int m, int n) // change dimensions
735 { DeleteRowPointer(); Matrix::resize_keep(m,n); MakeRowPointer(); }
736 void ReSize(int m, int n) // change dimensions
737 { DeleteRowPointer(); Matrix::resize(m,n); MakeRowPointer(); }
738 void resize(const GeneralMatrix& A);
739 void ReSize(const GeneralMatrix& A) { resize(A); }
740 ~nricMatrix() { DeleteRowPointer(); }
741 Real** nric() const { CheckStore(); return row_pointer-1; }
742 void cleanup(); // to clear store
743 void MiniCleanUp();
744 void operator+=(const Matrix& M) { PlusEqual(M); }
745 void operator-=(const Matrix& M) { MinusEqual(M); }
746 void operator+=(Real f) { GeneralMatrix::Add(f); }
747 void operator-=(Real f) { GeneralMatrix::Add(-f); }
748 void swap(nricMatrix& gm);
749 NEW_DELETE(nricMatrix)
750};
751
752/// Symmetric matrix.
753class SymmetricMatrix : public GeneralMatrix
754{
755 GeneralMatrix* Image() const; // copy of matrix
756public:
757 SymmetricMatrix() {}
758 ~SymmetricMatrix() {}
759 SymmetricMatrix(ArrayLengthSpecifier);
760 SymmetricMatrix(const BaseMatrix&);
761 void operator=(const BaseMatrix&);
762 void operator=(Real f) { GeneralMatrix::operator=(f); }
763 void operator=(const SymmetricMatrix& m) { Eq(m); }
764 Real& operator()(int, int); // access element
765 Real& element(int, int); // access element
766 Real operator()(int, int) const; // access element
767 Real element(int, int) const; // access element
768#ifdef SETUP_C_SUBSCRIPTS
769 Real* operator[](int m) { return store+(m*(m+1))/2; }
770 const Real* operator[](int m) const { return store+(m*(m+1))/2; }
771#endif
772 MatrixType type() const;
773 SymmetricMatrix(const SymmetricMatrix& gm)
774 : GeneralMatrix() { GetMatrix(&gm); }
775 Real sum_square() const;
776 Real sum_absolute_value() const;
777 Real sum() const;
778 Real trace() const;
779 void GetRow(MatrixRowCol&);
780 void GetCol(MatrixRowCol&);
781 void GetCol(MatrixColX&);
782 void RestoreCol(MatrixRowCol&) {}
783 void RestoreCol(MatrixColX&);
784 GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
785 void resize(int); // change dimensions
786 void ReSize(int m) { resize(m); }
787 void resize_keep(int);
788 void resize(const GeneralMatrix& A);
789 void ReSize(const GeneralMatrix& A) { resize(A); }
790 void operator+=(const SymmetricMatrix& M) { PlusEqual(M); }
791 void operator-=(const SymmetricMatrix& M) { MinusEqual(M); }
792 void operator+=(Real f) { GeneralMatrix::Add(f); }
793 void operator-=(Real f) { GeneralMatrix::Add(-f); }
794 void swap(SymmetricMatrix& gm) { GeneralMatrix::swap((GeneralMatrix&)gm); }
795 NEW_DELETE(SymmetricMatrix)
796};
797
798/// Upper triangular matrix.
799class UpperTriangularMatrix : public GeneralMatrix
800{
801 GeneralMatrix* Image() const; // copy of matrix
802public:
803 UpperTriangularMatrix() {}
804 ~UpperTriangularMatrix() {}
805 UpperTriangularMatrix(ArrayLengthSpecifier);
806 void operator=(const BaseMatrix&);
807 void operator=(const UpperTriangularMatrix& m) { Eq(m); }
808 UpperTriangularMatrix(const BaseMatrix&);
809 UpperTriangularMatrix(const UpperTriangularMatrix& gm)
810 : GeneralMatrix() { GetMatrix(&gm); }
811 void operator=(Real f) { GeneralMatrix::operator=(f); }
812 Real& operator()(int, int); // access element
813 Real& element(int, int); // access element
814 Real operator()(int, int) const; // access element
815 Real element(int, int) const; // access element
816#ifdef SETUP_C_SUBSCRIPTS
817 Real* operator[](int m) { return store+m*ncols_val-(m*(m+1))/2; }
818 const Real* operator[](int m) const
819 { return store+m*ncols_val-(m*(m+1))/2; }
820#endif
821 MatrixType type() const;
822 GeneralMatrix* MakeSolver() { return this; } // for solving
823 void Solver(MatrixColX&, const MatrixColX&);
824 LogAndSign log_determinant() const;
825 Real trace() const;
826 void GetRow(MatrixRowCol&);
827 void GetCol(MatrixRowCol&);
828 void GetCol(MatrixColX&);
829 void RestoreCol(MatrixRowCol&);
830 void RestoreCol(MatrixColX& c) { RestoreCol((MatrixRowCol&)c); }
831 void NextRow(MatrixRowCol&);
832 void resize(int); // change dimensions
833 void ReSize(int m) { resize(m); }
834 void resize(const GeneralMatrix& A);
835 void ReSize(const GeneralMatrix& A) { resize(A); }
836 void resize_keep(int);
837 MatrixBandWidth bandwidth() const;
838 void operator+=(const UpperTriangularMatrix& M) { PlusEqual(M); }
839 void operator-=(const UpperTriangularMatrix& M) { MinusEqual(M); }
840 void operator+=(Real f) { GeneralMatrix::operator+=(f); }
841 void operator-=(Real f) { GeneralMatrix::operator-=(f); }
842 void swap(UpperTriangularMatrix& gm)
843 { GeneralMatrix::swap((GeneralMatrix&)gm); }
844 NEW_DELETE(UpperTriangularMatrix)
845};
846
847/// Lower triangular matrix.
848class LowerTriangularMatrix : public GeneralMatrix
849{
850 GeneralMatrix* Image() const; // copy of matrix
851public:
852 LowerTriangularMatrix() {}
853 ~LowerTriangularMatrix() {}
854 LowerTriangularMatrix(ArrayLengthSpecifier);
855 LowerTriangularMatrix(const LowerTriangularMatrix& gm)
856 : GeneralMatrix() { GetMatrix(&gm); }
857 LowerTriangularMatrix(const BaseMatrix& M);
858 void operator=(const BaseMatrix&);
859 void operator=(Real f) { GeneralMatrix::operator=(f); }
860 void operator=(const LowerTriangularMatrix& m) { Eq(m); }
861 Real& operator()(int, int); // access element
862 Real& element(int, int); // access element
863 Real operator()(int, int) const; // access element
864 Real element(int, int) const; // access element
865#ifdef SETUP_C_SUBSCRIPTS
866 Real* operator[](int m) { return store+(m*(m+1))/2; }
867 const Real* operator[](int m) const { return store+(m*(m+1))/2; }
868#endif
869 MatrixType type() const;
870 GeneralMatrix* MakeSolver() { return this; } // for solving
871 void Solver(MatrixColX&, const MatrixColX&);
872 LogAndSign log_determinant() const;
873 Real trace() const;
874 void GetRow(MatrixRowCol&);
875 void GetCol(MatrixRowCol&);
876 void GetCol(MatrixColX&);
877 void RestoreCol(MatrixRowCol&);
878 void RestoreCol(MatrixColX& c) { RestoreCol((MatrixRowCol&)c); }
879 void NextRow(MatrixRowCol&);
880 void resize(int); // change dimensions
881 void ReSize(int m) { resize(m); }
882 void resize_keep(int);
883 void resize(const GeneralMatrix& A);
884 void ReSize(const GeneralMatrix& A) { resize(A); }
885 MatrixBandWidth bandwidth() const;
886 void operator+=(const LowerTriangularMatrix& M) { PlusEqual(M); }
887 void operator-=(const LowerTriangularMatrix& M) { MinusEqual(M); }
888 void operator+=(Real f) { GeneralMatrix::operator+=(f); }
889 void operator-=(Real f) { GeneralMatrix::operator-=(f); }
890 void swap(LowerTriangularMatrix& gm)
891 { GeneralMatrix::swap((GeneralMatrix&)gm); }
892 NEW_DELETE(LowerTriangularMatrix)
893};
894
895/// Diagonal matrix.
896class DiagonalMatrix : public GeneralMatrix
897{
898 GeneralMatrix* Image() const; // copy of matrix
899public:
900 DiagonalMatrix() {}
901 ~DiagonalMatrix() {}
902 DiagonalMatrix(ArrayLengthSpecifier);
903 DiagonalMatrix(const BaseMatrix&);
904 DiagonalMatrix(const DiagonalMatrix& gm)
905 : GeneralMatrix() { GetMatrix(&gm); }
906 void operator=(const BaseMatrix&);
907 void operator=(Real f) { GeneralMatrix::operator=(f); }
908 void operator=(const DiagonalMatrix& m) { Eq(m); }
909 Real& operator()(int, int); // access element
910 Real& operator()(int); // access element
911 Real operator()(int, int) const; // access element
912 Real operator()(int) const;
913 Real& element(int, int); // access element
914 Real& element(int); // access element
915 Real element(int, int) const; // access element
916 Real element(int) const; // access element
917#ifdef SETUP_C_SUBSCRIPTS
918 Real& operator[](int m) { return store[m]; }
919 const Real& operator[](int m) const { return store[m]; }
920#endif
921 MatrixType type() const;
922
923 LogAndSign log_determinant() const;
924 Real trace() const;
925 void GetRow(MatrixRowCol&);
926 void GetCol(MatrixRowCol&);
927 void GetCol(MatrixColX&);
928 void NextRow(MatrixRowCol&);
929 void NextCol(MatrixRowCol&);
930 void NextCol(MatrixColX&);
931 GeneralMatrix* MakeSolver() { return this; } // for solving
932 void Solver(MatrixColX&, const MatrixColX&);
933 GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
934 void resize(int); // change dimensions
935 void ReSize(int m) { resize(m); }
936 void resize_keep(int);
937 void resize(const GeneralMatrix& A);
938 void ReSize(const GeneralMatrix& A) { resize(A); }
939 Real* nric() const
940 { CheckStore(); return store-1; } // for use by NRIC
941 MatrixBandWidth bandwidth() const;
942// ReturnMatrix Reverse() const; // reverse order of elements
943 void operator+=(const DiagonalMatrix& M) { PlusEqual(M); }
944 void operator-=(const DiagonalMatrix& M) { MinusEqual(M); }
945 void operator+=(Real f) { GeneralMatrix::operator+=(f); }
946 void operator-=(Real f) { GeneralMatrix::operator-=(f); }
947 void swap(DiagonalMatrix& gm)
948 { GeneralMatrix::swap((GeneralMatrix&)gm); }
949 NEW_DELETE(DiagonalMatrix)
950};
951
952/// Row vector.
953class RowVector : public Matrix
954{
955 GeneralMatrix* Image() const; // copy of matrix
956public:
957 RowVector() { nrows_val = 1; }
958 ~RowVector() {}
959 RowVector(ArrayLengthSpecifier n) : Matrix(1,n.Value()) {}
960 RowVector(const BaseMatrix&);
961 RowVector(const RowVector& gm) : Matrix() { GetMatrix(&gm); }
962 void operator=(const BaseMatrix&);
963 void operator=(Real f) { GeneralMatrix::operator=(f); }
964 void operator=(const RowVector& m) { Eq(m); }
965 Real& operator()(int); // access element
966 Real& element(int); // access element
967 Real operator()(int) const; // access element
968 Real element(int) const; // access element
969#ifdef SETUP_C_SUBSCRIPTS
970 Real& operator[](int m) { return store[m]; }
971 const Real& operator[](int m) const { return store[m]; }
972 // following for Numerical Recipes in C++
973 RowVector(Real a, int n) : Matrix(a, 1, n) {}
974 RowVector(const Real* a, int n) : Matrix(a, 1, n) {}
975#endif
976 MatrixType type() const;
977 void GetCol(MatrixRowCol&);
978 void GetCol(MatrixColX&);
979 void NextCol(MatrixRowCol&);
980 void NextCol(MatrixColX&);
981 void RestoreCol(MatrixRowCol&) {}
982 void RestoreCol(MatrixColX& c);
983 GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
984 void resize(int); // change dimensions
985 void ReSize(int m) { resize(m); }
986 void resize_keep(int);
987 void resize_keep(int,int);
988 void resize(int,int); // in case access is matrix
989 void ReSize(int m,int n) { resize(m, n); }
990 void resize(const GeneralMatrix& A);
991 void ReSize(const GeneralMatrix& A) { resize(A); }
992 Real* nric() const
993 { CheckStore(); return store-1; } // for use by NRIC
994 void cleanup(); // to clear store
995 void MiniCleanUp()
996 { store = 0; storage = 0; nrows_val = 1; ncols_val = 0; tag_val = -1; }
997 // friend ReturnMatrix GetMatrixRow(Matrix& A, int row);
998 void operator+=(const Matrix& M) { PlusEqual(M); }
999 void operator-=(const Matrix& M) { MinusEqual(M); }
1000 void operator+=(Real f) { GeneralMatrix::Add(f); }
1001 void operator-=(Real f) { GeneralMatrix::Add(-f); }
1002 void swap(RowVector& gm)
1003 { GeneralMatrix::swap((GeneralMatrix&)gm); }
1004 NEW_DELETE(RowVector)
1005};
1006
1007/// Column vector.
1008class ColumnVector : public Matrix
1009{
1010 GeneralMatrix* Image() const; // copy of matrix
1011public:
1012 ColumnVector() { ncols_val = 1; }
1013 ~ColumnVector() {}
1014 ColumnVector(ArrayLengthSpecifier n) : Matrix(n.Value(),1) {}
1015 ColumnVector(const BaseMatrix&);
1016 ColumnVector(const ColumnVector& gm) : Matrix() { GetMatrix(&gm); }
1017 void operator=(const BaseMatrix&);
1018 void operator=(Real f) { GeneralMatrix::operator=(f); }
1019 void operator=(const ColumnVector& m) { Eq(m); }
1020 Real& operator()(int); // access element
1021 Real& element(int); // access element
1022 Real operator()(int) const; // access element
1023 Real element(int) const; // access element
1024#ifdef SETUP_C_SUBSCRIPTS
1025 Real& operator[](int m) { return store[m]; }
1026 const Real& operator[](int m) const { return store[m]; }
1027 // following for Numerical Recipes in C++
1028 ColumnVector(Real a, int m) : Matrix(a, m, 1) {}
1029 ColumnVector(const Real* a, int m) : Matrix(a, m, 1) {}
1030#endif
1031 MatrixType type() const;
1032 GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
1033 void resize(int); // change dimensions
1034 void ReSize(int m) { resize(m); }
1035 void resize_keep(int);
1036 void resize_keep(int,int);
1037 void resize(int,int); // in case access is matrix
1038 void ReSize(int m,int n) { resize(m, n); }
1039 void resize(const GeneralMatrix& A);
1040 void ReSize(const GeneralMatrix& A) { resize(A); }
1041 Real* nric() const
1042 { CheckStore(); return store-1; } // for use by NRIC
1043 void cleanup(); // to clear store
1044 void MiniCleanUp()
1045 { store = 0; storage = 0; nrows_val = 0; ncols_val = 1; tag_val = -1; }
1046// ReturnMatrix Reverse() const; // reverse order of elements
1047 void operator+=(const Matrix& M) { PlusEqual(M); }
1048 void operator-=(const Matrix& M) { MinusEqual(M); }
1049 void operator+=(Real f) { GeneralMatrix::Add(f); }
1050 void operator-=(Real f) { GeneralMatrix::Add(-f); }
1051 void swap(ColumnVector& gm)
1052 { GeneralMatrix::swap((GeneralMatrix&)gm); }
1053 NEW_DELETE(ColumnVector)
1054};
1055
1056/// LU matrix.
1057/// A square matrix decomposed into upper and lower triangular
1058/// in preparation for inverting or solving equations.
1059class CroutMatrix : public GeneralMatrix
1060{
1061 int* indx;
1062 bool d; // number of row swaps = even or odd
1063 bool sing;
1064 void ludcmp();
1065 void get_aux(CroutMatrix&); // for copying indx[] etc
1066 GeneralMatrix* Image() const; // copy of matrix
1067public:
1068 CroutMatrix(const BaseMatrix&);
1069 CroutMatrix() : indx(0), d(true), sing(true) {}
1070 CroutMatrix(const CroutMatrix&);
1071 void operator=(const CroutMatrix&);
1072 GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1073 MatrixType type() const;
1074 void lubksb(Real*, int=0);
1075 ~CroutMatrix();
1076 GeneralMatrix* MakeSolver() { return this; } // for solving
1077 LogAndSign log_determinant() const;
1078 void Solver(MatrixColX&, const MatrixColX&);
1079 void GetRow(MatrixRowCol&);
1080 void GetCol(MatrixRowCol&);
1081 void GetCol(MatrixColX& c) { GetCol((MatrixRowCol&)c); }
1082 void cleanup(); // to clear store
1083 void MiniCleanUp();
1084 bool IsEqual(const GeneralMatrix&) const;
1085 bool is_singular() const { return sing; }
1086 bool IsSingular() const { return sing; }
1087 const int* const_data_indx() const { return indx; }
1088 bool even_exchanges() const { return d; }
1089 void swap(CroutMatrix& gm);
1090 NEW_DELETE(CroutMatrix)
1091};
1092
1093// ***************************** band matrices ***************************/
1094
1095/// Band matrix.
1096class BandMatrix : public GeneralMatrix
1097{
1098 GeneralMatrix* Image() const; // copy of matrix
1099protected:
1100 void CornerClear() const; // set unused elements to zero
1101 short SimpleAddOK(const GeneralMatrix* gm);
1102public:
1103 int lower_val, upper_val; // band widths
1104 BandMatrix() { lower_val=0; upper_val=0; CornerClear(); }
1105 ~BandMatrix() {}
1106 BandMatrix(int n,int lb,int ub) { resize(n,lb,ub); CornerClear(); }
1107 // standard declaration
1108 BandMatrix(const BaseMatrix&); // evaluate BaseMatrix
1109 void operator=(const BaseMatrix&);
1110 void operator=(Real f) { GeneralMatrix::operator=(f); }
1111 void operator=(const BandMatrix& m) { Eq(m); }
1112 MatrixType type() const;
1113 Real& operator()(int, int); // access element
1114 Real& element(int, int); // access element
1115 Real operator()(int, int) const; // access element
1116 Real element(int, int) const; // access element
1117#ifdef SETUP_C_SUBSCRIPTS
1118 Real* operator[](int m) { return store+(upper_val+lower_val)*m+lower_val; }
1119 const Real* operator[](int m) const
1120 { return store+(upper_val+lower_val)*m+lower_val; }
1121#endif
1122 BandMatrix(const BandMatrix& gm) : GeneralMatrix() { GetMatrix(&gm); }
1123 LogAndSign log_determinant() const;
1124 GeneralMatrix* MakeSolver();
1125 Real trace() const;
1126 Real sum_square() const
1127 { CornerClear(); return GeneralMatrix::sum_square(); }
1128 Real sum_absolute_value() const
1129 { CornerClear(); return GeneralMatrix::sum_absolute_value(); }
1130 Real sum() const
1131 { CornerClear(); return GeneralMatrix::sum(); }
1132 Real maximum_absolute_value() const
1133 { CornerClear(); return GeneralMatrix::maximum_absolute_value(); }
1134 Real minimum_absolute_value() const
1135 { int i, j; return GeneralMatrix::minimum_absolute_value2(i, j); }
1136 Real maximum() const { int i, j; return GeneralMatrix::maximum2(i, j); }
1137 Real minimum() const { int i, j; return GeneralMatrix::minimum2(i, j); }
1138 void GetRow(MatrixRowCol&);
1139 void GetCol(MatrixRowCol&);
1140 void GetCol(MatrixColX&);
1141 void RestoreCol(MatrixRowCol&);
1142 void RestoreCol(MatrixColX& c) { RestoreCol((MatrixRowCol&)c); }
1143 void NextRow(MatrixRowCol&);
1144 virtual void resize(int, int, int); // change dimensions
1145 virtual void ReSize(int m, int n, int b) { resize(m, n, b); }
1146 void resize(const GeneralMatrix& A);
1147 void ReSize(const GeneralMatrix& A) { resize(A); }
1148 //bool SameStorageType(const GeneralMatrix& A) const;
1149 //void ReSizeForAdd(const GeneralMatrix& A, const GeneralMatrix& B);
1150 //void ReSizeForSP(const GeneralMatrix& A, const GeneralMatrix& B);
1151 MatrixBandWidth bandwidth() const;
1152 void SetParameters(const GeneralMatrix*);
1153 MatrixInput operator<<(double); // will give error
1154 MatrixInput operator<<(float); // will give error
1155 MatrixInput operator<<(int f);
1156 void operator<<(const double* r); // will give error
1157 void operator<<(const float* r); // will give error
1158 void operator<<(const int* r); // will give error
1159 // the next is included because Zortech and Borland
1160 // cannot find the copy in GeneralMatrix
1161 void operator<<(const BaseMatrix& X) { GeneralMatrix::operator<<(X); }
1162 void swap(BandMatrix& gm);
1163 NEW_DELETE(BandMatrix)
1164};
1165
1166/// Upper triangular band matrix.
1167class UpperBandMatrix : public BandMatrix
1168{
1169 GeneralMatrix* Image() const; // copy of matrix
1170public:
1171 UpperBandMatrix() {}
1172 ~UpperBandMatrix() {}
1173 UpperBandMatrix(int n, int ubw) // standard declaration
1174 : BandMatrix(n, 0, ubw) {}
1175 UpperBandMatrix(const BaseMatrix&); // evaluate BaseMatrix
1176 void operator=(const BaseMatrix&);
1177 void operator=(Real f) { GeneralMatrix::operator=(f); }
1178 void operator=(const UpperBandMatrix& m) { Eq(m); }
1179 MatrixType type() const;
1180 UpperBandMatrix(const UpperBandMatrix& gm) : BandMatrix() { GetMatrix(&gm); }
1181 GeneralMatrix* MakeSolver() { return this; }
1182 void Solver(MatrixColX&, const MatrixColX&);
1183 LogAndSign log_determinant() const;
1184 void resize(int, int, int); // change dimensions
1185 void ReSize(int m, int n, int b) { resize(m, n, b); }
1186 void resize(int n,int ubw) // change dimensions
1187 { BandMatrix::resize(n,0,ubw); }
1188 void ReSize(int n,int ubw) // change dimensions
1189 { BandMatrix::resize(n,0,ubw); }
1190 void resize(const GeneralMatrix& A) { BandMatrix::resize(A); }
1191 void ReSize(const GeneralMatrix& A) { BandMatrix::resize(A); }
1192 Real& operator()(int, int);
1193 Real operator()(int, int) const;
1194 Real& element(int, int);
1195 Real element(int, int) const;
1196#ifdef SETUP_C_SUBSCRIPTS
1197 Real* operator[](int m) { return store+upper_val*m; }
1198 const Real* operator[](int m) const { return store+upper_val*m; }
1199#endif
1200 void swap(UpperBandMatrix& gm)
1201 { BandMatrix::swap((BandMatrix&)gm); }
1202 NEW_DELETE(UpperBandMatrix)
1203};
1204
1205/// Lower triangular band matrix.
1206class LowerBandMatrix : public BandMatrix
1207{
1208 GeneralMatrix* Image() const; // copy of matrix
1209public:
1210 LowerBandMatrix() {}
1211 ~LowerBandMatrix() {}
1212 LowerBandMatrix(int n, int lbw) // standard declaration
1213 : BandMatrix(n, lbw, 0) {}
1214 LowerBandMatrix(const BaseMatrix&); // evaluate BaseMatrix
1215 void operator=(const BaseMatrix&);
1216 void operator=(Real f) { GeneralMatrix::operator=(f); }
1217 void operator=(const LowerBandMatrix& m) { Eq(m); }
1218 MatrixType type() const;
1219 LowerBandMatrix(const LowerBandMatrix& gm) : BandMatrix() { GetMatrix(&gm); }
1220 GeneralMatrix* MakeSolver() { return this; }
1221 void Solver(MatrixColX&, const MatrixColX&);
1222 LogAndSign log_determinant() const;
1223 void resize(int, int, int); // change dimensions
1224 void ReSize(int m, int n, int b) { resize(m, n, b); }
1225 void resize(int n,int lbw) // change dimensions
1226 { BandMatrix::resize(n,lbw,0); }
1227 void ReSize(int n,int lbw) // change dimensions
1228 { BandMatrix::resize(n,lbw,0); }
1229 void resize(const GeneralMatrix& A) { BandMatrix::resize(A); }
1230 void ReSize(const GeneralMatrix& A) { BandMatrix::resize(A); }
1231 Real& operator()(int, int);
1232 Real operator()(int, int) const;
1233 Real& element(int, int);
1234 Real element(int, int) const;
1235#ifdef SETUP_C_SUBSCRIPTS
1236 Real* operator[](int m) { return store+lower_val*(m+1); }
1237 const Real* operator[](int m) const { return store+lower_val*(m+1); }
1238#endif
1239 void swap(LowerBandMatrix& gm)
1240 { BandMatrix::swap((BandMatrix&)gm); }
1241 NEW_DELETE(LowerBandMatrix)
1242};
1243
1244/// Symmetric band matrix.
1245class SymmetricBandMatrix : public GeneralMatrix
1246{
1247 GeneralMatrix* Image() const; // copy of matrix
1248 void CornerClear() const; // set unused elements to zero
1249 short SimpleAddOK(const GeneralMatrix* gm);
1250public:
1251 int lower_val; // lower band width
1252 SymmetricBandMatrix() { lower_val=0; CornerClear(); }
1253 ~SymmetricBandMatrix() {}
1254 SymmetricBandMatrix(int n, int lb) { resize(n,lb); CornerClear(); }
1255 SymmetricBandMatrix(const BaseMatrix&);
1256 void operator=(const BaseMatrix&);
1257 void operator=(Real f) { GeneralMatrix::operator=(f); }
1258 void operator=(const SymmetricBandMatrix& m) { Eq(m); }
1259 Real& operator()(int, int); // access element
1260 Real& element(int, int); // access element
1261 Real operator()(int, int) const; // access element
1262 Real element(int, int) const; // access element
1263#ifdef SETUP_C_SUBSCRIPTS
1264 Real* operator[](int m) { return store+lower_val*(m+1); }
1265 const Real* operator[](int m) const { return store+lower_val*(m+1); }
1266#endif
1267 MatrixType type() const;
1268 SymmetricBandMatrix(const SymmetricBandMatrix& gm)
1269 : GeneralMatrix() { GetMatrix(&gm); }
1270 GeneralMatrix* MakeSolver();
1271 Real sum_square() const;
1272 Real sum_absolute_value() const;
1273 Real sum() const;
1274 Real maximum_absolute_value() const
1275 { CornerClear(); return GeneralMatrix::maximum_absolute_value(); }
1276 Real minimum_absolute_value() const
1277 { int i, j; return GeneralMatrix::minimum_absolute_value2(i, j); }
1278 Real maximum() const { int i, j; return GeneralMatrix::maximum2(i, j); }
1279 Real minimum() const { int i, j; return GeneralMatrix::minimum2(i, j); }
1280 Real trace() const;
1281 LogAndSign log_determinant() const;
1282 void GetRow(MatrixRowCol&);
1283 void GetCol(MatrixRowCol&);
1284 void GetCol(MatrixColX&);
1285 void RestoreCol(MatrixRowCol&) {}
1286 void RestoreCol(MatrixColX&);
1287 GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
1288 void resize(int,int); // change dimensions
1289 void ReSize(int m,int b) { resize(m, b); }
1290 void resize(const GeneralMatrix& A);
1291 void ReSize(const GeneralMatrix& A) { resize(A); }
1292 //bool SameStorageType(const GeneralMatrix& A) const;
1293 //void ReSizeForAdd(const GeneralMatrix& A, const GeneralMatrix& B);
1294 //void ReSizeForSP(const GeneralMatrix& A, const GeneralMatrix& B);
1295 MatrixBandWidth bandwidth() const;
1296 void SetParameters(const GeneralMatrix*);
1297 void operator<<(const double* r); // will give error
1298 void operator<<(const float* r); // will give error
1299 void operator<<(const int* r); // will give error
1300 void operator<<(const BaseMatrix& X) { GeneralMatrix::operator<<(X); }
1301 void swap(SymmetricBandMatrix& gm);
1302 NEW_DELETE(SymmetricBandMatrix)
1303};
1304
1305/// LU decomposition of a band matrix.
1306class BandLUMatrix : public GeneralMatrix
1307{
1308 int* indx;
1309 bool d;
1310 bool sing; // true if singular
1311 Real* store2;
1312 int storage2;
1313 int m1,m2; // lower and upper
1314 void ludcmp();
1315 void get_aux(BandLUMatrix&); // for copying indx[] etc
1316 GeneralMatrix* Image() const; // copy of matrix
1317public:
1318 BandLUMatrix(const BaseMatrix&);
1319 BandLUMatrix()
1320 : indx(0), d(true), sing(true), store2(0), storage2(0), m1(0), m2(0) {}
1321 BandLUMatrix(const BandLUMatrix&);
1322 void operator=(const BandLUMatrix&);
1323 GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1324 MatrixType type() const;
1325 void lubksb(Real*, int=0);
1326 ~BandLUMatrix();
1327 GeneralMatrix* MakeSolver() { return this; } // for solving
1328 LogAndSign log_determinant() const;
1329 void Solver(MatrixColX&, const MatrixColX&);
1330 void GetRow(MatrixRowCol&);
1331 void GetCol(MatrixRowCol&);
1332 void GetCol(MatrixColX& c) { GetCol((MatrixRowCol&)c); }
1333 void cleanup(); // to clear store
1334 void MiniCleanUp();
1335 bool IsEqual(const GeneralMatrix&) const;
1336 bool is_singular() const { return sing; }
1337 bool IsSingular() const { return sing; }
1338 const Real* const_data2() const { return store2; }
1339 int size2() const { return storage2; }
1340 const int* const_data_indx() const { return indx; }
1341 bool even_exchanges() const { return d; }
1342 MatrixBandWidth bandwidth() const;
1343 void swap(BandLUMatrix& gm);
1344 NEW_DELETE(BandLUMatrix)
1345};
1346
1347// ************************** special matrices ****************************
1348
1349/// Identity matrix.
1350class IdentityMatrix : public GeneralMatrix
1351{
1352 GeneralMatrix* Image() const; // copy of matrix
1353public:
1354 IdentityMatrix() {}
1355 ~IdentityMatrix() {}
1356 IdentityMatrix(ArrayLengthSpecifier n) : GeneralMatrix(1)
1357 { nrows_val = ncols_val = n.Value(); *store = 1; }
1358 IdentityMatrix(const IdentityMatrix& gm)
1359 : GeneralMatrix() { GetMatrix(&gm); }
1360 IdentityMatrix(const BaseMatrix&);
1361 void operator=(const BaseMatrix&);
1362 void operator=(const IdentityMatrix& m) { Eq(m); }
1363 void operator=(Real f) { GeneralMatrix::operator=(f); }
1364 MatrixType type() const;
1365
1366 LogAndSign log_determinant() const;
1367 Real trace() const;
1368 Real sum_square() const;
1369 Real sum_absolute_value() const;
1370 Real sum() const { return trace(); }
1371 void GetRow(MatrixRowCol&);
1372 void GetCol(MatrixRowCol&);
1373 void GetCol(MatrixColX&);
1374 void NextRow(MatrixRowCol&);
1375 void NextCol(MatrixRowCol&);
1376 void NextCol(MatrixColX&);
1377 GeneralMatrix* MakeSolver() { return this; } // for solving
1378 void Solver(MatrixColX&, const MatrixColX&);
1379 GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
1380 void resize(int n);
1381 void ReSize(int n) { resize(n); }
1382 void resize(const GeneralMatrix& A);
1383 void ReSize(const GeneralMatrix& A) { resize(A); }
1384 MatrixBandWidth bandwidth() const;
1385// ReturnMatrix Reverse() const; // reverse order of elements
1386 void swap(IdentityMatrix& gm)
1387 { GeneralMatrix::swap((GeneralMatrix&)gm); }
1388 NEW_DELETE(IdentityMatrix)
1389};
1390
1391
1392
1393
1394// ************************** GenericMatrix class ************************/
1395
1396/// A matrix which can be of any GeneralMatrix type.
1397class GenericMatrix : public BaseMatrix
1398{
1399 GeneralMatrix* gm;
1400 int search(const BaseMatrix* bm) const;
1401 friend class BaseMatrix;
1402public:
1403 GenericMatrix() : gm(0) {}
1404 GenericMatrix(const BaseMatrix& bm)
1405 { gm = ((BaseMatrix&)bm).Evaluate(); gm = gm->Image(); }
1406 GenericMatrix(const GenericMatrix& bm) : BaseMatrix()
1407 { gm = bm.gm->Image(); }
1408 void operator=(const GenericMatrix&);
1409 void operator=(const BaseMatrix&);
1410 void operator+=(const BaseMatrix&);
1411 void operator-=(const BaseMatrix&);
1412 void operator*=(const BaseMatrix&);
1413 void operator|=(const BaseMatrix&);
1414 void operator&=(const BaseMatrix&);
1415 void operator+=(Real);
1416 void operator-=(Real r) { operator+=(-r); }
1417 void operator*=(Real);
1418 void operator/=(Real r) { operator*=(1.0/r); }
1419 ~GenericMatrix() { delete gm; }
1420 void cleanup() { delete gm; gm = 0; }
1421 void Release() { gm->Release(); }
1422 void release() { gm->release(); }
1423 GeneralMatrix* Evaluate(MatrixType = MatrixTypeUnSp);
1424 MatrixBandWidth bandwidth() const;
1425 void swap(GenericMatrix& x);
1426 NEW_DELETE(GenericMatrix)
1427};
1428
1429// *************************** temporary classes *************************/
1430
1431/// Product of two matrices.
1432/// \internal
1433class MultipliedMatrix : public BaseMatrix
1434{
1435protected:
1436 // if these union statements cause problems, simply remove them
1437 // and declare the items individually
1438 union { const BaseMatrix* bm1; GeneralMatrix* gm1; };
1439 // pointers to summands
1440 union { const BaseMatrix* bm2; GeneralMatrix* gm2; };
1441 MultipliedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
1442 : bm1(bm1x),bm2(bm2x) {}
1443 int search(const BaseMatrix*) const;
1444 friend class BaseMatrix;
1445 friend class GeneralMatrix;
1446 friend class GenericMatrix;
1447public:
1448 ~MultipliedMatrix() {}
1449 GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1450 MatrixBandWidth bandwidth() const;
1451 NEW_DELETE(MultipliedMatrix)
1452};
1453
1454/// Sum of two matrices.
1455/// \internal
1456class AddedMatrix : public MultipliedMatrix
1457{
1458protected:
1459 AddedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
1460 : MultipliedMatrix(bm1x,bm2x) {}
1461
1462 friend class BaseMatrix;
1463 friend class GeneralMatrix;
1464 friend class GenericMatrix;
1465public:
1466 ~AddedMatrix() {}
1467 GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1468 MatrixBandWidth bandwidth() const;
1469 NEW_DELETE(AddedMatrix)
1470};
1471
1472/// Schur (elementwise) product of two matrices.
1473/// \internal
1474class SPMatrix : public AddedMatrix
1475{
1476protected:
1477 SPMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
1478 : AddedMatrix(bm1x,bm2x) {}
1479
1480 friend class BaseMatrix;
1481 friend class GeneralMatrix;
1482 friend class GenericMatrix;
1483public:
1484 ~SPMatrix() {}
1485 GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1486 MatrixBandWidth bandwidth() const;
1487
1488 friend SPMatrix SP(const BaseMatrix&, const BaseMatrix&);
1489
1490 NEW_DELETE(SPMatrix)
1491};
1492
1493/// Kronecker product of two matrices.
1494/// \internal
1495class KPMatrix : public MultipliedMatrix
1496{
1497protected:
1498 KPMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
1499 : MultipliedMatrix(bm1x,bm2x) {}
1500
1501 friend class BaseMatrix;
1502 friend class GeneralMatrix;
1503 friend class GenericMatrix;
1504public:
1505 ~KPMatrix() {}
1506 MatrixBandWidth bandwidth() const;
1507 GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1508 friend KPMatrix KP(const BaseMatrix&, const BaseMatrix&);
1509 NEW_DELETE(KPMatrix)
1510};
1511
1512/// Two matrices horizontally concatenated.
1513/// \internal
1514class ConcatenatedMatrix : public MultipliedMatrix
1515{
1516protected:
1517 ConcatenatedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
1518 : MultipliedMatrix(bm1x,bm2x) {}
1519
1520 friend class BaseMatrix;
1521 friend class GeneralMatrix;
1522 friend class GenericMatrix;
1523public:
1524 ~ConcatenatedMatrix() {}
1525 MatrixBandWidth bandwidth() const;
1526 GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1527 NEW_DELETE(ConcatenatedMatrix)
1528};
1529
1530/// Two matrices vertically concatenated.
1531/// \internal
1532class StackedMatrix : public ConcatenatedMatrix
1533{
1534protected:
1535 StackedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
1536 : ConcatenatedMatrix(bm1x,bm2x) {}
1537
1538 friend class BaseMatrix;
1539 friend class GeneralMatrix;
1540 friend class GenericMatrix;
1541public:
1542 ~StackedMatrix() {}
1543 GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1544 NEW_DELETE(StackedMatrix)
1545};
1546
1547/// Inverted matrix times matrix.
1548/// \internal
1549class SolvedMatrix : public MultipliedMatrix
1550{
1551 SolvedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
1552 : MultipliedMatrix(bm1x,bm2x) {}
1553 friend class BaseMatrix;
1554 friend class InvertedMatrix; // for operator*
1555public:
1556 ~SolvedMatrix() {}
1557 GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1558 MatrixBandWidth bandwidth() const;
1559 NEW_DELETE(SolvedMatrix)
1560};
1561
1562/// Difference between two matrices.
1563/// \internal
1564class SubtractedMatrix : public AddedMatrix
1565{
1566 SubtractedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
1567 : AddedMatrix(bm1x,bm2x) {}
1568 friend class BaseMatrix;
1569 friend class GeneralMatrix;
1570 friend class GenericMatrix;
1571public:
1572 ~SubtractedMatrix() {}
1573 GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1574 NEW_DELETE(SubtractedMatrix)
1575};
1576
1577/// Any type of matrix plus Real.
1578/// \internal
1579class ShiftedMatrix : public BaseMatrix
1580{
1581protected:
1582 union { const BaseMatrix* bm; GeneralMatrix* gm; };
1583 Real f;
1584 ShiftedMatrix(const BaseMatrix* bmx, Real fx) : bm(bmx),f(fx) {}
1585 int search(const BaseMatrix*) const;
1586 friend class BaseMatrix;
1587 friend class GeneralMatrix;
1588 friend class GenericMatrix;
1589public:
1590 ~ShiftedMatrix() {}
1591 GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1592 friend ShiftedMatrix operator+(Real f, const BaseMatrix& BM);
1593 NEW_DELETE(ShiftedMatrix)
1594};
1595
1596/// Real minus matrix.
1597/// \internal
1598class NegShiftedMatrix : public ShiftedMatrix
1599{
1600protected:
1601 NegShiftedMatrix(Real fx, const BaseMatrix* bmx) : ShiftedMatrix(bmx,fx) {}
1602 friend class BaseMatrix;
1603 friend class GeneralMatrix;
1604 friend class GenericMatrix;
1605public:
1606 ~NegShiftedMatrix() {}
1607 GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1608 friend NegShiftedMatrix operator-(Real, const BaseMatrix&);
1609 NEW_DELETE(NegShiftedMatrix)
1610};
1611
1612/// Any type of matrix times Real.
1613/// \internal
1614class ScaledMatrix : public ShiftedMatrix
1615{
1616 ScaledMatrix(const BaseMatrix* bmx, Real fx) : ShiftedMatrix(bmx,fx) {}
1617 friend class BaseMatrix;
1618 friend class GeneralMatrix;
1619 friend class GenericMatrix;
1620public:
1621 ~ScaledMatrix() {}
1622 GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1623 MatrixBandWidth bandwidth() const;
1624 friend ScaledMatrix operator*(Real f, const BaseMatrix& BM);
1625 NEW_DELETE(ScaledMatrix)
1626};
1627
1628/// Any type of matrix times -1.
1629/// \internal
1630class NegatedMatrix : public BaseMatrix
1631{
1632protected:
1633 union { const BaseMatrix* bm; GeneralMatrix* gm; };
1634 NegatedMatrix(const BaseMatrix* bmx) : bm(bmx) {}
1635 int search(const BaseMatrix*) const;
1636private:
1637 friend class BaseMatrix;
1638public:
1639 ~NegatedMatrix() {}
1640 GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1641 MatrixBandWidth bandwidth() const;
1642 NEW_DELETE(NegatedMatrix)
1643};
1644
1645/// Transposed matrix.
1646/// \internal
1647class TransposedMatrix : public NegatedMatrix
1648{
1649 TransposedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
1650 friend class BaseMatrix;
1651public:
1652 ~TransposedMatrix() {}
1653 GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1654 MatrixBandWidth bandwidth() const;
1655 NEW_DELETE(TransposedMatrix)
1656};
1657
1658/// Any type of matrix with order of elements reversed.
1659/// \internal
1660class ReversedMatrix : public NegatedMatrix
1661{
1662 ReversedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
1663 friend class BaseMatrix;
1664public:
1665 ~ReversedMatrix() {}
1666 GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1667 NEW_DELETE(ReversedMatrix)
1668};
1669
1670/// Inverse of matrix.
1671/// \internal
1672class InvertedMatrix : public NegatedMatrix
1673{
1674 InvertedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
1675public:
1676 ~InvertedMatrix() {}
1677 SolvedMatrix operator*(const BaseMatrix&) const; // inverse(A) * B
1678 ScaledMatrix operator*(Real t) const { return BaseMatrix::operator*(t); }
1679 friend class BaseMatrix;
1680 GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1681 MatrixBandWidth bandwidth() const;
1682 NEW_DELETE(InvertedMatrix)
1683};
1684
1685/// Any type of matrix interpreted as a RowVector.
1686/// \internal
1687class RowedMatrix : public NegatedMatrix
1688{
1689 RowedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
1690 friend class BaseMatrix;
1691public:
1692 ~RowedMatrix() {}
1693 GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1694 MatrixBandWidth bandwidth() const;
1695 NEW_DELETE(RowedMatrix)
1696};
1697
1698/// Any type of matrix interpreted as a ColumnVector.
1699/// \internal
1700class ColedMatrix : public NegatedMatrix
1701{
1702 ColedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
1703 friend class BaseMatrix;
1704public:
1705 ~ColedMatrix() {}
1706 GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1707 MatrixBandWidth bandwidth() const;
1708 NEW_DELETE(ColedMatrix)
1709};
1710
1711/// Any type of matrix interpreted as a DiagonalMatrix.
1712/// \internal
1713class DiagedMatrix : public NegatedMatrix
1714{
1715 DiagedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
1716 friend class BaseMatrix;
1717public:
1718 ~DiagedMatrix() {}
1719 GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1720 MatrixBandWidth bandwidth() const;
1721 NEW_DELETE(DiagedMatrix)
1722};
1723
1724/// Any type of matrix interpreted as a (rectangular) Matrix.
1725/// \internal
1726class MatedMatrix : public NegatedMatrix
1727{
1728 int nr, nc;
1729 MatedMatrix(const BaseMatrix* bmx, int nrx, int ncx)
1730 : NegatedMatrix(bmx), nr(nrx), nc(ncx) {}
1731 friend class BaseMatrix;
1732public:
1733 ~MatedMatrix() {}
1734 GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1735 MatrixBandWidth bandwidth() const;
1736 NEW_DELETE(MatedMatrix)
1737};
1738
1739/// A matrix in an "envelope' for return from a function.
1740/// \internal
1741class ReturnMatrix : public BaseMatrix
1742{
1743 GeneralMatrix* gm;
1744 int search(const BaseMatrix*) const;
1745public:
1746 ~ReturnMatrix() {}
1747 GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1748 friend class BaseMatrix;
1749 ReturnMatrix(const ReturnMatrix& tm) : BaseMatrix(), gm(tm.gm) {}
1750 ReturnMatrix(const GeneralMatrix* gmx) : gm((GeneralMatrix*&)gmx) {}
1751// ReturnMatrix(GeneralMatrix&);
1752 MatrixBandWidth bandwidth() const;
1753 NEW_DELETE(ReturnMatrix)
1754};
1755
1756
1757// ************************** submatrices ******************************/
1758
1759/// A submatrix of a matrix.
1760/// \internal
1761class GetSubMatrix : public NegatedMatrix
1762{
1763 int row_skip;
1764 int row_number;
1765 int col_skip;
1766 int col_number;
1767 bool IsSym;
1768
1769 GetSubMatrix
1770 (const BaseMatrix* bmx, int rs, int rn, int cs, int cn, bool is)
1771 : NegatedMatrix(bmx),
1772 row_skip(rs), row_number(rn), col_skip(cs), col_number(cn), IsSym(is) {}
1773 void SetUpLHS();
1774 friend class BaseMatrix;
1775public:
1776 GetSubMatrix(const GetSubMatrix& g)
1777 : NegatedMatrix(g.bm), row_skip(g.row_skip), row_number(g.row_number),
1778 col_skip(g.col_skip), col_number(g.col_number), IsSym(g.IsSym) {}
1779 ~GetSubMatrix() {}
1780 GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1781 void operator=(const BaseMatrix&);
1782 void operator+=(const BaseMatrix&);
1783 void operator-=(const BaseMatrix&);
1784 void operator=(const GetSubMatrix& m) { operator=((const BaseMatrix&)m); }
1785 void operator<<(const BaseMatrix&);
1786 void operator<<(const double*); // copy from array
1787 void operator<<(const float*); // copy from array
1788 void operator<<(const int*); // copy from array
1789 MatrixInput operator<<(double); // for loading a list
1790 MatrixInput operator<<(float); // for loading a list
1791 MatrixInput operator<<(int f);
1792 void operator=(Real); // copy from constant
1793 void operator+=(Real); // add constant
1794 void operator-=(Real r) { operator+=(-r); } // subtract constant
1795 void operator*=(Real); // multiply by constant
1796 void operator/=(Real r) { operator*=(1.0/r); } // divide by constant
1797 void inject(const GeneralMatrix&); // copy stored els only
1798 void Inject(const GeneralMatrix& GM) { inject(GM); }
1799 MatrixBandWidth bandwidth() const;
1800 NEW_DELETE(GetSubMatrix)
1801};
1802
1803// ******************** linear equation solving ****************************/
1804
1805/// A class for finding A.i() * B.
1806/// This is supposed to choose the appropriate method depending on the
1807/// type A. Not very satisfactory as it doesn't know about Cholesky for
1808/// for positive definite matrices.
1809class LinearEquationSolver : public BaseMatrix
1810{
1811 GeneralMatrix* gm;
1812 int search(const BaseMatrix*) const { return 0; }
1813 friend class BaseMatrix;
1814public:
1815 LinearEquationSolver(const BaseMatrix& bm);
1816 ~LinearEquationSolver() { delete gm; }
1817 void cleanup() { delete gm; }
1818 GeneralMatrix* Evaluate(MatrixType) { return gm; }
1819 // probably should have an error message if MatrixType != UnSp
1820 NEW_DELETE(LinearEquationSolver)
1821};
1822
1823// ************************** matrix input *******************************/
1824
1825/// Class for reading values into a (small) matrix within a program.
1826/// \internal
1827/// Is able to detect a mismatch in the number of elements.
1828
1829class MatrixInput
1830{
1831 int n; // number values still to be read
1832 Real* r; // pointer to next location to be read to
1833public:
1834 MatrixInput(const MatrixInput& mi) : n(mi.n), r(mi.r) {}
1835 MatrixInput(int nx, Real* rx) : n(nx), r(rx) {}
1836 ~MatrixInput();
1837 MatrixInput operator<<(double);
1838 MatrixInput operator<<(float);
1839 MatrixInput operator<<(int f);
1840 friend class GeneralMatrix;
1841};
1842
1843
1844
1845// **************** a very simple integer array class ********************/
1846
1847/// A very simple integer array class.
1848/// A minimal array class to imitate a C style array but giving dynamic storage
1849/// mostly intended for internal use by newmat.
1850/// Probably to be replaced by a templated class when I start using templates.
1851
1852class SimpleIntArray : public Janitor
1853{
1854protected:
1855 int* a; ///< pointer to the array
1856 int n; ///< length of the array
1857public:
1858 SimpleIntArray(int xn); ///< build an array length xn
1859 SimpleIntArray() : a(0), n(0) {} ///< build an array length 0
1860 ~SimpleIntArray(); ///< return the space to memory
1861 int& operator[](int i); ///< access element of the array - start at 0
1862 int operator[](int i) const;
1863 ///< access element of constant array
1864 void operator=(int ai); ///< set the array equal to a constant
1865 void operator=(const SimpleIntArray& b);
1866 ///< copy the elements of an array
1867 SimpleIntArray(const SimpleIntArray& b);
1868 ///< make a new array equal to an existing one
1869 int Size() const { return n; }
1870 ///< return the size of the array
1871 int size() const { return n; }
1872 ///< return the size of the array
1873 int* Data() { return a; } ///< pointer to the data
1874 const int* Data() const { return a; } ///< pointer to the data
1875 int* data() { return a; } ///< pointer to the data
1876 const int* data() const { return a; } ///< pointer to the data
1877 const int* const_data() const { return a; } ///< pointer to the data
1878 void resize(int i, bool keep = false);
1879 ///< change length, keep data if keep = true
1880 void ReSize(int i, bool keep = false) { resize(i, keep); }
1881 ///< change length, keep data if keep = true
1882 void resize_keep(int i) { resize(i, true); }
1883 ///< change length, keep data
1884 void cleanup() { resize(0); } ///< set length to zero
1885 void CleanUp() { resize(0); } ///< set length to zero
1886 NEW_DELETE(SimpleIntArray)
1887};
1888
1889// ********************** C subscript classes ****************************
1890
1891/// Let matrix simulate a C type two dimensional array
1892class RealStarStar
1893{
1894 Real** a;
1895public:
1896 RealStarStar(Matrix& A);
1897 ~RealStarStar() { delete [] a; }
1898 operator Real**() { return a; }
1899};
1900
1901/// Let matrix simulate a C type const two dimensional array
1902class ConstRealStarStar
1903{
1904 const Real** a;
1905public:
1906 ConstRealStarStar(const Matrix& A);
1907 ~ConstRealStarStar() { delete [] a; }
1908 operator const Real**() { return a; }
1909};
1910
1911// *************************** exceptions ********************************/
1912
1913/// Not positive definite exception.
1914class NPDException : public Runtime_error
1915{
1916public:
1917 static unsigned long Select;
1918 NPDException(const GeneralMatrix&);
1919};
1920
1921/// Covergence failure exception.
1922class ConvergenceException : public Runtime_error
1923{
1924public:
1925 static unsigned long Select;
1926 ConvergenceException(const GeneralMatrix& A);
1927 ConvergenceException(const char* c);
1928};
1929
1930/// Singular matrix exception.
1931class SingularException : public Runtime_error
1932{
1933public:
1934 static unsigned long Select;
1935 SingularException(const GeneralMatrix& A);
1936};
1937
1938/// Real overflow exception.
1939class OverflowException : public Runtime_error
1940{
1941public:
1942 static unsigned long Select;
1943 OverflowException(const char* c);
1944};
1945
1946/// Miscellaneous exception (details in character string).
1947class ProgramException : public Logic_error
1948{
1949protected:
1950 ProgramException();
1951public:
1952 static unsigned long Select;
1953 ProgramException(const char* c);
1954 ProgramException(const char* c, const GeneralMatrix&);
1955 ProgramException(const char* c, const GeneralMatrix&, const GeneralMatrix&);
1956 ProgramException(const char* c, MatrixType, MatrixType);
1957};
1958
1959/// Index exception.
1960class IndexException : public Logic_error
1961{
1962public:
1963 static unsigned long Select;
1964 IndexException(int i, const GeneralMatrix& A);
1965 IndexException(int i, int j, const GeneralMatrix& A);
1966 // next two are for access via element function
1967 IndexException(int i, const GeneralMatrix& A, bool);
1968 IndexException(int i, int j, const GeneralMatrix& A, bool);
1969};
1970
1971/// Cannot convert to vector exception.
1972class VectorException : public Logic_error
1973{
1974public:
1975 static unsigned long Select;
1976 VectorException();
1977 VectorException(const GeneralMatrix& A);
1978};
1979
1980/// A matrix is not square exception.
1981class NotSquareException : public Logic_error
1982{
1983public:
1984 static unsigned long Select;
1985 NotSquareException(const GeneralMatrix& A);
1986 NotSquareException();
1987};
1988
1989/// Submatrix dimension exception.
1990class SubMatrixDimensionException : public Logic_error
1991{
1992public:
1993 static unsigned long Select;
1994 SubMatrixDimensionException();
1995};
1996
1997/// Incompatible dimensions exception.
1998class IncompatibleDimensionsException : public Logic_error
1999{
2000public:
2001 static unsigned long Select;
2002 IncompatibleDimensionsException();
2003 IncompatibleDimensionsException(const GeneralMatrix&);
2004 IncompatibleDimensionsException(const GeneralMatrix&, const GeneralMatrix&);
2005};
2006
2007/// Not defined exception.
2008class NotDefinedException : public Logic_error
2009{
2010public:
2011 static unsigned long Select;
2012 NotDefinedException(const char* op, const char* matrix);
2013};
2014
2015/// Cannot build matrix with these properties exception.
2016class CannotBuildException : public Logic_error
2017{
2018public:
2019 static unsigned long Select;
2020 CannotBuildException(const char* matrix);
2021};
2022
2023
2024/// Internal newmat exception - shouldn't happen.
2025class InternalException : public Logic_error
2026{
2027public:
2028 static unsigned long Select; // for identifying exception
2029 InternalException(const char* c);
2030};
2031
2032// ************************ functions ************************************ //
2033
2034bool operator==(const GeneralMatrix& A, const GeneralMatrix& B);
2035bool operator==(const BaseMatrix& A, const BaseMatrix& B);
2036inline bool operator!=(const GeneralMatrix& A, const GeneralMatrix& B)
2037 { return ! (A==B); }
2038inline bool operator!=(const BaseMatrix& A, const BaseMatrix& B)
2039 { return ! (A==B); }
2040
2041 // inequality operators are dummies included for compatibility
2042 // with STL. They throw an exception if actually called.
2043inline bool operator<=(const BaseMatrix& A, const BaseMatrix&)
2044 { A.IEQND(); return true; }
2045inline bool operator>=(const BaseMatrix& A, const BaseMatrix&)
2046 { A.IEQND(); return true; }
2047inline bool operator<(const BaseMatrix& A, const BaseMatrix&)
2048 { A.IEQND(); return true; }
2049inline bool operator>(const BaseMatrix& A, const BaseMatrix&)
2050 { A.IEQND(); return true; }
2051
2052bool is_zero(const BaseMatrix& A);
2053inline bool IsZero(const BaseMatrix& A) { return is_zero(A); }
2054
2055Real dotproduct(const Matrix& A, const Matrix& B);
2056Matrix crossproduct(const Matrix& A, const Matrix& B);
2057ReturnMatrix crossproduct_rows(const Matrix& A, const Matrix& B);
2058ReturnMatrix crossproduct_columns(const Matrix& A, const Matrix& B);
2059
2060inline Real DotProduct(const Matrix& A, const Matrix& B)
2061 { return dotproduct(A, B); }
2062inline Matrix CrossProduct(const Matrix& A, const Matrix& B)
2063 { return crossproduct(A, B); }
2064inline ReturnMatrix CrossProductRows(const Matrix& A, const Matrix& B)
2065 { return crossproduct_rows(A, B); }
2066inline ReturnMatrix CrossProductColumns(const Matrix& A, const Matrix& B)
2067 { return crossproduct_columns(A, B); }
2068
2069void newmat_block_copy(int n, Real* from, Real* to);
2070
2071// ********************* friend functions ******************************** //
2072
2073// Functions declared as friends - G++ wants them declared externally as well
2074
2075bool Rectangular(MatrixType a, MatrixType b, MatrixType c);
2076bool Compare(const MatrixType&, MatrixType&);
2077Real dotproduct(const Matrix& A, const Matrix& B);
2078SPMatrix SP(const BaseMatrix&, const BaseMatrix&);
2079KPMatrix KP(const BaseMatrix&, const BaseMatrix&);
2080ShiftedMatrix operator+(Real f, const BaseMatrix& BM);
2081NegShiftedMatrix operator-(Real, const BaseMatrix&);
2082ScaledMatrix operator*(Real f, const BaseMatrix& BM);
2083
2084// ********************* inline functions ******************************** //
2085
2086inline LogAndSign log_determinant(const BaseMatrix& B)
2087 { return B.log_determinant(); }
2088inline LogAndSign LogDeterminant(const BaseMatrix& B)
2089 { return B.log_determinant(); }
2090inline Real determinant(const BaseMatrix& B)
2091 { return B.determinant(); }
2092inline Real Determinant(const BaseMatrix& B)
2093 { return B.determinant(); }
2094inline Real sum_square(const BaseMatrix& B) { return B.sum_square(); }
2095inline Real SumSquare(const BaseMatrix& B) { return B.sum_square(); }
2096inline Real norm_Frobenius(const BaseMatrix& B) { return B.norm_Frobenius(); }
2097inline Real norm_frobenius(const BaseMatrix& B) { return B.norm_Frobenius(); }
2098inline Real NormFrobenius(const BaseMatrix& B) { return B.norm_Frobenius(); }
2099inline Real trace(const BaseMatrix& B) { return B.trace(); }
2100inline Real Trace(const BaseMatrix& B) { return B.trace(); }
2101inline Real sum_absolute_value(const BaseMatrix& B)
2102 { return B.sum_absolute_value(); }
2103inline Real SumAbsoluteValue(const BaseMatrix& B)
2104 { return B.sum_absolute_value(); }
2105inline Real sum(const BaseMatrix& B)
2106 { return B.sum(); }
2107inline Real Sum(const BaseMatrix& B)
2108 { return B.sum(); }
2109inline Real maximum_absolute_value(const BaseMatrix& B)
2110 { return B.maximum_absolute_value(); }
2111inline Real MaximumAbsoluteValue(const BaseMatrix& B)
2112 { return B.maximum_absolute_value(); }
2113inline Real minimum_absolute_value(const BaseMatrix& B)
2114 { return B.minimum_absolute_value(); }
2115inline Real MinimumAbsoluteValue(const BaseMatrix& B)
2116 { return B.minimum_absolute_value(); }
2117inline Real maximum(const BaseMatrix& B) { return B.maximum(); }
2118inline Real Maximum(const BaseMatrix& B) { return B.maximum(); }
2119inline Real minimum(const BaseMatrix& B) { return B.minimum(); }
2120inline Real Minimum(const BaseMatrix& B) { return B.minimum(); }
2121inline Real norm1(const BaseMatrix& B) { return B.norm1(); }
2122inline Real Norm1(const BaseMatrix& B) { return B.norm1(); }
2123inline Real norm1(RowVector& RV) { return RV.maximum_absolute_value(); }
2124inline Real Norm1(RowVector& RV) { return RV.maximum_absolute_value(); }
2125inline Real norm_infinity(const BaseMatrix& B) { return B.norm_infinity(); }
2126inline Real NormInfinity(const BaseMatrix& B) { return B.norm_infinity(); }
2127inline Real norm_infinity(ColumnVector& CV)
2128 { return CV.maximum_absolute_value(); }
2129inline Real NormInfinity(ColumnVector& CV)
2130 { return CV.maximum_absolute_value(); }
2131inline bool IsZero(const GeneralMatrix& A) { return A.IsZero(); }
2132inline bool is_zero(const GeneralMatrix& A) { return A.is_zero(); }
2133
2134
2135inline MatrixInput MatrixInput::operator<<(int f) { return *this << (Real)f; }
2136inline MatrixInput GeneralMatrix::operator<<(int f) { return *this << (Real)f; }
2137inline MatrixInput BandMatrix::operator<<(int f) { return *this << (Real)f; }
2138inline MatrixInput GetSubMatrix::operator<<(int f) { return *this << (Real)f; }
2139
2140inline ReversedMatrix BaseMatrix::Reverse() const { return reverse(); }
2141inline RowedMatrix BaseMatrix::AsRow() const { return as_row(); }
2142inline ColedMatrix BaseMatrix::AsColumn() const { return as_column(); }
2143inline DiagedMatrix BaseMatrix::AsDiagonal() const { return as_diagonal(); }
2144inline MatedMatrix BaseMatrix::AsMatrix(int m, int n) const
2145 { return as_matrix(m, n); }
2146inline GetSubMatrix BaseMatrix::SubMatrix(int fr, int lr, int fc, int lc) const
2147 { return submatrix(fr, lr, fc, lc); }
2148inline GetSubMatrix BaseMatrix::SymSubMatrix(int f, int l) const
2149 { return sym_submatrix(f, l); }
2150inline GetSubMatrix BaseMatrix::Row(int f) const { return row(f); }
2151inline GetSubMatrix BaseMatrix::Rows(int f, int l) const { return rows(f, l); }
2152inline GetSubMatrix BaseMatrix::Column(int f) const { return column(f); }
2153inline GetSubMatrix BaseMatrix::Columns(int f, int l) const
2154 { return columns(f, l); }
2155inline Real BaseMatrix::AsScalar() const { return as_scalar(); }
2156
2157inline ReturnMatrix GeneralMatrix::ForReturn() const { return for_return(); }
2158
2159inline void swap(Matrix& A, Matrix& B) { A.swap(B); }
2160inline void swap(SquareMatrix& A, SquareMatrix& B) { A.swap(B); }
2161inline void swap(nricMatrix& A, nricMatrix& B) { A.swap(B); }
2162inline void swap(UpperTriangularMatrix& A, UpperTriangularMatrix& B)
2163 { A.swap(B); }
2164inline void swap(LowerTriangularMatrix& A, LowerTriangularMatrix& B)
2165 { A.swap(B); }
2166inline void swap(SymmetricMatrix& A, SymmetricMatrix& B) { A.swap(B); }
2167inline void swap(DiagonalMatrix& A, DiagonalMatrix& B) { A.swap(B); }
2168inline void swap(RowVector& A, RowVector& B) { A.swap(B); }
2169inline void swap(ColumnVector& A, ColumnVector& B) { A.swap(B); }
2170inline void swap(CroutMatrix& A, CroutMatrix& B) { A.swap(B); }
2171inline void swap(BandMatrix& A, BandMatrix& B) { A.swap(B); }
2172inline void swap(UpperBandMatrix& A, UpperBandMatrix& B) { A.swap(B); }
2173inline void swap(LowerBandMatrix& A, LowerBandMatrix& B) { A.swap(B); }
2174inline void swap(SymmetricBandMatrix& A, SymmetricBandMatrix& B) { A.swap(B); }
2175inline void swap(BandLUMatrix& A, BandLUMatrix& B) { A.swap(B); }
2176inline void swap(IdentityMatrix& A, IdentityMatrix& B) { A.swap(B); }
2177inline void swap(GenericMatrix& A, GenericMatrix& B) { A.swap(B); }
2178
2179#ifdef OPT_COMPATIBLE // for compatibility with opt++
2180
2181inline Real Norm2(const ColumnVector& CV) { return CV.norm_Frobenius(); }
2182inline Real Dot(ColumnVector& CV1, ColumnVector& CV2)
2183 { return dotproduct(CV1, CV2); }
2184
2185#endif
2186
2187
2188#ifdef use_namespace
2189}
2190#endif
2191
2192
2193#endif
2194
2195// body file: newmat1.cpp
2196// body file: newmat2.cpp
2197// body file: newmat3.cpp
2198// body file: newmat4.cpp
2199// body file: newmat5.cpp
2200// body file: newmat6.cpp
2201// body file: newmat7.cpp
2202// body file: newmat8.cpp
2203// body file: newmatex.cpp
2204// body file: bandmat.cpp
2205// body file: submat.cpp
2206
2207
2208
2209///@}
2210
2211
2212
2213
Note: See TracBrowser for help on using the repository browser.