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