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