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