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

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

* empty log message *

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