00001
00002
00003
00004
00005 #ifndef WANT_MATH
00006 #define WANT_MATH // include.h will get math fns
00007 #endif
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
00022 #ifdef DO_REPORT
00023 #define REPORT { static ExeCounter ExeCount(__LINE__,10); ++ExeCount; }
00024 #else
00025 #define REPORT {}
00026 #endif
00027
00028 static inline int my_min(int x, int y) { return x < y ? x : y; }
00029 static inline int my_max(int x, int y) { return x > y ? x : y; }
00030
00031
00032 BandMatrix::BandMatrix(const BaseMatrix& M)
00033 {
00034 REPORT
00035
00036 GeneralMatrix* gmx=(const_cast<BaseMatrix&>(M)).Evaluate(MatrixType::BM);
00037 GetMatrix(gmx); CornerClear();
00038 }
00039
00040 void BandMatrix::SetParameters(const GeneralMatrix* gmx)
00041 {
00042 REPORT
00043 MatrixBandWidth bw = gmx->BandWidth();
00044 lower = bw.lower; upper = bw.upper;
00045 }
00046
00047 void BandMatrix::ReSize(int n, int lb, int ub)
00048 {
00049 REPORT
00050 Tracer tr("BandMatrix::ReSize");
00051 if (lb<0 || ub<0) Throw(ProgramException("Undefined bandwidth"));
00052 lower = (lb<=n) ? lb : n-1; upper = (ub<=n) ? ub : n-1;
00053 GeneralMatrix::ReSize(n,n,n*(lower+1+upper)); CornerClear();
00054 }
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065 short BandMatrix::SimpleAddOK(const GeneralMatrix* gm)
00066 {
00067 const BandMatrix* bm = (const BandMatrix*)gm;
00068 if (bm->lower == lower && bm->upper == upper) { REPORT return 0; }
00069 else if (bm->lower >= lower && bm->upper >= upper) { REPORT return 1; }
00070 else if (bm->lower <= lower && bm->upper <= upper) { REPORT return 2; }
00071 else { REPORT return 3; }
00072 }
00073
00074 short SymmetricBandMatrix::SimpleAddOK(const GeneralMatrix* gm)
00075 {
00076 const SymmetricBandMatrix* bm = (const SymmetricBandMatrix*)gm;
00077 if (bm->lower == lower) { REPORT return 0; }
00078 else if (bm->lower > lower) { REPORT return 1; }
00079 else { REPORT return 2; }
00080 }
00081
00082 void UpperBandMatrix::ReSize(int n, int lb, int ub)
00083 {
00084 REPORT
00085 if (lb != 0)
00086 {
00087 Tracer tr("UpperBandMatrix::ReSize");
00088 Throw(ProgramException("UpperBandMatrix with non-zero lower band" ));
00089 }
00090 BandMatrix::ReSize(n, lb, ub);
00091 }
00092
00093 void LowerBandMatrix::ReSize(int n, int lb, int ub)
00094 {
00095 REPORT
00096 if (ub != 0)
00097 {
00098 Tracer tr("LowerBandMatrix::ReSize");
00099 Throw(ProgramException("LowerBandMatrix with non-zero upper band" ));
00100 }
00101 BandMatrix::ReSize(n, lb, ub);
00102 }
00103
00104 void BandMatrix::ReSize(const GeneralMatrix& A)
00105 {
00106 REPORT
00107 int n = A.Nrows();
00108 if (n != A.Ncols())
00109 {
00110 Tracer tr("BandMatrix::ReSize(GM)");
00111 Throw(NotSquareException(*this));
00112 }
00113 MatrixBandWidth mbw = A.BandWidth();
00114 ReSize(n, mbw.Lower(), mbw.Upper());
00115 }
00116
00117 bool BandMatrix::SameStorageType(const GeneralMatrix& A) const
00118 {
00119 if (Type() != A.Type()) { REPORT return false; }
00120 REPORT
00121 return BandWidth() == A.BandWidth();
00122 }
00123
00124 void BandMatrix::ReSizeForAdd(const GeneralMatrix& A, const GeneralMatrix& B)
00125 {
00126 REPORT
00127 Tracer tr("BandMatrix::ReSizeForAdd");
00128 MatrixBandWidth A_BW = A.BandWidth(); MatrixBandWidth B_BW = B.BandWidth();
00129 if ((A_BW.Lower() < 0) | (A_BW.Upper() < 0) | (B_BW.Lower() < 0)
00130 | (A_BW.Upper() < 0))
00131 Throw(ProgramException("Can't ReSize to BandMatrix" ));
00132
00133 ReSize(A.Nrows(), my_max(A_BW.Lower(), B_BW.Lower()),
00134 my_max(A_BW.Upper(), B_BW.Upper()));
00135 }
00136
00137 void BandMatrix::ReSizeForSP(const GeneralMatrix& A, const GeneralMatrix& B)
00138 {
00139 REPORT
00140 Tracer tr("BandMatrix::ReSizeForSP");
00141 MatrixBandWidth A_BW = A.BandWidth(); MatrixBandWidth B_BW = B.BandWidth();
00142 if ((A_BW.Lower() < 0) | (A_BW.Upper() < 0) | (B_BW.Lower() < 0)
00143 | (A_BW.Upper() < 0))
00144 Throw(ProgramException("Can't ReSize to BandMatrix" ));
00145
00146 ReSize(A.Nrows(), my_min(A_BW.Lower(), B_BW.Lower()),
00147 my_min(A_BW.Upper(), B_BW.Upper()));
00148 }
00149
00150
00151 void BandMatrix::operator=(const BaseMatrix& X)
00152 {
00153 REPORT
00154
00155 Eq(X,MatrixType::BM); CornerClear();
00156 }
00157
00158 void BandMatrix::CornerClear() const
00159 {
00160
00161 REPORT
00162 int i = lower; Real* s = store; int bw = lower + 1 + upper;
00163 while (i)
00164 { int j = i--; Real* sj = s; s += bw; while (j--) *sj++ = 0.0; }
00165 i = upper; s = store + storage;
00166 while (i)
00167 { int j = i--; Real* sj = s; s -= bw; while (j--) *(--sj) = 0.0; }
00168 }
00169
00170 MatrixBandWidth MatrixBandWidth::operator+(const MatrixBandWidth& bw) const
00171 {
00172 REPORT
00173 int l = bw.lower; int u = bw.upper;
00174 l = (lower < 0 || l < 0) ? -1 : (lower > l) ? lower : l;
00175 u = (upper < 0 || u < 0) ? -1 : (upper > u) ? upper : u;
00176 return MatrixBandWidth(l,u);
00177 }
00178
00179 MatrixBandWidth MatrixBandWidth::operator*(const MatrixBandWidth& bw) const
00180 {
00181 REPORT
00182 int l = bw.lower; int u = bw.upper;
00183 l = (lower < 0 || l < 0) ? -1 : lower+l;
00184 u = (upper < 0 || u < 0) ? -1 : upper+u;
00185 return MatrixBandWidth(l,u);
00186 }
00187
00188 MatrixBandWidth MatrixBandWidth::minimum(const MatrixBandWidth& bw) const
00189 {
00190 REPORT
00191 int l = bw.lower; int u = bw.upper;
00192 if ((lower >= 0) && ( (l < 0) || (l > lower) )) l = lower;
00193 if ((upper >= 0) && ( (u < 0) || (u > upper) )) u = upper;
00194 return MatrixBandWidth(l,u);
00195 }
00196
00197 UpperBandMatrix::UpperBandMatrix(const BaseMatrix& M)
00198 {
00199 REPORT
00200
00201 GeneralMatrix* gmx=(const_cast<BaseMatrix&>(M)).Evaluate(MatrixType::UB);
00202 GetMatrix(gmx); CornerClear();
00203 }
00204
00205 void UpperBandMatrix::operator=(const BaseMatrix& X)
00206 {
00207 REPORT
00208
00209 Eq(X,MatrixType::UB); CornerClear();
00210 }
00211
00212 LowerBandMatrix::LowerBandMatrix(const BaseMatrix& M)
00213 {
00214 REPORT
00215
00216 GeneralMatrix* gmx=(const_cast<BaseMatrix&>(M)).Evaluate(MatrixType::LB);
00217 GetMatrix(gmx); CornerClear();
00218 }
00219
00220 void LowerBandMatrix::operator=(const BaseMatrix& X)
00221 {
00222 REPORT
00223
00224 Eq(X,MatrixType::LB); CornerClear();
00225 }
00226
00227 BandLUMatrix::BandLUMatrix(const BaseMatrix& m)
00228 {
00229 REPORT
00230 Tracer tr("BandLUMatrix");
00231 storage2 = 0; store2 = 0;
00232 GeneralMatrix* gm = (const_cast<BaseMatrix&>(m)).Evaluate(MatrixType::BM);
00233 m1 = ((BandMatrix*)gm)->lower; m2 = ((BandMatrix*)gm)->upper;
00234 GetMatrix(gm);
00235 if (nrows_value != ncols_value) Throw(NotSquareException(*this));
00236 d = true; sing = false;
00237 indx = new int [nrows_value]; MatrixErrorNoSpace(indx);
00238 MONITOR_INT_NEW("Index (BndLUMat)",nrows_value,indx)
00239 storage2 = nrows_value * m1;
00240 store2 = new Real [storage2]; MatrixErrorNoSpace(store2);
00241 MONITOR_REAL_NEW("Make (BandLUMat)",storage2,store2)
00242 ludcmp();
00243 }
00244
00245 BandLUMatrix::~BandLUMatrix()
00246 {
00247 REPORT
00248 MONITOR_INT_DELETE("Index (BndLUMat)",nrows_value,indx)
00249 MONITOR_REAL_DELETE("Delete (BndLUMt)",storage2,store2)
00250 delete [] indx; delete [] store2;
00251 }
00252
00253 MatrixType BandLUMatrix::Type() const { REPORT return MatrixType::BC; }
00254
00255
00256 LogAndSign BandLUMatrix::LogDeterminant() const
00257 {
00258 REPORT
00259 if (sing) return 0.0;
00260 Real* a = store; int w = m1+1+m2; LogAndSign sum; int i = nrows_value;
00261
00262 if (i) for (;;) { sum *= *a; if (!(--i)) break; a += w; }
00263 if (!d) sum.ChangeSign(); return sum;
00264 }
00265
00266 GeneralMatrix* BandMatrix::MakeSolver()
00267 {
00268 REPORT
00269 GeneralMatrix* gm = new BandLUMatrix(*this);
00270 MatrixErrorNoSpace(gm); gm->ReleaseAndDelete(); return gm;
00271 }
00272
00273
00274 void BandLUMatrix::ludcmp()
00275 {
00276 REPORT
00277 Real* a = store2; int i = storage2;
00278
00279
00280 while (i--) *a++ = 0.0;
00281 a = store;
00282 i = m1; int j = m2; int k; int n = nrows_value; int w = m1 + 1 + m2;
00283 while (i)
00284 {
00285 Real* ai = a + i;
00286 k = ++j; while (k--) *a++ = *ai++;
00287 k = i--; while (k--) *a++ = 0.0;
00288 }
00289
00290 a = store; int l = m1;
00291 for (k=0; k<n; k++)
00292 {
00293 Real x = *a; i = k; Real* aj = a;
00294 if (l < n) l++;
00295 for (j=k+1; j<l; j++)
00296 { aj += w; if (fabs(x) < fabs(*aj)) { x = *aj; i = j; } }
00297 indx[k] = i;
00298 if (x==0) { sing = true; return; }
00299 if (i!=k)
00300 {
00301 d = !d; Real* ak = a; Real* ai = store + i * w; j = w;
00302 while (j--) { x = *ak; *ak++ = *ai; *ai++ = x; }
00303 }
00304 aj = a + w; Real* m = store2 + m1 * k;
00305 for (j=k+1; j<l; j++)
00306 {
00307 *m++ = x = *aj / *a; i = w; Real* ak = a;
00308 while (--i) { Real* aj1 = aj++; *aj1 = *aj - x * *(++ak); }
00309 *aj++ = 0.0;
00310 }
00311 a += w;
00312 }
00313 }
00314
00315 void BandLUMatrix::lubksb(Real* B, int mini)
00316 {
00317 REPORT
00318 Tracer tr("BandLUMatrix::lubksb");
00319 if (sing) Throw(SingularException(*this));
00320 int n = nrows_value; int l = m1; int w = m1 + 1 + m2;
00321
00322 for (int k=0; k<n; k++)
00323 {
00324 int i = indx[k];
00325 if (i!=k) { Real x=B[k]; B[k]=B[i]; B[i]=x; }
00326 if (l<n) l++;
00327 Real* m = store2 + k*m1; Real* b = B+k; Real* bi = b;
00328 for (i=k+1; i<l; i++) *(++bi) -= *m++ * *b;
00329 }
00330
00331 l = -m1;
00332 for (int i = n-1; i>=mini; i--)
00333 {
00334 Real* b = B + i; Real* bk = b; Real x = *bk;
00335 Real* a = store + w*i; Real y = *a;
00336 int k = l+m1; while (k--) x -= *(++a) * *(++bk);
00337 *b = x / y;
00338 if (l < m2) l++;
00339 }
00340 }
00341
00342 void BandLUMatrix::Solver(MatrixColX& mcout, const MatrixColX& mcin)
00343 {
00344 REPORT
00345 int i = mcin.skip; Real* el = mcin.data-i; Real* el1=el;
00346 while (i--) *el++ = 0.0;
00347 el += mcin.storage; i = nrows_value - mcin.skip - mcin.storage;
00348 while (i--) *el++ = 0.0;
00349 lubksb(el1, mcout.skip);
00350 }
00351
00352
00353
00354
00355 void UpperBandMatrix::Solver(MatrixColX& mcout,
00356 const MatrixColX& mcin)
00357 {
00358 REPORT
00359 int i = mcin.skip-mcout.skip; Real* elx = mcin.data-i;
00360 while (i-- > 0) *elx++ = 0.0;
00361 int nr = mcin.skip+mcin.storage;
00362 elx = mcin.data+mcin.storage; Real* el = elx;
00363 int j = mcout.skip+mcout.storage-nr; i = nr-mcout.skip;
00364 while (j-- > 0) *elx++ = 0.0;
00365
00366 Real* Ael = store + (upper+1)*(i-1)+1; j = 0;
00367 if (i > 0) for(;;)
00368 {
00369 elx = el; Real sum = 0.0; int jx = j;
00370 while (jx--) sum += *(--Ael) * *(--elx);
00371 elx--; *elx = (*elx - sum) / *(--Ael);
00372 if (--i <= 0) break;
00373 if (j<upper) Ael -= upper - (++j); else el--;
00374 }
00375 }
00376
00377 void LowerBandMatrix::Solver(MatrixColX& mcout,
00378 const MatrixColX& mcin)
00379 {
00380 REPORT
00381 int i = mcin.skip-mcout.skip; Real* elx = mcin.data-i;
00382 while (i-- > 0) *elx++ = 0.0;
00383 int nc = mcin.skip; i = nc+mcin.storage; elx = mcin.data+mcin.storage;
00384 int nr = mcout.skip+mcout.storage; int j = nr-i; i = nr-nc;
00385 while (j-- > 0) *elx++ = 0.0;
00386
00387 Real* el = mcin.data; Real* Ael = store + (lower+1)*nc + lower; j = 0;
00388 if (i > 0) for(;;)
00389 {
00390 elx = el; Real sum = 0.0; int jx = j;
00391 while (jx--) sum += *Ael++ * *elx++;
00392 *elx = (*elx - sum) / *Ael++;
00393 if (--i <= 0) break;
00394 if (j<lower) Ael += lower - (++j); else el++;
00395 }
00396 }
00397
00398
00399 LogAndSign BandMatrix::LogDeterminant() const
00400 {
00401 REPORT
00402 BandLUMatrix C(*this); return C.LogDeterminant();
00403 }
00404
00405 LogAndSign LowerBandMatrix::LogDeterminant() const
00406 {
00407 REPORT
00408 int i = nrows_value; LogAndSign sum;
00409 Real* s = store + lower; int j = lower + 1;
00410
00411 if (i) for (;;) { sum *= *s; if (!(--i)) break; s += j; }
00412 const_cast<GeneralMatrix&>((const GeneralMatrix&)*this).tDelete(); return sum;
00413 }
00414
00415 LogAndSign UpperBandMatrix::LogDeterminant() const
00416 {
00417 REPORT
00418 int i = nrows_value; LogAndSign sum; Real* s = store; int j = upper + 1;
00419
00420 if (i) for (;;) { sum *= *s; if (!(--i)) break; s += j; }
00421 const_cast<GeneralMatrix&>((const GeneralMatrix&)*this).tDelete(); return sum;
00422 }
00423
00424 GeneralMatrix* SymmetricBandMatrix::MakeSolver()
00425 {
00426 REPORT
00427 GeneralMatrix* gm = new BandLUMatrix(*this);
00428 MatrixErrorNoSpace(gm); gm->ReleaseAndDelete(); return gm;
00429 }
00430
00431 SymmetricBandMatrix::SymmetricBandMatrix(const BaseMatrix& M)
00432 {
00433 REPORT
00434
00435 GeneralMatrix* gmx=(const_cast<BaseMatrix&>(M)).Evaluate(MatrixType::SB);
00436 GetMatrix(gmx);
00437 }
00438
00439 GeneralMatrix* SymmetricBandMatrix::Transpose(TransposedMatrix*, MatrixType mt)
00440 { REPORT return Evaluate(mt); }
00441
00442 LogAndSign SymmetricBandMatrix::LogDeterminant() const
00443 {
00444 REPORT
00445 BandLUMatrix C(*this); return C.LogDeterminant();
00446 }
00447
00448 void SymmetricBandMatrix::SetParameters(const GeneralMatrix* gmx)
00449 { REPORT lower = gmx->BandWidth().lower; }
00450
00451 void SymmetricBandMatrix::ReSize(int n, int lb)
00452 {
00453 REPORT
00454 Tracer tr("SymmetricBandMatrix::ReSize");
00455 if (lb<0) Throw(ProgramException("Undefined bandwidth"));
00456 lower = (lb<=n) ? lb : n-1;
00457 GeneralMatrix::ReSize(n,n,n*(lower+1));
00458 }
00459
00460 void SymmetricBandMatrix::ReSize(const GeneralMatrix& A)
00461 {
00462 REPORT
00463 int n = A.Nrows();
00464 if (n != A.Ncols())
00465 {
00466 Tracer tr("SymmetricBandMatrix::ReSize(GM)");
00467 Throw(NotSquareException(*this));
00468 }
00469 MatrixBandWidth mbw = A.BandWidth(); int b = mbw.Lower();
00470 if (b != mbw.Upper())
00471 {
00472 Tracer tr("SymmetricBandMatrix::ReSize(GM)");
00473 Throw(ProgramException("Upper and lower band-widths not equal"));
00474 }
00475 ReSize(n, b);
00476 }
00477
00478 bool SymmetricBandMatrix::SameStorageType(const GeneralMatrix& A) const
00479 {
00480 if (Type() != A.Type()) { REPORT return false; }
00481 REPORT
00482 return BandWidth() == A.BandWidth();
00483 }
00484
00485 void SymmetricBandMatrix::ReSizeForAdd(const GeneralMatrix& A,
00486 const GeneralMatrix& B)
00487 {
00488 REPORT
00489 Tracer tr("SymmetricBandMatrix::ReSizeForAdd");
00490 MatrixBandWidth A_BW = A.BandWidth(); MatrixBandWidth B_BW = B.BandWidth();
00491 if ((A_BW.Lower() < 0) | (B_BW.Lower() < 0))
00492 Throw(ProgramException("Can't ReSize to SymmetricBandMatrix" ));
00493
00494 ReSize(A.Nrows(), my_max(A_BW.Lower(), B_BW.Lower()));
00495 }
00496
00497 void SymmetricBandMatrix::ReSizeForSP(const GeneralMatrix& A,
00498 const GeneralMatrix& B)
00499 {
00500 REPORT
00501 Tracer tr("SymmetricBandMatrix::ReSizeForSP");
00502 MatrixBandWidth A_BW = A.BandWidth(); MatrixBandWidth B_BW = B.BandWidth();
00503 if ((A_BW.Lower() < 0) | (B_BW.Lower() < 0))
00504 Throw(ProgramException("Can't ReSize to SymmetricBandMatrix" ));
00505
00506 ReSize(A.Nrows(), my_min(A_BW.Lower(), B_BW.Lower()));
00507 }
00508
00509
00510 void SymmetricBandMatrix::operator=(const BaseMatrix& X)
00511 {
00512 REPORT
00513
00514 Eq(X,MatrixType::SB);
00515 }
00516
00517 void SymmetricBandMatrix::CornerClear() const
00518 {
00519
00520 REPORT
00521 int i = lower; Real* s = store; int bw = lower + 1;
00522 if (i) for(;;)
00523 {
00524 int j = i;
00525 Real* sj = s;
00526 while (j--) *sj++ = 0.0;
00527 if (!(--i)) break;
00528 s += bw;
00529 }
00530 }
00531
00532 MatrixBandWidth SymmetricBandMatrix::BandWidth() const
00533 { REPORT return MatrixBandWidth(lower,lower); }
00534
00535 inline Real square(Real x) { return x*x; }
00536
00537
00538 Real SymmetricBandMatrix::SumSquare() const
00539 {
00540 REPORT
00541 CornerClear();
00542 Real sum1=0.0; Real sum2=0.0; Real* s=store; int i=nrows_value; int l=lower;
00543 while (i--)
00544 { int j = l; while (j--) sum2 += square(*s++); sum1 += square(*s++); }
00545 const_cast<GeneralMatrix&>((const GeneralMatrix&)*this).tDelete(); return sum1 + 2.0 * sum2;
00546 }
00547
00548 Real SymmetricBandMatrix::SumAbsoluteValue() const
00549 {
00550 REPORT
00551 CornerClear();
00552 Real sum1=0.0; Real sum2=0.0; Real* s=store; int i=nrows_value; int l=lower;
00553 while (i--)
00554 { int j = l; while (j--) sum2 += fabs(*s++); sum1 += fabs(*s++); }
00555 const_cast<GeneralMatrix&>((const GeneralMatrix&)*this).tDelete(); return sum1 + 2.0 * sum2;
00556 }
00557
00558 Real SymmetricBandMatrix::Sum() const
00559 {
00560 REPORT
00561 CornerClear();
00562 Real sum1=0.0; Real sum2=0.0; Real* s=store; int i=nrows_value; int l=lower;
00563 while (i--)
00564 { int j = l; while (j--) sum2 += *s++; sum1 += *s++; }
00565 const_cast<GeneralMatrix&>((const GeneralMatrix&)*this).tDelete(); return sum1 + 2.0 * sum2;
00566 }
00567
00568
00569 #ifdef use_namespace
00570 }
00571 #endif
00572