Planet
navi homePPSaboutscreenshotsdownloaddevelopmentforum

source: orxonox.OLD/orxonox/trunk/src/lib/newmat/newmat.h @ 4639

Last change on this file since 4639 was 4565, checked in by patrick, 19 years ago

orxonox/trunk: added the newmat library to the project. needs some translation in directory, temp under util/newmat. is needed by the collision detection engine to perform lin alg operations such as eigenvector decomposition. perhaps we will make our own library to do that later.
File size: 64.8 KB
Line 
1//$$ newmat.h           definition file for new version of matrix package
2
3// Copyright (C) 1991,2,3,4,7,2000,2002: R B Davies
4
5#ifndef NEWMAT_LIB
6#define NEWMAT_LIB 0
7
8#include "include.h"
9
10#include "boolean.h"
11#include "myexcept.h"
12
13
14#ifdef use_namespace
15namespace NEWMAT { using namespace RBD_COMMON; }
16namespace RBD_LIBRARIES { using namespace NEWMAT; }
17namespace NEWMAT {
18#endif
19
20//#define DO_REPORT                     // to activate REPORT
21
22#ifdef NO_LONG_NAMES
23#define UpperTriangularMatrix UTMatrix
24#define LowerTriangularMatrix LTMatrix
25#define SymmetricMatrix SMatrix
26#define DiagonalMatrix DMatrix
27#define BandMatrix BMatrix
28#define UpperBandMatrix UBMatrix
29#define LowerBandMatrix LBMatrix
30#define SymmetricBandMatrix SBMatrix
31#define BandLUMatrix BLUMatrix
32#endif
33
34#ifndef TEMPS_DESTROYED_QUICKLY_R
35#define ReturnMatrix ReturnMatrixX
36#else
37#define ReturnMatrix ReturnMatrixX&
38#endif
39
40// ************************** general utilities ****************************/
41
42class GeneralMatrix;
43
44void MatrixErrorNoSpace(void*);                 // no space handler
45
46class LogAndSign
47// Return from LogDeterminant function
48//    - value of the log plus the sign (+, - or 0)
49{
50   Real log_value;
51   int sign;
52public:
53   LogAndSign() { log_value=0.0; sign=1; }
54   LogAndSign(Real);
55   void operator*=(Real);
56   void PowEq(int k);  // raise to power of k
57   void ChangeSign() { sign = -sign; }
58   Real LogValue() const { return log_value; }
59   int Sign() const { return sign; }
60   Real Value() const;
61   FREE_CHECK(LogAndSign)
62};
63
64// the following class is for counting the number of times a piece of code
65// is executed. It is used for locating any code not executed by test
66// routines. Use turbo GREP locate all places this code is called and
67// check which ones are not accessed.
68// Somewhat implementation dependent as it relies on "cout" still being
69// present when ExeCounter objects are destructed.
70
71#ifdef DO_REPORT
72
73class ExeCounter
74{
75   int line;                                    // code line number
76   int fileid;                                  // file identifier
77   long nexe;                                   // number of executions
78   static int nreports;                         // number of reports
79public:
80   ExeCounter(int,int);
81   void operator++() { nexe++; }
82   ~ExeCounter();                               // prints out reports
83};
84
85#endif
86
87
88// ************************** class MatrixType *****************************/
89
90// Is used for finding the type of a matrix resulting from the binary operations
91// +, -, * and identifying what conversions are permissible.
92// This class must be updated when new matrix types are added.
93
94class GeneralMatrix;                            // defined later
95class BaseMatrix;                               // defined later
96class MatrixInput;                              // defined later
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                     LUDeco    = 64,
108                     Ones      = 128 };
109
110   enum            { US = 0,
111                     UT = Valid + Upper,
112                     LT = Valid + Lower,
113                     Rt = Valid,
114                     Sm = Valid + Symmetric,
115                     Dg = Valid + Diagonal + Band + Lower + Upper + Symmetric,
116                     Id = Valid + Diagonal + Band + Lower + Upper + Symmetric
117                        + Ones,
118                     RV = Valid,     //   do not separate out
119                     CV = Valid,     //   vectors
120                     BM = Valid + Band,
121                     UB = Valid + Band + Upper,
122                     LB = Valid + Band + Lower,
123                     SB = Valid + Band + Symmetric,
124                     Ct = Valid + LUDeco,
125                     BC = Valid + Band + LUDeco
126                   };
127
128
129   static int nTypes() { return 10; }           // number of different types
130                                               // exclude Ct, US, BC
131public:
132   int attribute;
133   bool DataLossOK;                            // true if data loss is OK when
134                                               // this represents a destination
135public:
136   MatrixType () : DataLossOK(false) {}
137   MatrixType (int i) : attribute(i), DataLossOK(false) {}
138   MatrixType (int i, bool dlok) : attribute(i), DataLossOK(dlok) {}
139   MatrixType (const MatrixType& mt)
140      : attribute(mt.attribute), DataLossOK(mt.DataLossOK) {}
141   void operator=(const MatrixType& mt)
142      { attribute = mt.attribute; DataLossOK = mt.DataLossOK; }
143   void SetDataLossOK() { DataLossOK = true; }
144   int operator+() const { return attribute; }
145   MatrixType operator+(MatrixType mt) const
146      { return MatrixType(attribute & mt.attribute); }
147   MatrixType operator*(const MatrixType&) const;
148   MatrixType SP(const MatrixType&) const;
149   MatrixType KP(const MatrixType&) const;
150   MatrixType operator|(const MatrixType& mt) const
151      { return MatrixType(attribute & mt.attribute & Valid); }
152   MatrixType operator&(const MatrixType& mt) const
153      { return MatrixType(attribute & mt.attribute & Valid); }
154   bool operator>=(MatrixType mt) const
155      { return ( attribute & mt.attribute ) == attribute; }
156   bool operator<(MatrixType mt) const         // for MS Visual C++ 4
157      { return ( attribute & mt.attribute ) != attribute; }
158   bool operator==(MatrixType t) const
159      { return (attribute == t.attribute); }
160   bool operator!=(MatrixType t) const
161      { return (attribute != t.attribute); }
162   bool operator!() const { return (attribute & Valid) == 0; }
163   MatrixType i() const;                       // type of inverse
164   MatrixType t() const;                       // type of transpose
165   MatrixType AddEqualEl() const               // Add constant to matrix
166      { return MatrixType(attribute & (Valid + Symmetric)); }
167   MatrixType MultRHS() const;                 // type for rhs of multiply
168   MatrixType sub() const                      // type of submatrix
169      { return MatrixType(attribute & Valid); }
170   MatrixType ssub() const                     // type of sym submatrix
171      { return MatrixType(attribute); }        // not for selection matrix
172   GeneralMatrix* New() const;                 // new matrix of given type
173   GeneralMatrix* New(int,int,BaseMatrix*) const;
174                                               // new matrix of given type
175   const char* Value() const;                  // to print type
176   friend bool Rectangular(MatrixType a, MatrixType b, MatrixType c);
177   friend bool Compare(const MatrixType&, MatrixType&);
178                                               // compare and check conv.
179   bool IsBand() const { return (attribute & Band) != 0; }
180   bool IsDiagonal() const { return (attribute & Diagonal) != 0; }
181   bool IsSymmetric() const { return (attribute & Symmetric) != 0; }
182   bool CannotConvert() const { return (attribute & LUDeco) != 0; }
183                                               // used by operator==
184   FREE_CHECK(MatrixType)
185};
186
187
188// *********************** class MatrixBandWidth ***********************/
189
190class MatrixBandWidth
191{
192public:
193   int lower;
194   int upper;
195   MatrixBandWidth(const int l, const int u) : lower(l), upper (u) {}
196   MatrixBandWidth(const int i) : lower(i), upper(i) {}
197   MatrixBandWidth operator+(const MatrixBandWidth&) const;
198   MatrixBandWidth operator*(const MatrixBandWidth&) const;
199   MatrixBandWidth minimum(const MatrixBandWidth&) const;
200   MatrixBandWidth t() const { return MatrixBandWidth(upper,lower); }
201   bool operator==(const MatrixBandWidth& bw) const
202      { return (lower == bw.lower) && (upper == bw.upper); }
203   bool operator!=(const MatrixBandWidth& bw) const { return !operator==(bw); }
204   int Upper() const { return upper; }
205   int Lower() const { return lower; }
206   FREE_CHECK(MatrixBandWidth)
207};
208
209
210// ********************* Array length specifier ************************/
211
212// This class is introduced to avoid constructors such as
213//   ColumnVector(int)
214// being used for conversions
215
216class ArrayLengthSpecifier
217{
218   int value;
219public:
220   int Value() const { return value; }
221   ArrayLengthSpecifier(int l) : value(l) {}
222};
223
224// ************************* Matrix routines ***************************/
225
226
227class MatrixRowCol;                             // defined later
228class MatrixRow;
229class MatrixCol;
230class MatrixColX;
231
232class GeneralMatrix;                            // defined later
233class AddedMatrix;
234class MultipliedMatrix;
235class SubtractedMatrix;
236class SPMatrix;
237class KPMatrix;
238class ConcatenatedMatrix;
239class StackedMatrix;
240class SolvedMatrix;
241class ShiftedMatrix;
242class NegShiftedMatrix;
243class ScaledMatrix;
244class TransposedMatrix;
245class ReversedMatrix;
246class NegatedMatrix;
247class InvertedMatrix;
248class RowedMatrix;
249class ColedMatrix;
250class DiagedMatrix;
251class MatedMatrix;
252class GetSubMatrix;
253class ReturnMatrixX;
254class Matrix;
255class nricMatrix;
256class RowVector;
257class ColumnVector;
258class SymmetricMatrix;
259class UpperTriangularMatrix;
260class LowerTriangularMatrix;
261class DiagonalMatrix;
262class CroutMatrix;
263class BandMatrix;
264class LowerBandMatrix;
265class UpperBandMatrix;
266class SymmetricBandMatrix;
267class LinearEquationSolver;
268class GenericMatrix;
269
270
271#define MatrixTypeUnSp 0
272//static MatrixType MatrixTypeUnSp(MatrixType::US);
273//                                              // AT&T needs this
274
275class BaseMatrix : public Janitor               // base of all matrix classes
276{
277protected:
278   virtual int search(const BaseMatrix*) const = 0;
279                                                // count number of times matrix
280                                                // is referred to
281
282public:
283   virtual GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp) = 0;
284                                                // evaluate temporary
285   // for old version of G++
286   //   virtual GeneralMatrix* Evaluate(MatrixType mt) = 0;
287   //   GeneralMatrix* Evaluate() { return Evaluate(MatrixTypeUnSp); }
288#ifndef TEMPS_DESTROYED_QUICKLY
289   AddedMatrix operator+(const BaseMatrix&) const;    // results of operations
290   MultipliedMatrix operator*(const BaseMatrix&) const;
291   SubtractedMatrix operator-(const BaseMatrix&) const;
292   ConcatenatedMatrix operator|(const BaseMatrix&) const;
293   StackedMatrix operator&(const BaseMatrix&) const;
294   ShiftedMatrix operator+(Real) const;
295   ScaledMatrix operator*(Real) const;
296   ScaledMatrix operator/(Real) const;
297   ShiftedMatrix operator-(Real) const;
298   TransposedMatrix t() const;
299//   TransposedMatrix t;
300   NegatedMatrix operator-() const;                   // change sign of elements
301   ReversedMatrix Reverse() const;
302   InvertedMatrix i() const;
303//   InvertedMatrix i;
304   RowedMatrix AsRow() const;
305   ColedMatrix AsColumn() const;
306   DiagedMatrix AsDiagonal() const;
307   MatedMatrix AsMatrix(int,int) const;
308   GetSubMatrix SubMatrix(int,int,int,int) const;
309   GetSubMatrix SymSubMatrix(int,int) const;
310   GetSubMatrix Row(int) const;
311   GetSubMatrix Rows(int,int) const;
312   GetSubMatrix Column(int) const;
313   GetSubMatrix Columns(int,int) const;
314#else
315   AddedMatrix& operator+(const BaseMatrix&) const;    // results of operations
316   MultipliedMatrix& operator*(const BaseMatrix&) const;
317   SubtractedMatrix& operator-(const BaseMatrix&) const;
318   ConcatenatedMatrix& operator|(const BaseMatrix&) const;
319   StackedMatrix& operator&(const BaseMatrix&) const;
320   ShiftedMatrix& operator+(Real) const;
321   ScaledMatrix& operator*(Real) const;
322   ScaledMatrix& operator/(Real) const;
323   ShiftedMatrix& operator-(Real) const;
324   TransposedMatrix& t() const;
325//   TransposedMatrix& t;
326   NegatedMatrix& operator-() const;                   // change sign of elements
327   ReversedMatrix& Reverse() const;
328   InvertedMatrix& i() const;
329//   InvertedMatrix& i;
330   RowedMatrix& AsRow() const;
331   ColedMatrix& AsColumn() const;
332   DiagedMatrix& AsDiagonal() const;
333   MatedMatrix& AsMatrix(int,int) const;
334   GetSubMatrix& SubMatrix(int,int,int,int) const;
335   GetSubMatrix& SymSubMatrix(int,int) const;
336   GetSubMatrix& Row(int) const;
337   GetSubMatrix& Rows(int,int) const;
338   GetSubMatrix& Column(int) const;
339   GetSubMatrix& Columns(int,int) const;
340#endif
341   Real AsScalar() const;                      // conversion of 1 x 1 matrix
342   virtual LogAndSign LogDeterminant() const;
343   Real Determinant() const;
344   virtual Real SumSquare() const;
345   Real NormFrobenius() const;
346   virtual Real SumAbsoluteValue() const;
347   virtual Real Sum() const;
348   virtual Real MaximumAbsoluteValue() const;
349   virtual Real MaximumAbsoluteValue1(int& i) const;
350   virtual Real MaximumAbsoluteValue2(int& i, int& j) const;
351   virtual Real MinimumAbsoluteValue() const;
352   virtual Real MinimumAbsoluteValue1(int& i) const;
353   virtual Real MinimumAbsoluteValue2(int& i, int& j) const;
354   virtual Real Maximum() const;
355   virtual Real Maximum1(int& i) const;
356   virtual Real Maximum2(int& i, int& j) const;
357   virtual Real Minimum() const;
358   virtual Real Minimum1(int& i) const;
359   virtual Real Minimum2(int& i, int& j) const;
360   virtual Real Trace() const;
361   Real Norm1() const;
362   Real NormInfinity() const;
363   virtual MatrixBandWidth BandWidth() const;  // bandwidths of band matrix
364   virtual void CleanUp() {}                   // to clear store
365   void IEQND() const;                         // called by ineq. ops
366//   virtual ReturnMatrix Reverse() const;       // reverse order of elements
367
368
369//protected:
370//   BaseMatrix() : t(this), i(this) {}
371
372   friend class GeneralMatrix;
373   friend class Matrix;
374   friend class nricMatrix;
375   friend class RowVector;
376   friend class ColumnVector;
377   friend class SymmetricMatrix;
378   friend class UpperTriangularMatrix;
379   friend class LowerTriangularMatrix;
380   friend class DiagonalMatrix;
381   friend class CroutMatrix;
382   friend class BandMatrix;
383   friend class LowerBandMatrix;
384   friend class UpperBandMatrix;
385   friend class SymmetricBandMatrix;
386   friend class AddedMatrix;
387   friend class MultipliedMatrix;
388   friend class SubtractedMatrix;
389   friend class SPMatrix;
390   friend class KPMatrix;
391   friend class ConcatenatedMatrix;
392   friend class StackedMatrix;
393   friend class SolvedMatrix;
394   friend class ShiftedMatrix;
395   friend class NegShiftedMatrix;
396   friend class ScaledMatrix;
397   friend class TransposedMatrix;
398   friend class ReversedMatrix;
399   friend class NegatedMatrix;
400   friend class InvertedMatrix;
401   friend class RowedMatrix;
402   friend class ColedMatrix;
403   friend class DiagedMatrix;
404   friend class MatedMatrix;
405   friend class GetSubMatrix;
406   friend class ReturnMatrixX;
407   friend class LinearEquationSolver;
408   friend class GenericMatrix;
409   NEW_DELETE(BaseMatrix)
410};
411
412
413// ***************************** working classes **************************/
414
415class GeneralMatrix : public BaseMatrix         // declarable matrix types
416{
417   virtual GeneralMatrix* Image() const;        // copy of matrix
418protected:
419   int tag;                                     // shows whether can reuse
420   int nrows, ncols;                            // dimensions
421   int storage;                                 // total store required
422   Real* store;                                 // point to store (0=not set)
423   GeneralMatrix();                             // initialise with no store
424   GeneralMatrix(ArrayLengthSpecifier);         // constructor getting store
425   void Add(GeneralMatrix*, Real);              // sum of GM and Real
426   void Add(Real);                              // add Real to this
427   void NegAdd(GeneralMatrix*, Real);           // Real - GM
428   void NegAdd(Real);                           // this = this - Real
429   void Multiply(GeneralMatrix*, Real);         // product of GM and Real
430   void Multiply(Real);                         // multiply this by Real
431   void Negate(GeneralMatrix*);                 // change sign
432   void Negate();                               // change sign
433   void ReverseElements();                      // internal reverse of elements
434   void ReverseElements(GeneralMatrix*);        // reverse order of elements
435   void operator=(Real);                        // set matrix to constant
436   Real* GetStore();                            // get store or copy
437   GeneralMatrix* BorrowStore(GeneralMatrix*, MatrixType);
438                                                // temporarily access store
439   void GetMatrix(const GeneralMatrix*);        // used by = and initialise
440   void Eq(const BaseMatrix&, MatrixType);      // used by =
441   void Eq(const BaseMatrix&, MatrixType, bool);// used by <<
442   void Eq2(const BaseMatrix&, MatrixType);     // cut down version of Eq
443   int search(const BaseMatrix*) const;
444   virtual GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
445   void CheckConversion(const BaseMatrix&);     // check conversion OK
446   void ReSize(int, int, int);                  // change dimensions
447   virtual short SimpleAddOK(const GeneralMatrix* gm) { return 0; }
448             // see bandmat.cpp for explanation
449public:
450   GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
451   virtual MatrixType Type() const = 0;         // type of a matrix
452   int Nrows() const { return nrows; }          // get dimensions
453   int Ncols() const { return ncols; }
454   int Storage() const { return storage; }
455   Real* Store() const { return store; }
456   virtual ~GeneralMatrix();                    // delete store if set
457   void tDelete();                              // delete if tag permits
458   bool reuse();                                // true if tag allows reuse
459   void Protect() { tag=-1; }                   // cannot delete or reuse
460   int Tag() const { return tag; }
461   bool IsZero() const;                         // test matrix has all zeros
462   void Release() { tag=1; }                    // del store after next use
463   void Release(int t) { tag=t; }               // del store after t accesses
464   void ReleaseAndDelete() { tag=0; }           // delete matrix after use
465   void operator<<(const Real*);                // assignment from an array
466   void operator<<(const BaseMatrix& X)
467      { Eq(X,this->Type(),true); }              // = without checking type
468   void Inject(const GeneralMatrix&);           // copy stored els only
469   void operator+=(const BaseMatrix&);
470   void operator-=(const BaseMatrix&);
471   void operator*=(const BaseMatrix&);
472   void operator|=(const BaseMatrix&);
473   void operator&=(const BaseMatrix&);
474   void operator+=(Real);
475   void operator-=(Real r) { operator+=(-r); }
476   void operator*=(Real);
477   void operator/=(Real r) { operator*=(1.0/r); }
478   virtual GeneralMatrix* MakeSolver();         // for solving
479   virtual void Solver(MatrixColX&, const MatrixColX&) {}
480   virtual void GetRow(MatrixRowCol&) = 0;      // Get matrix row
481   virtual void RestoreRow(MatrixRowCol&) {}    // Restore matrix row
482   virtual void NextRow(MatrixRowCol&);         // Go to next row
483   virtual void GetCol(MatrixRowCol&) = 0;      // Get matrix col
484   virtual void GetCol(MatrixColX&) = 0;        // Get matrix col
485   virtual void RestoreCol(MatrixRowCol&) {}    // Restore matrix col
486   virtual void RestoreCol(MatrixColX&) {}      // Restore matrix col
487   virtual void NextCol(MatrixRowCol&);         // Go to next col
488   virtual void NextCol(MatrixColX&);           // Go to next col
489   Real SumSquare() const;
490   Real SumAbsoluteValue() const;
491   Real Sum() const;
492   Real MaximumAbsoluteValue1(int& i) const;
493   Real MinimumAbsoluteValue1(int& i) const;
494   Real Maximum1(int& i) const;
495   Real Minimum1(int& i) const;
496   Real MaximumAbsoluteValue() const;
497   Real MaximumAbsoluteValue2(int& i, int& j) const;
498   Real MinimumAbsoluteValue() const;
499   Real MinimumAbsoluteValue2(int& i, int& j) const;
500   Real Maximum() const;
501   Real Maximum2(int& i, int& j) const;
502   Real Minimum() const;
503   Real Minimum2(int& i, int& j) const;
504   LogAndSign LogDeterminant() const;
505   virtual bool IsEqual(const GeneralMatrix&) const;
506                                                // same type, same values
507   void CheckStore() const;                     // check store is non-zero
508   virtual void SetParameters(const GeneralMatrix*) {}
509                                                // set parameters in GetMatrix
510   operator ReturnMatrix() const;               // for building a ReturnMatrix
511   ReturnMatrix ForReturn() const;
512   virtual bool SameStorageType(const GeneralMatrix& A) const;
513   virtual void ReSizeForAdd(const GeneralMatrix& A, const GeneralMatrix& B);
514   virtual void ReSizeForSP(const GeneralMatrix& A, const GeneralMatrix& B);
515   virtual void ReSize(const GeneralMatrix& A);
516   MatrixInput operator<<(Real);                // for loading a list
517   MatrixInput operator<<(int f);
518//   ReturnMatrix Reverse() const;                // reverse order of elements
519   void CleanUp();                              // to clear store
520
521   friend class Matrix;
522   friend class nricMatrix;
523   friend class SymmetricMatrix;
524   friend class UpperTriangularMatrix;
525   friend class LowerTriangularMatrix;
526   friend class DiagonalMatrix;
527   friend class CroutMatrix;
528   friend class RowVector;
529   friend class ColumnVector;
530   friend class BandMatrix;
531   friend class LowerBandMatrix;
532   friend class UpperBandMatrix;
533   friend class SymmetricBandMatrix;
534   friend class BaseMatrix;
535   friend class AddedMatrix;
536   friend class MultipliedMatrix;
537   friend class SubtractedMatrix;
538   friend class SPMatrix;
539   friend class KPMatrix;
540   friend class ConcatenatedMatrix;
541   friend class StackedMatrix;
542   friend class SolvedMatrix;
543   friend class ShiftedMatrix;
544   friend class NegShiftedMatrix;
545   friend class ScaledMatrix;
546   friend class TransposedMatrix;
547   friend class ReversedMatrix;
548   friend class NegatedMatrix;
549   friend class InvertedMatrix;
550   friend class RowedMatrix;
551   friend class ColedMatrix;
552   friend class DiagedMatrix;
553   friend class MatedMatrix;
554   friend class GetSubMatrix;
555   friend class ReturnMatrixX;
556   friend class LinearEquationSolver;
557   friend class GenericMatrix;
558   NEW_DELETE(GeneralMatrix)
559};
560
561
562
563class Matrix : public GeneralMatrix             // usual rectangular matrix
564{
565   GeneralMatrix* Image() const;                // copy of matrix
566public:
567   Matrix() {}
568   ~Matrix() {}
569   Matrix(int, int);                            // standard declaration
570   Matrix(const BaseMatrix&);                   // evaluate BaseMatrix
571   void operator=(const BaseMatrix&);
572   void operator=(Real f) { GeneralMatrix::operator=(f); }
573   void operator=(const Matrix& m) { operator=((const BaseMatrix&)m); }
574   MatrixType Type() const;
575   Real& operator()(int, int);                  // access element
576   Real& element(int, int);                     // access element
577   Real operator()(int, int) const;            // access element
578   Real element(int, int) const;               // access element
579#ifdef SETUP_C_SUBSCRIPTS
580   Real* operator[](int m) { return store+m*ncols; }
581   const Real* operator[](int m) const { return store+m*ncols; }
582#endif
583   Matrix(const Matrix& gm) { GetMatrix(&gm); }
584   GeneralMatrix* MakeSolver();
585   Real Trace() const;
586   void GetRow(MatrixRowCol&);
587   void GetCol(MatrixRowCol&);
588   void GetCol(MatrixColX&);
589   void RestoreCol(MatrixRowCol&);
590   void RestoreCol(MatrixColX&);
591   void NextRow(MatrixRowCol&);
592   void NextCol(MatrixRowCol&);
593   void NextCol(MatrixColX&);
594   virtual void ReSize(int,int);           // change dimensions
595      // virtual so we will catch it being used in a vector called as a matrix
596   void ReSize(const GeneralMatrix& A);
597   Real MaximumAbsoluteValue2(int& i, int& j) const;
598   Real MinimumAbsoluteValue2(int& i, int& j) const;
599   Real Maximum2(int& i, int& j) const;
600   Real Minimum2(int& i, int& j) const;
601   friend Real DotProduct(const Matrix& A, const Matrix& B);
602   NEW_DELETE(Matrix)
603};
604
605class nricMatrix : public Matrix                // for use with Numerical
606                                                // Recipes in C
607{
608   GeneralMatrix* Image() const;                // copy of matrix
609   Real** row_pointer;                          // points to rows
610   void MakeRowPointer();                       // build rowpointer
611   void DeleteRowPointer();
612public:
613   nricMatrix() {}
614   nricMatrix(int m, int n)                     // standard declaration
615      :  Matrix(m,n) { MakeRowPointer(); }
616   nricMatrix(const BaseMatrix& bm)             // evaluate BaseMatrix
617      :  Matrix(bm) { MakeRowPointer(); }
618   void operator=(const BaseMatrix& bm)
619      { DeleteRowPointer(); Matrix::operator=(bm); MakeRowPointer(); }
620   void operator=(Real f) { GeneralMatrix::operator=(f); }
621   void operator=(const nricMatrix& m) { operator=((const BaseMatrix&)m); }
622   void operator<<(const BaseMatrix& X)
623      { DeleteRowPointer(); Eq(X,this->Type(),true); MakeRowPointer(); }
624   nricMatrix(const nricMatrix& gm) { GetMatrix(&gm); MakeRowPointer(); }
625   void ReSize(int m, int n)               // change dimensions
626      { DeleteRowPointer(); Matrix::ReSize(m,n); MakeRowPointer(); }
627   void ReSize(const GeneralMatrix& A);
628   ~nricMatrix() { DeleteRowPointer(); }
629   Real** nric() const { CheckStore(); return row_pointer-1; }
630   void CleanUp();                                // to clear store
631   NEW_DELETE(nricMatrix)
632};
633
634class SymmetricMatrix : public GeneralMatrix
635{
636   GeneralMatrix* Image() const;                // copy of matrix
637public:
638   SymmetricMatrix() {}
639   ~SymmetricMatrix() {}
640   SymmetricMatrix(ArrayLengthSpecifier);
641   SymmetricMatrix(const BaseMatrix&);
642   void operator=(const BaseMatrix&);
643   void operator=(Real f) { GeneralMatrix::operator=(f); }
644   void operator=(const SymmetricMatrix& m) { operator=((const BaseMatrix&)m); }
645   Real& operator()(int, int);                  // access element
646   Real& element(int, int);                     // access element
647   Real operator()(int, int) const;             // access element
648   Real element(int, int) const;                // access element
649#ifdef SETUP_C_SUBSCRIPTS
650   Real* operator[](int m) { return store+(m*(m+1))/2; }
651   const Real* operator[](int m) const { return store+(m*(m+1))/2; }
652#endif
653   MatrixType Type() const;
654   SymmetricMatrix(const SymmetricMatrix& gm) { GetMatrix(&gm); }
655   Real SumSquare() const;
656   Real SumAbsoluteValue() const;
657   Real Sum() const;
658   Real Trace() const;
659   void GetRow(MatrixRowCol&);
660   void GetCol(MatrixRowCol&);
661   void GetCol(MatrixColX&);
662   void RestoreCol(MatrixRowCol&) {}
663   void RestoreCol(MatrixColX&);
664        GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
665   void ReSize(int);                       // change dimensions
666   void ReSize(const GeneralMatrix& A);
667   NEW_DELETE(SymmetricMatrix)
668};
669
670class UpperTriangularMatrix : public GeneralMatrix
671{
672   GeneralMatrix* Image() const;                // copy of matrix
673public:
674   UpperTriangularMatrix() {}
675   ~UpperTriangularMatrix() {}
676   UpperTriangularMatrix(ArrayLengthSpecifier);
677   void operator=(const BaseMatrix&);
678   void operator=(const UpperTriangularMatrix& m)
679      { operator=((const BaseMatrix&)m); }
680   UpperTriangularMatrix(const BaseMatrix&);
681   UpperTriangularMatrix(const UpperTriangularMatrix& gm) { GetMatrix(&gm); }
682   void operator=(Real f) { GeneralMatrix::operator=(f); }
683   Real& operator()(int, int);                  // access element
684   Real& element(int, int);                     // access element
685   Real operator()(int, int) const;             // access element
686   Real element(int, int) const;                // access element
687#ifdef SETUP_C_SUBSCRIPTS
688   Real* operator[](int m) { return store+m*ncols-(m*(m+1))/2; }
689   const Real* operator[](int m) const { return store+m*ncols-(m*(m+1))/2; }
690#endif
691   MatrixType Type() const;
692   GeneralMatrix* MakeSolver() { return this; } // for solving
693   void Solver(MatrixColX&, const MatrixColX&);
694   LogAndSign LogDeterminant() const;
695   Real Trace() const;
696   void GetRow(MatrixRowCol&);
697   void GetCol(MatrixRowCol&);
698   void GetCol(MatrixColX&);
699   void RestoreCol(MatrixRowCol&);
700   void RestoreCol(MatrixColX& c) { RestoreCol((MatrixRowCol&)c); }
701   void NextRow(MatrixRowCol&);
702   void ReSize(int);                       // change dimensions
703   void ReSize(const GeneralMatrix& A);
704   MatrixBandWidth BandWidth() const;
705   NEW_DELETE(UpperTriangularMatrix)
706};
707
708class LowerTriangularMatrix : public GeneralMatrix
709{
710   GeneralMatrix* Image() const;                // copy of matrix
711public:
712   LowerTriangularMatrix() {}
713   ~LowerTriangularMatrix() {}
714   LowerTriangularMatrix(ArrayLengthSpecifier);
715   LowerTriangularMatrix(const LowerTriangularMatrix& gm) { GetMatrix(&gm); }
716   LowerTriangularMatrix(const BaseMatrix& M);
717   void operator=(const BaseMatrix&);
718   void operator=(Real f) { GeneralMatrix::operator=(f); }
719   void operator=(const LowerTriangularMatrix& m)
720      { operator=((const BaseMatrix&)m); }
721   Real& operator()(int, int);                  // access element
722   Real& element(int, int);                     // access element
723   Real operator()(int, int) const;             // access element
724   Real element(int, int) const;                // access element
725#ifdef SETUP_C_SUBSCRIPTS
726   Real* operator[](int m) { return store+(m*(m+1))/2; }
727   const Real* operator[](int m) const { return store+(m*(m+1))/2; }
728#endif
729   MatrixType Type() const;
730   GeneralMatrix* MakeSolver() { return this; } // for solving
731   void Solver(MatrixColX&, const MatrixColX&);
732   LogAndSign LogDeterminant() const;
733   Real Trace() const;
734   void GetRow(MatrixRowCol&);
735   void GetCol(MatrixRowCol&);
736   void GetCol(MatrixColX&);
737   void RestoreCol(MatrixRowCol&);
738   void RestoreCol(MatrixColX& c) { RestoreCol((MatrixRowCol&)c); }
739   void NextRow(MatrixRowCol&);
740   void ReSize(int);                       // change dimensions
741   void ReSize(const GeneralMatrix& A);
742   MatrixBandWidth BandWidth() const;
743   NEW_DELETE(LowerTriangularMatrix)
744};
745
746class DiagonalMatrix : public GeneralMatrix
747{
748   GeneralMatrix* Image() const;                // copy of matrix
749public:
750   DiagonalMatrix() {}
751   ~DiagonalMatrix() {}
752   DiagonalMatrix(ArrayLengthSpecifier);
753   DiagonalMatrix(const BaseMatrix&);
754   DiagonalMatrix(const DiagonalMatrix& gm) { GetMatrix(&gm); }
755   void operator=(const BaseMatrix&);
756   void operator=(Real f) { GeneralMatrix::operator=(f); }
757   void operator=(const DiagonalMatrix& m) { operator=((const BaseMatrix&)m); }
758   Real& operator()(int, int);                  // access element
759   Real& operator()(int);                       // access element
760   Real operator()(int, int) const;             // access element
761   Real operator()(int) const;
762   Real& element(int, int);                     // access element
763   Real& element(int);                          // access element
764   Real element(int, int) const;                // access element
765   Real element(int) const;                     // access element
766#ifdef SETUP_C_SUBSCRIPTS
767   Real& operator[](int m) { return store[m]; }
768   const Real& operator[](int m) const { return store[m]; }
769#endif
770   MatrixType Type() const;
771
772   LogAndSign LogDeterminant() const;
773   Real Trace() const;
774   void GetRow(MatrixRowCol&);
775   void GetCol(MatrixRowCol&);
776   void GetCol(MatrixColX&);
777   void NextRow(MatrixRowCol&);
778   void NextCol(MatrixRowCol&);
779   void NextCol(MatrixColX&);
780   GeneralMatrix* MakeSolver() { return this; } // for solving
781   void Solver(MatrixColX&, const MatrixColX&);
782   GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
783   void ReSize(int);                       // change dimensions
784   void ReSize(const GeneralMatrix& A);
785   Real* nric() const
786      { CheckStore(); return store-1; }         // for use by NRIC
787   MatrixBandWidth BandWidth() const;
788//   ReturnMatrix Reverse() const;                // reverse order of elements
789   NEW_DELETE(DiagonalMatrix)
790};
791
792class RowVector : public Matrix
793{
794   GeneralMatrix* Image() const;                // copy of matrix
795public:
796   RowVector() { nrows = 1; }
797   ~RowVector() {}
798   RowVector(ArrayLengthSpecifier n) : Matrix(1,n.Value()) {}
799   RowVector(const BaseMatrix&);
800   RowVector(const RowVector& gm) { GetMatrix(&gm); }
801   void operator=(const BaseMatrix&);
802   void operator=(Real f) { GeneralMatrix::operator=(f); }
803   void operator=(const RowVector& m) { operator=((const BaseMatrix&)m); }
804   Real& operator()(int);                       // access element
805   Real& element(int);                          // access element
806   Real operator()(int) const;                  // access element
807   Real element(int) const;                     // access element
808#ifdef SETUP_C_SUBSCRIPTS
809   Real& operator[](int m) { return store[m]; }
810   const Real& operator[](int m) const { return store[m]; }
811#endif
812   MatrixType Type() const;
813   void GetCol(MatrixRowCol&);
814   void GetCol(MatrixColX&);
815   void NextCol(MatrixRowCol&);
816   void NextCol(MatrixColX&);
817   void RestoreCol(MatrixRowCol&) {}
818   void RestoreCol(MatrixColX& c);
819   GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
820   void ReSize(int);                       // change dimensions
821   void ReSize(int,int);                   // in case access is matrix
822   void ReSize(const GeneralMatrix& A);
823   Real* nric() const
824      { CheckStore(); return store-1; }         // for use by NRIC
825   void CleanUp();                              // to clear store
826   // friend ReturnMatrix GetMatrixRow(Matrix& A, int row);
827   NEW_DELETE(RowVector)
828};
829
830class ColumnVector : public Matrix
831{
832   GeneralMatrix* Image() const;                // copy of matrix
833public:
834   ColumnVector() { ncols = 1; }
835   ~ColumnVector() {}
836   ColumnVector(ArrayLengthSpecifier n) : Matrix(n.Value(),1) {}
837   ColumnVector(const BaseMatrix&);
838   ColumnVector(const ColumnVector& gm) { GetMatrix(&gm); }
839   void operator=(const BaseMatrix&);
840   void operator=(Real f) { GeneralMatrix::operator=(f); }
841   void operator=(const ColumnVector& m) { operator=((const BaseMatrix&)m); }
842   Real& operator()(int);                       // access element
843   Real& element(int);                          // access element
844   Real operator()(int) const;                  // access element
845   Real element(int) const;                     // access element
846#ifdef SETUP_C_SUBSCRIPTS
847   Real& operator[](int m) { return store[m]; }
848   const Real& operator[](int m) const { return store[m]; }
849#endif
850   MatrixType Type() const;
851   GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
852   void ReSize(int);                       // change dimensions
853   void ReSize(int,int);                   // in case access is matrix
854   void ReSize(const GeneralMatrix& A);
855   Real* nric() const
856      { CheckStore(); return store-1; }         // for use by NRIC
857   void CleanUp();                              // to clear store
858//   ReturnMatrix Reverse() const;                // reverse order of elements
859   NEW_DELETE(ColumnVector)
860};
861
862class CroutMatrix : public GeneralMatrix        // for LU decomposition
863{
864   int* indx;
865   bool d;
866   bool sing;
867   void ludcmp();
868public:
869   CroutMatrix(const BaseMatrix&);
870   MatrixType Type() const;
871   void lubksb(Real*, int=0);
872   ~CroutMatrix();
873   GeneralMatrix* MakeSolver() { return this; } // for solving
874   LogAndSign LogDeterminant() const;
875   void Solver(MatrixColX&, const MatrixColX&);
876   void GetRow(MatrixRowCol&);
877   void GetCol(MatrixRowCol&);
878   void GetCol(MatrixColX& c) { GetCol((MatrixRowCol&)c); }
879   void operator=(const BaseMatrix&);
880   void operator=(const CroutMatrix& m) { operator=((const BaseMatrix&)m); }
881   void CleanUp();                                // to clear store
882   bool IsEqual(const GeneralMatrix&) const;
883   bool IsSingular() const { return sing; }
884   NEW_DELETE(CroutMatrix)
885};
886
887// ***************************** band matrices ***************************/
888
889class BandMatrix : public GeneralMatrix         // band matrix
890{
891   GeneralMatrix* Image() const;                // copy of matrix
892protected:
893   void CornerClear() const;                    // set unused elements to zero
894   short SimpleAddOK(const GeneralMatrix* gm);
895public:
896   int lower, upper;                            // band widths
897   BandMatrix() { lower=0; upper=0; CornerClear(); }
898   ~BandMatrix() {}
899   BandMatrix(int n,int lb,int ub) { ReSize(n,lb,ub); CornerClear(); }
900                                                // standard declaration
901   BandMatrix(const BaseMatrix&);               // evaluate BaseMatrix
902   void operator=(const BaseMatrix&);
903   void operator=(Real f) { GeneralMatrix::operator=(f); }
904   void operator=(const BandMatrix& m) { operator=((const BaseMatrix&)m); }
905   MatrixType Type() const;
906   Real& operator()(int, int);                  // access element
907   Real& element(int, int);                     // access element
908   Real operator()(int, int) const;             // access element
909   Real element(int, int) const;                // access element
910#ifdef SETUP_C_SUBSCRIPTS
911   Real* operator[](int m) { return store+(upper+lower)*m+lower; }
912   const Real* operator[](int m) const { return store+(upper+lower)*m+lower; }
913#endif
914   BandMatrix(const BandMatrix& gm) { GetMatrix(&gm); }
915   LogAndSign LogDeterminant() const;
916   GeneralMatrix* MakeSolver();
917   Real Trace() const;
918   Real SumSquare() const { CornerClear(); return GeneralMatrix::SumSquare(); }
919   Real SumAbsoluteValue() const
920      { CornerClear(); return GeneralMatrix::SumAbsoluteValue(); }
921   Real Sum() const
922      { CornerClear(); return GeneralMatrix::Sum(); }
923   Real MaximumAbsoluteValue() const
924      { CornerClear(); return GeneralMatrix::MaximumAbsoluteValue(); }
925   Real MinimumAbsoluteValue() const
926      { int i, j; return GeneralMatrix::MinimumAbsoluteValue2(i, j); }
927   Real Maximum() const { int i, j; return GeneralMatrix::Maximum2(i, j); }
928   Real Minimum() const { int i, j; return GeneralMatrix::Minimum2(i, j); }
929   void GetRow(MatrixRowCol&);
930   void GetCol(MatrixRowCol&);
931   void GetCol(MatrixColX&);
932   void RestoreCol(MatrixRowCol&);
933   void RestoreCol(MatrixColX& c) { RestoreCol((MatrixRowCol&)c); }
934   void NextRow(MatrixRowCol&);
935   virtual void ReSize(int, int, int);             // change dimensions
936   void ReSize(const GeneralMatrix& A);
937   bool SameStorageType(const GeneralMatrix& A) const;
938   void ReSizeForAdd(const GeneralMatrix& A, const GeneralMatrix& B);
939   void ReSizeForSP(const GeneralMatrix& A, const GeneralMatrix& B);
940   MatrixBandWidth BandWidth() const;
941   void SetParameters(const GeneralMatrix*);
942   MatrixInput operator<<(Real);                // will give error
943   MatrixInput operator<<(int f);
944   void operator<<(const Real* r);              // will give error
945      // the next is included because Zortech and Borland
946      // cannot find the copy in GeneralMatrix
947   void operator<<(const BaseMatrix& X) { GeneralMatrix::operator<<(X); }
948   NEW_DELETE(BandMatrix)
949};
950
951class UpperBandMatrix : public BandMatrix       // upper band matrix
952{
953   GeneralMatrix* Image() const;                // copy of matrix
954public:
955   UpperBandMatrix() {}
956   ~UpperBandMatrix() {}
957   UpperBandMatrix(int n, int ubw)              // standard declaration
958      : BandMatrix(n, 0, ubw) {}
959   UpperBandMatrix(const BaseMatrix&);          // evaluate BaseMatrix
960   void operator=(const BaseMatrix&);
961   void operator=(Real f) { GeneralMatrix::operator=(f); }
962   void operator=(const UpperBandMatrix& m)
963      { operator=((const BaseMatrix&)m); }
964   MatrixType Type() const;
965   UpperBandMatrix(const UpperBandMatrix& gm) { GetMatrix(&gm); }
966   GeneralMatrix* MakeSolver() { return this; }
967   void Solver(MatrixColX&, const MatrixColX&);
968   LogAndSign LogDeterminant() const;
969   void ReSize(int, int, int);             // change dimensions
970   void ReSize(int n,int ubw)              // change dimensions
971      { BandMatrix::ReSize(n,0,ubw); }
972   void ReSize(const GeneralMatrix& A) { BandMatrix::ReSize(A); }
973   Real& operator()(int, int);
974   Real operator()(int, int) const;
975   Real& element(int, int);
976   Real element(int, int) const;
977#ifdef SETUP_C_SUBSCRIPTS
978   Real* operator[](int m) { return store+upper*m; }
979   const Real* operator[](int m) const { return store+upper*m; }
980#endif
981   NEW_DELETE(UpperBandMatrix)
982};
983
984class LowerBandMatrix : public BandMatrix       // upper band matrix
985{
986   GeneralMatrix* Image() const;                // copy of matrix
987public:
988   LowerBandMatrix() {}
989   ~LowerBandMatrix() {}
990   LowerBandMatrix(int n, int lbw)              // standard declaration
991      : BandMatrix(n, lbw, 0) {}
992   LowerBandMatrix(const BaseMatrix&);          // evaluate BaseMatrix
993   void operator=(const BaseMatrix&);
994   void operator=(Real f) { GeneralMatrix::operator=(f); }
995   void operator=(const LowerBandMatrix& m)
996      { operator=((const BaseMatrix&)m); }
997   MatrixType Type() const;
998   LowerBandMatrix(const LowerBandMatrix& gm) { GetMatrix(&gm); }
999   GeneralMatrix* MakeSolver() { return this; }
1000   void Solver(MatrixColX&, const MatrixColX&);
1001   LogAndSign LogDeterminant() const;
1002   void ReSize(int, int, int);             // change dimensions
1003   void ReSize(int n,int lbw)             // change dimensions
1004      { BandMatrix::ReSize(n,lbw,0); }
1005   void ReSize(const GeneralMatrix& A) { BandMatrix::ReSize(A); }
1006   Real& operator()(int, int);
1007   Real operator()(int, int) const;
1008   Real& element(int, int);
1009   Real element(int, int) const;
1010#ifdef SETUP_C_SUBSCRIPTS
1011   Real* operator[](int m) { return store+lower*(m+1); }
1012   const Real* operator[](int m) const { return store+lower*(m+1); }
1013#endif
1014   NEW_DELETE(LowerBandMatrix)
1015};
1016
1017class SymmetricBandMatrix : public GeneralMatrix
1018{
1019   GeneralMatrix* Image() const;                // copy of matrix
1020   void CornerClear() const;                    // set unused elements to zero
1021   short SimpleAddOK(const GeneralMatrix* gm);
1022public:
1023   int lower;                                   // lower band width
1024   SymmetricBandMatrix() { lower=0; CornerClear(); }
1025   ~SymmetricBandMatrix() {}
1026   SymmetricBandMatrix(int n, int lb) { ReSize(n,lb); CornerClear(); }
1027   SymmetricBandMatrix(const BaseMatrix&);
1028   void operator=(const BaseMatrix&);
1029   void operator=(Real f) { GeneralMatrix::operator=(f); }
1030   void operator=(const SymmetricBandMatrix& m)
1031      { operator=((const BaseMatrix&)m); }
1032   Real& operator()(int, int);                  // access element
1033   Real& element(int, int);                     // access element
1034   Real operator()(int, int) const;             // access element
1035   Real element(int, int) const;                // access element
1036#ifdef SETUP_C_SUBSCRIPTS
1037   Real* operator[](int m) { return store+lower*(m+1); }
1038   const Real* operator[](int m) const { return store+lower*(m+1); }
1039#endif
1040   MatrixType Type() const;
1041   SymmetricBandMatrix(const SymmetricBandMatrix& gm) { GetMatrix(&gm); }
1042   GeneralMatrix* MakeSolver();
1043   Real SumSquare() const;
1044   Real SumAbsoluteValue() const;
1045   Real Sum() const;
1046   Real MaximumAbsoluteValue() const
1047      { CornerClear(); return GeneralMatrix::MaximumAbsoluteValue(); }
1048   Real MinimumAbsoluteValue() const
1049      { int i, j; return GeneralMatrix::MinimumAbsoluteValue2(i, j); }
1050   Real Maximum() const { int i, j; return GeneralMatrix::Maximum2(i, j); }
1051   Real Minimum() const { int i, j; return GeneralMatrix::Minimum2(i, j); }
1052   Real Trace() const;
1053   LogAndSign LogDeterminant() const;
1054   void GetRow(MatrixRowCol&);
1055   void GetCol(MatrixRowCol&);
1056   void GetCol(MatrixColX&);
1057   void RestoreCol(MatrixRowCol&) {}
1058   void RestoreCol(MatrixColX&);
1059   GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
1060   void ReSize(int,int);                       // change dimensions
1061   void ReSize(const GeneralMatrix& A);
1062   bool SameStorageType(const GeneralMatrix& A) const;
1063   void ReSizeForAdd(const GeneralMatrix& A, const GeneralMatrix& B);
1064   void ReSizeForSP(const GeneralMatrix& A, const GeneralMatrix& B);
1065   MatrixBandWidth BandWidth() const;
1066   void SetParameters(const GeneralMatrix*);
1067   NEW_DELETE(SymmetricBandMatrix)
1068};
1069
1070class BandLUMatrix : public GeneralMatrix
1071// for LU decomposition of band matrix
1072{
1073   int* indx;
1074   bool d;
1075   bool sing;                                   // true if singular
1076   Real* store2;
1077   int storage2;
1078   void ludcmp();
1079   int m1,m2;                                   // lower and upper
1080public:
1081   BandLUMatrix(const BaseMatrix&);
1082   MatrixType Type() const;
1083   void lubksb(Real*, int=0);
1084   ~BandLUMatrix();
1085   GeneralMatrix* MakeSolver() { return this; } // for solving
1086   LogAndSign LogDeterminant() const;
1087   void Solver(MatrixColX&, const MatrixColX&);
1088   void GetRow(MatrixRowCol&);
1089   void GetCol(MatrixRowCol&);
1090   void GetCol(MatrixColX& c) { GetCol((MatrixRowCol&)c); }
1091   void operator=(const BaseMatrix&);
1092   void operator=(const BandLUMatrix& m) { operator=((const BaseMatrix&)m); }
1093   void CleanUp();                                // to clear store
1094   bool IsEqual(const GeneralMatrix&) const;
1095   bool IsSingular() const { return sing; }
1096   NEW_DELETE(BandLUMatrix)
1097};
1098
1099// ************************** special matrices ****************************
1100
1101class IdentityMatrix : public GeneralMatrix
1102{
1103   GeneralMatrix* Image() const;          // copy of matrix
1104public:
1105   IdentityMatrix() {}
1106   ~IdentityMatrix() {}
1107   IdentityMatrix(ArrayLengthSpecifier n) : GeneralMatrix(1)
1108      { nrows = ncols = n.Value(); *store = 1; }
1109   IdentityMatrix(const IdentityMatrix& gm) { GetMatrix(&gm); }
1110   IdentityMatrix(const BaseMatrix&);
1111   void operator=(const BaseMatrix&);
1112   void operator=(Real f) { GeneralMatrix::operator=(f); }
1113   MatrixType Type() const;
1114
1115   LogAndSign LogDeterminant() const;
1116   Real Trace() const;
1117   Real SumSquare() const;
1118   Real SumAbsoluteValue() const;
1119   Real Sum() const { return Trace(); }
1120   void GetRow(MatrixRowCol&);
1121   void GetCol(MatrixRowCol&);
1122   void GetCol(MatrixColX&);
1123   void NextRow(MatrixRowCol&);
1124   void NextCol(MatrixRowCol&);
1125   void NextCol(MatrixColX&);
1126   GeneralMatrix* MakeSolver() { return this; } // for solving
1127   void Solver(MatrixColX&, const MatrixColX&);
1128   GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
1129   void ReSize(int n);
1130   void ReSize(const GeneralMatrix& A);
1131   MatrixBandWidth BandWidth() const;
1132//   ReturnMatrix Reverse() const;                // reverse order of elements
1133   NEW_DELETE(IdentityMatrix)
1134};
1135
1136
1137
1138
1139// ************************** GenericMatrix class ************************/
1140
1141class GenericMatrix : public BaseMatrix
1142{
1143   GeneralMatrix* gm;
1144   int search(const BaseMatrix* bm) const;
1145   friend class BaseMatrix;
1146public:
1147   GenericMatrix() : gm(0) {}
1148   GenericMatrix(const BaseMatrix& bm)
1149      { gm = ((BaseMatrix&)bm).Evaluate(); gm = gm->Image(); }
1150   GenericMatrix(const GenericMatrix& bm)
1151      { gm = bm.gm->Image(); }
1152   void operator=(const GenericMatrix&);
1153   void operator=(const BaseMatrix&);
1154   void operator+=(const BaseMatrix&);
1155   void operator-=(const BaseMatrix&);
1156   void operator*=(const BaseMatrix&);
1157   void operator|=(const BaseMatrix&);
1158   void operator&=(const BaseMatrix&);
1159   void operator+=(Real);
1160   void operator-=(Real r) { operator+=(-r); }
1161   void operator*=(Real);
1162   void operator/=(Real r) { operator*=(1.0/r); }
1163   ~GenericMatrix() { delete gm; }
1164   void CleanUp() { delete gm; gm = 0; }
1165   void Release() { gm->Release(); }
1166   GeneralMatrix* Evaluate(MatrixType = MatrixTypeUnSp);
1167   MatrixBandWidth BandWidth() const;
1168   NEW_DELETE(GenericMatrix)
1169};
1170
1171// *************************** temporary classes *************************/
1172
1173class MultipliedMatrix : public BaseMatrix
1174{
1175protected:
1176   // if these union statements cause problems, simply remove them
1177   // and declare the items individually
1178   union { const BaseMatrix* bm1; GeneralMatrix* gm1; };
1179                                                  // pointers to summands
1180   union { const BaseMatrix* bm2; GeneralMatrix* gm2; };
1181   MultipliedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
1182      : bm1(bm1x),bm2(bm2x) {}
1183   int search(const BaseMatrix*) const;
1184   friend class BaseMatrix;
1185   friend class GeneralMatrix;
1186   friend class GenericMatrix;
1187public:
1188   ~MultipliedMatrix() {}
1189   GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1190   MatrixBandWidth BandWidth() const;
1191   NEW_DELETE(MultipliedMatrix)
1192};
1193
1194class AddedMatrix : public MultipliedMatrix
1195{
1196protected:
1197   AddedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
1198      : MultipliedMatrix(bm1x,bm2x) {}
1199
1200   friend class BaseMatrix;
1201   friend class GeneralMatrix;
1202   friend class GenericMatrix;
1203public:
1204   ~AddedMatrix() {}
1205   GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1206   MatrixBandWidth BandWidth() const;
1207   NEW_DELETE(AddedMatrix)
1208};
1209
1210class SPMatrix : public AddedMatrix
1211{
1212protected:
1213   SPMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
1214      : AddedMatrix(bm1x,bm2x) {}
1215
1216   friend class BaseMatrix;
1217   friend class GeneralMatrix;
1218   friend class GenericMatrix;
1219public:
1220   ~SPMatrix() {}
1221   GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1222   MatrixBandWidth BandWidth() const;
1223
1224#ifndef TEMPS_DESTROYED_QUICKLY
1225   friend SPMatrix SP(const BaseMatrix&, const BaseMatrix&);
1226#else
1227   friend SPMatrix& SP(const BaseMatrix&, const BaseMatrix&);
1228#endif
1229
1230   NEW_DELETE(SPMatrix)
1231};
1232
1233class KPMatrix : public MultipliedMatrix
1234{
1235protected:
1236   KPMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
1237      : MultipliedMatrix(bm1x,bm2x) {}
1238
1239   friend class BaseMatrix;
1240   friend class GeneralMatrix;
1241   friend class GenericMatrix;
1242public:
1243   ~KPMatrix() {}
1244   MatrixBandWidth BandWidth() const;
1245   GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1246#ifndef TEMPS_DESTROYED_QUICKLY
1247   friend KPMatrix KP(const BaseMatrix&, const BaseMatrix&);
1248#else
1249   friend KPMatrix& KP(const BaseMatrix&, const BaseMatrix&);
1250#endif
1251   NEW_DELETE(KPMatrix)
1252};
1253
1254class ConcatenatedMatrix : public MultipliedMatrix
1255{
1256protected:
1257   ConcatenatedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
1258      : MultipliedMatrix(bm1x,bm2x) {}
1259
1260   friend class BaseMatrix;
1261   friend class GeneralMatrix;
1262   friend class GenericMatrix;
1263public:
1264   ~ConcatenatedMatrix() {}
1265   MatrixBandWidth BandWidth() const;
1266   GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1267   NEW_DELETE(ConcatenatedMatrix)
1268};
1269
1270class StackedMatrix : public ConcatenatedMatrix
1271{
1272protected:
1273   StackedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
1274      : ConcatenatedMatrix(bm1x,bm2x) {}
1275
1276   friend class BaseMatrix;
1277   friend class GeneralMatrix;
1278   friend class GenericMatrix;
1279public:
1280   ~StackedMatrix() {}
1281   GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1282   NEW_DELETE(StackedMatrix)
1283};
1284
1285class SolvedMatrix : public MultipliedMatrix
1286{
1287   SolvedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
1288      : MultipliedMatrix(bm1x,bm2x) {}
1289   friend class BaseMatrix;
1290   friend class InvertedMatrix;                        // for operator*
1291public:
1292   ~SolvedMatrix() {}
1293   GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1294   MatrixBandWidth BandWidth() const;
1295   NEW_DELETE(SolvedMatrix)
1296};
1297
1298class SubtractedMatrix : public AddedMatrix
1299{
1300   SubtractedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
1301      : AddedMatrix(bm1x,bm2x) {}
1302   friend class BaseMatrix;
1303   friend class GeneralMatrix;
1304   friend class GenericMatrix;
1305public:
1306   ~SubtractedMatrix() {}
1307   GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1308   NEW_DELETE(SubtractedMatrix)
1309};
1310
1311class ShiftedMatrix : public BaseMatrix
1312{
1313protected:
1314   union { const BaseMatrix* bm; GeneralMatrix* gm; };
1315   Real f;
1316   ShiftedMatrix(const BaseMatrix* bmx, Real fx) : bm(bmx),f(fx) {}
1317   int search(const BaseMatrix*) const;
1318   friend class BaseMatrix;
1319   friend class GeneralMatrix;
1320   friend class GenericMatrix;
1321public:
1322   ~ShiftedMatrix() {}
1323   GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1324#ifndef TEMPS_DESTROYED_QUICKLY
1325   friend ShiftedMatrix operator+(Real f, const BaseMatrix& BM);
1326//      { return ShiftedMatrix(&BM, f); }
1327#endif
1328   NEW_DELETE(ShiftedMatrix)
1329};
1330
1331class NegShiftedMatrix : public ShiftedMatrix
1332{
1333protected:
1334   NegShiftedMatrix(Real fx, const BaseMatrix* bmx) : ShiftedMatrix(bmx,fx) {}
1335   friend class BaseMatrix;
1336   friend class GeneralMatrix;
1337   friend class GenericMatrix;
1338public:
1339   ~NegShiftedMatrix() {}
1340   GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1341#ifndef TEMPS_DESTROYED_QUICKLY
1342   friend NegShiftedMatrix operator-(Real, const BaseMatrix&);
1343#else
1344   friend NegShiftedMatrix& operator-(Real, const BaseMatrix&);
1345#endif
1346   NEW_DELETE(NegShiftedMatrix)
1347};
1348
1349class ScaledMatrix : public ShiftedMatrix
1350{
1351   ScaledMatrix(const BaseMatrix* bmx, Real fx) : ShiftedMatrix(bmx,fx) {}
1352   friend class BaseMatrix;
1353   friend class GeneralMatrix;
1354   friend class GenericMatrix;
1355public:
1356   ~ScaledMatrix() {}
1357   GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1358   MatrixBandWidth BandWidth() const;
1359#ifndef TEMPS_DESTROYED_QUICKLY
1360   friend ScaledMatrix operator*(Real f, const BaseMatrix& BM);
1361      //{ return ScaledMatrix(&BM, f); }
1362#endif
1363   NEW_DELETE(ScaledMatrix)
1364};
1365
1366class NegatedMatrix : public BaseMatrix
1367{
1368protected:
1369   union { const BaseMatrix* bm; GeneralMatrix* gm; };
1370   NegatedMatrix(const BaseMatrix* bmx) : bm(bmx) {}
1371   int search(const BaseMatrix*) const;
1372private:
1373   friend class BaseMatrix;
1374public:
1375   ~NegatedMatrix() {}
1376   GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1377   MatrixBandWidth BandWidth() const;
1378   NEW_DELETE(NegatedMatrix)
1379};
1380
1381class TransposedMatrix : public NegatedMatrix
1382{
1383   TransposedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
1384   friend class BaseMatrix;
1385public:
1386   ~TransposedMatrix() {}
1387   GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1388   MatrixBandWidth BandWidth() const;
1389   NEW_DELETE(TransposedMatrix)
1390};
1391
1392class ReversedMatrix : public NegatedMatrix
1393{
1394   ReversedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
1395   friend class BaseMatrix;
1396public:
1397   ~ReversedMatrix() {}
1398   GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1399   NEW_DELETE(ReversedMatrix)
1400};
1401
1402class InvertedMatrix : public NegatedMatrix
1403{
1404   InvertedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
1405public:
1406   ~InvertedMatrix() {}
1407#ifndef TEMPS_DESTROYED_QUICKLY
1408   SolvedMatrix operator*(const BaseMatrix&) const;       // inverse(A) * B
1409   ScaledMatrix operator*(Real t) const { return BaseMatrix::operator*(t); }
1410#else
1411   SolvedMatrix& operator*(const BaseMatrix&);            // inverse(A) * B
1412   ScaledMatrix& operator*(Real t) const { return BaseMatrix::operator*(t); }
1413#endif
1414   friend class BaseMatrix;
1415   GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1416   MatrixBandWidth BandWidth() const;
1417   NEW_DELETE(InvertedMatrix)
1418};
1419
1420class RowedMatrix : public NegatedMatrix
1421{
1422   RowedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
1423   friend class BaseMatrix;
1424public:
1425   ~RowedMatrix() {}
1426   GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1427   MatrixBandWidth BandWidth() const;
1428   NEW_DELETE(RowedMatrix)
1429};
1430
1431class ColedMatrix : public NegatedMatrix
1432{
1433   ColedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
1434   friend class BaseMatrix;
1435public:
1436   ~ColedMatrix() {}
1437   GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1438   MatrixBandWidth BandWidth() const;
1439   NEW_DELETE(ColedMatrix)
1440};
1441
1442class DiagedMatrix : public NegatedMatrix
1443{
1444   DiagedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
1445   friend class BaseMatrix;
1446public:
1447   ~DiagedMatrix() {}
1448   GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1449   MatrixBandWidth BandWidth() const;
1450   NEW_DELETE(DiagedMatrix)
1451};
1452
1453class MatedMatrix : public NegatedMatrix
1454{
1455   int nr, nc;
1456   MatedMatrix(const BaseMatrix* bmx, int nrx, int ncx)
1457      : NegatedMatrix(bmx), nr(nrx), nc(ncx) {}
1458   friend class BaseMatrix;
1459public:
1460   ~MatedMatrix() {}
1461   GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1462   MatrixBandWidth BandWidth() const;
1463   NEW_DELETE(MatedMatrix)
1464};
1465
1466class ReturnMatrixX : public BaseMatrix    // for matrix return
1467{
1468   GeneralMatrix* gm;
1469   int search(const BaseMatrix*) const;
1470public:
1471   ~ReturnMatrixX() {}
1472   GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1473   friend class BaseMatrix;
1474#ifdef TEMPS_DESTROYED_QUICKLY_R
1475   ReturnMatrixX(const ReturnMatrixX& tm);
1476#else
1477   ReturnMatrixX(const ReturnMatrixX& tm) : gm(tm.gm) {}
1478#endif
1479   ReturnMatrixX(const GeneralMatrix* gmx) : gm((GeneralMatrix*&)gmx) {}
1480//   ReturnMatrixX(GeneralMatrix&);
1481   MatrixBandWidth BandWidth() const;
1482   NEW_DELETE(ReturnMatrixX)
1483};
1484
1485
1486// ************************** submatrices ******************************/
1487
1488class GetSubMatrix : public NegatedMatrix
1489{
1490   int row_skip;
1491   int row_number;
1492   int col_skip;
1493   int col_number;
1494   bool IsSym;
1495
1496   GetSubMatrix
1497      (const BaseMatrix* bmx, int rs, int rn, int cs, int cn, bool is)
1498      : NegatedMatrix(bmx),
1499      row_skip(rs), row_number(rn), col_skip(cs), col_number(cn), IsSym(is) {}
1500   void SetUpLHS();
1501   friend class BaseMatrix;
1502public:
1503   GetSubMatrix(const GetSubMatrix& g)
1504      : NegatedMatrix(g.bm), row_skip(g.row_skip), row_number(g.row_number),
1505      col_skip(g.col_skip), col_number(g.col_number), IsSym(g.IsSym) {}
1506   ~GetSubMatrix() {}
1507   GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
1508   void operator=(const BaseMatrix&);
1509   void operator+=(const BaseMatrix&);
1510   void operator-=(const BaseMatrix&);
1511   void operator=(const GetSubMatrix& m) { operator=((const BaseMatrix&)m); }
1512   void operator<<(const BaseMatrix&);
1513   void operator<<(const Real*);                // copy from array
1514   MatrixInput operator<<(Real);                // for loading a list
1515   MatrixInput operator<<(int f);
1516   void operator=(Real);                        // copy from constant
1517   void operator+=(Real);                       // add constant
1518   void operator-=(Real r) { operator+=(-r); }  // subtract constant
1519   void operator*=(Real);                       // multiply by constant
1520   void operator/=(Real r) { operator*=(1.0/r); } // divide by constant
1521   void Inject(const GeneralMatrix&);           // copy stored els only
1522   MatrixBandWidth BandWidth() const;
1523   NEW_DELETE(GetSubMatrix)
1524};
1525
1526// ******************** linear equation solving ****************************/
1527
1528class LinearEquationSolver : public BaseMatrix
1529{
1530   GeneralMatrix* gm;
1531   int search(const BaseMatrix*) const { return 0; }
1532   friend class BaseMatrix;
1533public:
1534   LinearEquationSolver(const BaseMatrix& bm);
1535   ~LinearEquationSolver() { delete gm; }
1536   void CleanUp() { delete gm; } 
1537   GeneralMatrix* Evaluate(MatrixType) { return gm; }
1538   // probably should have an error message if MatrixType != UnSp
1539   NEW_DELETE(LinearEquationSolver)
1540};
1541
1542// ************************** matrix input *******************************/
1543
1544class MatrixInput          // for reading a list of values into a matrix
1545                           // the difficult part is detecting a mismatch
1546                           // in the number of elements
1547{
1548   int n;                  // number values still to be read
1549   Real* r;                // pointer to next location to be read to
1550public:
1551   MatrixInput(const MatrixInput& mi) : n(mi.n), r(mi.r) {}
1552   MatrixInput(int nx, Real* rx) : n(nx), r(rx) {}
1553   ~MatrixInput();
1554   MatrixInput operator<<(Real);
1555   MatrixInput operator<<(int f);
1556   friend class GeneralMatrix;
1557};
1558
1559
1560
1561// **************** a very simple integer array class ********************/
1562
1563// A minimal array class to imitate a C style array but giving dynamic storage
1564// mostly intended for internal use by newmat
1565
1566class SimpleIntArray : public Janitor
1567{
1568protected:
1569   int* a;                    // pointer to the array
1570   int n;                     // length of the array
1571public:
1572   SimpleIntArray(int xn);    // build an array length xn
1573   ~SimpleIntArray();         // return the space to memory
1574   int& operator[](int i);    // access element of the array - start at 0
1575   int operator[](int i) const;
1576                              // access element of constant array
1577   void operator=(int ai);    // set the array equal to a constant
1578   void operator=(const SimpleIntArray& b);
1579                              // copy the elements of an array
1580   SimpleIntArray(const SimpleIntArray& b);
1581                              // make a new array equal to an existing one
1582   int Size() const { return n; }
1583                              // return the size of the array
1584   int* Data() { return a; }  // pointer to the data
1585   const int* Data() const { return a; }
1586                              // pointer to the data
1587   void ReSize(int i, bool keep = false);
1588                              // change length, keep data if keep = true
1589   void CleanUp() { ReSize(0); }
1590   NEW_DELETE(SimpleIntArray)
1591};
1592
1593// *************************** exceptions ********************************/
1594
1595class NPDException : public Runtime_error     // Not positive definite
1596{
1597public:
1598   static unsigned long Select;          // for identifying exception
1599   NPDException(const GeneralMatrix&);
1600};
1601
1602class ConvergenceException : public Runtime_error
1603{
1604public:
1605   static unsigned long Select;          // for identifying exception
1606   ConvergenceException(const GeneralMatrix& A);
1607   ConvergenceException(const char* c);
1608};
1609
1610class SingularException : public Runtime_error
1611{
1612public:
1613   static unsigned long Select;          // for identifying exception
1614   SingularException(const GeneralMatrix& A);
1615};
1616
1617class OverflowException : public Runtime_error
1618{
1619public:
1620   static unsigned long Select;          // for identifying exception
1621   OverflowException(const char* c);
1622};
1623
1624class ProgramException : public Logic_error
1625{
1626protected:
1627   ProgramException();
1628public:
1629   static unsigned long Select;          // for identifying exception
1630   ProgramException(const char* c);
1631   ProgramException(const char* c, const GeneralMatrix&);
1632   ProgramException(const char* c, const GeneralMatrix&, const GeneralMatrix&);
1633   ProgramException(const char* c, MatrixType, MatrixType);
1634};
1635
1636class IndexException : public Logic_error
1637{
1638public:
1639   static unsigned long Select;          // for identifying exception
1640   IndexException(int i, const GeneralMatrix& A);
1641   IndexException(int i, int j, const GeneralMatrix& A);
1642   // next two are for access via element function
1643   IndexException(int i, const GeneralMatrix& A, bool);
1644   IndexException(int i, int j, const GeneralMatrix& A, bool);
1645};
1646
1647class VectorException : public Logic_error    // cannot convert to vector
1648{
1649public:
1650   static unsigned long Select;          // for identifying exception
1651   VectorException();
1652   VectorException(const GeneralMatrix& A);
1653};
1654
1655class NotSquareException : public Logic_error
1656{
1657public:
1658   static unsigned long Select;          // for identifying exception
1659   NotSquareException(const GeneralMatrix& A);
1660};
1661
1662class SubMatrixDimensionException : public Logic_error
1663{
1664public:
1665   static unsigned long Select;          // for identifying exception
1666   SubMatrixDimensionException();
1667};
1668
1669class IncompatibleDimensionsException : public Logic_error
1670{
1671public:
1672   static unsigned long Select;          // for identifying exception
1673   IncompatibleDimensionsException();
1674   IncompatibleDimensionsException(const GeneralMatrix&, const GeneralMatrix&);
1675};
1676
1677class NotDefinedException : public Logic_error
1678{
1679public:
1680   static unsigned long Select;          // for identifying exception
1681   NotDefinedException(const char* op, const char* matrix);
1682};
1683
1684class CannotBuildException : public Logic_error
1685{
1686public:
1687   static unsigned long Select;          // for identifying exception
1688   CannotBuildException(const char* matrix);
1689};
1690
1691
1692class InternalException : public Logic_error
1693{
1694public:
1695   static unsigned long Select;          // for identifying exception
1696   InternalException(const char* c);
1697};
1698
1699// ************************ functions ************************************ //
1700
1701bool operator==(const GeneralMatrix& A, const GeneralMatrix& B);
1702bool operator==(const BaseMatrix& A, const BaseMatrix& B);
1703inline bool operator!=(const GeneralMatrix& A, const GeneralMatrix& B)
1704   { return ! (A==B); }
1705inline bool operator!=(const BaseMatrix& A, const BaseMatrix& B)
1706   { return ! (A==B); }
1707
1708   // inequality operators are dummies included for compatibility
1709   // with STL. They throw an exception if actually called.
1710inline bool operator<=(const BaseMatrix& A, const BaseMatrix&)
1711   { A.IEQND(); return true; }
1712inline bool operator>=(const BaseMatrix& A, const BaseMatrix&)
1713   { A.IEQND(); return true; }
1714inline bool operator<(const BaseMatrix& A, const BaseMatrix&)
1715   { A.IEQND(); return true; }
1716inline bool operator>(const BaseMatrix& A, const BaseMatrix&)
1717   { A.IEQND(); return true; }
1718
1719bool IsZero(const BaseMatrix& A);
1720
1721
1722// ********************* inline functions ******************************** //
1723
1724
1725inline LogAndSign LogDeterminant(const BaseMatrix& B)
1726   { return B.LogDeterminant(); }
1727inline Real Determinant(const BaseMatrix& B)
1728   { return B.Determinant(); }
1729inline Real SumSquare(const BaseMatrix& B) { return B.SumSquare(); }
1730inline Real NormFrobenius(const BaseMatrix& B) { return B.NormFrobenius(); }
1731inline Real Trace(const BaseMatrix& B) { return B.Trace(); }
1732inline Real SumAbsoluteValue(const BaseMatrix& B)
1733   { return B.SumAbsoluteValue(); }
1734inline Real Sum(const BaseMatrix& B)
1735   { return B.Sum(); }
1736inline Real MaximumAbsoluteValue(const BaseMatrix& B)
1737   { return B.MaximumAbsoluteValue(); }
1738inline Real MinimumAbsoluteValue(const BaseMatrix& B)
1739   { return B.MinimumAbsoluteValue(); }
1740inline Real Maximum(const BaseMatrix& B) { return B.Maximum(); }
1741inline Real Minimum(const BaseMatrix& B) { return B.Minimum(); }
1742inline Real Norm1(const BaseMatrix& B) { return B.Norm1(); }
1743inline Real Norm1(RowVector& RV) { return RV.MaximumAbsoluteValue(); }
1744inline Real NormInfinity(const BaseMatrix& B) { return B.NormInfinity(); }
1745inline Real NormInfinity(ColumnVector& CV)
1746   { return CV.MaximumAbsoluteValue(); }
1747inline bool IsZero(const GeneralMatrix& A) { return A.IsZero(); }
1748
1749#ifdef TEMPS_DESTROYED_QUICKLY
1750inline ShiftedMatrix& operator+(Real f, const BaseMatrix& BM)
1751   { return BM + f; }
1752inline ScaledMatrix& operator*(Real f, const BaseMatrix& BM)
1753   { return BM * f; }
1754#endif
1755
1756// these are moved out of the class definitions because of a problem
1757// with the Intel 8.1 compiler
1758#ifndef TEMPS_DESTROYED_QUICKLY
1759   inline ShiftedMatrix operator+(Real f, const BaseMatrix& BM)
1760      { return ShiftedMatrix(&BM, f); }
1761   inline ScaledMatrix operator*(Real f, const BaseMatrix& BM)
1762      { return ScaledMatrix(&BM, f); }
1763#endif
1764
1765
1766inline MatrixInput MatrixInput::operator<<(int f) { return *this << (Real)f; }
1767inline MatrixInput GeneralMatrix::operator<<(int f) { return *this << (Real)f; }
1768inline MatrixInput BandMatrix::operator<<(int f) { return *this << (Real)f; }
1769inline MatrixInput GetSubMatrix::operator<<(int f) { return *this << (Real)f; }
1770
1771
1772
1773#ifdef use_namespace
1774}
1775#endif
1776
1777
1778#endif
1779
1780// body file: newmat1.cpp
1781// body file: newmat2.cpp
1782// body file: newmat3.cpp
1783// body file: newmat4.cpp
1784// body file: newmat5.cpp
1785// body file: newmat6.cpp
1786// body file: newmat7.cpp
1787// body file: newmat8.cpp
1788// body file: newmatex.cpp
1789// body file: bandmat.cpp
1790// body file: submat.cpp
1791
1792
1793
1794
1795
1796
1797
Note: See TracBrowser for help on using the repository browser.