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


kni
Author(s): Martin Günther
autogenerated on Mon Aug 14 2017 02:44:12