newmat1.cc
Go to the documentation of this file.
00001 //$$ newmat1.cpp   Matrix type functions
00002 
00003 // Copyright (C) 1991,2,3,4: R B Davies
00004 
00005 //#define WANT_STREAM
00006 
00007 #include "newmat.h"
00008 
00009 #ifdef use_namespace
00010 namespace NEWMAT {
00011 #endif
00012 
00013 #ifdef DO_REPORT
00014 #define REPORT { static ExeCounter ExeCount(__LINE__,1); ++ExeCount; }
00015 #else
00016 #define REPORT {}
00017 #endif
00018 
00019 
00020 /************************* MatrixType functions *****************************/
00021 
00022 
00023 // all operations to return MatrixTypes which correspond to valid types
00024 // of matrices.
00025 // Eg: if it has the Diagonal attribute, then it must also have
00026 // Upper, Lower, Band and Symmetric.
00027 
00028 
00029 MatrixType MatrixType::operator*(const MatrixType& mt) const
00030 {
00031    REPORT
00032    int a = attribute & mt.attribute & ~Symmetric;
00033    a |= (a & Diagonal) * 31;                   // recognise diagonal
00034    return MatrixType(a);
00035 }
00036 
00037 MatrixType MatrixType::SP(const MatrixType& mt) const
00038 // elementwise product
00039 // Lower, Upper, Diag, Band if only one is
00040 // Symmetric, Ones, Valid (and Real) if both are
00041 // Need to include Lower & Upper => Diagonal
00042 // Will need to include both Skew => Symmetric
00043 {
00044    REPORT
00045    int a = ((attribute | mt.attribute) & ~(Symmetric + Valid + Ones))
00046       | (attribute & mt.attribute);
00047    if ((a & Lower) != 0  &&  (a & Upper) != 0) a |= Diagonal;
00048    a |= (a & Diagonal) * 31;                   // recognise diagonal
00049    return MatrixType(a);
00050 }
00051 
00052 MatrixType MatrixType::KP(const MatrixType& mt) const
00053 // Kronecker product
00054 // Lower, Upper, Diag, Symmetric, Band, Valid if both are
00055 // Do not treat Band separately
00056 // Ones is complicated so leave this out
00057 {
00058    REPORT
00059    int a = (attribute & mt.attribute) & ~Ones;
00060    return MatrixType(a);
00061 }
00062 
00063 MatrixType MatrixType::i() const               // type of inverse
00064 {
00065    REPORT
00066    int a = attribute & ~(Band+LUDeco);
00067    a |= (a & Diagonal) * 31;                   // recognise diagonal
00068    return MatrixType(a);
00069 }
00070 
00071 MatrixType MatrixType::t() const
00072 // swap lower and upper attributes
00073 // assume Upper is in bit above Lower
00074 {
00075    REPORT
00076    int a = attribute;
00077    a ^= (((a >> 1) ^ a) & Lower) * 3;
00078    return MatrixType(a);
00079 }
00080 
00081 MatrixType MatrixType::MultRHS() const
00082 {
00083    REPORT
00084    // remove symmetric attribute unless diagonal
00085    return (attribute >= Dg) ? attribute : (attribute & ~Symmetric);
00086 }
00087 
00088 bool Rectangular(MatrixType a, MatrixType b, MatrixType c)
00089 {
00090    REPORT
00091    return
00092       ((a.attribute | b.attribute | c.attribute) & ~MatrixType::Valid) == 0;
00093 }
00094 
00095 const char* MatrixType::Value() const
00096 {
00097 // make a string with the name of matrix with the given attributes
00098    switch (attribute)
00099    {
00100    case Valid:                              REPORT return "Rect ";
00101    case Valid+Symmetric:                    REPORT return "Sym  ";
00102    case Valid+Band:                         REPORT return "Band ";
00103    case Valid+Symmetric+Band:               REPORT return "SmBnd";
00104    case Valid+Upper:                        REPORT return "UT   ";
00105    case Valid+Diagonal+Symmetric+Band+Upper+Lower:
00106                                             REPORT return "Diag ";
00107    case Valid+Diagonal+Symmetric+Band+Upper+Lower+Ones:
00108                                             REPORT return "Ident";
00109    case Valid+Band+Upper:                   REPORT return "UpBnd";
00110    case Valid+Lower:                        REPORT return "LT   ";
00111    case Valid+Band+Lower:                   REPORT return "LwBnd";
00112    default:
00113       REPORT
00114       if (!(attribute & Valid))             return "UnSp ";
00115       if (attribute & LUDeco)
00116          return (attribute & Band) ?     "BndLU" : "Crout";
00117                                             return "?????";
00118    }
00119 }
00120 
00121 
00122 GeneralMatrix* MatrixType::New(int nr, int nc, BaseMatrix* bm) const
00123 {
00124 // make a new matrix with the given attributes
00125 
00126    Tracer tr("New"); GeneralMatrix* gm=0;   // initialised to keep gnu happy
00127    switch (attribute)
00128    {
00129    case Valid:
00130       REPORT
00131       if (nc==1) { gm = new ColumnVector(nr); break; }
00132       if (nr==1) { gm = new RowVector(nc); break; }
00133       gm = new Matrix(nr, nc); break;
00134 
00135    case Valid+Symmetric:
00136       REPORT gm = new SymmetricMatrix(nr); break;
00137 
00138    case Valid+Band:
00139       {
00140          REPORT
00141          MatrixBandWidth bw = bm->BandWidth();
00142          gm = new BandMatrix(nr,bw.lower,bw.upper); break;
00143       }
00144 
00145    case Valid+Symmetric+Band:
00146       REPORT gm = new SymmetricBandMatrix(nr,bm->BandWidth().lower); break;
00147 
00148    case Valid+Upper:
00149       REPORT gm = new UpperTriangularMatrix(nr); break;
00150 
00151    case Valid+Diagonal+Symmetric+Band+Upper+Lower:
00152       REPORT gm = new DiagonalMatrix(nr); break;
00153 
00154    case Valid+Band+Upper:
00155       REPORT gm = new UpperBandMatrix(nr,bm->BandWidth().upper); break;
00156 
00157    case Valid+Lower:
00158       REPORT gm = new LowerTriangularMatrix(nr); break;
00159 
00160    case Valid+Band+Lower:
00161       REPORT gm = new LowerBandMatrix(nr,bm->BandWidth().lower); break;
00162 
00163    case Valid+Diagonal+Symmetric+Band+Upper+Lower+Ones:
00164       REPORT gm = new IdentityMatrix(nr); break;
00165 
00166    default:
00167       Throw(ProgramException("Invalid matrix type"));
00168    }
00169    
00170    MatrixErrorNoSpace(gm); gm->Protect(); return gm;
00171 }
00172 
00173 
00174 #ifdef use_namespace
00175 }
00176 #endif
00177 


rl_agent
Author(s): Todd Hester
autogenerated on Thu Jun 6 2019 22:00:13