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

Last change on this file since 8901 was 8901, checked in by stuerze, 7 months ago

upgrade to newmat11 library

File size: 88.8 KB
Line 
1/// \defgroup newmat Newmat matrix manipulation library
2///@{
3
4/// \file newmat.h
5/// Definition file for matrix library.
6
7// Copyright (C) 2004: R B Davies
8
9#ifndef NEWMAT_LIB
10#define NEWMAT_LIB 0
11
12#include "include.h"
13
14#include "myexcept.h"
15
16
17#ifdef use_namespace
18namespace NEWMAT { using namespace RBD_COMMON; }
19namespace RBD_LIBRARIES { using namespace NEWMAT; }
20namespace NEWMAT {
21#endif
22
23//#define DO_REPORT                     // to activate REPORT
24
25#ifdef NO_LONG_NAMES
26#define UpperTriangularMatrix UTMatrix
27#define LowerTriangularMatrix LTMatrix
28#define SymmetricMatrix SMatrix
29#define DiagonalMatrix DMatrix
30#define BandMatrix BMatrix
31#define UpperBandMatrix UBMatrix
32#define LowerBandMatrix LBMatrix
33#define SymmetricBandMatrix SBMatrix
34#define BandLUMatrix BLUMatrix
35#endif
36
37// ************************** general utilities ****************************/
38
39class GeneralMatrix;                            // defined later
40class BaseMatrix;                               // defined later
41class MatrixInput;                              // defined later
42
43void MatrixErrorNoSpace(const void*);           ///< test for allocation fails
44
45/// Return from LogDeterminant function.
46/// Members are the log of the absolute value and the sign (+1, -1 or 0)
47class LogAndSign
48{
49   Real log_val;
50   int sign_val;
51public:
52   LogAndSign() { log_val=0.0; sign_val=1; }
53   LogAndSign(Real);
54   void operator*=(Real);                       ///< multiply by a real
55   void pow_eq(int k);                          ///< raise to power of k
56   void PowEq(int k) { pow_eq(k); }
57   void ChangeSign() { sign_val = -sign_val; }
58   void change_sign() { sign_val = -sign_val; } ///< change sign
59   Real LogValue() const { return log_val; }
60   Real log_value() const { return log_val; }   ///< log of the absolute value
61   int Sign() const { return sign_val; }
62   int sign() const { return sign_val; }        ///< sign of the value
63   Real value() const;                          ///< the value
64   Real Value() const { return value(); }
65   FREE_CHECK(LogAndSign)
66};
67
68// the following class is for counting the number of times a piece of code
69// is executed. It is used for locating any code not executed by test
70// routines. Use turbo GREP locate all places this code is called and
71// check which ones are not accessed.
72// Somewhat implementation dependent as it relies on "cout" still being
73// present when ExeCounter objects are destructed.
74
75#ifdef DO_REPORT
76
77class ExeCounter
78{
79   int line;                                    // code line number
80   int fileid;                                  // file identifier
81   long nexe;                                   // number of executions
82   static int nreports;                         // number of reports
83public:
84   ExeCounter(int,int);
85   void operator++() { nexe++; }
86   ~ExeCounter();                               // prints out reports
87};
88
89#endif
90
91
92// ************************** class MatrixType *****************************
93
94/// Find the type of a matrix resulting from matrix operations.
95/// Also identify what conversions are permissible.
96/// This class must be updated when new matrix types are added.
97
98class MatrixType
99{
100public:
101   enum Attribute {  Valid     = 1,
102                     Diagonal  = 2,             // order of these is important
103                     Symmetric = 4,
104                     Band      = 8,
105                     Lower     = 16,
106                     Upper     = 32,
107                     Square    = 64,
108                     Skew      = 128,
109                     LUDeco    = 256,
110                     Ones      = 512 };
111
112   enum            { US = 0,
113                     UT = Valid + Upper + Square,
114                     LT = Valid + Lower + Square,
115                     Rt = Valid,
116                     Sq = Valid + Square,
117                     Sm = Valid + Symmetric + Square,
118                     Sk = Valid + Skew + Square,
119                     Dg = Valid + Diagonal + Band + Lower + Upper + Symmetric
120                        + Square,
121                     Id = Valid + Diagonal + Band + Lower + Upper + Symmetric
122                        + Square + Ones,
123                     RV = Valid,     //   do not separate out
124                     CV = Valid,     //   vectors
125                     BM = Valid + Band + Square,
126                     UB = Valid + Band + Upper + Square,
127                     LB = Valid + Band + Lower + Square,
128                     SB = Valid + Band + Symmetric + Square,
129                     KB = Valid + Band + Skew + Square,
130                     Ct = Valid + LUDeco + Square,
131                     BC = Valid + Band + LUDeco + Square,
132                     Mask = ~Square
133                   };
134
135
136   static int nTypes() { return 13; }          // number of different types
137                                               // exclude Ct, US, BC
138public:
139   int attribute;
140   bool DataLossOK;                            // true if data loss is OK when
141                                               // this represents a destination
142public:
143   MatrixType () : DataLossOK(false) {}
144   MatrixType (int i) : attribute(i), DataLossOK(false) {}
145   MatrixType (int i, bool dlok) : attribute(i), DataLossOK(dlok) {}
146   MatrixType (const MatrixType& mt)
147      : attribute(mt.attribute), DataLossOK(mt.DataLossOK) {}
148   void operator=(const MatrixType& mt)
149      { attribute = mt.attribute; DataLossOK = mt.DataLossOK; }
150   void SetDataLossOK() { DataLossOK = true; }
151   int operator+() const { return attribute; }
152   MatrixType operator+(MatrixType mt) const
153      { return MatrixType(attribute & mt.attribute); }
154   MatrixType operator*(const MatrixType&) const;
155   MatrixType SP(const MatrixType&) const;
156   MatrixType KP(const MatrixType&) const;
157   MatrixType operator|(const MatrixType& mt) const
158      { return MatrixType(attribute & mt.attribute & Valid); }
159   MatrixType operator&(const MatrixType& mt) const
160      { return MatrixType(attribute & mt.attribute & Valid); }
161   bool operator>=(MatrixType mt) const
162      { return ( attribute & ~mt.attribute & Mask ) == 0; }
163   bool operator<(MatrixType mt) const         // for MS Visual C++ 4
164      { return ( attribute & ~mt.attribute & Mask ) != 0; }
165   bool operator==(MatrixType t) const
166      { return (attribute == t.attribute); }
167   bool operator!=(MatrixType t) const
168      { return (attribute != t.attribute); }
169   bool operator!() const { return (attribute & Valid) == 0; }
170   MatrixType i() const;                       ///< type of inverse
171   MatrixType t() const;                       ///< type of transpose
172   MatrixType AddEqualEl() const               ///< add constant to matrix
173      { return MatrixType(attribute & (Valid + Symmetric + Square)); }
174   MatrixType MultRHS() const;                 ///< type for rhs of multiply
175   MatrixType sub() const                      ///< type of submatrix
176      { return MatrixType(attribute & Valid); }
177   MatrixType ssub() const                     ///< type of sym submatrix
178      { return MatrixType(attribute); }        // not for selection matrix
179   GeneralMatrix* New() const;                 ///< new matrix of given type
180   GeneralMatrix* New(int,int,BaseMatrix*) const;
181                                               ///< new matrix of given type
182   const char* value() const;                  ///< type as char string
183   const char* Value() const { return value(); }
184   friend bool Rectangular(MatrixType a, MatrixType b, MatrixType c);
185   friend bool Compare(const MatrixType&, MatrixType&);
186                                               ///< compare and check conversion
187   bool is_band() const { return (attribute & Band) != 0; }
188   bool is_diagonal() const { return (attribute & Diagonal) != 0; }
189   bool is_symmetric() const { return (attribute & Symmetric) != 0; }
190   bool CannotConvert() const { return (attribute & LUDeco) != 0; }
191                                               // used by operator==
192   FREE_CHECK(MatrixType)
193};
194
195
196// *********************** class MatrixBandWidth ***********************/
197
198///Upper and lower bandwidths of a matrix.
199///That is number of diagonals strictly above or below main diagonal,
200///e.g. diagonal matrix has 0 upper and lower bandwiths.
201///-1 means the matrix may have the maximum bandwidth.
202class MatrixBandWidth
203{
204public:
205   int lower_val;
206   int upper_val;
207   MatrixBandWidth(const int l, const int u) : lower_val(l), upper_val(u) {}
208   MatrixBandWidth(const int i) : lower_val(i), upper_val(i) {}
209   MatrixBandWidth operator+(const MatrixBandWidth&) const;
210   MatrixBandWidth operator*(const MatrixBandWidth&) const;
211   MatrixBandWidth minimum(const MatrixBandWidth&) const;
212   MatrixBandWidth t() const { return MatrixBandWidth(upper_val,lower_val); }
213   bool operator==(const MatrixBandWidth& bw) const
214      { return (lower_val == bw.lower_val) && (upper_val == bw.upper_val); }
215   bool operator!=(const MatrixBandWidth& bw) const { return !operator==(bw); }
216   int Upper() const { return upper_val; }
217   int upper() const { return upper_val; }
218   int Lower() const { return lower_val; }
219   int lower() const { return lower_val; }
220   FREE_CHECK(MatrixBandWidth)
221};
222
223
224// ********************* Array length specifier ************************/
225
226/// This class is used to avoid constructors such as
227/// ColumnVector(int) being used for conversions.
228/// Eventually this should be replaced by the use of the keyword "explicit".
229
230class ArrayLengthSpecifier
231{
232   int v;
233public:
234   int Value() const { return v; }
235   int value() const { return v; }
236   ArrayLengthSpecifier(int l) : v(l) {}
237};
238
239// ************************* Matrix routines ***************************/
240
241
242class MatrixRowCol;                             // defined later
243class MatrixRow;
244class MatrixCol;
245class MatrixColX;
246
247class GeneralMatrix;                            // defined later
248class AddedMatrix;
249class MultipliedMatrix;
250class SubtractedMatrix;
251class SPMatrix;
252class KPMatrix;
253class ConcatenatedMatrix;
254class StackedMatrix;
255class SolvedMatrix;
256class ShiftedMatrix;
257class NegShiftedMatrix;
258class ScaledMatrix;
259class TransposedMatrix;
260class ReversedMatrix;
261class NegatedMatrix;
262class InvertedMatrix;
263class RowedMatrix;
264class ColedMatrix;
265class DiagedMatrix;
266class MatedMatrix;
267class GetSubMatrix;
268class ReturnMatrix;
269class Matrix;
270class SquareMatrix;
271class nricMatrix;
272class RowVector;
273class ColumnVector;
274class SymmetricMatrix;
275class UpperTriangularMatrix;
276class LowerTriangularMatrix;
277class DiagonalMatrix;
278class CroutMatrix;
279class BandMatrix;
280class LowerBandMatrix;
281class UpperBandMatrix;
282class SymmetricBandMatrix;
283class LinearEquationSolver;
284class GenericMatrix;
285
286
287#define MatrixTypeUnSp 0
288//static MatrixType MatrixTypeUnSp(MatrixType::US);
289//                                              // AT&T needs this
290
291/// Base of the matrix classes.
292class BaseMatrix : public Janitor
293{
294protected:
295   virtual int search(const BaseMatrix*) const = 0;
296                                                // count number of times matrix is referred to
297public:
298   virtual GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp) = 0;
299                                                // evaluate temporary
300   // for old version of G++
301   //   virtual GeneralMatrix* Evaluate(MatrixType mt) = 0;
302   //   GeneralMatrix* Evaluate() { return Evaluate(MatrixTypeUnSp); }
303   AddedMatrix operator+(const BaseMatrix&) const;    // results of operations
304   MultipliedMatrix operator*(const BaseMatrix&) const;
305   SubtractedMatrix operator-(const BaseMatrix&) const;
306   ConcatenatedMatrix operator|(const BaseMatrix&) const;
307   StackedMatrix operator&(const BaseMatrix&) const;
308   ShiftedMatrix operator+(Real) const;
309   ScaledMatrix operator*(Real) const;
310   ScaledMatrix operator/(Real) const;
311   ShiftedMatrix operator-(Real) const;
312   TransposedMatrix t() const;
313//   TransposedMatrix t;
314   NegatedMatrix operator-() const;                   // change sign of elements
315   ReversedMatrix reverse() const;
316   ReversedMatrix Reverse() const;
317   InvertedMatrix i() const;
318//   InvertedMatrix i;
319   RowedMatrix as_row() const;
320   RowedMatrix AsRow() const;
321   ColedMatrix as_column() const;
322   ColedMatrix AsColumn() const;
323   DiagedMatrix as_diagonal() const;
324   DiagedMatrix AsDiagonal() const;
325   MatedMatrix as_matrix(int,int) const;
326   MatedMatrix AsMatrix(int m, int n) const;
327   GetSubMatrix submatrix(int,int,int,int) const;
328   GetSubMatrix SubMatrix(int fr, int lr, int fc, int lc) const;
329   GetSubMatrix sym_submatrix(int,int) const;
330   GetSubMatrix SymSubMatrix(int f, int l) const;
331   GetSubMatrix row(int) const;
332   GetSubMatrix rows(int,int) const;
333   GetSubMatrix column(int) const;
334   GetSubMatrix columns(int,int) const;
335   GetSubMatrix Row(int f) const;
336   GetSubMatrix Rows(int f, int l) const;
337   GetSubMatrix Column(int f) const;
338   GetSubMatrix Columns(int f, int l) const;
339   Real as_scalar() const;                      // conversion of 1 x 1 matrix
340   Real AsScalar() const;
341   virtual LogAndSign log_determinant() const;
342   LogAndSign LogDeterminant() const { return log_determinant(); }
343   Real determinant() const;
344   Real Determinant() const { return determinant(); }
345   virtual Real sum_square() const;
346   Real SumSquare() const { return sum_square(); }
347   Real norm_Frobenius() const;
348   Real norm_frobenius() const { return norm_Frobenius(); }
349   Real NormFrobenius() const { return norm_Frobenius(); }
350   virtual Real sum_absolute_value() const;
351   Real SumAbsoluteValue() const { return sum_absolute_value(); }
352   virtual Real sum() const;
353   virtual Real Sum() const { return sum(); }
354   virtual Real maximum_absolute_value() const;
355   Real MaximumAbsoluteValue() const { return maximum_absolute_value(); }
356   virtual Real maximum_absolute_value1(int& i) const;
357   Real MaximumAbsoluteValue1(int& i) const
358      { return maximum_absolute_value1(i); }
359   virtual Real maximum_absolute_value2(int& i, int& j) const;
360   Real MaximumAbsoluteValue2(int& i, int& j) const
361      { return maximum_absolute_value2(i,j); }
362   virtual Real minimum_absolute_value() const;
363   Real MinimumAbsoluteValue() const { return minimum_absolute_value(); }
364   virtual Real minimum_absolute_value1(int& i) const;
365   Real MinimumAbsoluteValue1(int& i) const
366      { return minimum_absolute_value1(i); }
367   virtual Real minimum_absolute_value2(int& i, int& j) const;
368   Real MinimumAbsoluteValue2(int& i, int& j) const
369      { return minimum_absolute_value2(i,j); }
370   virtual Real maximum() const;
371   Real Maximum() const { return maximum(); }
372   virtual Real maximum1(int& i) const;
373   Real Maximum1(int& i) const { return maximum1(i); }
374   virtual Real maximum2(int& i, int& j) const;
375   Real Maximum2(int& i, int& j) const { return maximum2(i,j); }
376   virtual Real minimum() const;
377   Real Minimum() const { return minimum(); }
378   virtual Real minimum1(int& i) const;
379   Real Minimum1(int& i) const { return minimum1(i); }
380   virtual Real minimum2(int& i, int& j) const;
381   Real Minimum2(int& i, int& j) const { return minimum2(i,j); }
382   virtual Real trace() const;
383   Real Trace() const { return trace(); }
384   Real norm1() const;
385   Real Norm1() const { return norm1(); }
386   Real norm_infinity() const;
387   Real NormInfinity() const { return norm_infinity(); }
388   virtual MatrixBandWidth bandwidth() const;  // bandwidths of band matrix
389   virtual MatrixBandWidth BandWidth() const { return bandwidth(); }
390   void IEQND() const;                         // called by ineq. ops
391   ReturnMatrix sum_square_columns() const;
392   ReturnMatrix sum_square_rows() const;
393   ReturnMatrix sum_columns() const;
394   ReturnMatrix sum_rows() const;
395   virtual void cleanup() {}
396   void CleanUp() { cleanup(); }
397
398//   virtual ReturnMatrix Reverse() const;       // reverse order of elements
399//protected:
400//   BaseMatrix() : t(this), i(this) {}
401
402   friend class GeneralMatrix;
403   friend class Matrix;
404   friend class SquareMatrix;
405   friend class nricMatrix;
406   friend class RowVector;
407   friend class ColumnVector;
408   friend class SymmetricMatrix;
409   friend class UpperTriangularMatrix;
410   friend class LowerTriangularMatrix;
411   friend class DiagonalMatrix;
412   friend class CroutMatrix;
413   friend class BandMatrix;
414   friend class LowerBandMatrix;
415   friend class UpperBandMatrix;
416   friend class SymmetricBandMatrix;
417   friend class AddedMatrix;
418   friend class MultipliedMatrix;
419   friend class SubtractedMatrix;
420   friend class SPMatrix;
421   friend class KPMatrix;
422   friend class ConcatenatedMatrix;
423   friend class StackedMatrix;
424   friend class SolvedMatrix;
425   friend class ShiftedMatrix;
426   friend class NegShiftedMatrix;
427   friend class ScaledMatrix;
428   friend class TransposedMatrix;
429   friend class ReversedMatrix;
430   friend class NegatedMatrix;
431   friend class InvertedMatrix;
432   friend class RowedMatrix;
433   friend class ColedMatrix;
434   friend class DiagedMatrix;
435   friend class MatedMatrix;
436   friend class GetSubMatrix;
437   friend class ReturnMatrix;
438   friend class LinearEquationSolver;
439   friend class GenericMatrix;
440   NEW_DELETE(BaseMatrix)
441};
442
443
444// ***************************** working classes **************************/
445
446/// The classes for matrices that can contain data are derived from this.
447class GeneralMatrix : public BaseMatrix         // declarable matrix types
448{
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.