Go to the documentation of this file.00001
00002
00003
00008
00009
00010
00011
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
00027
00028
00029
00030
00031
00032
00033
00034
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;
00042 return MatrixType(a);
00043 }
00044
00045 MatrixType MatrixType::SP(const MatrixType& mt) const
00046
00047
00048
00049
00050
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;
00064 return MatrixType(a);
00065 }
00066
00067 MatrixType MatrixType::KP(const MatrixType& mt) const
00068
00069
00070
00071
00072 {
00073 REPORT
00074 int a = (attribute & mt.attribute) & ~Ones;
00075 if ((attribute & Band) != 0 && (mt.attribute & Square) != 0)
00076 a |= Band;
00077
00078
00079 return MatrixType(a);
00080 }
00081
00082 MatrixType MatrixType::i() const
00083 {
00084 REPORT
00085 int a = attribute & ~(Band+LUDeco);
00086 a |= (a & Diagonal) * 63;
00087 return MatrixType(a);
00088 }
00089
00090 MatrixType MatrixType::t() const
00091
00092
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
00104 return (attribute >= Dg) ? attribute : (attribute & ~Symmetric);
00105 }
00106
00107
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
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
00149
00150 Tracer tr("New"); GeneralMatrix* gm=0;
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