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

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

upgrade to newmat11 library
File size: 88.8 KB
RevLine 
[8901]1/// \defgroup newmat Newmat matrix manipulation library
2///@{
3
4/// \file newmat.h
5/// Definition file for matrix library.
6
7// Copyright (C) 2004: R B Davies
8
9#ifndef NEWMAT_LIB
10#define NEWMAT_LIB 0
11
12#include "include.h"
13
14#include "myexcept.h"
15
16
17#ifdef use_namespace
18namespace NEWMAT { using namespace RBD_COMMON; }
19namespace RBD_LIBRARIES { using namespace NEWMAT; }
20namespace NEWMAT {
21#endif
22
23//#define DO_REPORT // to activate REPORT
24
25#ifdef NO_LONG_NAMES
26#define UpperTriangularMatrix UTMatrix
27#define LowerTriangularMatrix LTMatrix
28#define SymmetricMatrix SMatrix
29#define DiagonalMatrix DMatrix
30#define BandMatrix BMatrix
31#define UpperBandMatrix UBMatrix
32#define LowerBandMatrix LBMatrix
33#define SymmetricBandMatrix SBMatrix
34#define BandLUMatrix BLUMatrix
35#endif
36
37// ************************** general utilities ****************************/
38
39class GeneralMatrix; // defined later
40class BaseMatrix; // defined later
41class MatrixInput; // defined later
42
43void MatrixErrorNoSpace(const void*); ///< test for allocation fails
44
45/// Return from LogDeterminant function.
46/// Members are the log of the absolute value and the sign (+1, -1 or 0)
47class LogAndSign
48{
49 Real log_val;
50 int sign_val;
51public:
52 LogAndSign() { log_val=0.0; sign_val=1; }
53 LogAndSign(Real);
54 void operator*=(Real); ///< multiply by a real
55 void pow_eq(int k); ///< raise to power of k
56 void PowEq(int k) { pow_eq(k); }
57 void ChangeSign() { sign_val = -sign_val; }
58 void change_sign() { sign_val = -sign_val; } ///< change sign
59 Real LogValue() const { return log_val; }
60 Real log_value() const { return log_val; } ///< log of the absolute value
61 int Sign() const { return sign_val; }
62 int sign() const { return sign_val; } ///< sign of the value
63 Real value() const; ///< the value
64 Real Value() const { return value(); }
65 FREE_CHECK(LogAndSign)
66};
67
68// the following class is for counting the number of times a piece of code
69// is executed. It is used for locating any code not executed by test
70// routines. Use turbo GREP locate all places this code is called and
71// check which ones are not accessed.
72// Somewhat implementation dependent as it relies on "cout" still being
73// present when ExeCounter objects are destructed.
74
75#ifdef DO_REPORT
76
77class ExeCounter
78{
79 int line; // code line number
80 int fileid; // file identifier
81 long nexe; // number of executions
82 static int nreports; // number of reports
83public:
84 ExeCounter(int,int);
85 void operator++() { nexe++; }
86 ~ExeCounter(); // prints out reports
87};
88
89#endif
90
91
92// ************************** class MatrixType *****************************
93
94/// Find the type of a matrix resulting from matrix operations.
95/// Also identify what conversions are permissible.
96/// This class must be updated when new matrix types are added.
97
98class MatrixType
99{
100public:
101 enum Attribute { Valid = 1,
102 Diagonal = 2, // order of these is important
103 Symmetric = 4,
104 Band = 8,
105 Lower = 16,
106 Upper = 32,
107 Square = 64,
108 Skew = 128,
109 LUDeco = 256,
110 Ones = 512 };
111
112 enum { US = 0,
113 UT = Valid + Upper + Square,
114 LT = Valid + Lower + Square,
115 Rt = Valid,
116 Sq = Valid + Square,
117 Sm = Valid + Symmetric + Square,
118 Sk = Valid + Skew + Square,
119 Dg = Valid + Diagonal + Band + Lower + Upper + Symmetric
120 + Square,
121 Id = Valid + Diagonal + Band + Lower + Upper + Symmetric
122 + Square + Ones,
123 RV = Valid, // do not separate out
124 CV = Valid, // vectors
125 BM = Valid + Band + Square,
126 UB = Valid + Band + Upper + Square,
127 LB = Valid + Band + Lower + Square,
128 SB = Valid + Band + Symmetric + Square,
129 KB = Valid + Band + Skew + Square,
130 Ct = Valid + LUDeco + Square,
131 BC = Valid + Band + LUDeco + Square,
132 Mask = ~Square
133 };
134
135
136 static int nTypes() { return 13; } // number of different types
137 // exclude Ct, US, BC
138public:
139 int attribute;
140 bool DataLossOK; // true if data loss is OK when
141 // this represents a destination
142public:
143 MatrixType () : DataLossOK(false) {}
144 MatrixType (int i) : attribute(i), DataLossOK(false) {}
145 MatrixType (int i, bool dlok) : attribute(i), DataLossOK(dlok) {}
146 MatrixType (const MatrixType& mt)
147 : attribute(mt.attribute), DataLossOK(mt.DataLossOK) {}
148 void operator=(const MatrixType& mt)
149 { attribute = mt.attribute; DataLossOK = mt.DataLossOK; }
150 void SetDataLossOK() { DataLossOK = true; }
151 int operator+() const { return attribute; }
152 MatrixType operator+(MatrixType mt) const
153 { return MatrixType(attribute & mt.attribute); }
154 MatrixType operator*(const MatrixType&) const;
155 MatrixType SP(const MatrixType&) const;
156 MatrixType KP(const MatrixType&) const;
157 MatrixType operator|(const MatrixType& mt) const
158 { return MatrixType(attribute & mt.attribute & Valid); }
159 MatrixType operator&(const MatrixType& mt) const
160 { return MatrixType(attribute & mt.attribute & Valid); }
161 bool operator>=(MatrixType mt) const
162 { return ( attribute & ~mt.attribute & Mask ) == 0; }
163 bool operator<(MatrixType mt) const // for MS Visual C++ 4
164 { return ( attribute & ~mt.attribute & Mask ) != 0; }
165 bool operator==(MatrixType t) const
166 { return (attribute == t.attribute); }
167 bool operator!=(MatrixType t) const
168 { return (attribute != t.attribute); }
169 bool operator!() const { return (attribute & Valid) == 0; }
170 MatrixType i() const; ///< type of inverse
171 MatrixType t() const; ///< type of transpose
172 MatrixType AddEqualEl() const ///< add constant to matrix
173 { return MatrixType(attribute & (Valid + Symmetric + Square)); }
174 MatrixType MultRHS() const; ///< type for rhs of multiply
175 MatrixType sub() const ///< type of submatrix
176 { return MatrixType(attribute & Valid); }
177 MatrixType ssub() const ///< type of sym submatrix
178 { return MatrixType(attribute); } // not for selection matrix
179 GeneralMatrix* New() const; ///< new matrix of given type
180 GeneralMatrix* New(int,int,BaseMatrix*) const;
181 ///< new matrix of given type
182 const char* value() const; ///< type as char string
183 const char* Value() const { return value(); }
184 friend bool Rectangular(MatrixType a, MatrixType b, MatrixType c);
185 friend bool Compare(const MatrixType&, MatrixType&);
186 ///< compare and check conversion
187 bool is_band() const { return (attribute & Band) != 0; }
188 bool is_diagonal() const { return (attribute & Diagonal) != 0; }
189 bool is_symmetric() const { return (attribute & Symmetric) != 0; }
190 bool CannotConvert() const { return (attribute & LUDeco) != 0; }
191 // used by operator==
192 FREE_CHECK(MatrixType)
193};
194
195
196// *********************** class MatrixBandWidth ***********************/
197
198///Upper and lower bandwidths of a matrix.
199///That is number of diagonals strictly above or below main diagonal,
200///e.g. diagonal matrix has 0 upper and lower bandwiths.
201///-1 means the matrix may have the maximum bandwidth.
202class MatrixBandWidth
203{
204public:
205 int lower_val;
206 int upper_val;
207 MatrixBandWidth(const int l, const int u) : lower_val(l), upper_val(u) {}
208 MatrixBandWidth(const int i) : lower_val(i), upper_val(i) {}
209 MatrixBandWidth operator+(const MatrixBandWidth&) const;
210 MatrixBandWidth operator*(const MatrixBandWidth&) const;
211 MatrixBandWidth minimum(const MatrixBandWidth&) const;
212 MatrixBandWidth t() const { return MatrixBandWidth(upper_val,lower_val); }
213 bool operator==(const MatrixBandWidth& bw) const
214 { return (lower_val == bw.lower_val) && (upper_val == bw.upper_val); }
215 bool operator!=(const MatrixBandWidth& bw) const { return !operator==(bw); }
216 int Upper() const { return upper_val; }
217 int upper() const { return upper_val; }
218 int Lower() const { return lower_val; }
219 int lower() const { return lower_val; }
220 FREE_CHECK(MatrixBandWidth)
221};
222
223
224// ********************* Array length specifier ************************/
225
226/// This class is used to avoid constructors such as
227/// ColumnVector(int) being used for conversions.
228/// Eventually this should be replaced by the use of the keyword "explicit".
229
230class ArrayLengthSpecifier
231{
232 int v;
233public:
234 int Value() const { return v; }
235 int value() const { return v; }
236 ArrayLengthSpecifier(int l) : v(l) {}
237};
238
239// ************************* Matrix routines ***************************/
240
241
242class MatrixRowCol; // defined later
243class MatrixRow;
244class MatrixCol;
245class MatrixColX;
246
247class GeneralMatrix; // defined later
248class AddedMatrix;
249class MultipliedMatrix;
250class SubtractedMatrix;
251class SPMatrix;
252class KPMatrix;
253class ConcatenatedMatrix;
254class StackedMatrix;
255class SolvedMatrix;
256class ShiftedMatrix;
257class NegShiftedMatrix;
258class ScaledMatrix;
259class TransposedMatrix;
260class ReversedMatrix;
261class NegatedMatrix;
262class InvertedMatrix;
263class RowedMatrix;
264class ColedMatrix;
265class DiagedMatrix;
266class MatedMatrix;
267class GetSubMatrix;
268class ReturnMatrix;
269class Matrix;
270class SquareMatrix;
271class nricMatrix;
272class RowVector;
273class ColumnVector;
274class SymmetricMatrix;
275class UpperTriangularMatrix;
276class LowerTriangularMatrix;
277class DiagonalMatrix;
278class CroutMatrix;
279class BandMatrix;
280class LowerBandMatrix;
281class UpperBandMatrix;
282class SymmetricBandMatrix;
283class LinearEquationSolver;
284class GenericMatrix;
285
286
287#define MatrixTypeUnSp 0
288//static MatrixType MatrixTypeUnSp(MatrixType::US);
289// // AT&T needs this
290
291/// Base of the matrix classes.
292class BaseMatrix : public Janitor
293{
294protected:
295 virtual int search(const BaseMatrix*) const = 0;
296 // count number of times matrix is referred to
297public:
298 virtual GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp) = 0;
299 // evaluate temporary
300 // for old version of G++
301 // virtual GeneralMatrix* Evaluate(MatrixType mt) = 0;
302 // GeneralMatrix* Evaluate() { return Evaluate(MatrixTypeUnSp); }
303 AddedMatrix operator+(const BaseMatrix&) const; // results of operations
304 MultipliedMatrix operator*(const BaseMatrix&) const;
305 SubtractedMatrix operator-(const BaseMatrix&) const;
306 ConcatenatedMatrix operator|(const BaseMatrix&) const;
307 StackedMatrix operator&(const BaseMatrix&) const;
308 ShiftedMatrix operator+(Real) const;
309 ScaledMatrix operator*(Real) const;
310 ScaledMatrix operator/(Real) const;
311 ShiftedMatrix operator-(Real) const;
312 TransposedMatrix t() const;
313// TransposedMatrix t;
314 NegatedMatrix operator-() const; // change sign of elements
315 ReversedMatrix reverse() const;
316 ReversedMatrix Reverse() const;
317 InvertedMatrix i() const;
318// InvertedMatrix i;
319 RowedMatrix as_row() const;
320 RowedMatrix AsRow() const;
321 ColedMatrix as_column() const;
322 ColedMatrix AsColumn() const;
323 DiagedMatrix as_diagonal() const;
324 DiagedMatrix AsDiagonal() const;
325 MatedMatrix as_matrix(int,int) const;
326 MatedMatrix AsMatrix(int m, int n) const;
327 GetSubMatrix submatrix(int,int,int,int) const;
328 GetSubMatrix SubMatrix(int fr, int lr, int fc, int lc) const;
329 GetSubMatrix sym_submatrix(int,int) const;
330 GetSubMatrix SymSubMatrix(int f, int l) const;
331 GetSubMatrix row(int) const;
332 GetSubMatrix rows(int,int) const;
333 GetSubMatrix column(int) const;
334 GetSubMatrix columns(int,int) const;
335 GetSubMatrix Row(int f) const;
336 GetSubMatrix Rows(int f, int l) const;
337 GetSubMatrix Column(int f) const;
338 GetSubMatrix Columns(int f, int l) const;
339 Real as_scalar() const; // conversion of 1 x 1 matrix
340 Real AsScalar() const;
341 virtual LogAndSign log_determinant() const;
342 LogAndSign LogDeterminant() const { return log_determinant(); }
343 Real determinant() const;
344 Real Determinant() const { return determinant(); }
345 virtual Real sum_square() const;
346 Real SumSquare() const { return sum_square(); }
347 Real norm_Frobenius() const;
348 Real norm_frobenius() const { return norm_Frobenius(); }
349 Real NormFrobenius() const { return norm_Frobenius(); }
350 virtual Real sum_absolute_value() const;
351 Real SumAbsoluteValue() const { return sum_absolute_value(); }
352 virtual Real sum() const;
353 virtual Real Sum() const { return sum(); }
354 virtual Real maximum_absolute_value() const;
355 Real MaximumAbsoluteValue() const { return maximum_absolute_value(); }
356 virtual Real maximum_absolute_value1(int& i) const;
357 Real MaximumAbsoluteValue1(int& i) const
358 { return maximum_absolute_value1(i); }
359 virtual Real maximum_absolute_value2(int& i, int& j) const;
360 Real MaximumAbsoluteValue2(int& i, int& j) const
361 { return maximum_absolute_value2(i,j); }
362 virtual Real minimum_absolute_value() const;
363 Real MinimumAbsoluteValue() const { return minimum_absolute_value(); }
364 virtual Real minimum_absolute_value1(int& i) const;
365 Real MinimumAbsoluteValue1(int& i) const
366 { return minimum_absolute_value1(i); }
367 virtual Real minimum_absolute_value2(int& i, int& j) const;
368 Real MinimumAbsoluteValue2(int& i, int& j) const
369 { return minimum_absolute_value2(i,j); }
370 virtual Real maximum() const;
371 Real Maximum() const { return maximum(); }
372 virtual Real maximum1(int& i) const;
373 Real Maximum1(int& i) const { return maximum1(i); }
374 virtual Real maximum2(int& i, int& j) const;
375 Real Maximum2(int& i, int& j) const { return maximum2(i,j); }
376 virtual Real minimum() const;
377 Real Minimum() const { return minimum(); }
378 virtual Real minimum1(int& i) const;
379 Real Minimum1(int& i) const { return minimum1(i); }
380 virtual Real minimum2(int& i, int& j) const;
381 Real Minimum2(int& i, int& j) const { return minimum2(i,j); }
382 virtual Real trace() const;
383 Real Trace() const { return trace(); }
384 Real norm1() const;
385 Real Norm1() const { return norm1(); }
386 Real norm_infinity() const;
387 Real NormInfinity() const { return norm_infinity(); }
388 virtual MatrixBandWidth bandwidth() const; // bandwidths of band matrix
389 virtual MatrixBandWidth BandWidth() const { return bandwidth(); }
390 void IEQND() const; // called by ineq. ops
391 ReturnMatrix sum_square_columns() const;
392 ReturnMatrix sum_square_rows() const;
393 ReturnMatrix sum_columns() const;
394 ReturnMatrix sum_rows() const;
395 virtual void cleanup() {}
396 void CleanUp() { cleanup(); }
397
398// virtual ReturnMatrix Reverse() const; // reverse order of elements
399//protected:
400// BaseMatrix() : t(this), i(this) {}
401
402 friend class GeneralMatrix;
403 friend class Matrix;
404 friend class SquareMatrix;
405 friend class nricMatrix;
406 friend class RowVector;
407 friend class ColumnVector;
408 friend class SymmetricMatrix;
409 friend class UpperTriangularMatrix;
410 friend class LowerTriangularMatrix;
411 friend class DiagonalMatrix;
412 friend class CroutMatrix;
413 friend class BandMatrix;
414 friend class LowerBandMatrix;
415 friend class UpperBandMatrix;
416 friend class SymmetricBandMatrix;
417 friend class AddedMatrix;
418 friend class MultipliedMatrix;
419 friend class SubtractedMatrix;
420 friend class SPMatrix;
421 friend class KPMatrix;
422 friend class ConcatenatedMatrix;
423 friend class StackedMatrix;
424 friend class SolvedMatrix;
425 friend class ShiftedMatrix;
426 friend class NegShiftedMatrix;
427 friend class ScaledMatrix;
428 friend class TransposedMatrix;
429 friend class ReversedMatrix;
430 friend class NegatedMatrix;
431 friend class InvertedMatrix;
432 friend class RowedMatrix;
433 friend class ColedMatrix;
434 friend class DiagedMatrix;
435 friend class MatedMatrix;
436 friend class GetSubMatrix;
437 friend class ReturnMatrix;
438 friend class LinearEquationSolver;
439 friend class GenericMatrix;
440 NEW_DELETE(BaseMatrix)
441};
442
443
444// ***************************** working classes **************************/
445
446/// The classes for matrices that can contain data are derived from this.
447class GeneralMatrix : public BaseMatrix // declarable matrix types
448{
449protected:
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
474 virtual short SimpleAddOK(const GeneralMatrix*) { return 0; }
475 // see bandmat.cpp for explanation
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
485public:
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
492 MatrixType Type() const { return type(); }
493 int Nrows() const { return nrows_val; } // get dimensions
494 int Ncols() const { return ncols_val; }
495 int Storage() const { return storage; }
496 Real* Store() const { return store; }
497 // updated names
498 int nrows() const { return nrows_val; } // get dimensions
499 int ncols() const { return ncols_val; }
500 int size() const { return storage; }
501 Real* data() { return store; }
502 const Real* data() const { return store; }
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
508 void protect() { tag_val=-1; } // cannot delete or reuse
509 void Protect() { tag_val=-1; } // cannot delete or reuse
510 int tag() const { return tag_val; }
511 int Tag() const { return tag_val; }
512 bool is_zero() const; // test matrix has all zeros
513 bool IsZero() const { return is_zero(); } // test matrix has all zeros
514 void Release() { tag_val=1; } // del store after next use
515 void Release(int t) { tag_val=t; } // del store after t accesses
516 void ReleaseAndDelete() { tag_val=0; } // delete matrix after use
517 void release() { tag_val=1; } // del store after next use
518 void release(int t) { tag_val=t; } // del store after t accesses
519 void release_and_delete() { tag_val=0; } // delete matrix after use
520 void operator<<(const double*); // assignment from an array
521 void operator<<(const float*); // assignment from an array
522 void operator<<(const int*); // assignment from an array
523 void operator<<(const BaseMatrix& X)
524 { Eq(X,this->type(),true); } // = without checking type
525 void inject(const GeneralMatrix&); // copy stored els only
526 void Inject(const GeneralMatrix& GM) { inject(GM); }
527 void operator+=(const BaseMatrix&);
528 void 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
538 virtual void Solver(MatrixColX&, const MatrixColX&) {}
539 virtual void GetRow(MatrixRowCol&) = 0; // Get matrix row
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
544 virtual void RestoreCol(MatrixRowCol&) {} // Restore matrix col
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
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);
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.
627class Matrix : public GeneralMatrix
628{
629public:
630 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
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);
663 virtual void ReSize(int m, int n) { resize(m, n); }
664 void 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); }
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); }
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); }
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.
685class SquareMatrix : public Matrix
686{
687public:
688 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;
697 SquareMatrix(const SquareMatrix& gm) : Matrix() { GetMatrix(&gm); }
698 SquareMatrix(const Matrix& gm);
699 void resize(int); // change dimensions
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
704 void ReSize(int m, int n) { resize(m, n); }
705 void resize(const GeneralMatrix& A);
706 void ReSize(const GeneralMatrix& A) { resize(A); }
707 void operator+=(const Matrix& M) { PlusEqual(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); }
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); }
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.
721class nricMatrix : public Matrix
722{
723 Real** row_pointer; // points to rows
724 void MakeRowPointer(); // build rowpointer
725 void DeleteRowPointer();
726public:
727 nricMatrix() {}
728 nricMatrix(int m, int n) // standard declaration
729 : Matrix(m,n) { MakeRowPointer(); }
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(); }
739 nricMatrix(const nricMatrix& gm) : Matrix()
740 { GetMatrix(&gm); MakeRowPointer(); }
741 void resize(int m, int n) // change dimensions
742 { DeleteRowPointer(); Matrix::resize(m,n); MakeRowPointer(); }
743 void resize_keep(int m, int n) // change dimensions
744 { DeleteRowPointer(); Matrix::resize_keep(m,n); MakeRowPointer(); }
745 void ReSize(int m, int n) // change dimensions
746 { DeleteRowPointer(); Matrix::resize(m,n); MakeRowPointer(); }
747 void resize(const GeneralMatrix& A);
748 void ReSize(const GeneralMatrix& A) { resize(A); }
749 ~nricMatrix() { DeleteRowPointer(); }
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); }
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); }
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.
767class SymmetricMatrix : public GeneralMatrix
768{
769public:
770 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;
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&);
795 void RestoreCol(MatrixRowCol&) {}
796 void RestoreCol(MatrixColX&);
797 GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
798 void resize(int); // change dimensions
799 void ReSize(int m) { resize(m); }
800 void resize_keep(int);
801 void resize(const GeneralMatrix& A);
802 void ReSize(const GeneralMatrix& A) { resize(A); }
803 void operator+=(const SymmetricMatrix& M) { PlusEqual(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); }
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); }
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.
817class UpperTriangularMatrix : public GeneralMatrix
818{
819public:
820 UpperTriangularMatrix() {}
821 ~UpperTriangularMatrix() {}
822 UpperTriangularMatrix(ArrayLengthSpecifier);
823 void operator=(const BaseMatrix&);
824 void operator=(const UpperTriangularMatrix& m) { Eq(m); }
825 UpperTriangularMatrix(const BaseMatrix&);
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;
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&);
847 void RestoreCol(MatrixColX& c) { RestoreCol((MatrixRowCol&)c); }
848 void NextRow(MatrixRowCol&);
849 void resize(int); // change dimensions
850 void ReSize(int m) { resize(m); }
851 void 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); }
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); }
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); }
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.
870class LowerTriangularMatrix : public GeneralMatrix
871{
872public:
873 LowerTriangularMatrix() {}
874 ~LowerTriangularMatrix() {}
875 LowerTriangularMatrix(ArrayLengthSpecifier);
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;
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&);
899 void RestoreCol(MatrixColX& c) { RestoreCol((MatrixRowCol&)c); }
900 void NextRow(MatrixRowCol&);
901 void resize(int); // change dimensions
902 void ReSize(int m) { resize(m); }
903 void resize_keep(int);
904 void resize(const GeneralMatrix& A);
905 void ReSize(const GeneralMatrix& A) { resize(A); }
906 MatrixBandWidth bandwidth() const;
907 void operator+=(const LowerTriangularMatrix& M) { PlusEqual(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); }
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); }
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.
922class DiagonalMatrix : public GeneralMatrix
923{
924public:
925 DiagonalMatrix() {}
926 ~DiagonalMatrix() {}
927 DiagonalMatrix(ArrayLengthSpecifier);
928 DiagonalMatrix(const BaseMatrix&);
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&);
956 GeneralMatrix* MakeSolver() { return this; } // for solving
957 void Solver(MatrixColX&, const MatrixColX&);
958 GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
959 void resize(int); // change dimensions
960 void ReSize(int m) { resize(m); }
961 void resize_keep(int);
962 void resize(const GeneralMatrix& A);
963 void ReSize(const GeneralMatrix& A) { resize(A); }
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); }
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); }
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); }
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.
983class RowVector : public Matrix
984{
985public:
986 RowVector() { nrows_val = 1; }
987 ~RowVector() {}
988 RowVector(ArrayLengthSpecifier n) : Matrix(1,n.Value()) {}
989 RowVector(const BaseMatrix&);
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++
1002 RowVector(Real a, int n) : Matrix(a, 1, 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&);
1010 void RestoreCol(MatrixRowCol&) {}
1011 void RestoreCol(MatrixColX& c);
1012 GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
1013 void resize(int); // change dimensions
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
1018 void ReSize(int m,int n) { resize(m, n); }
1019 void resize(const GeneralMatrix& A);
1020 void ReSize(const GeneralMatrix& A) { resize(A); }
1021 Real* nric() const
1022 { CheckStore(); return store-1; } // for use by NRIC
1023 void cleanup(); // to clear store
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); }
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); }
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); }
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.
1042class ColumnVector : public Matrix
1043{
1044public:
1045 ColumnVector() { ncols_val = 1; }
1046 ~ColumnVector() {}
1047 ColumnVector(ArrayLengthSpecifier n) : Matrix(n.Value(),1) {}
1048 ColumnVector(const BaseMatrix&);
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++
1061 ColumnVector(Real a, int m) : Matrix(a, m, 1) {}
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
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
1071 void ReSize(int m,int n) { resize(m, n); }
1072 void resize(const GeneralMatrix& A);
1073 void ReSize(const GeneralMatrix& A) { resize(A); }
1074 Real* nric() const
1075 { CheckStore(); return store-1; } // for use by NRIC
1076 void cleanup(); // to clear store
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); }
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); }
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); }
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.
1097class 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
1104public:
1105 CroutMatrix(const BaseMatrix&);
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();
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&);
1118 void GetCol(MatrixColX& c) { GetCol((MatrixRowCol&)c); }
1119 void cleanup(); // to clear store
1120 void MiniCleanUp();
1121 bool IsEqual(const GeneralMatrix&) const;
1122 bool is_singular() const { return sing; }
1123 bool IsSingular() const { return sing; }
1124 const int* const_data_indx() const { return indx; }
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.
1134class BandMatrix : public GeneralMatrix
1135{
1136protected:
1137 void CornerClear() const; // set unused elements to zero
1138 short SimpleAddOK(const GeneralMatrix* gm);
1139public:
1140 int lower_val, upper_val; // band widths
1141 BandMatrix() { lower_val=0; upper_val=0; CornerClear(); }
1142 ~BandMatrix() {}
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
1159 BandMatrix(const BandMatrix& gm) : GeneralMatrix() { GetMatrix(&gm); }
1160 LogAndSign log_determinant() const;
1161 GeneralMatrix* MakeSolver();
1162 Real trace() const;
1163 Real sum_square() const
1164 { CornerClear(); return GeneralMatrix::sum_square(); }
1165 Real sum_absolute_value() const
1166 { CornerClear(); return GeneralMatrix::sum_absolute_value(); }
1167 Real sum() const
1168 { CornerClear(); return GeneralMatrix::sum(); }
1169 Real maximum_absolute_value() const
1170 { CornerClear(); return GeneralMatrix::maximum_absolute_value(); }
1171 Real minimum_absolute_value() const
1172 { int i, j; return GeneralMatrix::minimum_absolute_value2(i, j); }
1173 Real maximum() const { int i, j; return GeneralMatrix::maximum2(i, j); }
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&);
1179 void RestoreCol(MatrixColX& c) { RestoreCol((MatrixRowCol&)c); }
1180 void NextRow(MatrixRowCol&);
1181 virtual void resize(int, int, int); // change dimensions
1182 virtual void ReSize(int m, int n, int b) { resize(m, n, b); }
1183 void 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.
1205class UpperBandMatrix : public BandMatrix
1206{
1207public:
1208 UpperBandMatrix() {}
1209 ~UpperBandMatrix() {}
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;
1217 UpperBandMatrix(const UpperBandMatrix& gm) : BandMatrix() { GetMatrix(&gm); }
1218 GeneralMatrix* MakeSolver() { return this; }
1219 void Solver(MatrixColX&, const MatrixColX&);
1220 LogAndSign log_determinant() const;
1221 void resize(int, int, int); // change dimensions
1222 void ReSize(int m, int n, int b) { resize(m, n, b); }
1223 void resize(int n,int ubw) // change dimensions
1224 { BandMatrix::resize(n,0,ubw); }
1225 void ReSize(int n,int ubw) // change dimensions
1226 { BandMatrix::resize(n,0,ubw); }
1227 void resize(const GeneralMatrix& A) { BandMatrix::resize(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
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.
1244class LowerBandMatrix : public BandMatrix
1245{
1246public:
1247 LowerBandMatrix() {}
1248 ~LowerBandMatrix() {}
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;
1256 LowerBandMatrix(const LowerBandMatrix& gm) : BandMatrix() { GetMatrix(&gm); }
1257 GeneralMatrix* MakeSolver() { return this; }
1258 void Solver(MatrixColX&, const MatrixColX&);
1259 LogAndSign log_determinant() const;
1260 void resize(int, int, int); // change dimensions
1261 void ReSize(int m, int n, int b) { resize(m, n, b); }
1262 void resize(int n,int lbw) // change dimensions
1263 { BandMatrix::resize(n,lbw,0); }
1264 void ReSize(int n,int lbw) // change dimensions
1265 { BandMatrix::resize(n,lbw,0); }
1266 void resize(const GeneralMatrix& A) { BandMatrix::resize(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
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.
1283class SymmetricBandMatrix : public GeneralMatrix
1284{
1285 void CornerClear() const; // set unused elements to zero
1286 short SimpleAddOK(const GeneralMatrix* gm);
1287public:
1288 int lower_val; // lower band width
1289 SymmetricBandMatrix() { lower_val=0; CornerClear(); }
1290 ~SymmetricBandMatrix() {}
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;
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;
1311 Real maximum_absolute_value() const
1312 { CornerClear(); return GeneralMatrix::maximum_absolute_value(); }
1313 Real minimum_absolute_value() const
1314 { int i, j; return GeneralMatrix::minimum_absolute_value2(i, j); }
1315 Real maximum() const { int i, j; return GeneralMatrix::maximum2(i, j); }
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&);
1322 void RestoreCol(MatrixRowCol&) {}
1323 void RestoreCol(MatrixColX&);
1324 GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
1325 void resize(int,int); // change dimensions
1326 void ReSize(int m,int b) { resize(m, b); }
1327 void 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.
1344class 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
1354public:
1355 BandLUMatrix(const BaseMatrix&);
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();
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&);
1369 void GetCol(MatrixColX& c) { GetCol((MatrixRowCol&)c); }
1370 void cleanup(); // to clear store
1371 void MiniCleanUp();
1372 bool IsEqual(const GeneralMatrix&) const;
1373 bool is_singular() const { return sing; }
1374 bool IsSingular() const { return sing; }
1375 const Real* const_data2() const { return store2; }
1376 int size2() const { return storage2; }
1377 const int* const_data_indx() const { return indx; }
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.
1388class IdentityMatrix : public GeneralMatrix
1389{
1390public:
1391 IdentityMatrix() {}
1392 ~IdentityMatrix() {}
1393 IdentityMatrix(ArrayLengthSpecifier n) : GeneralMatrix(1)
1394 { nrows_val = ncols_val = n.Value(); *store = 1; }
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;
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&);
1414 GeneralMatrix* MakeSolver() { return this; } // for solving
1415 void Solver(MatrixColX&, const MatrixColX&);
1416 GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
1417 void resize(int n);
1418 void ReSize(int n) { resize(n); }
1419 void resize(const GeneralMatrix& A);
1420 void ReSize(const GeneralMatrix& A) { resize(A); }
1421 MatrixBandWidth bandwidth() const;
1422// ReturnMatrix Reverse() const; // reverse order of elements
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.
1435class GenericMatrix : public BaseMatrix
1436{
1437 GeneralMatrix* gm;
1438 int search(const BaseMatrix* bm) const;
1439 friend class BaseMatrix;
1440public:
1441 GenericMatrix() : gm(0) {}
1442 GenericMatrix(const BaseMatrix& bm)
1443 { gm = ((BaseMatrix&)bm).Evaluate(); gm = gm->Image(); }
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); }
1458 ~GenericMatrix() { delete gm; }
1459 void cleanup() { delete gm; gm = 0; }
1460 void Release() { gm->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
1472class MultipliedMatrix : public BaseMatrix
1473{
1474protected:
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; };
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;
1486public:
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
1495class AddedMatrix : public MultipliedMatrix
1496{
1497protected:
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;
1504public:
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
1513class SPMatrix : public AddedMatrix
1514{
1515protected:
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;
1522public:
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
1534class KPMatrix : public MultipliedMatrix
1535{
1536protected:
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;
1543public:
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
1553class ConcatenatedMatrix : public MultipliedMatrix
1554{
1555protected:
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;
1562public:
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
1571class StackedMatrix : public ConcatenatedMatrix
1572{
1573protected:
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;
1580public:
1581 ~StackedMatrix() {}
1582 GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1583 NEW_DELETE(StackedMatrix)
1584};
1585
1586/// Inverted matrix times matrix.
1587/// \internal
1588class SolvedMatrix : public MultipliedMatrix
1589{
1590 SolvedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
1591 : MultipliedMatrix(bm1x,bm2x) {}
1592 friend class BaseMatrix;
1593 friend class InvertedMatrix; // for operator*
1594public:
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
1603class SubtractedMatrix : public AddedMatrix
1604{
1605 SubtractedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
1606 : AddedMatrix(bm1x,bm2x) {}
1607 friend class BaseMatrix;
1608 friend class GeneralMatrix;
1609 friend class GenericMatrix;
1610public:
1611 ~SubtractedMatrix() {}
1612 GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1613 NEW_DELETE(SubtractedMatrix)
1614};
1615
1616/// Any type of matrix plus Real.
1617/// \internal
1618class ShiftedMatrix : public BaseMatrix
1619{
1620protected:
1621 union { const BaseMatrix* bm; GeneralMatrix* gm; };
1622 Real f;
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;
1628public:
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
1637class NegShiftedMatrix : public ShiftedMatrix
1638{
1639protected:
1640 NegShiftedMatrix(Real fx, const BaseMatrix* bmx) : ShiftedMatrix(bmx,fx) {}
1641 friend class BaseMatrix;
1642 friend class GeneralMatrix;
1643 friend class GenericMatrix;
1644public:
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
1653class ScaledMatrix : public ShiftedMatrix
1654{
1655 ScaledMatrix(const BaseMatrix* bmx, Real fx) : ShiftedMatrix(bmx,fx) {}
1656 friend class BaseMatrix;
1657 friend class GeneralMatrix;
1658 friend class GenericMatrix;
1659public:
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
1669class NegatedMatrix : public BaseMatrix
1670{
1671protected:
1672 union { const BaseMatrix* bm; GeneralMatrix* gm; };
1673 NegatedMatrix(const BaseMatrix* bmx) : bm(bmx) {}
1674 int search(const BaseMatrix*) const;
1675private:
1676 friend class BaseMatrix;
1677public:
1678 ~NegatedMatrix() {}
1679 GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1680 MatrixBandWidth bandwidth() const;
1681 NEW_DELETE(NegatedMatrix)
1682};
1683
1684/// Transposed matrix.
1685/// \internal
1686class TransposedMatrix : public NegatedMatrix
1687{
1688 TransposedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
1689 friend class BaseMatrix;
1690public:
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
1699class ReversedMatrix : public NegatedMatrix
1700{
1701 ReversedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
1702 friend class BaseMatrix;
1703public:
1704 ~ReversedMatrix() {}
1705 GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1706 NEW_DELETE(ReversedMatrix)
1707};
1708
1709/// Inverse of matrix.
1710/// \internal
1711class InvertedMatrix : public NegatedMatrix
1712{
1713 InvertedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
1714public:
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
1726class RowedMatrix : public NegatedMatrix
1727{
1728 RowedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
1729 friend class BaseMatrix;
1730public:
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
1739class ColedMatrix : public NegatedMatrix
1740{
1741 ColedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
1742 friend class BaseMatrix;
1743public:
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
1752class DiagedMatrix : public NegatedMatrix
1753{
1754 DiagedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
1755 friend class BaseMatrix;
1756public:
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
1765class MatedMatrix : public NegatedMatrix
1766{
1767 int nr, nc;
1768 MatedMatrix(const BaseMatrix* bmx, int nrx, int ncx)
1769 : NegatedMatrix(bmx), nr(nrx), nc(ncx) {}
1770 friend class BaseMatrix;
1771public:
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
1780class ReturnMatrix : public BaseMatrix
1781{
1782 GeneralMatrix* gm;
1783 int search(const BaseMatrix*) const;
1784public:
1785 ~ReturnMatrix() {}
1786 GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1787 friend class BaseMatrix;
1788 ReturnMatrix(const ReturnMatrix& tm) : BaseMatrix(), gm(tm.gm) {}
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
1800class GetSubMatrix : public NegatedMatrix
1801{
1802 int row_skip;
1803 int row_number;
1804 int col_skip;
1805 int col_number;
1806 bool IsSym;
1807
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;
1814public:
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) {}
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
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.
1849class LinearEquationSolver : public BaseMatrix
1850{
1851 GeneralMatrix* gm;
1852 int search(const BaseMatrix*) const { return 0; }
1853 friend class BaseMatrix;
1854public:
1855 LinearEquationSolver(const BaseMatrix& bm);
1856 ~LinearEquationSolver() { delete gm; }
1857 void cleanup() { delete gm; }
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
1869class MatrixInput
1870{
1871 int n; // number values still to be read
1872 Real* r; // pointer to next location to be read to
1873public:
1874 MatrixInput(const MatrixInput& mi) : n(mi.n), r(mi.r) {}
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
1892class SimpleIntArray : public Janitor
1893{
1894protected:
1895 int* a; ///< pointer to the array
1896 int n; ///< length of the array
1897public:
1898 SimpleIntArray(int xn); ///< build an array length xn
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
1909 int Size() const { return n; }
1910 ///< return the size of the array
1911 int size() const { return n; }
1912 ///< return the size of the array
1913 int* Data() { return a; } ///< pointer to the data
1914 const int* Data() const { return a; } ///< pointer to the data
1915 int* data() { return a; } ///< pointer to the data
1916 const int* data() const { return a; } ///< pointer to the 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
1922 void resize_keep(int i) { resize(i, true); }
1923 ///< change length, keep data
1924 void cleanup() { resize(0); } ///< set length to zero
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
1932class RealStarStar
1933{
1934 Real** a;
1935public:
1936 RealStarStar(Matrix& A);
1937 ~RealStarStar() { delete [] a; }
1938 operator Real**() { return a; }
1939};
1940
1941/// Let matrix simulate a C type const two dimensional array
1942class ConstRealStarStar
1943{
1944 const Real** a;
1945public:
1946 ConstRealStarStar(const Matrix& A);
1947 ~ConstRealStarStar() { delete [] a; }
1948 operator const Real**() { return a; }
1949};
1950
1951// *************************** exceptions ********************************/
1952
1953/// Not positive definite exception.
1954class NPDException : public Runtime_error
1955{
1956public:
1957 static unsigned long Select;
1958 NPDException(const GeneralMatrix&);
1959};
1960
1961/// Covergence failure exception.
1962class ConvergenceException : public Runtime_error
1963{
1964public:
1965 static unsigned long Select;
1966 ConvergenceException(const GeneralMatrix& A);
1967 ConvergenceException(const char* c);
1968};
1969
1970/// Singular matrix exception.
1971class SingularException : public Runtime_error
1972{
1973public:
1974 static unsigned long Select;
1975 SingularException(const GeneralMatrix& A);
1976};
1977
1978/// Real overflow exception.
1979class OverflowException : public Runtime_error
1980{
1981public:
1982 static unsigned long Select;
1983 OverflowException(const char* c);
1984};
1985
1986/// Miscellaneous exception (details in character string).
1987class ProgramException : public Logic_error
1988{
1989protected:
1990 ProgramException();
1991public:
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.
2000class IndexException : public Logic_error
2001{
2002public:
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.
2012class VectorException : public Logic_error
2013{
2014public:
2015 static unsigned long Select;
2016 VectorException();
2017 VectorException(const GeneralMatrix& A);
2018};
2019
2020/// A matrix is not square exception.
2021class NotSquareException : public Logic_error
2022{
2023public:
2024 static unsigned long Select;
2025 NotSquareException(const GeneralMatrix& A);
2026 NotSquareException();
2027};
2028
2029/// Submatrix dimension exception.
2030class SubMatrixDimensionException : public Logic_error
2031{
2032public:
2033 static unsigned long Select;
2034 SubMatrixDimensionException();
2035};
2036
2037/// Incompatible dimensions exception.
2038class IncompatibleDimensionsException : public Logic_error
2039{
2040public:
2041 static unsigned long Select;
2042 IncompatibleDimensionsException();
2043 IncompatibleDimensionsException(const GeneralMatrix&);
2044 IncompatibleDimensionsException(const GeneralMatrix&, const GeneralMatrix&);
2045};
2046
2047/// Not defined exception.
2048class NotDefinedException : public Logic_error
2049{
2050public:
2051 static unsigned long Select;
2052 NotDefinedException(const char* op, const char* matrix);
2053};
2054
2055/// Cannot build matrix with these properties exception.
2056class CannotBuildException : public Logic_error
2057{
2058public:
2059 static unsigned long Select;
2060 CannotBuildException(const char* matrix);
2061};
2062
2063
2064/// Internal newmat exception - shouldn't happen.
2065class InternalException : public Logic_error
2066{
2067public:
2068 static unsigned long Select; // for identifying exception
2069 InternalException(const char* c);
2070};
2071
2072// ************************ functions ************************************ //
2073
2074bool operator==(const GeneralMatrix& A, const GeneralMatrix& B);
2075bool operator==(const BaseMatrix& A, const BaseMatrix& B);
2076inline bool operator!=(const GeneralMatrix& A, const GeneralMatrix& B)
2077 { return ! (A==B); }
2078inline 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.
2083inline bool operator<=(const BaseMatrix& A, const BaseMatrix&)
2084 { A.IEQND(); return true; }
2085inline bool operator>=(const BaseMatrix& A, const BaseMatrix&)
2086 { A.IEQND(); return true; }
2087inline bool operator<(const BaseMatrix& A, const BaseMatrix&)
2088 { A.IEQND(); return true; }
2089inline bool operator>(const BaseMatrix& A, const BaseMatrix&)
2090 { A.IEQND(); return true; }
2091
2092bool is_zero(const BaseMatrix& A);
2093inline bool IsZero(const BaseMatrix& A) { return is_zero(A); }
2094
2095Real dotproduct(const Matrix& A, const Matrix& B);
2096Matrix crossproduct(const Matrix& A, const Matrix& B);
2097ReturnMatrix crossproduct_rows(const Matrix& A, const Matrix& B);
2098ReturnMatrix crossproduct_columns(const Matrix& A, const Matrix& B);
2099
2100inline Real DotProduct(const Matrix& A, const Matrix& B)
2101 { return dotproduct(A, B); }
2102inline Matrix CrossProduct(const Matrix& A, const Matrix& B)
2103 { return crossproduct(A, B); }
2104inline ReturnMatrix CrossProductRows(const Matrix& A, const Matrix& B)
2105 { return crossproduct_rows(A, B); }
2106inline ReturnMatrix CrossProductColumns(const Matrix& A, const Matrix& B)
2107 { return crossproduct_columns(A, B); }
2108
2109void 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
2115bool Rectangular(MatrixType a, MatrixType b, MatrixType c);
2116bool Compare(const MatrixType&, MatrixType&);
2117Real dotproduct(const Matrix& A, const Matrix& B);
2118SPMatrix SP(const BaseMatrix&, const BaseMatrix&);
2119KPMatrix KP(const BaseMatrix&, const BaseMatrix&);
2120ShiftedMatrix operator+(Real f, const BaseMatrix& BM);
2121NegShiftedMatrix operator-(Real, const BaseMatrix&);
2122ScaledMatrix operator*(Real f, const BaseMatrix& BM);
2123
2124// ********************* inline functions ******************************** //
2125
2126inline LogAndSign log_determinant(const BaseMatrix& B)
2127 { return B.log_determinant(); }
2128inline LogAndSign LogDeterminant(const BaseMatrix& B)
2129 { return B.log_determinant(); }
2130inline Real determinant(const BaseMatrix& B)
2131 { return B.determinant(); }
2132inline Real Determinant(const BaseMatrix& B)
2133 { return B.determinant(); }
2134inline Real sum_square(const BaseMatrix& B) { return B.sum_square(); }
2135inline Real SumSquare(const BaseMatrix& B) { return B.sum_square(); }
2136inline Real norm_Frobenius(const BaseMatrix& B) { return B.norm_Frobenius(); }
2137inline Real norm_frobenius(const BaseMatrix& B) { return B.norm_Frobenius(); }
2138inline Real NormFrobenius(const BaseMatrix& B) { return B.norm_Frobenius(); }
2139inline Real trace(const BaseMatrix& B) { return B.trace(); }
2140inline Real Trace(const BaseMatrix& B) { return B.trace(); }
2141inline Real sum_absolute_value(const BaseMatrix& B)
2142 { return B.sum_absolute_value(); }
2143inline Real SumAbsoluteValue(const BaseMatrix& B)
2144 { return B.sum_absolute_value(); }
2145inline Real sum(const BaseMatrix& B)
2146 { return B.sum(); }
2147inline Real Sum(const BaseMatrix& B)
2148 { return B.sum(); }
2149inline Real maximum_absolute_value(const BaseMatrix& B)
2150 { return B.maximum_absolute_value(); }
2151inline Real MaximumAbsoluteValue(const BaseMatrix& B)
2152 { return B.maximum_absolute_value(); }
2153inline Real minimum_absolute_value(const BaseMatrix& B)
2154 { return B.minimum_absolute_value(); }
2155inline Real MinimumAbsoluteValue(const BaseMatrix& B)
2156 { return B.minimum_absolute_value(); }
2157inline Real maximum(const BaseMatrix& B) { return B.maximum(); }
2158inline Real Maximum(const BaseMatrix& B) { return B.maximum(); }
2159inline Real minimum(const BaseMatrix& B) { return B.minimum(); }
2160inline Real Minimum(const BaseMatrix& B) { return B.minimum(); }
2161inline Real norm1(const BaseMatrix& B) { return B.norm1(); }
2162inline Real Norm1(const BaseMatrix& B) { return B.norm1(); }
2163inline Real norm1(RowVector& RV) { return RV.maximum_absolute_value(); }
2164inline Real Norm1(RowVector& RV) { return RV.maximum_absolute_value(); }
2165inline Real norm_infinity(const BaseMatrix& B) { return B.norm_infinity(); }
2166inline Real NormInfinity(const BaseMatrix& B) { return B.norm_infinity(); }
2167inline Real norm_infinity(ColumnVector& CV)
2168 { return CV.maximum_absolute_value(); }
2169inline Real NormInfinity(ColumnVector& CV)
2170 { return CV.maximum_absolute_value(); }
2171inline bool IsZero(const GeneralMatrix& A) { return A.IsZero(); }
2172inline bool is_zero(const GeneralMatrix& A) { return A.is_zero(); }
2173
2174
2175inline MatrixInput MatrixInput::operator<<(int f) { return *this << (Real)f; }
2176inline MatrixInput GeneralMatrix::operator<<(int f) { return *this << (Real)f; }
2177inline MatrixInput BandMatrix::operator<<(int f) { return *this << (Real)f; }
2178inline MatrixInput GetSubMatrix::operator<<(int f) { return *this << (Real)f; }
2179
2180inline ReversedMatrix BaseMatrix::Reverse() const { return reverse(); }
2181inline RowedMatrix BaseMatrix::AsRow() const { return as_row(); }
2182inline ColedMatrix BaseMatrix::AsColumn() const { return as_column(); }
2183inline DiagedMatrix BaseMatrix::AsDiagonal() const { return as_diagonal(); }
2184inline MatedMatrix BaseMatrix::AsMatrix(int m, int n) const
2185 { return as_matrix(m, n); }
2186inline GetSubMatrix BaseMatrix::SubMatrix(int fr, int lr, int fc, int lc) const
2187 { return submatrix(fr, lr, fc, lc); }
2188inline GetSubMatrix BaseMatrix::SymSubMatrix(int f, int l) const
2189 { return sym_submatrix(f, l); }
2190inline GetSubMatrix BaseMatrix::Row(int f) const { return row(f); }
2191inline GetSubMatrix BaseMatrix::Rows(int f, int l) const { return rows(f, l); }
2192inline GetSubMatrix BaseMatrix::Column(int f) const { return column(f); }
2193inline GetSubMatrix BaseMatrix::Columns(int f, int l) const
2194 { return columns(f, l); }
2195inline Real BaseMatrix::AsScalar() const { return as_scalar(); }
2196
2197inline ReturnMatrix GeneralMatrix::ForReturn() const { return for_return(); }
2198
2199inline void swap(Matrix& A, Matrix& B) { A.swap(B); }
2200inline void swap(SquareMatrix& A, SquareMatrix& B) { A.swap(B); }
2201inline void swap(nricMatrix& A, nricMatrix& B) { A.swap(B); }
2202inline void swap(UpperTriangularMatrix& A, UpperTriangularMatrix& B)
2203 { A.swap(B); }
2204inline void swap(LowerTriangularMatrix& A, LowerTriangularMatrix& B)
2205 { A.swap(B); }
2206inline void swap(SymmetricMatrix& A, SymmetricMatrix& B) { A.swap(B); }
2207inline void swap(DiagonalMatrix& A, DiagonalMatrix& B) { A.swap(B); }
2208inline void swap(RowVector& A, RowVector& B) { A.swap(B); }
2209inline void swap(ColumnVector& A, ColumnVector& B) { A.swap(B); }
2210inline void swap(CroutMatrix& A, CroutMatrix& B) { A.swap(B); }
2211inline void swap(BandMatrix& A, BandMatrix& B) { A.swap(B); }
2212inline void swap(UpperBandMatrix& A, UpperBandMatrix& B) { A.swap(B); }
2213inline void swap(LowerBandMatrix& A, LowerBandMatrix& B) { A.swap(B); }
2214inline void swap(SymmetricBandMatrix& A, SymmetricBandMatrix& B) { A.swap(B); }
2215inline void swap(BandLUMatrix& A, BandLUMatrix& B) { A.swap(B); }
2216inline void swap(IdentityMatrix& A, IdentityMatrix& B) { A.swap(B); }
2217inline void swap(GenericMatrix& A, GenericMatrix& B) { A.swap(B); }
2218
2219#ifdef OPT_COMPATIBLE // for compatibility with opt++
2220
2221inline Real Norm2(const ColumnVector& CV) { return CV.norm_Frobenius(); }
2222inline 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
Note: See TracBrowser for help on using the repository browser.