00001
00002
00003
00004
00005
00006
00007 #include "include.h"
00008
00009 #include "newmat.h"
00010 #include "newmatrc.h"
00011
00012 #ifdef use_namespace
00013 namespace NEWMAT {
00014 #endif
00015
00016
00017 #ifdef DO_REPORT
00018 #define REPORT { static ExeCounter ExeCount(__LINE__,5); ++ExeCount; }
00019 #else
00020 #define REPORT {}
00021 #endif
00022
00023
00024
00025
00026
00027 GeneralMatrix* GeneralMatrix::Transpose(TransposedMatrix* tm, MatrixType mt)
00028 {
00029 GeneralMatrix* gm1;
00030
00031 if (Compare(Type().t(),mt))
00032 {
00033 REPORT
00034 gm1 = mt.New(ncols_value,nrows_value,tm);
00035 for (int i=0; i<ncols_value; i++)
00036 {
00037 MatrixRow mr(gm1, StoreOnExit+DirectPart, i);
00038 MatrixCol mc(this, mr.Data(), LoadOnEntry, i);
00039 }
00040 }
00041 else
00042 {
00043 REPORT
00044 gm1 = mt.New(ncols_value,nrows_value,tm);
00045 MatrixRow mr(this, LoadOnEntry);
00046 MatrixCol mc(gm1, StoreOnExit+DirectPart);
00047 int i = nrows_value;
00048 while (i--) { mc.Copy(mr); mr.Next(); mc.Next(); }
00049 }
00050 tDelete(); gm1->ReleaseAndDelete(); return gm1;
00051 }
00052
00053 GeneralMatrix* SymmetricMatrix::Transpose(TransposedMatrix*, MatrixType mt)
00054 { REPORT return Evaluate(mt); }
00055
00056
00057 GeneralMatrix* DiagonalMatrix::Transpose(TransposedMatrix*, MatrixType mt)
00058 { REPORT return Evaluate(mt); }
00059
00060 GeneralMatrix* ColumnVector::Transpose(TransposedMatrix*, MatrixType mt)
00061 {
00062 REPORT
00063 GeneralMatrix* gmx = new RowVector; MatrixErrorNoSpace(gmx);
00064 gmx->nrows_value = 1; gmx->ncols_value = gmx->storage = storage;
00065 return BorrowStore(gmx,mt);
00066 }
00067
00068 GeneralMatrix* RowVector::Transpose(TransposedMatrix*, MatrixType mt)
00069 {
00070 REPORT
00071 GeneralMatrix* gmx = new ColumnVector; MatrixErrorNoSpace(gmx);
00072 gmx->ncols_value = 1; gmx->nrows_value = gmx->storage = storage;
00073 return BorrowStore(gmx,mt);
00074 }
00075
00076 GeneralMatrix* IdentityMatrix::Transpose(TransposedMatrix*, MatrixType mt)
00077 { REPORT return Evaluate(mt); }
00078
00079 GeneralMatrix* GeneralMatrix::Evaluate(MatrixType mt)
00080 {
00081 if (Compare(this->Type(),mt)) { REPORT return this; }
00082 REPORT
00083 GeneralMatrix* gmx = mt.New(nrows_value,ncols_value,this);
00084 MatrixRow mr(this, LoadOnEntry);
00085 MatrixRow mrx(gmx, StoreOnExit+DirectPart);
00086 int i=nrows_value;
00087 while (i--) { mrx.Copy(mr); mrx.Next(); mr.Next(); }
00088 tDelete(); gmx->ReleaseAndDelete(); return gmx;
00089 }
00090
00091 GeneralMatrix* GenericMatrix::Evaluate(MatrixType mt)
00092 { REPORT return gm->Evaluate(mt); }
00093
00094 GeneralMatrix* ShiftedMatrix::Evaluate(MatrixType mt)
00095 {
00096 gm=(const_cast<BaseMatrix*&>(bm))->Evaluate();
00097 int nr=gm->Nrows(); int nc=gm->Ncols();
00098 Compare(gm->Type().AddEqualEl(),mt);
00099 if (!(mt==gm->Type()))
00100 {
00101 REPORT
00102 GeneralMatrix* gmx = mt.New(nr,nc,this);
00103 MatrixRow mr(gm, LoadOnEntry);
00104 MatrixRow mrx(gmx, StoreOnExit+DirectPart);
00105 while (nr--) { mrx.Add(mr,f); mrx.Next(); mr.Next(); }
00106 gmx->ReleaseAndDelete(); gm->tDelete();
00107 return gmx;
00108 }
00109 else if (gm->reuse())
00110 {
00111 REPORT gm->Add(f);
00112 return gm;
00113 }
00114 else
00115 {
00116 REPORT GeneralMatrix* gmy = gm->Type().New(nr,nc,this);
00117 gmy->ReleaseAndDelete(); gmy->Add(gm,f);
00118 return gmy;
00119 }
00120 }
00121
00122 GeneralMatrix* NegShiftedMatrix::Evaluate(MatrixType mt)
00123 {
00124 gm=(const_cast<BaseMatrix*&>(bm))->Evaluate();
00125 int nr=gm->Nrows(); int nc=gm->Ncols();
00126 Compare(gm->Type().AddEqualEl(),mt);
00127 if (!(mt==gm->Type()))
00128 {
00129 REPORT
00130 GeneralMatrix* gmx = mt.New(nr,nc,this);
00131 MatrixRow mr(gm, LoadOnEntry);
00132 MatrixRow mrx(gmx, StoreOnExit+DirectPart);
00133 while (nr--) { mrx.NegAdd(mr,f); mrx.Next(); mr.Next(); }
00134 gmx->ReleaseAndDelete(); gm->tDelete();
00135 return gmx;
00136 }
00137 else if (gm->reuse())
00138 {
00139 REPORT gm->NegAdd(f);
00140 return gm;
00141 }
00142 else
00143 {
00144 REPORT GeneralMatrix* gmy = gm->Type().New(nr,nc,this);
00145 gmy->ReleaseAndDelete(); gmy->NegAdd(gm,f);
00146 return gmy;
00147 }
00148 }
00149
00150 GeneralMatrix* ScaledMatrix::Evaluate(MatrixType mt)
00151 {
00152 gm=(const_cast<BaseMatrix*&>(bm))->Evaluate();
00153 int nr=gm->Nrows(); int nc=gm->Ncols();
00154 if (Compare(gm->Type(),mt))
00155 {
00156 if (gm->reuse())
00157 {
00158 REPORT gm->Multiply(f);
00159 return gm;
00160 }
00161 else
00162 {
00163 REPORT GeneralMatrix* gmx = gm->Type().New(nr,nc,this);
00164 gmx->ReleaseAndDelete(); gmx->Multiply(gm,f);
00165 return gmx;
00166 }
00167 }
00168 else
00169 {
00170 REPORT
00171 GeneralMatrix* gmx = mt.New(nr,nc,this);
00172 MatrixRow mr(gm, LoadOnEntry);
00173 MatrixRow mrx(gmx, StoreOnExit+DirectPart);
00174 while (nr--) { mrx.Multiply(mr,f); mrx.Next(); mr.Next(); }
00175 gmx->ReleaseAndDelete(); gm->tDelete();
00176 return gmx;
00177 }
00178 }
00179
00180 GeneralMatrix* NegatedMatrix::Evaluate(MatrixType mt)
00181 {
00182 gm=(const_cast<BaseMatrix*&>(bm))->Evaluate();
00183 int nr=gm->Nrows(); int nc=gm->Ncols();
00184 if (Compare(gm->Type(),mt))
00185 {
00186 if (gm->reuse())
00187 {
00188 REPORT gm->Negate();
00189 return gm;
00190 }
00191 else
00192 {
00193 REPORT
00194 GeneralMatrix* gmx = gm->Type().New(nr,nc,this);
00195 gmx->ReleaseAndDelete(); gmx->Negate(gm);
00196 return gmx;
00197 }
00198 }
00199 else
00200 {
00201 REPORT
00202 GeneralMatrix* gmx = mt.New(nr,nc,this);
00203 MatrixRow mr(gm, LoadOnEntry);
00204 MatrixRow mrx(gmx, StoreOnExit+DirectPart);
00205 while (nr--) { mrx.Negate(mr); mrx.Next(); mr.Next(); }
00206 gmx->ReleaseAndDelete(); gm->tDelete();
00207 return gmx;
00208 }
00209 }
00210
00211 GeneralMatrix* ReversedMatrix::Evaluate(MatrixType mt)
00212 {
00213 gm=(const_cast<BaseMatrix*&>(bm))->Evaluate(); GeneralMatrix* gmx;
00214
00215 if ((gm->Type()).IsBand() && ! (gm->Type()).IsDiagonal())
00216 {
00217 gm->tDelete();
00218 Throw(NotDefinedException("Reverse", "band matrices"));
00219 }
00220
00221 if (gm->reuse()) { REPORT gm->ReverseElements(); gmx = gm; }
00222 else
00223 {
00224 REPORT
00225 gmx = gm->Type().New(gm->Nrows(), gm->Ncols(), this);
00226 gmx->ReverseElements(gm); gmx->ReleaseAndDelete();
00227 }
00228 return gmx->Evaluate(mt);
00229
00230 }
00231
00232 GeneralMatrix* TransposedMatrix::Evaluate(MatrixType mt)
00233 {
00234 REPORT
00235 gm=(const_cast<BaseMatrix*&>(bm))->Evaluate();
00236 Compare(gm->Type().t(),mt);
00237 GeneralMatrix* gmx=gm->Transpose(this, mt);
00238 return gmx;
00239 }
00240
00241 GeneralMatrix* RowedMatrix::Evaluate(MatrixType mt)
00242 {
00243 gm = (const_cast<BaseMatrix*&>(bm))->Evaluate();
00244 GeneralMatrix* gmx = new RowVector; MatrixErrorNoSpace(gmx);
00245 gmx->nrows_value = 1; gmx->ncols_value = gmx->storage = gm->storage;
00246 return gm->BorrowStore(gmx,mt);
00247 }
00248
00249 GeneralMatrix* ColedMatrix::Evaluate(MatrixType mt)
00250 {
00251 gm = (const_cast<BaseMatrix*&>(bm))->Evaluate();
00252 GeneralMatrix* gmx = new ColumnVector; MatrixErrorNoSpace(gmx);
00253 gmx->ncols_value = 1; gmx->nrows_value = gmx->storage = gm->storage;
00254 return gm->BorrowStore(gmx,mt);
00255 }
00256
00257 GeneralMatrix* DiagedMatrix::Evaluate(MatrixType mt)
00258 {
00259 gm = (const_cast<BaseMatrix*&>(bm))->Evaluate();
00260 GeneralMatrix* gmx = new DiagonalMatrix; MatrixErrorNoSpace(gmx);
00261 gmx->nrows_value = gmx->ncols_value = gmx->storage = gm->storage;
00262 return gm->BorrowStore(gmx,mt);
00263 }
00264
00265 GeneralMatrix* MatedMatrix::Evaluate(MatrixType mt)
00266 {
00267 Tracer tr("MatedMatrix::Evaluate");
00268 gm = (const_cast<BaseMatrix*&>(bm))->Evaluate();
00269 GeneralMatrix* gmx = new Matrix; MatrixErrorNoSpace(gmx);
00270 gmx->nrows_value = nr; gmx->ncols_value = nc; gmx->storage = gm->storage;
00271 if (nr*nc != gmx->storage)
00272 Throw(IncompatibleDimensionsException());
00273 return gm->BorrowStore(gmx,mt);
00274 }
00275
00276 GeneralMatrix* GetSubMatrix::Evaluate(MatrixType mt)
00277 {
00278 REPORT
00279 Tracer tr("SubMatrix(evaluate)");
00280 gm = (const_cast<BaseMatrix*&>(bm))->Evaluate();
00281 if (row_number < 0) row_number = gm->Nrows();
00282 if (col_number < 0) col_number = gm->Ncols();
00283 if (row_skip+row_number > gm->Nrows() || col_skip+col_number > gm->Ncols())
00284 {
00285 gm->tDelete();
00286 Throw(SubMatrixDimensionException());
00287 }
00288 if (IsSym) Compare(gm->Type().ssub(), mt);
00289 else Compare(gm->Type().sub(), mt);
00290 GeneralMatrix* gmx = mt.New(row_number, col_number, this);
00291 int i = row_number;
00292 MatrixRow mr(gm, LoadOnEntry, row_skip);
00293 MatrixRow mrx(gmx, StoreOnExit+DirectPart);
00294 MatrixRowCol sub;
00295 while (i--)
00296 {
00297 mr.SubRowCol(sub, col_skip, col_number);
00298 mrx.Copy(sub); mrx.Next(); mr.Next();
00299 }
00300 gmx->ReleaseAndDelete(); gm->tDelete();
00301 return gmx;
00302 }
00303
00304
00305 GeneralMatrix* ReturnMatrix::Evaluate(MatrixType mt)
00306 {
00307 return gm->Evaluate(mt);
00308 }
00309
00310
00311 void GeneralMatrix::Add(GeneralMatrix* gm1, Real f)
00312 {
00313 REPORT
00314 Real* s1=gm1->store; Real* s=store; int i=(storage >> 2);
00315 while (i--)
00316 { *s++ = *s1++ + f; *s++ = *s1++ + f; *s++ = *s1++ + f; *s++ = *s1++ + f; }
00317 i = storage & 3; while (i--) *s++ = *s1++ + f;
00318 }
00319
00320 void GeneralMatrix::Add(Real f)
00321 {
00322 REPORT
00323 Real* s=store; int i=(storage >> 2);
00324 while (i--) { *s++ += f; *s++ += f; *s++ += f; *s++ += f; }
00325 i = storage & 3; while (i--) *s++ += f;
00326 }
00327
00328 void GeneralMatrix::NegAdd(GeneralMatrix* gm1, Real f)
00329 {
00330 REPORT
00331 Real* s1=gm1->store; Real* s=store; int i=(storage >> 2);
00332 while (i--)
00333 { *s++ = f - *s1++; *s++ = f - *s1++; *s++ = f - *s1++; *s++ = f - *s1++; }
00334 i = storage & 3; while (i--) *s++ = f - *s1++;
00335 }
00336
00337 void GeneralMatrix::NegAdd(Real f)
00338 {
00339 REPORT
00340 Real* s=store; int i=(storage >> 2);
00341 while (i--)
00342 {
00343 *s = f - *s; s++; *s = f - *s; s++;
00344 *s = f - *s; s++; *s = f - *s; s++;
00345 }
00346 i = storage & 3; while (i--) { *s = f - *s; s++; }
00347 }
00348
00349 void GeneralMatrix::Negate(GeneralMatrix* gm1)
00350 {
00351
00352 REPORT
00353 Real* s1=gm1->store; Real* s=store; int i=(storage >> 2);
00354 while (i--)
00355 { *s++ = -(*s1++); *s++ = -(*s1++); *s++ = -(*s1++); *s++ = -(*s1++); }
00356 i = storage & 3; while(i--) *s++ = -(*s1++);
00357 }
00358
00359 void GeneralMatrix::Negate()
00360 {
00361 REPORT
00362 Real* s=store; int i=(storage >> 2);
00363 while (i--)
00364 { *s = -(*s); s++; *s = -(*s); s++; *s = -(*s); s++; *s = -(*s); s++; }
00365 i = storage & 3; while(i--) { *s = -(*s); s++; }
00366 }
00367
00368 void GeneralMatrix::Multiply(GeneralMatrix* gm1, Real f)
00369 {
00370 REPORT
00371 Real* s1=gm1->store; Real* s=store; int i=(storage >> 2);
00372 while (i--)
00373 { *s++ = *s1++ * f; *s++ = *s1++ * f; *s++ = *s1++ * f; *s++ = *s1++ * f; }
00374 i = storage & 3; while (i--) *s++ = *s1++ * f;
00375 }
00376
00377 void GeneralMatrix::Multiply(Real f)
00378 {
00379 REPORT
00380 Real* s=store; int i=(storage >> 2);
00381 while (i--) { *s++ *= f; *s++ *= f; *s++ *= f; *s++ *= f; }
00382 i = storage & 3; while (i--) *s++ *= f;
00383 }
00384
00385
00386
00387
00388
00389
00390
00391 MatrixInput MatrixInput::operator<<(Real f)
00392 {
00393 REPORT
00394 Tracer et("MatrixInput");
00395 if (n<=0) Throw(ProgramException("List of values too long"));
00396 *r = f; int n1 = n-1; n=0;
00397 return MatrixInput(n1, r+1);
00398 }
00399
00400
00401 MatrixInput GeneralMatrix::operator<<(Real f)
00402 {
00403 REPORT
00404 Tracer et("MatrixInput");
00405 int n = Storage();
00406 if (n<=0) Throw(ProgramException("Loading data to zero length matrix"));
00407 Real* r; r = Store(); *r = f; n--;
00408 return MatrixInput(n, r+1);
00409 }
00410
00411 MatrixInput GetSubMatrix::operator<<(Real f)
00412 {
00413 REPORT
00414 Tracer et("MatrixInput (GetSubMatrix)");
00415 SetUpLHS();
00416 if (row_number != 1 || col_skip != 0 || col_number != gm->Ncols())
00417 {
00418 Throw(ProgramException("MatrixInput requires complete rows"));
00419 }
00420 MatrixRow mr(gm, DirectPart, row_skip);
00421 int n = mr.Storage();
00422 if (n<=0)
00423 {
00424 Throw(ProgramException("Loading data to zero length row"));
00425 }
00426 Real* r; r = mr.Data(); *r = f; n--;
00427 if (+(mr.cw*HaveStore))
00428 {
00429 Throw(ProgramException("Fails with this matrix type"));
00430 }
00431 return MatrixInput(n, r+1);
00432 }
00433
00434 MatrixInput::~MatrixInput()
00435 {
00436 REPORT
00437 Tracer et("MatrixInput");
00438 if (n!=0) Throw(ProgramException("A list of values was too short"));
00439 }
00440
00441 MatrixInput BandMatrix::operator<<(Real)
00442 {
00443 Tracer et("MatrixInput");
00444 bool dummy = true;
00445 if (dummy)
00446 Throw(ProgramException("Cannot use list read with a BandMatrix"));
00447 return MatrixInput(0, 0);
00448 }
00449
00450 void BandMatrix::operator<<(const Real*)
00451 { Throw(ProgramException("Cannot use array read with a BandMatrix")); }
00452
00453 void BandMatrix::operator<<(const int*)
00454 { Throw(ProgramException("Cannot use array read with a BandMatrix")); }
00455
00456 void SymmetricBandMatrix::operator<<(const Real*)
00457 { Throw(ProgramException("Cannot use array read with a BandMatrix")); }
00458
00459 void SymmetricBandMatrix::operator<<(const int*)
00460 { Throw(ProgramException("Cannot use array read with a BandMatrix")); }
00461
00462
00463
00464 void GeneralMatrix::ReverseElements(GeneralMatrix* gm)
00465 {
00466
00467 REPORT
00468 int n = Storage(); Real* rx = Store() + n; Real* x = gm->Store();
00469 while (n--) *(--rx) = *(x++);
00470 }
00471
00472 void GeneralMatrix::ReverseElements()
00473 {
00474
00475 REPORT
00476 int n = Storage(); Real* x = Store(); Real* rx = x + n;
00477 n /= 2;
00478 while (n--) { Real t = *(--rx); *rx = *x; *(x++) = t; }
00479 }
00480
00481
00482 #ifdef use_namespace
00483 }
00484 #endif
00485