/// \defgroup newmat Newmat matrix manipulation library
///@{
/// \file newmat.h
/// Definition file for matrix library.
// Copyright (C) 2004: R B Davies
#ifndef NEWMAT_LIB
#define NEWMAT_LIB 0
#include "include.h"
#include "myexcept.h"
#ifdef use_namespace
namespace NEWMAT { using namespace RBD_COMMON; }
namespace RBD_LIBRARIES { using namespace NEWMAT; }
namespace NEWMAT {
#endif
//#define DO_REPORT // to activate REPORT
#ifdef NO_LONG_NAMES
#define UpperTriangularMatrix UTMatrix
#define LowerTriangularMatrix LTMatrix
#define SymmetricMatrix SMatrix
#define DiagonalMatrix DMatrix
#define BandMatrix BMatrix
#define UpperBandMatrix UBMatrix
#define LowerBandMatrix LBMatrix
#define SymmetricBandMatrix SBMatrix
#define BandLUMatrix BLUMatrix
#endif
// ************************** general utilities ****************************/
class GeneralMatrix; // defined later
class BaseMatrix; // defined later
class MatrixInput; // defined later
void MatrixErrorNoSpace(const void*); ///< test for allocation fails
/// Return from LogDeterminant function.
/// Members are the log of the absolute value and the sign (+1, -1 or 0)
class LogAndSign
{
Real log_val;
int sign_val;
public:
LogAndSign() { log_val=0.0; sign_val=1; }
LogAndSign(Real);
void operator*=(Real); ///< multiply by a real
void pow_eq(int k); ///< raise to power of k
void PowEq(int k) { pow_eq(k); }
void ChangeSign() { sign_val = -sign_val; }
void change_sign() { sign_val = -sign_val; } ///< change sign
Real LogValue() const { return log_val; }
Real log_value() const { return log_val; } ///< log of the absolute value
int Sign() const { return sign_val; }
int sign() const { return sign_val; } ///< sign of the value
Real value() const; ///< the value
Real Value() const { return value(); }
FREE_CHECK(LogAndSign)
};
// the following class is for counting the number of times a piece of code
// is executed. It is used for locating any code not executed by test
// routines. Use turbo GREP locate all places this code is called and
// check which ones are not accessed.
// Somewhat implementation dependent as it relies on "cout" still being
// present when ExeCounter objects are destructed.
#ifdef DO_REPORT
class ExeCounter
{
int line; // code line number
int fileid; // file identifier
long nexe; // number of executions
static int nreports; // number of reports
public:
ExeCounter(int,int);
void operator++() { nexe++; }
~ExeCounter(); // prints out reports
};
#endif
// ************************** class MatrixType *****************************
/// Find the type of a matrix resulting from matrix operations.
/// Also identify what conversions are permissible.
/// This class must be updated when new matrix types are added.
class MatrixType
{
public:
enum Attribute { Valid = 1,
Diagonal = 2, // order of these is important
Symmetric = 4,
Band = 8,
Lower = 16,
Upper = 32,
Square = 64,
Skew = 128,
LUDeco = 256,
Ones = 512 };
enum { US = 0,
UT = Valid + Upper + Square,
LT = Valid + Lower + Square,
Rt = Valid,
Sq = Valid + Square,
Sm = Valid + Symmetric + Square,
Sk = Valid + Skew + Square,
Dg = Valid + Diagonal + Band + Lower + Upper + Symmetric
+ Square,
Id = Valid + Diagonal + Band + Lower + Upper + Symmetric
+ Square + Ones,
RV = Valid, // do not separate out
CV = Valid, // vectors
BM = Valid + Band + Square,
UB = Valid + Band + Upper + Square,
LB = Valid + Band + Lower + Square,
SB = Valid + Band + Symmetric + Square,
KB = Valid + Band + Skew + Square,
Ct = Valid + LUDeco + Square,
BC = Valid + Band + LUDeco + Square,
Mask = ~Square
};
static int nTypes() { return 13; } // number of different types
// exclude Ct, US, BC
public:
int attribute;
bool DataLossOK; // true if data loss is OK when
// this represents a destination
public:
MatrixType () : DataLossOK(false) {}
MatrixType (int i) : attribute(i), DataLossOK(false) {}
MatrixType (int i, bool dlok) : attribute(i), DataLossOK(dlok) {}
MatrixType (const MatrixType& mt)
: attribute(mt.attribute), DataLossOK(mt.DataLossOK) {}
void operator=(const MatrixType& mt)
{ attribute = mt.attribute; DataLossOK = mt.DataLossOK; }
void SetDataLossOK() { DataLossOK = true; }
int operator+() const { return attribute; }
MatrixType operator+(MatrixType mt) const
{ return MatrixType(attribute & mt.attribute); }
MatrixType operator*(const MatrixType&) const;
MatrixType SP(const MatrixType&) const;
MatrixType KP(const MatrixType&) const;
MatrixType operator|(const MatrixType& mt) const
{ return MatrixType(attribute & mt.attribute & Valid); }
MatrixType operator&(const MatrixType& mt) const
{ return MatrixType(attribute & mt.attribute & Valid); }
bool operator>=(MatrixType mt) const
{ return ( attribute & ~mt.attribute & Mask ) == 0; }
bool operator<(MatrixType mt) const // for MS Visual C++ 4
{ return ( attribute & ~mt.attribute & Mask ) != 0; }
bool operator==(MatrixType t) const
{ return (attribute == t.attribute); }
bool operator!=(MatrixType t) const
{ return (attribute != t.attribute); }
bool operator!() const { return (attribute & Valid) == 0; }
MatrixType i() const; ///< type of inverse
MatrixType t() const; ///< type of transpose
MatrixType AddEqualEl() const ///< add constant to matrix
{ return MatrixType(attribute & (Valid + Symmetric + Square)); }
MatrixType MultRHS() const; ///< type for rhs of multiply
MatrixType sub() const ///< type of submatrix
{ return MatrixType(attribute & Valid); }
MatrixType ssub() const ///< type of sym submatrix
{ return MatrixType(attribute); } // not for selection matrix
GeneralMatrix* New() const; ///< new matrix of given type
GeneralMatrix* New(int,int,BaseMatrix*) const;
///< new matrix of given type
const char* value() const; ///< type as char string
const char* Value() const { return value(); }
friend bool Rectangular(MatrixType a, MatrixType b, MatrixType c);
friend bool Compare(const MatrixType&, MatrixType&);
///< compare and check conversion
bool is_band() const { return (attribute & Band) != 0; }
bool is_diagonal() const { return (attribute & Diagonal) != 0; }
bool is_symmetric() const { return (attribute & Symmetric) != 0; }
bool CannotConvert() const { return (attribute & LUDeco) != 0; }
// used by operator==
FREE_CHECK(MatrixType)
};
// *********************** class MatrixBandWidth ***********************/
///Upper and lower bandwidths of a matrix.
///That is number of diagonals strictly above or below main diagonal,
///e.g. diagonal matrix has 0 upper and lower bandwiths.
///-1 means the matrix may have the maximum bandwidth.
class MatrixBandWidth
{
public:
int lower_val;
int upper_val;
MatrixBandWidth(const int l, const int u) : lower_val(l), upper_val(u) {}
MatrixBandWidth(const int i) : lower_val(i), upper_val(i) {}
MatrixBandWidth operator+(const MatrixBandWidth&) const;
MatrixBandWidth operator*(const MatrixBandWidth&) const;
MatrixBandWidth minimum(const MatrixBandWidth&) const;
MatrixBandWidth t() const { return MatrixBandWidth(upper_val,lower_val); }
bool operator==(const MatrixBandWidth& bw) const
{ return (lower_val == bw.lower_val) && (upper_val == bw.upper_val); }
bool operator!=(const MatrixBandWidth& bw) const { return !operator==(bw); }
int Upper() const { return upper_val; }
int upper() const { return upper_val; }
int Lower() const { return lower_val; }
int lower() const { return lower_val; }
FREE_CHECK(MatrixBandWidth)
};
// ********************* Array length specifier ************************/
/// This class is used to avoid constructors such as
/// ColumnVector(int) being used for conversions.
/// Eventually this should be replaced by the use of the keyword "explicit".
class ArrayLengthSpecifier
{
int v;
public:
int Value() const { return v; }
int value() const { return v; }
ArrayLengthSpecifier(int l) : v(l) {}
};
// ************************* Matrix routines ***************************/
class MatrixRowCol; // defined later
class MatrixRow;
class MatrixCol;
class MatrixColX;
class GeneralMatrix; // defined later
class AddedMatrix;
class MultipliedMatrix;
class SubtractedMatrix;
class SPMatrix;
class KPMatrix;
class ConcatenatedMatrix;
class StackedMatrix;
class SolvedMatrix;
class ShiftedMatrix;
class NegShiftedMatrix;
class ScaledMatrix;
class TransposedMatrix;
class ReversedMatrix;
class NegatedMatrix;
class InvertedMatrix;
class RowedMatrix;
class ColedMatrix;
class DiagedMatrix;
class MatedMatrix;
class GetSubMatrix;
class ReturnMatrix;
class Matrix;
class SquareMatrix;
class nricMatrix;
class RowVector;
class ColumnVector;
class SymmetricMatrix;
class UpperTriangularMatrix;
class LowerTriangularMatrix;
class DiagonalMatrix;
class CroutMatrix;
class BandMatrix;
class LowerBandMatrix;
class UpperBandMatrix;
class SymmetricBandMatrix;
class LinearEquationSolver;
class GenericMatrix;
#define MatrixTypeUnSp 0
//static MatrixType MatrixTypeUnSp(MatrixType::US);
// // AT&T needs this
/// Base of the matrix classes.
class BaseMatrix : public Janitor
{
protected:
virtual int search(const BaseMatrix*) const = 0;
// count number of times matrix is referred to
public:
virtual GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp) = 0;
// evaluate temporary
// for old version of G++
// virtual GeneralMatrix* Evaluate(MatrixType mt) = 0;
// GeneralMatrix* Evaluate() { return Evaluate(MatrixTypeUnSp); }
AddedMatrix operator+(const BaseMatrix&) const; // results of operations
MultipliedMatrix operator*(const BaseMatrix&) const;
SubtractedMatrix operator-(const BaseMatrix&) const;
ConcatenatedMatrix operator|(const BaseMatrix&) const;
StackedMatrix operator&(const BaseMatrix&) const;
ShiftedMatrix operator+(Real) const;
ScaledMatrix operator*(Real) const;
ScaledMatrix operator/(Real) const;
ShiftedMatrix operator-(Real) const;
TransposedMatrix t() const;
// TransposedMatrix t;
NegatedMatrix operator-() const; // change sign of elements
ReversedMatrix reverse() const;
ReversedMatrix Reverse() const;
InvertedMatrix i() const;
// InvertedMatrix i;
RowedMatrix as_row() const;
RowedMatrix AsRow() const;
ColedMatrix as_column() const;
ColedMatrix AsColumn() const;
DiagedMatrix as_diagonal() const;
DiagedMatrix AsDiagonal() const;
MatedMatrix as_matrix(int,int) const;
MatedMatrix AsMatrix(int m, int n) const;
GetSubMatrix submatrix(int,int,int,int) const;
GetSubMatrix SubMatrix(int fr, int lr, int fc, int lc) const;
GetSubMatrix sym_submatrix(int,int) const;
GetSubMatrix SymSubMatrix(int f, int l) const;
GetSubMatrix row(int) const;
GetSubMatrix rows(int,int) const;
GetSubMatrix column(int) const;
GetSubMatrix columns(int,int) const;
GetSubMatrix Row(int f) const;
GetSubMatrix Rows(int f, int l) const;
GetSubMatrix Column(int f) const;
GetSubMatrix Columns(int f, int l) const;
Real as_scalar() const; // conversion of 1 x 1 matrix
Real AsScalar() const;
virtual LogAndSign log_determinant() const;
LogAndSign LogDeterminant() const { return log_determinant(); }
Real determinant() const;
Real Determinant() const { return determinant(); }
virtual Real sum_square() const;
Real SumSquare() const { return sum_square(); }
Real norm_Frobenius() const;
Real norm_frobenius() const { return norm_Frobenius(); }
Real NormFrobenius() const { return norm_Frobenius(); }
virtual Real sum_absolute_value() const;
Real SumAbsoluteValue() const { return sum_absolute_value(); }
virtual Real sum() const;
virtual Real Sum() const { return sum(); }
virtual Real maximum_absolute_value() const;
Real MaximumAbsoluteValue() const { return maximum_absolute_value(); }
virtual Real maximum_absolute_value1(int& i) const;
Real MaximumAbsoluteValue1(int& i) const
{ return maximum_absolute_value1(i); }
virtual Real maximum_absolute_value2(int& i, int& j) const;
Real MaximumAbsoluteValue2(int& i, int& j) const
{ return maximum_absolute_value2(i,j); }
virtual Real minimum_absolute_value() const;
Real MinimumAbsoluteValue() const { return minimum_absolute_value(); }
virtual Real minimum_absolute_value1(int& i) const;
Real MinimumAbsoluteValue1(int& i) const
{ return minimum_absolute_value1(i); }
virtual Real minimum_absolute_value2(int& i, int& j) const;
Real MinimumAbsoluteValue2(int& i, int& j) const
{ return minimum_absolute_value2(i,j); }
virtual Real maximum() const;
Real Maximum() const { return maximum(); }
virtual Real maximum1(int& i) const;
Real Maximum1(int& i) const { return maximum1(i); }
virtual Real maximum2(int& i, int& j) const;
Real Maximum2(int& i, int& j) const { return maximum2(i,j); }
virtual Real minimum() const;
Real Minimum() const { return minimum(); }
virtual Real minimum1(int& i) const;
Real Minimum1(int& i) const { return minimum1(i); }
virtual Real minimum2(int& i, int& j) const;
Real Minimum2(int& i, int& j) const { return minimum2(i,j); }
virtual Real trace() const;
Real Trace() const { return trace(); }
Real norm1() const;
Real Norm1() const { return norm1(); }
Real norm_infinity() const;
Real NormInfinity() const { return norm_infinity(); }
virtual MatrixBandWidth bandwidth() const; // bandwidths of band matrix
virtual MatrixBandWidth BandWidth() const { return bandwidth(); }
void IEQND() const; // called by ineq. ops
ReturnMatrix sum_square_columns() const;
ReturnMatrix sum_square_rows() const;
ReturnMatrix sum_columns() const;
ReturnMatrix sum_rows() const;
virtual void cleanup() {}
void CleanUp() { cleanup(); }
// virtual ReturnMatrix Reverse() const; // reverse order of elements
//protected:
// BaseMatrix() : t(this), i(this) {}
friend class GeneralMatrix;
friend class Matrix;
friend class SquareMatrix;
friend class nricMatrix;
friend class RowVector;
friend class ColumnVector;
friend class SymmetricMatrix;
friend class UpperTriangularMatrix;
friend class LowerTriangularMatrix;
friend class DiagonalMatrix;
friend class CroutMatrix;
friend class BandMatrix;
friend class LowerBandMatrix;
friend class UpperBandMatrix;
friend class SymmetricBandMatrix;
friend class AddedMatrix;
friend class MultipliedMatrix;
friend class SubtractedMatrix;
friend class SPMatrix;
friend class KPMatrix;
friend class ConcatenatedMatrix;
friend class StackedMatrix;
friend class SolvedMatrix;
friend class ShiftedMatrix;
friend class NegShiftedMatrix;
friend class ScaledMatrix;
friend class TransposedMatrix;
friend class ReversedMatrix;
friend class NegatedMatrix;
friend class InvertedMatrix;
friend class RowedMatrix;
friend class ColedMatrix;
friend class DiagedMatrix;
friend class MatedMatrix;
friend class GetSubMatrix;
friend class ReturnMatrix;
friend class LinearEquationSolver;
friend class GenericMatrix;
NEW_DELETE(BaseMatrix)
};
// ***************************** working classes **************************/
/// The classes for matrices that can contain data are derived from this.
class GeneralMatrix : public BaseMatrix // declarable matrix types
{
virtual GeneralMatrix* Image() const; // copy of matrix
protected:
int tag_val; // shows whether can reuse
int nrows_val, ncols_val; // dimensions
int storage; // total store required
Real* store; // point to store (0=not set)
GeneralMatrix(); // initialise with no store
GeneralMatrix(ArrayLengthSpecifier); // constructor getting store
void Add(GeneralMatrix*, Real); // sum of GM and Real
void Add(Real); // add Real to this
void NegAdd(GeneralMatrix*, Real); // Real - GM
void NegAdd(Real); // this = this - Real
void Multiply(GeneralMatrix*, Real); // product of GM and Real
void Multiply(Real); // multiply this by Real
void Negate(GeneralMatrix*); // change sign
void Negate(); // change sign
void ReverseElements(); // internal reverse of elements
void ReverseElements(GeneralMatrix*); // reverse order of elements
void operator=(Real); // set matrix to constant
Real* GetStore(); // get store or copy
GeneralMatrix* BorrowStore(GeneralMatrix*, MatrixType);
// temporarily access store
void GetMatrix(const GeneralMatrix*); // used by = and initialise
void Eq(const BaseMatrix&, MatrixType); // used by =
void Eq(const GeneralMatrix&); // version with no conversion
void Eq(const BaseMatrix&, MatrixType, bool);// used by <<
void Eq2(const BaseMatrix&, MatrixType); // cut down version of Eq
int search(const BaseMatrix*) const;
virtual GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
void CheckConversion(const BaseMatrix&); // check conversion OK
void resize(int, int, int); // change dimensions
virtual short SimpleAddOK(const GeneralMatrix*) { return 0; }
// see bandmat.cpp for explanation
virtual void MiniCleanUp()
{ store = 0; storage = 0; nrows_val = 0; ncols_val = 0; tag_val = -1;}
// CleanUp when the data array has already been deleted
void PlusEqual(const GeneralMatrix& gm);
void MinusEqual(const GeneralMatrix& gm);
void PlusEqual(Real f);
void MinusEqual(Real f);
void swap(GeneralMatrix& gm); // swap values
public:
GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
virtual MatrixType type() const = 0; // type of a matrix
MatrixType Type() const { return type(); }
int Nrows() const { return nrows_val; } // get dimensions
int Ncols() const { return ncols_val; }
int Storage() const { return storage; }
Real* Store() const { return store; }
// updated names
int nrows() const { return nrows_val; } // get dimensions
int ncols() const { return ncols_val; }
int size() const { return storage; }
Real* data() { return store; }
const Real* data() const { return store; }
const Real* const_data() const { return store; }
virtual ~GeneralMatrix(); // delete store if set
void tDelete(); // delete if tag_val permits
bool reuse(); // true if tag_val allows reuse
void protect() { tag_val=-1; } // cannot delete or reuse
void Protect() { tag_val=-1; } // cannot delete or reuse
int tag() const { return tag_val; }
int Tag() const { return tag_val; }
bool is_zero() const; // test matrix has all zeros
bool IsZero() const { return is_zero(); } // test matrix has all zeros
void Release() { tag_val=1; } // del store after next use
void Release(int t) { tag_val=t; } // del store after t accesses
void ReleaseAndDelete() { tag_val=0; } // delete matrix after use
void release() { tag_val=1; } // del store after next use
void release(int t) { tag_val=t; } // del store after t accesses
void release_and_delete() { tag_val=0; } // delete matrix after use
void operator<<(const double*); // assignment from an array
void operator<<(const float*); // assignment from an array
void operator<<(const int*); // assignment from an array
void operator<<(const BaseMatrix& X)
{ Eq(X,this->type(),true); } // = without checking type
void inject(const GeneralMatrix&); // copy stored els only
void Inject(const GeneralMatrix& GM) { inject(GM); }
void operator+=(const BaseMatrix&);
void operator-=(const BaseMatrix&);
void operator*=(const BaseMatrix&);
void operator|=(const BaseMatrix&);
void operator&=(const BaseMatrix&);
void operator+=(Real);
void operator-=(Real r) { operator+=(-r); }
void operator*=(Real);
void operator/=(Real r) { operator*=(1.0/r); }
virtual GeneralMatrix* MakeSolver(); // for solving
virtual void Solver(MatrixColX&, const MatrixColX&) {}
virtual void GetRow(MatrixRowCol&) = 0; // Get matrix row
virtual void RestoreRow(MatrixRowCol&) {} // Restore matrix row
virtual void NextRow(MatrixRowCol&); // Go to next row
virtual void GetCol(MatrixRowCol&) = 0; // Get matrix col
virtual void GetCol(MatrixColX&) = 0; // Get matrix col
virtual void RestoreCol(MatrixRowCol&) {} // Restore matrix col
virtual void RestoreCol(MatrixColX&) {} // Restore matrix col
virtual void NextCol(MatrixRowCol&); // Go to next col
virtual void NextCol(MatrixColX&); // Go to next col
Real sum_square() const;
Real sum_absolute_value() const;
Real sum() const;
Real maximum_absolute_value1(int& i) const;
Real minimum_absolute_value1(int& i) const;
Real maximum1(int& i) const;
Real minimum1(int& i) const;
Real maximum_absolute_value() const;
Real maximum_absolute_value2(int& i, int& j) const;
Real minimum_absolute_value() const;
Real minimum_absolute_value2(int& i, int& j) const;
Real maximum() const;
Real maximum2(int& i, int& j) const;
Real minimum() const;
Real minimum2(int& i, int& j) const;
LogAndSign log_determinant() const;
virtual bool IsEqual(const GeneralMatrix&) const;
// same type, same values
void CheckStore() const; // check store is non-zero
virtual void SetParameters(const GeneralMatrix*) {}
// set parameters in GetMatrix
operator ReturnMatrix() const; // for building a ReturnMatrix
ReturnMatrix for_return() const;
ReturnMatrix ForReturn() const;
//virtual bool SameStorageType(const GeneralMatrix& A) const;
//virtual void ReSizeForAdd(const GeneralMatrix& A, const GeneralMatrix& B);
//virtual void ReSizeForSP(const GeneralMatrix& A, const GeneralMatrix& B);
virtual void resize(const GeneralMatrix& A);
virtual void ReSize(const GeneralMatrix& A) { resize(A); }
MatrixInput operator<<(double); // for loading a list
MatrixInput operator<<(float); // for loading a list
MatrixInput operator<<(int f);
// ReturnMatrix Reverse() const; // reverse order of elements
void cleanup(); // to clear store
friend class Matrix;
friend class SquareMatrix;
friend class nricMatrix;
friend class SymmetricMatrix;
friend class UpperTriangularMatrix;
friend class LowerTriangularMatrix;
friend class DiagonalMatrix;
friend class CroutMatrix;
friend class RowVector;
friend class ColumnVector;
friend class BandMatrix;
friend class LowerBandMatrix;
friend class UpperBandMatrix;
friend class SymmetricBandMatrix;
friend class BaseMatrix;
friend class AddedMatrix;
friend class MultipliedMatrix;
friend class SubtractedMatrix;
friend class SPMatrix;
friend class KPMatrix;
friend class ConcatenatedMatrix;
friend class StackedMatrix;
friend class SolvedMatrix;
friend class ShiftedMatrix;
friend class NegShiftedMatrix;
friend class ScaledMatrix;
friend class TransposedMatrix;
friend class ReversedMatrix;
friend class NegatedMatrix;
friend class InvertedMatrix;
friend class RowedMatrix;
friend class ColedMatrix;
friend class DiagedMatrix;
friend class MatedMatrix;
friend class GetSubMatrix;
friend class ReturnMatrix;
friend class LinearEquationSolver;
friend class GenericMatrix;
NEW_DELETE(GeneralMatrix)
};
/// The usual rectangular matrix.
class Matrix : public GeneralMatrix
{
GeneralMatrix* Image() const; // copy of matrix
public:
Matrix() {}
~Matrix() {}
Matrix(int, int); // standard declaration
Matrix(const BaseMatrix&); // evaluate BaseMatrix
void operator=(const BaseMatrix&);
void operator=(Real f) { GeneralMatrix::operator=(f); }
void operator=(const Matrix& m) { Eq(m); }
MatrixType type() const;
Real& operator()(int, int); // access element
Real& element(int, int); // access element
Real operator()(int, int) const; // access element
Real element(int, int) const; // access element
#ifdef SETUP_C_SUBSCRIPTS
Real* operator[](int m) { return store+m*ncols_val; }
const Real* operator[](int m) const { return store+m*ncols_val; }
// following for Numerical Recipes in C++
Matrix(Real, int, int);
Matrix(const Real*, int, int);
#endif
Matrix(const Matrix& gm) : GeneralMatrix() { GetMatrix(&gm); }
GeneralMatrix* MakeSolver();
Real trace() const;
void GetRow(MatrixRowCol&);
void GetCol(MatrixRowCol&);
void GetCol(MatrixColX&);
void RestoreCol(MatrixRowCol&);
void RestoreCol(MatrixColX&);
void NextRow(MatrixRowCol&);
void NextCol(MatrixRowCol&);
void NextCol(MatrixColX&);
virtual void resize(int,int); // change dimensions
// virtual so we will catch it being used in a vector called as a matrix
virtual void resize_keep(int,int);
virtual void ReSize(int m, int n) { resize(m, n); }
void resize(const GeneralMatrix& A);
void ReSize(const GeneralMatrix& A) { resize(A); }
Real maximum_absolute_value2(int& i, int& j) const;
Real minimum_absolute_value2(int& i, int& j) const;
Real maximum2(int& i, int& j) const;
Real minimum2(int& i, int& j) const;
void operator+=(const Matrix& M) { PlusEqual(M); }
void operator-=(const Matrix& M) { MinusEqual(M); }
void operator+=(Real f) { GeneralMatrix::Add(f); }
void operator-=(Real f) { GeneralMatrix::Add(-f); }
void swap(Matrix& gm) { GeneralMatrix::swap((GeneralMatrix&)gm); }
friend Real dotproduct(const Matrix& A, const Matrix& B);
NEW_DELETE(Matrix)
};
/// Square matrix.
class SquareMatrix : public Matrix
{
GeneralMatrix* Image() const; // copy of matrix
public:
SquareMatrix() {}
~SquareMatrix() {}
SquareMatrix(ArrayLengthSpecifier); // standard declaration
SquareMatrix(const BaseMatrix&); // evaluate BaseMatrix
void operator=(const BaseMatrix&);
void operator=(Real f) { GeneralMatrix::operator=(f); }
void operator=(const SquareMatrix& m) { Eq(m); }
void operator=(const Matrix& m);
MatrixType type() const;
SquareMatrix(const SquareMatrix& gm) : Matrix() { GetMatrix(&gm); }
SquareMatrix(const Matrix& gm);
void resize(int); // change dimensions
void ReSize(int m) { resize(m); }
void resize_keep(int);
void resize_keep(int,int);
void resize(int,int); // change dimensions
void ReSize(int m, int n) { resize(m, n); }
void resize(const GeneralMatrix& A);
void ReSize(const GeneralMatrix& A) { resize(A); }
void operator+=(const Matrix& M) { PlusEqual(M); }
void operator-=(const Matrix& M) { MinusEqual(M); }
void operator+=(Real f) { GeneralMatrix::Add(f); }
void operator-=(Real f) { GeneralMatrix::Add(-f); }
void swap(SquareMatrix& gm) { GeneralMatrix::swap((GeneralMatrix&)gm); }
NEW_DELETE(SquareMatrix)
};
/// Rectangular matrix for use with Numerical Recipes in C.
class nricMatrix : public Matrix
{
GeneralMatrix* Image() const; // copy of matrix
Real** row_pointer; // points to rows
void MakeRowPointer(); // build rowpointer
void DeleteRowPointer();
public:
nricMatrix() {}
nricMatrix(int m, int n) // standard declaration
: Matrix(m,n) { MakeRowPointer(); }
nricMatrix(const BaseMatrix& bm) // evaluate BaseMatrix
: Matrix(bm) { MakeRowPointer(); }
void operator=(const BaseMatrix& bm)
{ DeleteRowPointer(); Matrix::operator=(bm); MakeRowPointer(); }
void operator=(Real f) { GeneralMatrix::operator=(f); }
void operator=(const nricMatrix& m)
{ DeleteRowPointer(); Eq(m); MakeRowPointer(); }
void operator<<(const BaseMatrix& X)
{ DeleteRowPointer(); Eq(X,this->type(),true); MakeRowPointer(); }
nricMatrix(const nricMatrix& gm) : Matrix()
{ GetMatrix(&gm); MakeRowPointer(); }
void resize(int m, int n) // change dimensions
{ DeleteRowPointer(); Matrix::resize(m,n); MakeRowPointer(); }
void resize_keep(int m, int n) // change dimensions
{ DeleteRowPointer(); Matrix::resize_keep(m,n); MakeRowPointer(); }
void ReSize(int m, int n) // change dimensions
{ DeleteRowPointer(); Matrix::resize(m,n); MakeRowPointer(); }
void resize(const GeneralMatrix& A);
void ReSize(const GeneralMatrix& A) { resize(A); }
~nricMatrix() { DeleteRowPointer(); }
Real** nric() const { CheckStore(); return row_pointer-1; }
void cleanup(); // to clear store
void MiniCleanUp();
void operator+=(const Matrix& M) { PlusEqual(M); }
void operator-=(const Matrix& M) { MinusEqual(M); }
void operator+=(Real f) { GeneralMatrix::Add(f); }
void operator-=(Real f) { GeneralMatrix::Add(-f); }
void swap(nricMatrix& gm);
NEW_DELETE(nricMatrix)
};
/// Symmetric matrix.
class SymmetricMatrix : public GeneralMatrix
{
GeneralMatrix* Image() const; // copy of matrix
public:
SymmetricMatrix() {}
~SymmetricMatrix() {}
SymmetricMatrix(ArrayLengthSpecifier);
SymmetricMatrix(const BaseMatrix&);
void operator=(const BaseMatrix&);
void operator=(Real f) { GeneralMatrix::operator=(f); }
void operator=(const SymmetricMatrix& m) { Eq(m); }
Real& operator()(int, int); // access element
Real& element(int, int); // access element
Real operator()(int, int) const; // access element
Real element(int, int) const; // access element
#ifdef SETUP_C_SUBSCRIPTS
Real* operator[](int m) { return store+(m*(m+1))/2; }
const Real* operator[](int m) const { return store+(m*(m+1))/2; }
#endif
MatrixType type() const;
SymmetricMatrix(const SymmetricMatrix& gm)
: GeneralMatrix() { GetMatrix(&gm); }
Real sum_square() const;
Real sum_absolute_value() const;
Real sum() const;
Real trace() const;
void GetRow(MatrixRowCol&);
void GetCol(MatrixRowCol&);
void GetCol(MatrixColX&);
void RestoreCol(MatrixRowCol&) {}
void RestoreCol(MatrixColX&);
GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
void resize(int); // change dimensions
void ReSize(int m) { resize(m); }
void resize_keep(int);
void resize(const GeneralMatrix& A);
void ReSize(const GeneralMatrix& A) { resize(A); }
void operator+=(const SymmetricMatrix& M) { PlusEqual(M); }
void operator-=(const SymmetricMatrix& M) { MinusEqual(M); }
void operator+=(Real f) { GeneralMatrix::Add(f); }
void operator-=(Real f) { GeneralMatrix::Add(-f); }
void swap(SymmetricMatrix& gm) { GeneralMatrix::swap((GeneralMatrix&)gm); }
NEW_DELETE(SymmetricMatrix)
};
/// Upper triangular matrix.
class UpperTriangularMatrix : public GeneralMatrix
{
GeneralMatrix* Image() const; // copy of matrix
public:
UpperTriangularMatrix() {}
~UpperTriangularMatrix() {}
UpperTriangularMatrix(ArrayLengthSpecifier);
void operator=(const BaseMatrix&);
void operator=(const UpperTriangularMatrix& m) { Eq(m); }
UpperTriangularMatrix(const BaseMatrix&);
UpperTriangularMatrix(const UpperTriangularMatrix& gm)
: GeneralMatrix() { GetMatrix(&gm); }
void operator=(Real f) { GeneralMatrix::operator=(f); }
Real& operator()(int, int); // access element
Real& element(int, int); // access element
Real operator()(int, int) const; // access element
Real element(int, int) const; // access element
#ifdef SETUP_C_SUBSCRIPTS
Real* operator[](int m) { return store+m*ncols_val-(m*(m+1))/2; }
const Real* operator[](int m) const
{ return store+m*ncols_val-(m*(m+1))/2; }
#endif
MatrixType type() const;
GeneralMatrix* MakeSolver() { return this; } // for solving
void Solver(MatrixColX&, const MatrixColX&);
LogAndSign log_determinant() const;
Real trace() const;
void GetRow(MatrixRowCol&);
void GetCol(MatrixRowCol&);
void GetCol(MatrixColX&);
void RestoreCol(MatrixRowCol&);
void RestoreCol(MatrixColX& c) { RestoreCol((MatrixRowCol&)c); }
void NextRow(MatrixRowCol&);
void resize(int); // change dimensions
void ReSize(int m) { resize(m); }
void resize(const GeneralMatrix& A);
void ReSize(const GeneralMatrix& A) { resize(A); }
void resize_keep(int);
MatrixBandWidth bandwidth() const;
void operator+=(const UpperTriangularMatrix& M) { PlusEqual(M); }
void operator-=(const UpperTriangularMatrix& M) { MinusEqual(M); }
void operator+=(Real f) { GeneralMatrix::operator+=(f); }
void operator-=(Real f) { GeneralMatrix::operator-=(f); }
void swap(UpperTriangularMatrix& gm)
{ GeneralMatrix::swap((GeneralMatrix&)gm); }
NEW_DELETE(UpperTriangularMatrix)
};
/// Lower triangular matrix.
class LowerTriangularMatrix : public GeneralMatrix
{
GeneralMatrix* Image() const; // copy of matrix
public:
LowerTriangularMatrix() {}
~LowerTriangularMatrix() {}
LowerTriangularMatrix(ArrayLengthSpecifier);
LowerTriangularMatrix(const LowerTriangularMatrix& gm)
: GeneralMatrix() { GetMatrix(&gm); }
LowerTriangularMatrix(const BaseMatrix& M);
void operator=(const BaseMatrix&);
void operator=(Real f) { GeneralMatrix::operator=(f); }
void operator=(const LowerTriangularMatrix& m) { Eq(m); }
Real& operator()(int, int); // access element
Real& element(int, int); // access element
Real operator()(int, int) const; // access element
Real element(int, int) const; // access element
#ifdef SETUP_C_SUBSCRIPTS
Real* operator[](int m) { return store+(m*(m+1))/2; }
const Real* operator[](int m) const { return store+(m*(m+1))/2; }
#endif
MatrixType type() const;
GeneralMatrix* MakeSolver() { return this; } // for solving
void Solver(MatrixColX&, const MatrixColX&);
LogAndSign log_determinant() const;
Real trace() const;
void GetRow(MatrixRowCol&);
void GetCol(MatrixRowCol&);
void GetCol(MatrixColX&);
void RestoreCol(MatrixRowCol&);
void RestoreCol(MatrixColX& c) { RestoreCol((MatrixRowCol&)c); }
void NextRow(MatrixRowCol&);
void resize(int); // change dimensions
void ReSize(int m) { resize(m); }
void resize_keep(int);
void resize(const GeneralMatrix& A);
void ReSize(const GeneralMatrix& A) { resize(A); }
MatrixBandWidth bandwidth() const;
void operator+=(const LowerTriangularMatrix& M) { PlusEqual(M); }
void operator-=(const LowerTriangularMatrix& M) { MinusEqual(M); }
void operator+=(Real f) { GeneralMatrix::operator+=(f); }
void operator-=(Real f) { GeneralMatrix::operator-=(f); }
void swap(LowerTriangularMatrix& gm)
{ GeneralMatrix::swap((GeneralMatrix&)gm); }
NEW_DELETE(LowerTriangularMatrix)
};
/// Diagonal matrix.
class DiagonalMatrix : public GeneralMatrix
{
GeneralMatrix* Image() const; // copy of matrix
public:
DiagonalMatrix() {}
~DiagonalMatrix() {}
DiagonalMatrix(ArrayLengthSpecifier);
DiagonalMatrix(const BaseMatrix&);
DiagonalMatrix(const DiagonalMatrix& gm)
: GeneralMatrix() { GetMatrix(&gm); }
void operator=(const BaseMatrix&);
void operator=(Real f) { GeneralMatrix::operator=(f); }
void operator=(const DiagonalMatrix& m) { Eq(m); }
Real& operator()(int, int); // access element
Real& operator()(int); // access element
Real operator()(int, int) const; // access element
Real operator()(int) const;
Real& element(int, int); // access element
Real& element(int); // access element
Real element(int, int) const; // access element
Real element(int) const; // access element
#ifdef SETUP_C_SUBSCRIPTS
Real& operator[](int m) { return store[m]; }
const Real& operator[](int m) const { return store[m]; }
#endif
MatrixType type() const;
LogAndSign log_determinant() const;
Real trace() const;
void GetRow(MatrixRowCol&);
void GetCol(MatrixRowCol&);
void GetCol(MatrixColX&);
void NextRow(MatrixRowCol&);
void NextCol(MatrixRowCol&);
void NextCol(MatrixColX&);
GeneralMatrix* MakeSolver() { return this; } // for solving
void Solver(MatrixColX&, const MatrixColX&);
GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
void resize(int); // change dimensions
void ReSize(int m) { resize(m); }
void resize_keep(int);
void resize(const GeneralMatrix& A);
void ReSize(const GeneralMatrix& A) { resize(A); }
Real* nric() const
{ CheckStore(); return store-1; } // for use by NRIC
MatrixBandWidth bandwidth() const;
// ReturnMatrix Reverse() const; // reverse order of elements
void operator+=(const DiagonalMatrix& M) { PlusEqual(M); }
void operator-=(const DiagonalMatrix& M) { MinusEqual(M); }
void operator+=(Real f) { GeneralMatrix::operator+=(f); }
void operator-=(Real f) { GeneralMatrix::operator-=(f); }
void swap(DiagonalMatrix& gm)
{ GeneralMatrix::swap((GeneralMatrix&)gm); }
NEW_DELETE(DiagonalMatrix)
};
/// Row vector.
class RowVector : public Matrix
{
GeneralMatrix* Image() const; // copy of matrix
public:
RowVector() { nrows_val = 1; }
~RowVector() {}
RowVector(ArrayLengthSpecifier n) : Matrix(1,n.Value()) {}
RowVector(const BaseMatrix&);
RowVector(const RowVector& gm) : Matrix() { GetMatrix(&gm); }
void operator=(const BaseMatrix&);
void operator=(Real f) { GeneralMatrix::operator=(f); }
void operator=(const RowVector& m) { Eq(m); }
Real& operator()(int); // access element
Real& element(int); // access element
Real operator()(int) const; // access element
Real element(int) const; // access element
#ifdef SETUP_C_SUBSCRIPTS
Real& operator[](int m) { return store[m]; }
const Real& operator[](int m) const { return store[m]; }
// following for Numerical Recipes in C++
RowVector(Real a, int n) : Matrix(a, 1, n) {}
RowVector(const Real* a, int n) : Matrix(a, 1, n) {}
#endif
MatrixType type() const;
void GetCol(MatrixRowCol&);
void GetCol(MatrixColX&);
void NextCol(MatrixRowCol&);
void NextCol(MatrixColX&);
void RestoreCol(MatrixRowCol&) {}
void RestoreCol(MatrixColX& c);
GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
void resize(int); // change dimensions
void ReSize(int m) { resize(m); }
void resize_keep(int);
void resize_keep(int,int);
void resize(int,int); // in case access is matrix
void ReSize(int m,int n) { resize(m, n); }
void resize(const GeneralMatrix& A);
void ReSize(const GeneralMatrix& A) { resize(A); }
Real* nric() const
{ CheckStore(); return store-1; } // for use by NRIC
void cleanup(); // to clear store
void MiniCleanUp()
{ store = 0; storage = 0; nrows_val = 1; ncols_val = 0; tag_val = -1; }
// friend ReturnMatrix GetMatrixRow(Matrix& A, int row);
void operator+=(const Matrix& M) { PlusEqual(M); }
void operator-=(const Matrix& M) { MinusEqual(M); }
void operator+=(Real f) { GeneralMatrix::Add(f); }
void operator-=(Real f) { GeneralMatrix::Add(-f); }
void swap(RowVector& gm)
{ GeneralMatrix::swap((GeneralMatrix&)gm); }
NEW_DELETE(RowVector)
};
/// Column vector.
class ColumnVector : public Matrix
{
GeneralMatrix* Image() const; // copy of matrix
public:
ColumnVector() { ncols_val = 1; }
~ColumnVector() {}
ColumnVector(ArrayLengthSpecifier n) : Matrix(n.Value(),1) {}
ColumnVector(const BaseMatrix&);
ColumnVector(const ColumnVector& gm) : Matrix() { GetMatrix(&gm); }
void operator=(const BaseMatrix&);
void operator=(Real f) { GeneralMatrix::operator=(f); }
void operator=(const ColumnVector& m) { Eq(m); }
Real& operator()(int); // access element
Real& element(int); // access element
Real operator()(int) const; // access element
Real element(int) const; // access element
#ifdef SETUP_C_SUBSCRIPTS
Real& operator[](int m) { return store[m]; }
const Real& operator[](int m) const { return store[m]; }
// following for Numerical Recipes in C++
ColumnVector(Real a, int m) : Matrix(a, m, 1) {}
ColumnVector(const Real* a, int m) : Matrix(a, m, 1) {}
#endif
MatrixType type() const;
GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
void resize(int); // change dimensions
void ReSize(int m) { resize(m); }
void resize_keep(int);
void resize_keep(int,int);
void resize(int,int); // in case access is matrix
void ReSize(int m,int n) { resize(m, n); }
void resize(const GeneralMatrix& A);
void ReSize(const GeneralMatrix& A) { resize(A); }
Real* nric() const
{ CheckStore(); return store-1; } // for use by NRIC
void cleanup(); // to clear store
void MiniCleanUp()
{ store = 0; storage = 0; nrows_val = 0; ncols_val = 1; tag_val = -1; }
// ReturnMatrix Reverse() const; // reverse order of elements
void operator+=(const Matrix& M) { PlusEqual(M); }
void operator-=(const Matrix& M) { MinusEqual(M); }
void operator+=(Real f) { GeneralMatrix::Add(f); }
void operator-=(Real f) { GeneralMatrix::Add(-f); }
void swap(ColumnVector& gm)
{ GeneralMatrix::swap((GeneralMatrix&)gm); }
NEW_DELETE(ColumnVector)
};
/// LU matrix.
/// A square matrix decomposed into upper and lower triangular
/// in preparation for inverting or solving equations.
class CroutMatrix : public GeneralMatrix
{
int* indx;
bool d; // number of row swaps = even or odd
bool sing;
void ludcmp();
void get_aux(CroutMatrix&); // for copying indx[] etc
GeneralMatrix* Image() const; // copy of matrix
public:
CroutMatrix(const BaseMatrix&);
CroutMatrix() : indx(0), d(true), sing(true) {}
CroutMatrix(const CroutMatrix&);
void operator=(const CroutMatrix&);
GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
MatrixType type() const;
void lubksb(Real*, int=0);
~CroutMatrix();
GeneralMatrix* MakeSolver() { return this; } // for solving
LogAndSign log_determinant() const;
void Solver(MatrixColX&, const MatrixColX&);
void GetRow(MatrixRowCol&);
void GetCol(MatrixRowCol&);
void GetCol(MatrixColX& c) { GetCol((MatrixRowCol&)c); }
void cleanup(); // to clear store
void MiniCleanUp();
bool IsEqual(const GeneralMatrix&) const;
bool is_singular() const { return sing; }
bool IsSingular() const { return sing; }
const int* const_data_indx() const { return indx; }
bool even_exchanges() const { return d; }
void swap(CroutMatrix& gm);
NEW_DELETE(CroutMatrix)
};
// ***************************** band matrices ***************************/
/// Band matrix.
class BandMatrix : public GeneralMatrix
{
GeneralMatrix* Image() const; // copy of matrix
protected:
void CornerClear() const; // set unused elements to zero
short SimpleAddOK(const GeneralMatrix* gm);
public:
int lower_val, upper_val; // band widths
BandMatrix() { lower_val=0; upper_val=0; CornerClear(); }
~BandMatrix() {}
BandMatrix(int n,int lb,int ub) { resize(n,lb,ub); CornerClear(); }
// standard declaration
BandMatrix(const BaseMatrix&); // evaluate BaseMatrix
void operator=(const BaseMatrix&);
void operator=(Real f) { GeneralMatrix::operator=(f); }
void operator=(const BandMatrix& m) { Eq(m); }
MatrixType type() const;
Real& operator()(int, int); // access element
Real& element(int, int); // access element
Real operator()(int, int) const; // access element
Real element(int, int) const; // access element
#ifdef SETUP_C_SUBSCRIPTS
Real* operator[](int m) { return store+(upper_val+lower_val)*m+lower_val; }
const Real* operator[](int m) const
{ return store+(upper_val+lower_val)*m+lower_val; }
#endif
BandMatrix(const BandMatrix& gm) : GeneralMatrix() { GetMatrix(&gm); }
LogAndSign log_determinant() const;
GeneralMatrix* MakeSolver();
Real trace() const;
Real sum_square() const
{ CornerClear(); return GeneralMatrix::sum_square(); }
Real sum_absolute_value() const
{ CornerClear(); return GeneralMatrix::sum_absolute_value(); }
Real sum() const
{ CornerClear(); return GeneralMatrix::sum(); }
Real maximum_absolute_value() const
{ CornerClear(); return GeneralMatrix::maximum_absolute_value(); }
Real minimum_absolute_value() const
{ int i, j; return GeneralMatrix::minimum_absolute_value2(i, j); }
Real maximum() const { int i, j; return GeneralMatrix::maximum2(i, j); }
Real minimum() const { int i, j; return GeneralMatrix::minimum2(i, j); }
void GetRow(MatrixRowCol&);
void GetCol(MatrixRowCol&);
void GetCol(MatrixColX&);
void RestoreCol(MatrixRowCol&);
void RestoreCol(MatrixColX& c) { RestoreCol((MatrixRowCol&)c); }
void NextRow(MatrixRowCol&);
virtual void resize(int, int, int); // change dimensions
virtual void ReSize(int m, int n, int b) { resize(m, n, b); }
void resize(const GeneralMatrix& A);
void ReSize(const GeneralMatrix& A) { resize(A); }
//bool SameStorageType(const GeneralMatrix& A) const;
//void ReSizeForAdd(const GeneralMatrix& A, const GeneralMatrix& B);
//void ReSizeForSP(const GeneralMatrix& A, const GeneralMatrix& B);
MatrixBandWidth bandwidth() const;
void SetParameters(const GeneralMatrix*);
MatrixInput operator<<(double); // will give error
MatrixInput operator<<(float); // will give error
MatrixInput operator<<(int f);
void operator<<(const double* r); // will give error
void operator<<(const float* r); // will give error
void operator<<(const int* r); // will give error
// the next is included because Zortech and Borland
// cannot find the copy in GeneralMatrix
void operator<<(const BaseMatrix& X) { GeneralMatrix::operator<<(X); }
void swap(BandMatrix& gm);
NEW_DELETE(BandMatrix)
};
/// Upper triangular band matrix.
class UpperBandMatrix : public BandMatrix
{
GeneralMatrix* Image() const; // copy of matrix
public:
UpperBandMatrix() {}
~UpperBandMatrix() {}
UpperBandMatrix(int n, int ubw) // standard declaration
: BandMatrix(n, 0, ubw) {}
UpperBandMatrix(const BaseMatrix&); // evaluate BaseMatrix
void operator=(const BaseMatrix&);
void operator=(Real f) { GeneralMatrix::operator=(f); }
void operator=(const UpperBandMatrix& m) { Eq(m); }
MatrixType type() const;
UpperBandMatrix(const UpperBandMatrix& gm) : BandMatrix() { GetMatrix(&gm); }
GeneralMatrix* MakeSolver() { return this; }
void Solver(MatrixColX&, const MatrixColX&);
LogAndSign log_determinant() const;
void resize(int, int, int); // change dimensions
void ReSize(int m, int n, int b) { resize(m, n, b); }
void resize(int n,int ubw) // change dimensions
{ BandMatrix::resize(n,0,ubw); }
void ReSize(int n,int ubw) // change dimensions
{ BandMatrix::resize(n,0,ubw); }
void resize(const GeneralMatrix& A) { BandMatrix::resize(A); }
void ReSize(const GeneralMatrix& A) { BandMatrix::resize(A); }
Real& operator()(int, int);
Real operator()(int, int) const;
Real& element(int, int);
Real element(int, int) const;
#ifdef SETUP_C_SUBSCRIPTS
Real* operator[](int m) { return store+upper_val*m; }
const Real* operator[](int m) const { return store+upper_val*m; }
#endif
void swap(UpperBandMatrix& gm)
{ BandMatrix::swap((BandMatrix&)gm); }
NEW_DELETE(UpperBandMatrix)
};
/// Lower triangular band matrix.
class LowerBandMatrix : public BandMatrix
{
GeneralMatrix* Image() const; // copy of matrix
public:
LowerBandMatrix() {}
~LowerBandMatrix() {}
LowerBandMatrix(int n, int lbw) // standard declaration
: BandMatrix(n, lbw, 0) {}
LowerBandMatrix(const BaseMatrix&); // evaluate BaseMatrix
void operator=(const BaseMatrix&);
void operator=(Real f) { GeneralMatrix::operator=(f); }
void operator=(const LowerBandMatrix& m) { Eq(m); }
MatrixType type() const;
LowerBandMatrix(const LowerBandMatrix& gm) : BandMatrix() { GetMatrix(&gm); }
GeneralMatrix* MakeSolver() { return this; }
void Solver(MatrixColX&, const MatrixColX&);
LogAndSign log_determinant() const;
void resize(int, int, int); // change dimensions
void ReSize(int m, int n, int b) { resize(m, n, b); }
void resize(int n,int lbw) // change dimensions
{ BandMatrix::resize(n,lbw,0); }
void ReSize(int n,int lbw) // change dimensions
{ BandMatrix::resize(n,lbw,0); }
void resize(const GeneralMatrix& A) { BandMatrix::resize(A); }
void ReSize(const GeneralMatrix& A) { BandMatrix::resize(A); }
Real& operator()(int, int);
Real operator()(int, int) const;
Real& element(int, int);
Real element(int, int) const;
#ifdef SETUP_C_SUBSCRIPTS
Real* operator[](int m) { return store+lower_val*(m+1); }
const Real* operator[](int m) const { return store+lower_val*(m+1); }
#endif
void swap(LowerBandMatrix& gm)
{ BandMatrix::swap((BandMatrix&)gm); }
NEW_DELETE(LowerBandMatrix)
};
/// Symmetric band matrix.
class SymmetricBandMatrix : public GeneralMatrix
{
GeneralMatrix* Image() const; // copy of matrix
void CornerClear() const; // set unused elements to zero
short SimpleAddOK(const GeneralMatrix* gm);
public:
int lower_val; // lower band width
SymmetricBandMatrix() { lower_val=0; CornerClear(); }
~SymmetricBandMatrix() {}
SymmetricBandMatrix(int n, int lb) { resize(n,lb); CornerClear(); }
SymmetricBandMatrix(const BaseMatrix&);
void operator=(const BaseMatrix&);
void operator=(Real f) { GeneralMatrix::operator=(f); }
void operator=(const SymmetricBandMatrix& m) { Eq(m); }
Real& operator()(int, int); // access element
Real& element(int, int); // access element
Real operator()(int, int) const; // access element
Real element(int, int) const; // access element
#ifdef SETUP_C_SUBSCRIPTS
Real* operator[](int m) { return store+lower_val*(m+1); }
const Real* operator[](int m) const { return store+lower_val*(m+1); }
#endif
MatrixType type() const;
SymmetricBandMatrix(const SymmetricBandMatrix& gm)
: GeneralMatrix() { GetMatrix(&gm); }
GeneralMatrix* MakeSolver();
Real sum_square() const;
Real sum_absolute_value() const;
Real sum() const;
Real maximum_absolute_value() const
{ CornerClear(); return GeneralMatrix::maximum_absolute_value(); }
Real minimum_absolute_value() const
{ int i, j; return GeneralMatrix::minimum_absolute_value2(i, j); }
Real maximum() const { int i, j; return GeneralMatrix::maximum2(i, j); }
Real minimum() const { int i, j; return GeneralMatrix::minimum2(i, j); }
Real trace() const;
LogAndSign log_determinant() const;
void GetRow(MatrixRowCol&);
void GetCol(MatrixRowCol&);
void GetCol(MatrixColX&);
void RestoreCol(MatrixRowCol&) {}
void RestoreCol(MatrixColX&);
GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
void resize(int,int); // change dimensions
void ReSize(int m,int b) { resize(m, b); }
void resize(const GeneralMatrix& A);
void ReSize(const GeneralMatrix& A) { resize(A); }
//bool SameStorageType(const GeneralMatrix& A) const;
//void ReSizeForAdd(const GeneralMatrix& A, const GeneralMatrix& B);
//void ReSizeForSP(const GeneralMatrix& A, const GeneralMatrix& B);
MatrixBandWidth bandwidth() const;
void SetParameters(const GeneralMatrix*);
void operator<<(const double* r); // will give error
void operator<<(const float* r); // will give error
void operator<<(const int* r); // will give error
void operator<<(const BaseMatrix& X) { GeneralMatrix::operator<<(X); }
void swap(SymmetricBandMatrix& gm);
NEW_DELETE(SymmetricBandMatrix)
};
/// LU decomposition of a band matrix.
class BandLUMatrix : public GeneralMatrix
{
int* indx;
bool d;
bool sing; // true if singular
Real* store2;
int storage2;
int m1,m2; // lower and upper
void ludcmp();
void get_aux(BandLUMatrix&); // for copying indx[] etc
GeneralMatrix* Image() const; // copy of matrix
public:
BandLUMatrix(const BaseMatrix&);
BandLUMatrix()
: indx(0), d(true), sing(true), store2(0), storage2(0), m1(0), m2(0) {}
BandLUMatrix(const BandLUMatrix&);
void operator=(const BandLUMatrix&);
GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
MatrixType type() const;
void lubksb(Real*, int=0);
~BandLUMatrix();
GeneralMatrix* MakeSolver() { return this; } // for solving
LogAndSign log_determinant() const;
void Solver(MatrixColX&, const MatrixColX&);
void GetRow(MatrixRowCol&);
void GetCol(MatrixRowCol&);
void GetCol(MatrixColX& c) { GetCol((MatrixRowCol&)c); }
void cleanup(); // to clear store
void MiniCleanUp();
bool IsEqual(const GeneralMatrix&) const;
bool is_singular() const { return sing; }
bool IsSingular() const { return sing; }
const Real* const_data2() const { return store2; }
int size2() const { return storage2; }
const int* const_data_indx() const { return indx; }
bool even_exchanges() const { return d; }
MatrixBandWidth bandwidth() const;
void swap(BandLUMatrix& gm);
NEW_DELETE(BandLUMatrix)
};
// ************************** special matrices ****************************
/// Identity matrix.
class IdentityMatrix : public GeneralMatrix
{
GeneralMatrix* Image() const; // copy of matrix
public:
IdentityMatrix() {}
~IdentityMatrix() {}
IdentityMatrix(ArrayLengthSpecifier n) : GeneralMatrix(1)
{ nrows_val = ncols_val = n.Value(); *store = 1; }
IdentityMatrix(const IdentityMatrix& gm)
: GeneralMatrix() { GetMatrix(&gm); }
IdentityMatrix(const BaseMatrix&);
void operator=(const BaseMatrix&);
void operator=(const IdentityMatrix& m) { Eq(m); }
void operator=(Real f) { GeneralMatrix::operator=(f); }
MatrixType type() const;
LogAndSign log_determinant() const;
Real trace() const;
Real sum_square() const;
Real sum_absolute_value() const;
Real sum() const { return trace(); }
void GetRow(MatrixRowCol&);
void GetCol(MatrixRowCol&);
void GetCol(MatrixColX&);
void NextRow(MatrixRowCol&);
void NextCol(MatrixRowCol&);
void NextCol(MatrixColX&);
GeneralMatrix* MakeSolver() { return this; } // for solving
void Solver(MatrixColX&, const MatrixColX&);
GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
void resize(int n);
void ReSize(int n) { resize(n); }
void resize(const GeneralMatrix& A);
void ReSize(const GeneralMatrix& A) { resize(A); }
MatrixBandWidth bandwidth() const;
// ReturnMatrix Reverse() const; // reverse order of elements
void swap(IdentityMatrix& gm)
{ GeneralMatrix::swap((GeneralMatrix&)gm); }
NEW_DELETE(IdentityMatrix)
};
// ************************** GenericMatrix class ************************/
/// A matrix which can be of any GeneralMatrix type.
class GenericMatrix : public BaseMatrix
{
GeneralMatrix* gm;
int search(const BaseMatrix* bm) const;
friend class BaseMatrix;
public:
GenericMatrix() : gm(0) {}
GenericMatrix(const BaseMatrix& bm)
{ gm = ((BaseMatrix&)bm).Evaluate(); gm = gm->Image(); }
GenericMatrix(const GenericMatrix& bm) : BaseMatrix()
{ gm = bm.gm->Image(); }
void operator=(const GenericMatrix&);
void operator=(const BaseMatrix&);
void operator+=(const BaseMatrix&);
void operator-=(const BaseMatrix&);
void operator*=(const BaseMatrix&);
void operator|=(const BaseMatrix&);
void operator&=(const BaseMatrix&);
void operator+=(Real);
void operator-=(Real r) { operator+=(-r); }
void operator*=(Real);
void operator/=(Real r) { operator*=(1.0/r); }
~GenericMatrix() { delete gm; }
void cleanup() { delete gm; gm = 0; }
void Release() { gm->Release(); }
void release() { gm->release(); }
GeneralMatrix* Evaluate(MatrixType = MatrixTypeUnSp);
MatrixBandWidth bandwidth() const;
void swap(GenericMatrix& x);
NEW_DELETE(GenericMatrix)
};
// *************************** temporary classes *************************/
/// Product of two matrices.
/// \internal
class MultipliedMatrix : public BaseMatrix
{
protected:
// if these union statements cause problems, simply remove them
// and declare the items individually
union { const BaseMatrix* bm1; GeneralMatrix* gm1; };
// pointers to summands
union { const BaseMatrix* bm2; GeneralMatrix* gm2; };
MultipliedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
: bm1(bm1x),bm2(bm2x) {}
int search(const BaseMatrix*) const;
friend class BaseMatrix;
friend class GeneralMatrix;
friend class GenericMatrix;
public:
~MultipliedMatrix() {}
GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
MatrixBandWidth bandwidth() const;
NEW_DELETE(MultipliedMatrix)
};
/// Sum of two matrices.
/// \internal
class AddedMatrix : public MultipliedMatrix
{
protected:
AddedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
: MultipliedMatrix(bm1x,bm2x) {}
friend class BaseMatrix;
friend class GeneralMatrix;
friend class GenericMatrix;
public:
~AddedMatrix() {}
GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
MatrixBandWidth bandwidth() const;
NEW_DELETE(AddedMatrix)
};
/// Schur (elementwise) product of two matrices.
/// \internal
class SPMatrix : public AddedMatrix
{
protected:
SPMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
: AddedMatrix(bm1x,bm2x) {}
friend class BaseMatrix;
friend class GeneralMatrix;
friend class GenericMatrix;
public:
~SPMatrix() {}
GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
MatrixBandWidth bandwidth() const;
friend SPMatrix SP(const BaseMatrix&, const BaseMatrix&);
NEW_DELETE(SPMatrix)
};
/// Kronecker product of two matrices.
/// \internal
class KPMatrix : public MultipliedMatrix
{
protected:
KPMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
: MultipliedMatrix(bm1x,bm2x) {}
friend class BaseMatrix;
friend class GeneralMatrix;
friend class GenericMatrix;
public:
~KPMatrix() {}
MatrixBandWidth bandwidth() const;
GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
friend KPMatrix KP(const BaseMatrix&, const BaseMatrix&);
NEW_DELETE(KPMatrix)
};
/// Two matrices horizontally concatenated.
/// \internal
class ConcatenatedMatrix : public MultipliedMatrix
{
protected:
ConcatenatedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
: MultipliedMatrix(bm1x,bm2x) {}
friend class BaseMatrix;
friend class GeneralMatrix;
friend class GenericMatrix;
public:
~ConcatenatedMatrix() {}
MatrixBandWidth bandwidth() const;
GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
NEW_DELETE(ConcatenatedMatrix)
};
/// Two matrices vertically concatenated.
/// \internal
class StackedMatrix : public ConcatenatedMatrix
{
protected:
StackedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
: ConcatenatedMatrix(bm1x,bm2x) {}
friend class BaseMatrix;
friend class GeneralMatrix;
friend class GenericMatrix;
public:
~StackedMatrix() {}
GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
NEW_DELETE(StackedMatrix)
};
/// Inverted matrix times matrix.
/// \internal
class SolvedMatrix : public MultipliedMatrix
{
SolvedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
: MultipliedMatrix(bm1x,bm2x) {}
friend class BaseMatrix;
friend class InvertedMatrix; // for operator*
public:
~SolvedMatrix() {}
GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
MatrixBandWidth bandwidth() const;
NEW_DELETE(SolvedMatrix)
};
/// Difference between two matrices.
/// \internal
class SubtractedMatrix : public AddedMatrix
{
SubtractedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
: AddedMatrix(bm1x,bm2x) {}
friend class BaseMatrix;
friend class GeneralMatrix;
friend class GenericMatrix;
public:
~SubtractedMatrix() {}
GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
NEW_DELETE(SubtractedMatrix)
};
/// Any type of matrix plus Real.
/// \internal
class ShiftedMatrix : public BaseMatrix
{
protected:
union { const BaseMatrix* bm; GeneralMatrix* gm; };
Real f;
ShiftedMatrix(const BaseMatrix* bmx, Real fx) : bm(bmx),f(fx) {}
int search(const BaseMatrix*) const;
friend class BaseMatrix;
friend class GeneralMatrix;
friend class GenericMatrix;
public:
~ShiftedMatrix() {}
GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
friend ShiftedMatrix operator+(Real f, const BaseMatrix& BM);
NEW_DELETE(ShiftedMatrix)
};
/// Real minus matrix.
/// \internal
class NegShiftedMatrix : public ShiftedMatrix
{
protected:
NegShiftedMatrix(Real fx, const BaseMatrix* bmx) : ShiftedMatrix(bmx,fx) {}
friend class BaseMatrix;
friend class GeneralMatrix;
friend class GenericMatrix;
public:
~NegShiftedMatrix() {}
GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
friend NegShiftedMatrix operator-(Real, const BaseMatrix&);
NEW_DELETE(NegShiftedMatrix)
};
/// Any type of matrix times Real.
/// \internal
class ScaledMatrix : public ShiftedMatrix
{
ScaledMatrix(const BaseMatrix* bmx, Real fx) : ShiftedMatrix(bmx,fx) {}
friend class BaseMatrix;
friend class GeneralMatrix;
friend class GenericMatrix;
public:
~ScaledMatrix() {}
GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
MatrixBandWidth bandwidth() const;
friend ScaledMatrix operator*(Real f, const BaseMatrix& BM);
NEW_DELETE(ScaledMatrix)
};
/// Any type of matrix times -1.
/// \internal
class NegatedMatrix : public BaseMatrix
{
protected:
union { const BaseMatrix* bm; GeneralMatrix* gm; };
NegatedMatrix(const BaseMatrix* bmx) : bm(bmx) {}
int search(const BaseMatrix*) const;
private:
friend class BaseMatrix;
public:
~NegatedMatrix() {}
GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
MatrixBandWidth bandwidth() const;
NEW_DELETE(NegatedMatrix)
};
/// Transposed matrix.
/// \internal
class TransposedMatrix : public NegatedMatrix
{
TransposedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
friend class BaseMatrix;
public:
~TransposedMatrix() {}
GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
MatrixBandWidth bandwidth() const;
NEW_DELETE(TransposedMatrix)
};
/// Any type of matrix with order of elements reversed.
/// \internal
class ReversedMatrix : public NegatedMatrix
{
ReversedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
friend class BaseMatrix;
public:
~ReversedMatrix() {}
GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
NEW_DELETE(ReversedMatrix)
};
/// Inverse of matrix.
/// \internal
class InvertedMatrix : public NegatedMatrix
{
InvertedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
public:
~InvertedMatrix() {}
SolvedMatrix operator*(const BaseMatrix&) const; // inverse(A) * B
ScaledMatrix operator*(Real t) const { return BaseMatrix::operator*(t); }
friend class BaseMatrix;
GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
MatrixBandWidth bandwidth() const;
NEW_DELETE(InvertedMatrix)
};
/// Any type of matrix interpreted as a RowVector.
/// \internal
class RowedMatrix : public NegatedMatrix
{
RowedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
friend class BaseMatrix;
public:
~RowedMatrix() {}
GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
MatrixBandWidth bandwidth() const;
NEW_DELETE(RowedMatrix)
};
/// Any type of matrix interpreted as a ColumnVector.
/// \internal
class ColedMatrix : public NegatedMatrix
{
ColedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
friend class BaseMatrix;
public:
~ColedMatrix() {}
GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
MatrixBandWidth bandwidth() const;
NEW_DELETE(ColedMatrix)
};
/// Any type of matrix interpreted as a DiagonalMatrix.
/// \internal
class DiagedMatrix : public NegatedMatrix
{
DiagedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
friend class BaseMatrix;
public:
~DiagedMatrix() {}
GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
MatrixBandWidth bandwidth() const;
NEW_DELETE(DiagedMatrix)
};
/// Any type of matrix interpreted as a (rectangular) Matrix.
/// \internal
class MatedMatrix : public NegatedMatrix
{
int nr, nc;
MatedMatrix(const BaseMatrix* bmx, int nrx, int ncx)
: NegatedMatrix(bmx), nr(nrx), nc(ncx) {}
friend class BaseMatrix;
public:
~MatedMatrix() {}
GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
MatrixBandWidth bandwidth() const;
NEW_DELETE(MatedMatrix)
};
/// A matrix in an "envelope' for return from a function.
/// \internal
class ReturnMatrix : public BaseMatrix
{
GeneralMatrix* gm;
int search(const BaseMatrix*) const;
public:
~ReturnMatrix() {}
GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
friend class BaseMatrix;
ReturnMatrix(const ReturnMatrix& tm) : BaseMatrix(), gm(tm.gm) {}
ReturnMatrix(const GeneralMatrix* gmx) : gm((GeneralMatrix*&)gmx) {}
// ReturnMatrix(GeneralMatrix&);
MatrixBandWidth bandwidth() const;
NEW_DELETE(ReturnMatrix)
};
// ************************** submatrices ******************************/
/// A submatrix of a matrix.
/// \internal
class GetSubMatrix : public NegatedMatrix
{
int row_skip;
int row_number;
int col_skip;
int col_number;
bool IsSym;
GetSubMatrix
(const BaseMatrix* bmx, int rs, int rn, int cs, int cn, bool is)
: NegatedMatrix(bmx),
row_skip(rs), row_number(rn), col_skip(cs), col_number(cn), IsSym(is) {}
void SetUpLHS();
friend class BaseMatrix;
public:
GetSubMatrix(const GetSubMatrix& g)
: NegatedMatrix(g.bm), row_skip(g.row_skip), row_number(g.row_number),
col_skip(g.col_skip), col_number(g.col_number), IsSym(g.IsSym) {}
~GetSubMatrix() {}
GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
void operator=(const BaseMatrix&);
void operator+=(const BaseMatrix&);
void operator-=(const BaseMatrix&);
void operator=(const GetSubMatrix& m) { operator=((const BaseMatrix&)m); }
void operator<<(const BaseMatrix&);
void operator<<(const double*); // copy from array
void operator<<(const float*); // copy from array
void operator<<(const int*); // copy from array
MatrixInput operator<<(double); // for loading a list
MatrixInput operator<<(float); // for loading a list
MatrixInput operator<<(int f);
void operator=(Real); // copy from constant
void operator+=(Real); // add constant
void operator-=(Real r) { operator+=(-r); } // subtract constant
void operator*=(Real); // multiply by constant
void operator/=(Real r) { operator*=(1.0/r); } // divide by constant
void inject(const GeneralMatrix&); // copy stored els only
void Inject(const GeneralMatrix& GM) { inject(GM); }
MatrixBandWidth bandwidth() const;
NEW_DELETE(GetSubMatrix)
};
// ******************** linear equation solving ****************************/
/// A class for finding A.i() * B.
/// This is supposed to choose the appropriate method depending on the
/// type A. Not very satisfactory as it doesn't know about Cholesky for
/// for positive definite matrices.
class LinearEquationSolver : public BaseMatrix
{
GeneralMatrix* gm;
int search(const BaseMatrix*) const { return 0; }
friend class BaseMatrix;
public:
LinearEquationSolver(const BaseMatrix& bm);
~LinearEquationSolver() { delete gm; }
void cleanup() { delete gm; }
GeneralMatrix* Evaluate(MatrixType) { return gm; }
// probably should have an error message if MatrixType != UnSp
NEW_DELETE(LinearEquationSolver)
};
// ************************** matrix input *******************************/
/// Class for reading values into a (small) matrix within a program.
/// \internal
/// Is able to detect a mismatch in the number of elements.
class MatrixInput
{
int n; // number values still to be read
Real* r; // pointer to next location to be read to
public:
MatrixInput(const MatrixInput& mi) : n(mi.n), r(mi.r) {}
MatrixInput(int nx, Real* rx) : n(nx), r(rx) {}
~MatrixInput();
MatrixInput operator<<(double);
MatrixInput operator<<(float);
MatrixInput operator<<(int f);
friend class GeneralMatrix;
};
// **************** a very simple integer array class ********************/
/// A very simple integer array class.
/// A minimal array class to imitate a C style array but giving dynamic storage
/// mostly intended for internal use by newmat.
/// Probably to be replaced by a templated class when I start using templates.
class SimpleIntArray : public Janitor
{
protected:
int* a; ///< pointer to the array
int n; ///< length of the array
public:
SimpleIntArray(int xn); ///< build an array length xn
SimpleIntArray() : a(0), n(0) {} ///< build an array length 0
~SimpleIntArray(); ///< return the space to memory
int& operator[](int i); ///< access element of the array - start at 0
int operator[](int i) const;
///< access element of constant array
void operator=(int ai); ///< set the array equal to a constant
void operator=(const SimpleIntArray& b);
///< copy the elements of an array
SimpleIntArray(const SimpleIntArray& b);
///< make a new array equal to an existing one
int Size() const { return n; }
///< return the size of the array
int size() const { return n; }
///< return the size of the array
int* Data() { return a; } ///< pointer to the data
const int* Data() const { return a; } ///< pointer to the data
int* data() { return a; } ///< pointer to the data
const int* data() const { return a; } ///< pointer to the data
const int* const_data() const { return a; } ///< pointer to the data
void resize(int i, bool keep = false);
///< change length, keep data if keep = true
void ReSize(int i, bool keep = false) { resize(i, keep); }
///< change length, keep data if keep = true
void resize_keep(int i) { resize(i, true); }
///< change length, keep data
void cleanup() { resize(0); } ///< set length to zero
void CleanUp() { resize(0); } ///< set length to zero
NEW_DELETE(SimpleIntArray)
};
// ********************** C subscript classes ****************************
/// Let matrix simulate a C type two dimensional array
class RealStarStar
{
Real** a;
public:
RealStarStar(Matrix& A);
~RealStarStar() { delete [] a; }
operator Real**() { return a; }
};
/// Let matrix simulate a C type const two dimensional array
class ConstRealStarStar
{
const Real** a;
public:
ConstRealStarStar(const Matrix& A);
~ConstRealStarStar() { delete [] a; }
operator const Real**() { return a; }
};
// *************************** exceptions ********************************/
/// Not positive definite exception.
class NPDException : public Runtime_error
{
public:
static unsigned long Select;
NPDException(const GeneralMatrix&);
};
/// Covergence failure exception.
class ConvergenceException : public Runtime_error
{
public:
static unsigned long Select;
ConvergenceException(const GeneralMatrix& A);
ConvergenceException(const char* c);
};
/// Singular matrix exception.
class SingularException : public Runtime_error
{
public:
static unsigned long Select;
SingularException(const GeneralMatrix& A);
};
/// Real overflow exception.
class OverflowException : public Runtime_error
{
public:
static unsigned long Select;
OverflowException(const char* c);
};
/// Miscellaneous exception (details in character string).
class ProgramException : public Logic_error
{
protected:
ProgramException();
public:
static unsigned long Select;
ProgramException(const char* c);
ProgramException(const char* c, const GeneralMatrix&);
ProgramException(const char* c, const GeneralMatrix&, const GeneralMatrix&);
ProgramException(const char* c, MatrixType, MatrixType);
};
/// Index exception.
class IndexException : public Logic_error
{
public:
static unsigned long Select;
IndexException(int i, const GeneralMatrix& A);
IndexException(int i, int j, const GeneralMatrix& A);
// next two are for access via element function
IndexException(int i, const GeneralMatrix& A, bool);
IndexException(int i, int j, const GeneralMatrix& A, bool);
};
/// Cannot convert to vector exception.
class VectorException : public Logic_error
{
public:
static unsigned long Select;
VectorException();
VectorException(const GeneralMatrix& A);
};
/// A matrix is not square exception.
class NotSquareException : public Logic_error
{
public:
static unsigned long Select;
NotSquareException(const GeneralMatrix& A);
NotSquareException();
};
/// Submatrix dimension exception.
class SubMatrixDimensionException : public Logic_error
{
public:
static unsigned long Select;
SubMatrixDimensionException();
};
/// Incompatible dimensions exception.
class IncompatibleDimensionsException : public Logic_error
{
public:
static unsigned long Select;
IncompatibleDimensionsException();
IncompatibleDimensionsException(const GeneralMatrix&);
IncompatibleDimensionsException(const GeneralMatrix&, const GeneralMatrix&);
};
/// Not defined exception.
class NotDefinedException : public Logic_error
{
public:
static unsigned long Select;
NotDefinedException(const char* op, const char* matrix);
};
/// Cannot build matrix with these properties exception.
class CannotBuildException : public Logic_error
{
public:
static unsigned long Select;
CannotBuildException(const char* matrix);
};
/// Internal newmat exception - shouldn't happen.
class InternalException : public Logic_error
{
public:
static unsigned long Select; // for identifying exception
InternalException(const char* c);
};
// ************************ functions ************************************ //
bool operator==(const GeneralMatrix& A, const GeneralMatrix& B);
bool operator==(const BaseMatrix& A, const BaseMatrix& B);
inline bool operator!=(const GeneralMatrix& A, const GeneralMatrix& B)
{ return ! (A==B); }
inline bool operator!=(const BaseMatrix& A, const BaseMatrix& B)
{ return ! (A==B); }
// inequality operators are dummies included for compatibility
// with STL. They throw an exception if actually called.
inline bool operator<=(const BaseMatrix& A, const BaseMatrix&)
{ A.IEQND(); return true; }
inline bool operator>=(const BaseMatrix& A, const BaseMatrix&)
{ A.IEQND(); return true; }
inline bool operator<(const BaseMatrix& A, const BaseMatrix&)
{ A.IEQND(); return true; }
inline bool operator>(const BaseMatrix& A, const BaseMatrix&)
{ A.IEQND(); return true; }
bool is_zero(const BaseMatrix& A);
inline bool IsZero(const BaseMatrix& A) { return is_zero(A); }
Real dotproduct(const Matrix& A, const Matrix& B);
Matrix crossproduct(const Matrix& A, const Matrix& B);
ReturnMatrix crossproduct_rows(const Matrix& A, const Matrix& B);
ReturnMatrix crossproduct_columns(const Matrix& A, const Matrix& B);
inline Real DotProduct(const Matrix& A, const Matrix& B)
{ return dotproduct(A, B); }
inline Matrix CrossProduct(const Matrix& A, const Matrix& B)
{ return crossproduct(A, B); }
inline ReturnMatrix CrossProductRows(const Matrix& A, const Matrix& B)
{ return crossproduct_rows(A, B); }
inline ReturnMatrix CrossProductColumns(const Matrix& A, const Matrix& B)
{ return crossproduct_columns(A, B); }
void newmat_block_copy(int n, Real* from, Real* to);
// ********************* friend functions ******************************** //
// Functions declared as friends - G++ wants them declared externally as well
bool Rectangular(MatrixType a, MatrixType b, MatrixType c);
bool Compare(const MatrixType&, MatrixType&);
Real dotproduct(const Matrix& A, const Matrix& B);
SPMatrix SP(const BaseMatrix&, const BaseMatrix&);
KPMatrix KP(const BaseMatrix&, const BaseMatrix&);
ShiftedMatrix operator+(Real f, const BaseMatrix& BM);
NegShiftedMatrix operator-(Real, const BaseMatrix&);
ScaledMatrix operator*(Real f, const BaseMatrix& BM);
// ********************* inline functions ******************************** //
inline LogAndSign log_determinant(const BaseMatrix& B)
{ return B.log_determinant(); }
inline LogAndSign LogDeterminant(const BaseMatrix& B)
{ return B.log_determinant(); }
inline Real determinant(const BaseMatrix& B)
{ return B.determinant(); }
inline Real Determinant(const BaseMatrix& B)
{ return B.determinant(); }
inline Real sum_square(const BaseMatrix& B) { return B.sum_square(); }
inline Real SumSquare(const BaseMatrix& B) { return B.sum_square(); }
inline Real norm_Frobenius(const BaseMatrix& B) { return B.norm_Frobenius(); }
inline Real norm_frobenius(const BaseMatrix& B) { return B.norm_Frobenius(); }
inline Real NormFrobenius(const BaseMatrix& B) { return B.norm_Frobenius(); }
inline Real trace(const BaseMatrix& B) { return B.trace(); }
inline Real Trace(const BaseMatrix& B) { return B.trace(); }
inline Real sum_absolute_value(const BaseMatrix& B)
{ return B.sum_absolute_value(); }
inline Real SumAbsoluteValue(const BaseMatrix& B)
{ return B.sum_absolute_value(); }
inline Real sum(const BaseMatrix& B)
{ return B.sum(); }
inline Real Sum(const BaseMatrix& B)
{ return B.sum(); }
inline Real maximum_absolute_value(const BaseMatrix& B)
{ return B.maximum_absolute_value(); }
inline Real MaximumAbsoluteValue(const BaseMatrix& B)
{ return B.maximum_absolute_value(); }
inline Real minimum_absolute_value(const BaseMatrix& B)
{ return B.minimum_absolute_value(); }
inline Real MinimumAbsoluteValue(const BaseMatrix& B)
{ return B.minimum_absolute_value(); }
inline Real maximum(const BaseMatrix& B) { return B.maximum(); }
inline Real Maximum(const BaseMatrix& B) { return B.maximum(); }
inline Real minimum(const BaseMatrix& B) { return B.minimum(); }
inline Real Minimum(const BaseMatrix& B) { return B.minimum(); }
inline Real norm1(const BaseMatrix& B) { return B.norm1(); }
inline Real Norm1(const BaseMatrix& B) { return B.norm1(); }
inline Real norm1(RowVector& RV) { return RV.maximum_absolute_value(); }
inline Real Norm1(RowVector& RV) { return RV.maximum_absolute_value(); }
inline Real norm_infinity(const BaseMatrix& B) { return B.norm_infinity(); }
inline Real NormInfinity(const BaseMatrix& B) { return B.norm_infinity(); }
inline Real norm_infinity(ColumnVector& CV)
{ return CV.maximum_absolute_value(); }
inline Real NormInfinity(ColumnVector& CV)
{ return CV.maximum_absolute_value(); }
inline bool IsZero(const GeneralMatrix& A) { return A.IsZero(); }
inline bool is_zero(const GeneralMatrix& A) { return A.is_zero(); }
inline MatrixInput MatrixInput::operator<<(int f) { return *this << (Real)f; }
inline MatrixInput GeneralMatrix::operator<<(int f) { return *this << (Real)f; }
inline MatrixInput BandMatrix::operator<<(int f) { return *this << (Real)f; }
inline MatrixInput GetSubMatrix::operator<<(int f) { return *this << (Real)f; }
inline ReversedMatrix BaseMatrix::Reverse() const { return reverse(); }
inline RowedMatrix BaseMatrix::AsRow() const { return as_row(); }
inline ColedMatrix BaseMatrix::AsColumn() const { return as_column(); }
inline DiagedMatrix BaseMatrix::AsDiagonal() const { return as_diagonal(); }
inline MatedMatrix BaseMatrix::AsMatrix(int m, int n) const
{ return as_matrix(m, n); }
inline GetSubMatrix BaseMatrix::SubMatrix(int fr, int lr, int fc, int lc) const
{ return submatrix(fr, lr, fc, lc); }
inline GetSubMatrix BaseMatrix::SymSubMatrix(int f, int l) const
{ return sym_submatrix(f, l); }
inline GetSubMatrix BaseMatrix::Row(int f) const { return row(f); }
inline GetSubMatrix BaseMatrix::Rows(int f, int l) const { return rows(f, l); }
inline GetSubMatrix BaseMatrix::Column(int f) const { return column(f); }
inline GetSubMatrix BaseMatrix::Columns(int f, int l) const
{ return columns(f, l); }
inline Real BaseMatrix::AsScalar() const { return as_scalar(); }
inline ReturnMatrix GeneralMatrix::ForReturn() const { return for_return(); }
inline void swap(Matrix& A, Matrix& B) { A.swap(B); }
inline void swap(SquareMatrix& A, SquareMatrix& B) { A.swap(B); }
inline void swap(nricMatrix& A, nricMatrix& B) { A.swap(B); }
inline void swap(UpperTriangularMatrix& A, UpperTriangularMatrix& B)
{ A.swap(B); }
inline void swap(LowerTriangularMatrix& A, LowerTriangularMatrix& B)
{ A.swap(B); }
inline void swap(SymmetricMatrix& A, SymmetricMatrix& B) { A.swap(B); }
inline void swap(DiagonalMatrix& A, DiagonalMatrix& B) { A.swap(B); }
inline void swap(RowVector& A, RowVector& B) { A.swap(B); }
inline void swap(ColumnVector& A, ColumnVector& B) { A.swap(B); }
inline void swap(CroutMatrix& A, CroutMatrix& B) { A.swap(B); }
inline void swap(BandMatrix& A, BandMatrix& B) { A.swap(B); }
inline void swap(UpperBandMatrix& A, UpperBandMatrix& B) { A.swap(B); }
inline void swap(LowerBandMatrix& A, LowerBandMatrix& B) { A.swap(B); }
inline void swap(SymmetricBandMatrix& A, SymmetricBandMatrix& B) { A.swap(B); }
inline void swap(BandLUMatrix& A, BandLUMatrix& B) { A.swap(B); }
inline void swap(IdentityMatrix& A, IdentityMatrix& B) { A.swap(B); }
inline void swap(GenericMatrix& A, GenericMatrix& B) { A.swap(B); }
#ifdef OPT_COMPATIBLE // for compatibility with opt++
inline Real Norm2(const ColumnVector& CV) { return CV.norm_Frobenius(); }
inline Real Dot(ColumnVector& CV1, ColumnVector& CV2)
{ return dotproduct(CV1, CV2); }
#endif
#ifdef use_namespace
}
#endif
#endif
// body file: newmat1.cpp
// body file: newmat2.cpp
// body file: newmat3.cpp
// body file: newmat4.cpp
// body file: newmat5.cpp
// body file: newmat6.cpp
// body file: newmat7.cpp
// body file: newmat8.cpp
// body file: newmatex.cpp
// body file: bandmat.cpp
// body file: submat.cpp
///@}
syntax highlighted by Code2HTML, v. 0.9.1