bandmat.cpp
Go to the documentation of this file.
00001 
00002 
00003 
00006 
00007 
00008 // Copyright (C) 1991,2,3,4,9: R B Davies
00009 
00010 #define WANT_MATH                    // include.h will get math fns
00011 
00012 //#define WANT_STREAM
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 // CheckConversion(M);
00038    // MatrixConversionCheck mcc;
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 // SimpleAddOK shows when we can add etc two matrices by a simple vector add
00060 // and when we can add one matrix into another
00061 //
00062 // *gm must be the same type as *this
00063 // - return 0 if simple add is OK
00064 // - return 1 if we can add into *gm only
00065 // - return 2 if we can add into *this only
00066 // - return 3 if we can't add either way
00067 //
00068 // For SP this will still be valid if we swap 1 and 2
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 bool BandMatrix::SameStorageType(const GeneralMatrix& A) const
00138 {
00139    if (type() != A.type()) { REPORT return false; }
00140    REPORT
00141    return bandwidth() == A.bandwidth();
00142 }
00143 
00144 void BandMatrix::resizeForAdd(const GeneralMatrix& A, const GeneralMatrix& B)
00145 {
00146    REPORT
00147    Tracer tr("BandMatrix::resizeForAdd");
00148    MatrixBandWidth A_BW = A.bandwidth(); MatrixBandWidth B_BW = B.bandwidth();
00149    if ((A_BW.Lower() < 0) | (A_BW.Upper() < 0) | (B_BW.Lower() < 0)
00150       | (A_BW.Upper() < 0))
00151          Throw(ProgramException("Can't resize to BandMatrix" ));
00152    // already know A and B are square
00153    resize(A.Nrows(), my_max(A_BW.Lower(), B_BW.Lower()),
00154       my_max(A_BW.Upper(), B_BW.Upper()));
00155 }
00156 
00157 void BandMatrix::resizeForSP(const GeneralMatrix& A, const GeneralMatrix& B)
00158 {
00159    REPORT
00160    Tracer tr("BandMatrix::resizeForSP");
00161    MatrixBandWidth A_BW = A.bandwidth(); MatrixBandWidth B_BW = B.bandwidth();
00162    if ((A_BW.Lower() < 0) | (A_BW.Upper() < 0) | (B_BW.Lower() < 0)
00163       | (A_BW.Upper() < 0))
00164          Throw(ProgramException("Can't resize to BandMatrix" ));
00165    // already know A and B are square
00166    resize(A.Nrows(), my_min(A_BW.Lower(), B_BW.Lower()),
00167       my_min(A_BW.Upper(), B_BW.Upper()));
00168 }
00169 */
00170 
00172 void BandMatrix::operator=(const BaseMatrix& X)
00173 {
00174    REPORT // CheckConversion(X);
00175    // MatrixConversionCheck mcc;
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 // CheckConversion(M);
00221    // MatrixConversionCheck mcc;
00222    GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::UB);
00223    GetMatrix(gmx); CornerClear();
00224 }
00225 
00226 void UpperBandMatrix::operator=(const BaseMatrix& X)
00227 {
00228    REPORT // CheckConversion(X);
00229    // MatrixConversionCheck mcc;
00230    Eq(X,MatrixType::UB); CornerClear();
00231 }
00232 
00233 LowerBandMatrix::LowerBandMatrix(const BaseMatrix& M)
00234 {
00235    REPORT // CheckConversion(M);
00236    // MatrixConversionCheck mcc;
00237    GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::LB);
00238    GetMatrix(gmx); CornerClear();
00239 }
00240 
00241 void LowerBandMatrix::operator=(const BaseMatrix& X)
00242 {
00243    REPORT // CheckConversion(X);
00244    // MatrixConversionCheck mcc;
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; // in event of exception during build
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 // could we use SetParameters instead of this
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) // reuse the array 
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                              // copy the array
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    // while (i--) { sum *= *a; a += w; }
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    // clear store2 - so unused locations are always zero -
00378    // required by operator==
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 // Do we need check for entirely zero output?
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 //   while (i--) { sum *= *s; s += j; }
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 //   while (i--) { sum *= *s; s += j; }
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  // CheckConversion(M);
00535    // MatrixConversionCheck mcc;
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 bool SymmetricBandMatrix::SameStorageType(const GeneralMatrix& A) const
00580 {
00581    if (type() != A.type()) { REPORT return false; }
00582    REPORT
00583    return bandwidth() == A.bandwidth();
00584 }
00585 
00586 void SymmetricBandMatrix::resizeForAdd(const GeneralMatrix& A,
00587    const GeneralMatrix& B)
00588 {
00589    REPORT
00590    Tracer tr("SymmetricBandMatrix::resizeForAdd");
00591    MatrixBandWidth A_BW = A.bandwidth(); MatrixBandWidth B_BW = B.bandwidth();
00592    if ((A_BW.Lower() < 0) | (B_BW.Lower() < 0))
00593          Throw(ProgramException("Can't resize to SymmetricBandMatrix" ));
00594    // already know A and B are square
00595    resize(A.Nrows(), my_max(A_BW.Lower(), B_BW.Lower()));
00596 }
00597 
00598 void SymmetricBandMatrix::resizeForSP(const GeneralMatrix& A,
00599    const GeneralMatrix& B)
00600 {
00601    REPORT
00602    Tracer tr("SymmetricBandMatrix::resizeForSP");
00603    MatrixBandWidth A_BW = A.bandwidth(); MatrixBandWidth B_BW = B.bandwidth();
00604    if ((A_BW.Lower() < 0) | (B_BW.Lower() < 0))
00605          Throw(ProgramException("Can't resize to SymmetricBandMatrix" ));
00606    // already know A and B are square
00607    resize(A.Nrows(), my_min(A_BW.Lower(), B_BW.Lower()));
00608 }
00609 */
00610 
00611 void SymmetricBandMatrix::operator=(const BaseMatrix& X)
00612 {
00613    REPORT // CheckConversion(X);
00614    // MatrixConversionCheck mcc;
00615    Eq(X,MatrixType::SB);
00616 }
00617 
00618 void SymmetricBandMatrix::CornerClear() const
00619 {
00620    // set unused parts of BandMatrix to zero
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 


kni
Author(s): Martin Günther
autogenerated on Mon Aug 14 2017 02:44:12