Tekkotsu Homepage
Demos
Overview
Downloads
Dev. Resources
Reference
Credits

newmat1.cpp

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 // Skew needs more work <<<<<<<<<
00024 
00025 // all operations to return MatrixTypes which correspond to valid types
00026 // of matrices.
00027 // Eg: if it has the Diagonal attribute, then it must also have
00028 // Upper, Lower, Band, Square and Symmetric.
00029 
00030 
00031 MatrixType MatrixType::operator*(const MatrixType& mt) const
00032 {
00033    REPORT
00034    int a = attribute & mt.attribute & ~(Symmetric | Skew);
00035    a |= (a & Diagonal) * 63;                   // recognise diagonal
00036    return MatrixType(a);
00037 }
00038 
00039 MatrixType MatrixType::SP(const MatrixType& mt) const
00040 // elementwise product
00041 // Lower, Upper, Diag, Band if only one is
00042 // Symmetric, Ones, Valid (and Real) if both are
00043 // Need to include Lower & Upper => Diagonal
00044 // Will need to include both Skew => Symmetric
00045 {
00046    REPORT
00047    int a = ((attribute | mt.attribute) & ~(Symmetric + Skew + Valid + Ones))
00048       | (attribute & mt.attribute);
00049    if ((a & Lower) != 0  &&  (a & Upper) != 0) a |= Diagonal;
00050    if ((attribute & Skew) != 0)
00051    {
00052       if ((mt.attribute & Symmetric) != 0) a |= Skew;  
00053       if ((mt.attribute & Skew) != 0) { a &= ~Skew; a |= Symmetric; }
00054    }
00055    else if ((mt.attribute & Skew) != 0 && (attribute & Symmetric) != 0)
00056       a |= Skew;  
00057    a |= (a & Diagonal) * 63;                   // recognise diagonal
00058    return MatrixType(a);
00059 }
00060 
00061 MatrixType MatrixType::KP(const MatrixType& mt) const
00062 // Kronecker product
00063 // Lower, Upper, Diag, Symmetric, Band, Valid if both are
00064 // Band if LHS is band & other is square 
00065 // Ones is complicated so leave this out
00066 {
00067    REPORT
00068    int a = (attribute & mt.attribute)  & ~Ones;
00069    if ((attribute & Band) != 0 && (mt.attribute & Square) != 0)
00070       a |= Band;
00071    //int a = ((attribute & mt.attribute) | (attribute & Band)) & ~Ones;
00072 
00073    return MatrixType(a);
00074 }
00075 
00076 MatrixType MatrixType::i() const               // type of inverse
00077 {
00078    REPORT
00079    int a = attribute & ~(Band+LUDeco);
00080    a |= (a & Diagonal) * 63;                   // recognise diagonal
00081    return MatrixType(a);
00082 }
00083 
00084 MatrixType MatrixType::t() const
00085 // swap lower and upper attributes
00086 // assume Upper is in bit above Lower
00087 {
00088    REPORT
00089    int a = attribute;
00090    a ^= (((a >> 1) ^ a) & Lower) * 3;
00091    return MatrixType(a);
00092 }
00093 
00094 MatrixType MatrixType::MultRHS() const
00095 {
00096    REPORT
00097    // remove symmetric attribute unless diagonal
00098    return (attribute >= Dg) ? attribute : (attribute & ~Symmetric);
00099 }
00100 
00101 // this is used for deciding type of multiplication
00102 bool Rectangular(MatrixType a, MatrixType b, MatrixType c)
00103 {
00104    REPORT
00105    return
00106       ((a.attribute | b.attribute | c.attribute)
00107       & ~(MatrixType::Valid | MatrixType::Square)) == 0;
00108 }
00109 
00110 const char* MatrixType::Value() const
00111 {
00112 // make a string with the name of matrix with the given attributes
00113    switch (attribute)
00114    {
00115    case Valid:                              REPORT return "Rect ";
00116    case Valid+Square:                       REPORT return "Squ  ";
00117    case Valid+Symmetric+Square:             REPORT return "Sym  ";
00118    case Valid+Skew+Square:                  REPORT return "Skew ";
00119    case Valid+Band+Square:                  REPORT return "Band ";
00120    case Valid+Symmetric+Band+Square:        REPORT return "SmBnd";
00121    case Valid+Upper+Square:                 REPORT return "UT   ";
00122    case Valid+Diagonal+Symmetric+Band+Upper+Lower+Square:
00123                                             REPORT return "Diag ";
00124    case Valid+Diagonal+Symmetric+Band+Upper+Lower+Ones+Square:
00125                                             REPORT return "Ident";
00126    case Valid+Band+Upper+Square:            REPORT return "UpBnd";
00127    case Valid+Lower+Square:                 REPORT return "LT   ";
00128    case Valid+Band+Lower+Square:            REPORT return "LwBnd";
00129    default:
00130       REPORT
00131       if (!(attribute & Valid))             return "UnSp ";
00132       if (attribute & LUDeco)
00133          return (attribute & Band) ?     "BndLU" : "Crout";
00134                                             return "?????";
00135    }
00136 }
00137 
00138 
00139 GeneralMatrix* MatrixType::New(int nr, int nc, BaseMatrix* bm) const
00140 {
00141 // make a new matrix with the given attributes
00142 
00143    Tracer tr("New"); GeneralMatrix* gm=0;   // initialised to keep gnu happy
00144    switch (attribute)
00145    {
00146    case Valid:
00147       REPORT
00148       if (nc==1) { gm = new ColumnVector(nr); break; }
00149       if (nr==1) { gm = new RowVector(nc); break; }
00150       gm = new Matrix(nr, nc); break;
00151 
00152    case Valid+Square:
00153       REPORT
00154       if (nc!=nr) { Throw(NotSquareException()); }
00155       gm = new SquareMatrix(nr); break;
00156 
00157    case Valid+Symmetric+Square:
00158       REPORT gm = new SymmetricMatrix(nr); break;
00159 
00160    case Valid+Band+Square:
00161       {
00162          REPORT
00163          MatrixBandWidth bw = bm->BandWidth();
00164          gm = new BandMatrix(nr,bw.lower,bw.upper); break;
00165       }
00166 
00167    case Valid+Symmetric+Band+Square:
00168       REPORT gm = new SymmetricBandMatrix(nr,bm->BandWidth().lower); break;
00169 
00170    case Valid+Upper+Square:
00171       REPORT gm = new UpperTriangularMatrix(nr); break;
00172 
00173    case Valid+Diagonal+Symmetric+Band+Upper+Lower+Square:
00174       REPORT gm = new DiagonalMatrix(nr); break;
00175 
00176    case Valid+Band+Upper+Square:
00177       REPORT gm = new UpperBandMatrix(nr,bm->BandWidth().upper); break;
00178 
00179    case Valid+Lower+Square:
00180       REPORT gm = new LowerTriangularMatrix(nr); break;
00181 
00182    case Valid+Band+Lower+Square:
00183       REPORT gm = new LowerBandMatrix(nr,bm->BandWidth().lower); break;
00184 
00185    case Valid+Diagonal+Symmetric+Band+Upper+Lower+Ones+Square:
00186       REPORT gm = new IdentityMatrix(nr); break;
00187 
00188    default:
00189       Throw(ProgramException("Invalid matrix type"));
00190    }
00191    
00192    MatrixErrorNoSpace(gm); gm->Protect(); return gm;
00193 }
00194 
00195 
00196 #ifdef use_namespace
00197 }
00198 #endif
00199 

newmat11b
Generated Mon May 9 04:54:18 2016 by Doxygen 1.6.3