00001
00002
00003
00006
00007
00008
00009
00010 #define WANT_MATH // include.h will get math fns
00011
00012
00013
00014 #include "include.h"
00015
00016 #include "newmat.h"
00017 #include "newmatrc.h"
00018
00019 #ifdef use_namespace
00020 namespace NEWMAT {
00021 #endif
00022
00023
00024
00025 #ifdef DO_REPORT
00026 #define REPORT { static ExeCounter ExeCount(__LINE__,10); ++ExeCount; }
00027 #else
00028 #define REPORT {}
00029 #endif
00030
00031 static inline int my_min(int x, int y) { return x < y ? x : y; }
00032 static inline int my_max(int x, int y) { return x > y ? x : y; }
00033
00034
00035 BandMatrix::BandMatrix(const BaseMatrix& M)
00036 {
00037 REPORT
00038
00039 GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::BM);
00040 GetMatrix(gmx); CornerClear();
00041 }
00042
00043 void BandMatrix::SetParameters(const GeneralMatrix* gmx)
00044 {
00045 REPORT
00046 MatrixBandWidth bw = gmx->bandwidth();
00047 lower_val = bw.lower_val; upper_val = bw.upper_val;
00048 }
00049
00050 void BandMatrix::resize(int n, int lb, int ub)
00051 {
00052 REPORT
00053 Tracer tr("BandMatrix::resize");
00054 if (lb<0 || ub<0) Throw(ProgramException("Undefined bandwidth"));
00055 lower_val = (lb<=n) ? lb : n-1; upper_val = (ub<=n) ? ub : n-1;
00056 GeneralMatrix::resize(n,n,n*(lower_val+1+upper_val)); CornerClear();
00057 }
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00073
00074 short BandMatrix::SimpleAddOK(const GeneralMatrix* gm)
00075 {
00076 const BandMatrix* bm = (const BandMatrix*)gm;
00077 if (bm->lower_val == lower_val && bm->upper_val == upper_val)
00078 { REPORT return 0; }
00079 else if (bm->lower_val >= lower_val && bm->upper_val >= upper_val)
00080 { REPORT return 1; }
00081 else if (bm->lower_val <= lower_val && bm->upper_val <= upper_val)
00082 { REPORT return 2; }
00083 else { REPORT return 3; }
00084 }
00085
00089
00090 short SymmetricBandMatrix::SimpleAddOK(const GeneralMatrix* gm)
00091 {
00092 const SymmetricBandMatrix* bm = (const SymmetricBandMatrix*)gm;
00093 if (bm->lower_val == lower_val) { REPORT return 0; }
00094 else if (bm->lower_val > lower_val) { REPORT return 1; }
00095 else { REPORT return 2; }
00096 }
00097
00099 void UpperBandMatrix::resize(int n, int lb, int ub)
00100 {
00101 REPORT
00102 if (lb != 0)
00103 {
00104 Tracer tr("UpperBandMatrix::resize");
00105 Throw(ProgramException("UpperBandMatrix with non-zero lower band" ));
00106 }
00107 BandMatrix::resize(n, lb, ub);
00108 }
00109
00111 void LowerBandMatrix::resize(int n, int lb, int ub)
00112 {
00113 REPORT
00114 if (ub != 0)
00115 {
00116 Tracer tr("LowerBandMatrix::resize");
00117 Throw(ProgramException("LowerBandMatrix with non-zero upper band" ));
00118 }
00119 BandMatrix::resize(n, lb, ub);
00120 }
00121
00123 void BandMatrix::resize(const GeneralMatrix& A)
00124 {
00125 REPORT
00126 int n = A.Nrows();
00127 if (n != A.Ncols())
00128 {
00129 Tracer tr("BandMatrix::resize(GM)");
00130 Throw(NotSquareException(*this));
00131 }
00132 MatrixBandWidth mbw = A.bandwidth();
00133 resize(n, mbw.Lower(), mbw.Upper());
00134 }
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00172 void BandMatrix::operator=(const BaseMatrix& X)
00173 {
00174 REPORT
00175
00176 Eq(X,MatrixType::BM); CornerClear();
00177 }
00178
00180 void BandMatrix::CornerClear() const
00181 {
00182 REPORT
00183 int i = lower_val; Real* s = store; int bw = lower_val + 1 + upper_val;
00184 while (i)
00185 { int j = i--; Real* sj = s; s += bw; while (j--) *sj++ = 0.0; }
00186 i = upper_val; s = store + storage;
00187 while (i)
00188 { int j = i--; Real* sj = s; s -= bw; while (j--) *(--sj) = 0.0; }
00189 }
00190
00191 MatrixBandWidth MatrixBandWidth::operator+(const MatrixBandWidth& bw) const
00192 {
00193 REPORT
00194 int l = bw.lower_val; int u = bw.upper_val;
00195 l = (lower_val < 0 || l < 0) ? -1 : (lower_val > l) ? lower_val : l;
00196 u = (upper_val < 0 || u < 0) ? -1 : (upper_val > u) ? upper_val : u;
00197 return MatrixBandWidth(l,u);
00198 }
00199
00200 MatrixBandWidth MatrixBandWidth::operator*(const MatrixBandWidth& bw) const
00201 {
00202 REPORT
00203 int l = bw.lower_val; int u = bw.upper_val;
00204 l = (lower_val < 0 || l < 0) ? -1 : lower_val+l;
00205 u = (upper_val < 0 || u < 0) ? -1 : upper_val+u;
00206 return MatrixBandWidth(l,u);
00207 }
00208
00209 MatrixBandWidth MatrixBandWidth::minimum(const MatrixBandWidth& bw) const
00210 {
00211 REPORT
00212 int l = bw.lower_val; int u = bw.upper_val;
00213 if ((lower_val >= 0) && ( (l < 0) || (l > lower_val) )) l = lower_val;
00214 if ((upper_val >= 0) && ( (u < 0) || (u > upper_val) )) u = upper_val;
00215 return MatrixBandWidth(l,u);
00216 }
00217
00218 UpperBandMatrix::UpperBandMatrix(const BaseMatrix& M)
00219 {
00220 REPORT
00221
00222 GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::UB);
00223 GetMatrix(gmx); CornerClear();
00224 }
00225
00226 void UpperBandMatrix::operator=(const BaseMatrix& X)
00227 {
00228 REPORT
00229
00230 Eq(X,MatrixType::UB); CornerClear();
00231 }
00232
00233 LowerBandMatrix::LowerBandMatrix(const BaseMatrix& M)
00234 {
00235 REPORT
00236
00237 GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::LB);
00238 GetMatrix(gmx); CornerClear();
00239 }
00240
00241 void LowerBandMatrix::operator=(const BaseMatrix& X)
00242 {
00243 REPORT
00244
00245 Eq(X,MatrixType::LB); CornerClear();
00246 }
00247
00248 BandLUMatrix::BandLUMatrix(const BaseMatrix& m)
00249 {
00250 REPORT
00251 Tracer tr("BandLUMatrix");
00252 storage2 = 0; store2 = 0; indx = 0;
00253 GeneralMatrix* gm = ((BaseMatrix&)m).Evaluate();
00254 if (gm->nrows() != gm->ncols())
00255 { gm->tDelete(); Throw(NotSquareException(*this)); }
00256 if (gm->type() == MatrixType::BC)
00257 { REPORT ((BandLUMatrix*)gm)->get_aux(*this); GetMatrix(gm); }
00258 else
00259 {
00260 REPORT
00261 BandMatrix* gm1 = (BandMatrix*)(gm->Evaluate(MatrixType::BM));
00262 m1 = gm1->lower_val; m2 = gm1->upper_val;
00263 GetMatrix(gm1);
00264 d = true; sing = false;
00265 indx = new int [nrows_val]; MatrixErrorNoSpace(indx);
00266 MONITOR_INT_NEW("Index (BndLUMat)",nrows_val,indx)
00267 storage2 = nrows_val * m1;
00268 store2 = new Real [storage2]; MatrixErrorNoSpace(store2);
00269 MONITOR_REAL_NEW("Make (BandLUMat)",storage2,store2)
00270 ludcmp();
00271 }
00272 }
00273
00274 GeneralMatrix* BandLUMatrix::Evaluate(MatrixType mt)
00275 {
00276 if (Compare(this->Type(),mt)) { REPORT return this; }
00277 REPORT
00278 Tracer et("BandLUMatrix::Evaluate");
00279 bool dummy = true;
00280 if (dummy) Throw(ProgramException("Illegal use of BandLUMatrix", *this));
00281 return this;
00282 }
00283
00284
00285 void BandLUMatrix::get_aux(BandLUMatrix& X)
00286 {
00287 X.d = d; X.sing = sing; X.storage2 = storage2; X.m1 = m1; X.m2 = m2;
00288 if (tag_val == 0 || tag_val == 1)
00289 {
00290 REPORT
00291 X.indx = indx; indx = 0;
00292 X.store2 = store2; store2 = 0;
00293 d = true; sing = true; storage2 = 0; m1 = 0; m2 = 0;
00294 return;
00295 }
00296 else if (nrows_val == 0)
00297 {
00298 REPORT
00299 indx = 0; store2 = 0; storage2 = 0;
00300 d = true; sing = true; m1 = m2 = 0;
00301 return;
00302 }
00303 else
00304 {
00305 REPORT
00306 Tracer tr("BandLUMatrix::get_aux");
00307 int *ix = new int [nrows_val]; MatrixErrorNoSpace(ix);
00308 MONITOR_INT_NEW("Index (BLUM::get_aux)", nrows_val, ix)
00309 int n = nrows_val; int* i = ix; int* j = indx;
00310 while(n--) *i++ = *j++;
00311 X.indx = ix;
00312 Real *rx = new Real [storage2]; MatrixErrorNoSpace(indx);
00313 MONITOR_REAL_NEW("Index (BLUM::get_aux)", storage2, rx)
00314 newmat_block_copy(storage2, store2, rx);
00315 X.store2 = rx;
00316 }
00317 }
00318
00319 BandLUMatrix::BandLUMatrix(const BandLUMatrix& gm) : GeneralMatrix()
00320 {
00321 REPORT
00322 Tracer tr("BandLUMatrix(const BandLUMatrix&)");
00323 ((BandLUMatrix&)gm).get_aux(*this);
00324 GetMatrix(&gm);
00325 }
00326
00327 void BandLUMatrix::operator=(const BandLUMatrix& gm)
00328 {
00329 if (&gm == this) { REPORT tag_val = -1; return; }
00330 REPORT
00331 delete [] indx; indx = 0;
00332 delete [] store2; store2 = 0; storage2 = 0;
00333 ((BandLUMatrix&)gm).get_aux(*this);
00334 Eq(gm);
00335 }
00336
00337
00338
00339
00340
00341
00342
00343
00344 BandLUMatrix::~BandLUMatrix()
00345 {
00346 REPORT
00347 MONITOR_INT_DELETE("Index (BndLUMat)",nrows_val,indx)
00348 MONITOR_REAL_DELETE("Delete (BndLUMt)",storage2,store2)
00349 delete [] indx; delete [] store2;
00350 }
00351
00352 MatrixType BandLUMatrix::type() const { REPORT return MatrixType::BC; }
00353
00354
00355 LogAndSign BandLUMatrix::log_determinant() const
00356 {
00357 REPORT
00358 if (sing) return 0.0;
00359 Real* a = store; int w = m1+1+m2; LogAndSign sum; int i = nrows_val;
00360
00361 if (i) for (;;) { sum *= *a; if (!(--i)) break; a += w; }
00362 if (!d) sum.ChangeSign(); return sum;
00363 }
00364
00365 GeneralMatrix* BandMatrix::MakeSolver()
00366 {
00367 REPORT
00368 GeneralMatrix* gm = new BandLUMatrix(*this);
00369 MatrixErrorNoSpace(gm); gm->ReleaseAndDelete(); return gm;
00370 }
00371
00372
00373 void BandLUMatrix::ludcmp()
00374 {
00375 REPORT
00376 Real* a = store2; int i = storage2;
00377
00378
00379 while (i--) *a++ = 0.0;
00380 a = store;
00381 i = m1; int j = m2; int k; int n = nrows_val; int w = m1 + 1 + m2;
00382 while (i)
00383 {
00384 Real* ai = a + i;
00385 k = ++j; while (k--) *a++ = *ai++;
00386 k = i--; while (k--) *a++ = 0.0;
00387 }
00388
00389 a = store; int l = m1;
00390 for (k=0; k<n; k++)
00391 {
00392 Real x = *a; i = k; Real* aj = a;
00393 if (l < n) l++;
00394 for (j=k+1; j<l; j++)
00395 { aj += w; if (fabs(x) < fabs(*aj)) { x = *aj; i = j; } }
00396 indx[k] = i;
00397 if (x==0) { sing = true; return; }
00398 if (i!=k)
00399 {
00400 d = !d; Real* ak = a; Real* ai = store + i * w; j = w;
00401 while (j--) { x = *ak; *ak++ = *ai; *ai++ = x; }
00402 }
00403 aj = a + w; Real* m = store2 + m1 * k;
00404 for (j=k+1; j<l; j++)
00405 {
00406 *m++ = x = *aj / *a; i = w; Real* ak = a;
00407 while (--i) { Real* aj1 = aj++; *aj1 = *aj - x * *(++ak); }
00408 *aj++ = 0.0;
00409 }
00410 a += w;
00411 }
00412 }
00413
00414 void BandLUMatrix::lubksb(Real* B, int mini)
00415 {
00416 REPORT
00417 Tracer tr("BandLUMatrix::lubksb");
00418 if (sing) Throw(SingularException(*this));
00419 int n = nrows_val; int l = m1; int w = m1 + 1 + m2;
00420
00421 for (int k=0; k<n; k++)
00422 {
00423 int i = indx[k];
00424 if (i!=k) { Real x=B[k]; B[k]=B[i]; B[i]=x; }
00425 if (l<n) l++;
00426 Real* m = store2 + k*m1; Real* b = B+k; Real* bi = b;
00427 for (i=k+1; i<l; i++) *(++bi) -= *m++ * *b;
00428 }
00429
00430 l = -m1;
00431 for (int i = n-1; i>=mini; i--)
00432 {
00433 Real* b = B + i; Real* bk = b; Real x = *bk;
00434 Real* a = store + w*i; Real y = *a;
00435 int k = l+m1; while (k--) x -= *(++a) * *(++bk);
00436 *b = x / y;
00437 if (l < m2) l++;
00438 }
00439 }
00440
00441 void BandLUMatrix::Solver(MatrixColX& mcout, const MatrixColX& mcin)
00442 {
00443 REPORT
00444 int i = mcin.skip; Real* el = mcin.data-i; Real* el1=el;
00445 while (i--) *el++ = 0.0;
00446 el += mcin.storage; i = nrows_val - mcin.skip - mcin.storage;
00447 while (i--) *el++ = 0.0;
00448 lubksb(el1, mcout.skip);
00449 }
00450
00451
00452
00453
00454 void UpperBandMatrix::Solver(MatrixColX& mcout,
00455 const MatrixColX& mcin)
00456 {
00457 REPORT
00458 int i = mcin.skip-mcout.skip; Real* elx = mcin.data-i;
00459 while (i-- > 0) *elx++ = 0.0;
00460 int nr = mcin.skip+mcin.storage;
00461 elx = mcin.data+mcin.storage; Real* el = elx;
00462 int j = mcout.skip+mcout.storage-nr; i = nr-mcout.skip;
00463 while (j-- > 0) *elx++ = 0.0;
00464
00465 Real* Ael = store + (upper_val+1)*(i-1)+1; j = 0;
00466 if (i > 0) for(;;)
00467 {
00468 elx = el; Real sum = 0.0; int jx = j;
00469 while (jx--) sum += *(--Ael) * *(--elx);
00470 elx--; *elx = (*elx - sum) / *(--Ael);
00471 if (--i <= 0) break;
00472 if (j<upper_val) Ael -= upper_val - (++j); else el--;
00473 }
00474 }
00475
00476 void LowerBandMatrix::Solver(MatrixColX& mcout,
00477 const MatrixColX& mcin)
00478 {
00479 REPORT
00480 int i = mcin.skip-mcout.skip; Real* elx = mcin.data-i;
00481 while (i-- > 0) *elx++ = 0.0;
00482 int nc = mcin.skip; i = nc+mcin.storage; elx = mcin.data+mcin.storage;
00483 int nr = mcout.skip+mcout.storage; int j = nr-i; i = nr-nc;
00484 while (j-- > 0) *elx++ = 0.0;
00485
00486 Real* el = mcin.data;
00487 Real* Ael = store + (lower_val+1)*nc + lower_val;
00488 j = 0;
00489 if (i > 0) for(;;)
00490 {
00491 elx = el; Real sum = 0.0; int jx = j;
00492 while (jx--) sum += *Ael++ * *elx++;
00493 *elx = (*elx - sum) / *Ael++;
00494 if (--i <= 0) break;
00495 if (j<lower_val) Ael += lower_val - (++j); else el++;
00496 }
00497 }
00498
00499
00500 LogAndSign BandMatrix::log_determinant() const
00501 {
00502 REPORT
00503 BandLUMatrix C(*this); return C.log_determinant();
00504 }
00505
00506 LogAndSign LowerBandMatrix::log_determinant() const
00507 {
00508 REPORT
00509 int i = nrows_val; LogAndSign sum;
00510 Real* s = store + lower_val; int j = lower_val + 1;
00511
00512 if (i) for (;;) { sum *= *s; if (!(--i)) break; s += j; }
00513 ((GeneralMatrix&)*this).tDelete(); return sum;
00514 }
00515
00516 LogAndSign UpperBandMatrix::log_determinant() const
00517 {
00518 REPORT
00519 int i = nrows_val; LogAndSign sum; Real* s = store; int j = upper_val + 1;
00520
00521 if (i) for (;;) { sum *= *s; if (!(--i)) break; s += j; }
00522 ((GeneralMatrix&)*this).tDelete(); return sum;
00523 }
00524
00525 GeneralMatrix* SymmetricBandMatrix::MakeSolver()
00526 {
00527 REPORT
00528 GeneralMatrix* gm = new BandLUMatrix(*this);
00529 MatrixErrorNoSpace(gm); gm->ReleaseAndDelete(); return gm;
00530 }
00531
00532 SymmetricBandMatrix::SymmetricBandMatrix(const BaseMatrix& M)
00533 {
00534 REPORT
00535
00536 GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::SB);
00537 GetMatrix(gmx);
00538 }
00539
00540 GeneralMatrix* SymmetricBandMatrix::Transpose(TransposedMatrix*, MatrixType mt)
00541 { REPORT return Evaluate(mt); }
00542
00543 LogAndSign SymmetricBandMatrix::log_determinant() const
00544 {
00545 REPORT
00546 BandLUMatrix C(*this); return C.log_determinant();
00547 }
00548
00549 void SymmetricBandMatrix::SetParameters(const GeneralMatrix* gmx)
00550 { REPORT lower_val = gmx->bandwidth().lower_val; }
00551
00552 void SymmetricBandMatrix::resize(int n, int lb)
00553 {
00554 REPORT
00555 Tracer tr("SymmetricBandMatrix::resize");
00556 if (lb<0) Throw(ProgramException("Undefined bandwidth"));
00557 lower_val = (lb<=n) ? lb : n-1;
00558 GeneralMatrix::resize(n,n,n*(lower_val+1));
00559 }
00560
00561 void SymmetricBandMatrix::resize(const GeneralMatrix& A)
00562 {
00563 REPORT
00564 int n = A.Nrows();
00565 if (n != A.Ncols())
00566 {
00567 Tracer tr("SymmetricBandMatrix::resize(GM)");
00568 Throw(NotSquareException(*this));
00569 }
00570 MatrixBandWidth mbw = A.bandwidth(); int b = mbw.Lower();
00571 if (b != mbw.Upper())
00572 {
00573 Tracer tr("SymmetricBandMatrix::resize(GM)");
00574 Throw(ProgramException("Upper and lower band-widths not equal"));
00575 }
00576 resize(n, b);
00577 }
00578
00579
00580
00581
00582
00583
00584
00585
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611 void SymmetricBandMatrix::operator=(const BaseMatrix& X)
00612 {
00613 REPORT
00614
00615 Eq(X,MatrixType::SB);
00616 }
00617
00618 void SymmetricBandMatrix::CornerClear() const
00619 {
00620
00621 REPORT
00622 int i = lower_val; Real* s = store; int bw = lower_val + 1;
00623 if (i) for(;;)
00624 {
00625 int j = i;
00626 Real* sj = s;
00627 while (j--) *sj++ = 0.0;
00628 if (!(--i)) break;
00629 s += bw;
00630 }
00631 }
00632
00633 MatrixBandWidth SymmetricBandMatrix::bandwidth() const
00634 { REPORT return MatrixBandWidth(lower_val,lower_val); }
00635
00636 GeneralMatrix* BandMatrix::Image() const
00637 {
00638 REPORT
00639 GeneralMatrix* gm = new BandMatrix(*this); MatrixErrorNoSpace(gm);
00640 return gm;
00641 }
00642
00643 GeneralMatrix* UpperBandMatrix::Image() const
00644 {
00645 REPORT
00646 GeneralMatrix* gm = new UpperBandMatrix(*this); MatrixErrorNoSpace(gm);
00647 return gm;
00648 }
00649
00650 GeneralMatrix* LowerBandMatrix::Image() const
00651 {
00652 REPORT
00653 GeneralMatrix* gm = new LowerBandMatrix(*this); MatrixErrorNoSpace(gm);
00654 return gm;
00655 }
00656
00657 GeneralMatrix* SymmetricBandMatrix::Image() const
00658 {
00659 REPORT
00660 GeneralMatrix* gm = new SymmetricBandMatrix(*this); MatrixErrorNoSpace(gm);
00661 return gm;
00662 }
00663
00664 GeneralMatrix* BandLUMatrix::Image() const
00665 {
00666 REPORT
00667 GeneralMatrix* gm = new BandLUMatrix(*this); MatrixErrorNoSpace(gm);
00668 return gm;
00669 }
00670
00671
00672 inline Real square(Real x) { return x*x; }
00673
00674 Real SymmetricBandMatrix::sum_square() const
00675 {
00676 REPORT
00677 CornerClear();
00678 Real sum1=0.0; Real sum2=0.0; Real* s=store; int i=nrows_val;
00679 int l=lower_val;
00680 while (i--)
00681 { int j = l; while (j--) sum2 += square(*s++); sum1 += square(*s++); }
00682 ((GeneralMatrix&)*this).tDelete(); return sum1 + 2.0 * sum2;
00683 }
00684
00685 Real SymmetricBandMatrix::sum_absolute_value() const
00686 {
00687 REPORT
00688 CornerClear();
00689 Real sum1=0.0; Real sum2=0.0; Real* s=store; int i=nrows_val;
00690 int l=lower_val;
00691 while (i--)
00692 { int j = l; while (j--) sum2 += fabs(*s++); sum1 += fabs(*s++); }
00693 ((GeneralMatrix&)*this).tDelete(); return sum1 + 2.0 * sum2;
00694 }
00695
00696 Real SymmetricBandMatrix::sum() const
00697 {
00698 REPORT
00699 CornerClear();
00700 Real sum1=0.0; Real sum2=0.0; Real* s=store; int i=nrows_val;
00701 int l=lower_val;
00702 while (i--)
00703 { int j = l; while (j--) sum2 += *s++; sum1 += *s++; }
00704 ((GeneralMatrix&)*this).tDelete(); return sum1 + 2.0 * sum2;
00705 }
00706
00707
00708
00709
00710
00711 #ifdef use_namespace
00712 }
00713 #endif
00714
00716
00717