BANDMAT.CPP Source File Main Page Namespace List Class Hierarchy Alphabetical List Compound List File List Namespace Members Compound Members File Members
BANDMAT.CPP

Go to the documentation of this file.
00001 //Owner: Fred
00002 //$$ bandmat.cpp                     Band matrix definitions
00003 
00004 // Copyright (C) 1991,2,3,4: R B Davies
00005 
00006 #define WANT_MATH                                 // include.h will get math fns
00007 
00008 #include "include.h"
00009 
00010 #include "newmat.h"
00011 #include "newmatrc.h"
00012 
00013 #ifdef use_namespace
00014 namespace NEWMAT {
00015 #endif
00016 
00017 #ifdef DO_REPORT
00018 #define REPORT { static ExeCounter ExeCount(__LINE__,10); ++ExeCount; }
00019 #else
00020 #define REPORT {}
00021 #endif
00022 
00023     BandMatrix::BandMatrix(const BaseMatrix& M) {
00024         REPORT                                        // CheckConversion(M);
00025             MatrixConversionCheck mcc;
00026         GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::BM);
00027         GetMatrix(gmx); CornerClear();
00028     }
00029 
00030     void BandMatrix::SetParameters(const GeneralMatrix* gmx) {
00031         MatrixBandWidth bw = gmx->BandWidth();
00032         lower = bw.lower; upper = bw.upper;
00033     }
00034 
00035     void BandMatrix::ReSize(int n, int lb, int ub) {
00036         REPORT
00037             Tracer tr("BandMatrix::ReSize");
00038         if (lb<0 || ub<0) Throw(ProgramException("Undefined bandwidth"));
00039         lower = (lb<=n) ? lb : n-1; upper = (ub<=n) ? ub : n-1;
00040         GeneralMatrix::ReSize(n,n,n*(lower+1+upper)); CornerClear();
00041     }
00042 
00043     void BandMatrix::operator=(const BaseMatrix& X) {
00044         REPORT                                        // CheckConversion(X);
00045             MatrixConversionCheck mcc;
00046         Eq(X,MatrixType::BM); CornerClear();
00047     }
00048 
00049     void BandMatrix::CornerClear() const {
00050         // set unused parts of BandMatrix to zero
00051         REPORT
00052             int i = lower; Real* s = store; int bw = lower + 1 + upper;
00053             while (i)
00054             { int j = i--; Real* sj = s; s += bw; while (j--) *sj++ = 0.0; }
00055             i = upper; s = store + storage;
00056             while (i)
00057             { int j = i--; Real* sj = s; s -= bw; while (j--) *(--sj) = 0.0; }
00058     }
00059 
00060     MatrixBandWidth MatrixBandWidth::operator+(const MatrixBandWidth& bw) const {
00061         int l = bw.lower; int u = bw.upper;
00062         l = (lower < 0 || l < 0) ? -1 : (lower > l) ? lower : l;
00063         u = (upper < 0 || u < 0) ? -1 : (upper > u) ? upper : u;
00064         return MatrixBandWidth(l,u);
00065     }
00066 
00067     MatrixBandWidth MatrixBandWidth::operator*(const MatrixBandWidth& bw) const {
00068         int l = bw.lower; int u = bw.upper;
00069         l = (lower < 0 || l < 0) ? -1 : lower+l;
00070         u = (upper < 0 || u < 0) ? -1 : upper+u;
00071         return MatrixBandWidth(l,u);
00072     }
00073 
00074     MatrixBandWidth MatrixBandWidth::minimum(const MatrixBandWidth& bw) const {
00075         int l = bw.lower; int u = bw.upper;
00076         if ((lower >= 0) && ( (l < 0) || (l > lower) )) l = lower;
00077         if ((upper >= 0) && ( (u < 0) || (u > upper) )) u = upper;
00078         return MatrixBandWidth(l,u);
00079     }
00080 
00081     UpperBandMatrix::UpperBandMatrix(const BaseMatrix& M) {
00082         REPORT                                        // CheckConversion(M);
00083             MatrixConversionCheck mcc;
00084         GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::UB);
00085         GetMatrix(gmx); CornerClear();
00086     }
00087 
00088     void UpperBandMatrix::operator=(const BaseMatrix& X) {
00089         REPORT                                        // CheckConversion(X);
00090             MatrixConversionCheck mcc;
00091         Eq(X,MatrixType::UB); CornerClear();
00092     }
00093 
00094     LowerBandMatrix::LowerBandMatrix(const BaseMatrix& M) {
00095         REPORT                                        // CheckConversion(M);
00096             MatrixConversionCheck mcc;
00097         GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::LB);
00098         GetMatrix(gmx); CornerClear();
00099     }
00100 
00101     void LowerBandMatrix::operator=(const BaseMatrix& X) {
00102         REPORT                                        // CheckConversion(X);
00103             MatrixConversionCheck mcc;
00104         Eq(X,MatrixType::LB); CornerClear();
00105     }
00106 
00107     BandLUMatrix::BandLUMatrix(const BaseMatrix& m) {
00108         REPORT
00109             Tracer tr("BandLUMatrix");
00110         GeneralMatrix* gm = ((BaseMatrix&)m).Evaluate(MatrixType::BM);
00111         m1 = ((BandMatrix*)gm)->lower; m2 = ((BandMatrix*)gm)->upper;
00112         GetMatrix(gm);
00113         if (nrows!=ncols) Throw(NotSquareException(*this));
00114         d = true; sing = false;
00115         indx = new int [nrows]; MatrixErrorNoSpace(indx);
00116         MONITOR_INT_NEW("Index (BndLUMat)",nrows,indx)
00117             storage2 = nrows * m1;
00118         store2 = new Real [storage2]; MatrixErrorNoSpace(store2);
00119         MONITOR_REAL_NEW("Make (BandLUMat)",storage2,store2)
00120             ludcmp();
00121     }
00122 
00123     BandLUMatrix::~BandLUMatrix() {
00124         MONITOR_INT_DELETE("Index (BndLUMat)",nrows,indx)
00125             MONITOR_REAL_DELETE("Delete (BndLUMt)",storage2,store2)
00126             delete [] indx; delete [] store2;
00127     }
00128 
00129     MatrixType BandLUMatrix::Type() const { return MatrixType::BC; }
00130 
00131     LogAndSign BandLUMatrix::LogDeterminant() const {
00132         if (sing) return 0.0;
00133         Real* a = store; int w = m1+1+m2; LogAndSign sum; int i = nrows;
00134         while (i--) { sum *= *a; a += w; }
00135         if (!d) sum.ChangeSign(); return sum;
00136     }
00137 
00138     GeneralMatrix* BandMatrix::MakeSolver() {
00139         REPORT
00140             GeneralMatrix* gm = new BandLUMatrix(*this);
00141         MatrixErrorNoSpace(gm); gm->ReleaseAndDelete(); return gm;
00142     }
00143 
00144     void BandLUMatrix::ludcmp() {
00145         REPORT
00146             Real* a = store2; int i = storage2;
00147         // clear store2 - so unused locations are always zero -
00148         // required by operator==
00149         while (i--) *a++ = 0.0;
00150         a = store;
00151         i = m1; int j = m2; int k; int n = nrows; int w = m1 + 1 + m2;
00152         while (i) {
00153             Real* ai = a + i;
00154             k = ++j; while (k--) *a++ = *ai++;
00155             k = i--; while (k--) *a++ = 0.0;
00156         }
00157 
00158         a = store; int l = m1;
00159         for (k=0; k<n; k++) {
00160             Real x = *a; i = k; Real* aj = a;
00161             if (l < n) l++;
00162             for (j=k+1; j<l; j++) {
00163                 aj += w; if (fabs(x) < fabs(*aj)) {
00164                     x = *aj; i = j;
00165                 }
00166             }
00167             indx[k] = i;
00168             if (x==0) { sing = true; return; }
00169             if (i!=k) {
00170                 d = !d; Real* ak = a; Real* ai = store + i * w; j = w;
00171                 while (j--) { x = *ak; *ak++ = *ai; *ai++ = x; }
00172             }
00173             aj = a + w; Real* m = store2 + m1 * k;
00174             for (j=k+1; j<l; j++) {
00175                 *m++ = x = *aj / *a; i = w; Real* ak = a;
00176                 while (--i) { Real* aj1 = aj++; *aj1 = *aj - x * *(++ak); }
00177                 *aj++ = 0.0;
00178             }
00179             a += w;
00180         }
00181     }
00182 
00183     void BandLUMatrix::lubksb(Real* B, int mini) {
00184         REPORT
00185             Tracer tr("BandLUMatrix::lubksb");
00186         if (sing) Throw(SingularException(*this));
00187         int n = nrows; int l = m1; int w = m1 + 1 + m2;
00188 
00189         for (int k=0; k<n; k++) {
00190             int i = indx[k];
00191             if (i!=k) { Real x=B[k]; B[k]=B[i]; B[i]=x; }
00192             if (l<n) l++;
00193             Real* m = store2 + k*m1; Real* b = B+k; Real* bi = b;
00194             for (i=k+1; i<l; i++)  *(++bi) -= *m++ * *b;
00195         }
00196 
00197         l = -m1;
00198         for (int i = n-1; i>=mini; i--) {
00199             Real* b = B + i; Real* bk = b; Real x = *bk;
00200             Real* a = store + w*i; Real y = *a;
00201             int k = l+m1; while (k--) x -=  *(++a) * *(++bk);
00202             *b = x / y;
00203             if (l < m2) l++;
00204         }
00205     }
00206 
00207     void BandLUMatrix::Solver(MatrixColX& mcout, const MatrixColX& mcin) {
00208         REPORT
00209             int i = mcin.skip; Real* el = mcin.data-i; Real* el1=el;
00210             while (i--) *el++ = 0.0;
00211             el += mcin.storage; i = nrows - mcin.skip - mcin.storage;
00212             while (i--) *el++ = 0.0;
00213             lubksb(el1, mcout.skip);
00214     }
00215 
00216     // Do we need check for entirely zero output?
00217 
00218     void UpperBandMatrix::Solver(MatrixColX& mcout,
00219                                  const MatrixColX& mcin) {
00220         REPORT
00221             int i = mcin.skip-mcout.skip; Real* elx = mcin.data-i;
00222         while (i-- > 0) *elx++ = 0.0;
00223         int nr = mcin.skip+mcin.storage;
00224         elx = mcin.data+mcin.storage; Real* el = elx;
00225         int j = mcout.skip+mcout.storage-nr; i = nr-mcout.skip;
00226         while (j-- > 0) *elx++ = 0.0;
00227 
00228         Real* Ael = store + (upper+1)*(i-1)+1; j = 0;
00229         while (i-- > 0) {
00230             elx = el; Real sum = 0.0; int jx = j;
00231             while (jx--) sum += *(--Ael) * *(--elx);
00232             elx--; *elx = (*elx - sum) / *(--Ael);
00233             if (j<upper) Ael -= upper - (++j); else el--;
00234         }
00235     }
00236 
00237     void LowerBandMatrix::Solver(MatrixColX& mcout,
00238                                  const MatrixColX& mcin) {
00239         REPORT
00240             int i = mcin.skip-mcout.skip; Real* elx = mcin.data-i;
00241         while (i-- > 0) *elx++ = 0.0;
00242         int nc = mcin.skip; i = nc+mcin.storage; elx = mcin.data+mcin.storage;
00243         int nr = mcout.skip+mcout.storage; int j = nr-i; i = nr-nc;
00244         while (j-- > 0) *elx++ = 0.0;
00245 
00246         Real* el = mcin.data; Real* Ael = store + (lower+1)*nc + lower; j = 0;
00247         while (i-- > 0) {
00248             elx = el; Real sum = 0.0; int jx = j;
00249             while (jx--) sum += *Ael++ * *elx++;
00250             *elx = (*elx - sum) / *Ael++;
00251             if (j<lower) Ael += lower - (++j); else el++;
00252         }
00253     }
00254 
00255     LogAndSign BandMatrix::LogDeterminant() const {
00256         REPORT
00257             BandLUMatrix C(*this); return C.LogDeterminant();
00258     }
00259 
00260     LogAndSign LowerBandMatrix::LogDeterminant() const {
00261         REPORT
00262             int i = nrows; LogAndSign sum; Real* s = store + lower; int j = lower + 1;
00263             while (i--) { sum *= *s; s += j; }
00264             ((GeneralMatrix&)*this).tDelete(); return sum;
00265     }
00266 
00267     LogAndSign UpperBandMatrix::LogDeterminant() const {
00268         REPORT
00269             int i = nrows; LogAndSign sum; Real* s = store; int j = upper + 1;
00270             while (i--) { sum *= *s; s += j; }
00271             ((GeneralMatrix&)*this).tDelete(); return sum;
00272     }
00273 
00274     GeneralMatrix* SymmetricBandMatrix::MakeSolver() {
00275         REPORT
00276             GeneralMatrix* gm = new BandLUMatrix(*this);
00277         MatrixErrorNoSpace(gm); gm->ReleaseAndDelete(); return gm;
00278     }
00279 
00280     SymmetricBandMatrix::SymmetricBandMatrix(const BaseMatrix& M) {
00281         REPORT                                        // CheckConversion(M);
00282             MatrixConversionCheck mcc;
00283         GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::SB);
00284         GetMatrix(gmx);
00285     }
00286 
00287     GeneralMatrix* SymmetricBandMatrix::Transpose(TransposedMatrix*, MatrixType mt)
00288     { REPORT  return Evaluate(mt); }
00289 
00290     LogAndSign SymmetricBandMatrix::LogDeterminant() const {
00291         REPORT
00292             BandLUMatrix C(*this); return C.LogDeterminant();
00293     }
00294 
00295     void SymmetricBandMatrix::SetParameters(const GeneralMatrix* gmx)
00296     { lower = gmx->BandWidth().lower; }
00297 
00298     void SymmetricBandMatrix::ReSize(int n, int lb) {
00299         REPORT
00300             Tracer tr("SymmetricBandMatrix::ReSize");
00301         if (lb<0) Throw(ProgramException("Undefined bandwidth"));
00302         lower = (lb<=n) ? lb : n-1;
00303         GeneralMatrix::ReSize(n,n,n*(lower+1));
00304     }
00305 
00306     void SymmetricBandMatrix::operator=(const BaseMatrix& X) {
00307         REPORT                                        // CheckConversion(X);
00308             MatrixConversionCheck mcc;
00309         Eq(X,MatrixType::SB);
00310     }
00311 
00312     void SymmetricBandMatrix::CornerClear() const {
00313         // set unused parts of BandMatrix to zero
00314         REPORT
00315             int i = lower; Real* s = store; int bw = lower + 1;
00316             while (i)
00317             { int j = i--; Real* sj = s; s += bw; while (j--) *sj++ = 0.0; }
00318     }
00319 
00320     MatrixBandWidth SymmetricBandMatrix::BandWidth() const
00321     { return MatrixBandWidth(lower,lower); }
00322 
00323     inline Real square(Real x) { return x*x; }
00324 
00325     Real SymmetricBandMatrix::SumSquare() const {
00326         REPORT
00327             CornerClear();
00328         Real sum1=0.0; Real sum2=0.0; Real* s=store; int i=nrows; int l=lower;
00329         while (i--)
00330         { int j = l; while (j--) sum2 += square(*s++); sum1 += square(*s++); }
00331         ((GeneralMatrix&)*this).tDelete(); return sum1 + 2.0 * sum2;
00332     }
00333 
00334     Real SymmetricBandMatrix::SumAbsoluteValue() const {
00335         REPORT
00336             CornerClear();
00337         Real sum1=0.0; Real sum2=0.0; Real* s=store; int i=nrows; int l=lower;
00338         while (i--)
00339         { int j = l; while (j--) sum2 += fabs(*s++); sum1 += fabs(*s++); }
00340         ((GeneralMatrix&)*this).tDelete(); return sum1 + 2.0 * sum2;
00341     }
00342 
00343     Real SymmetricBandMatrix::Sum() const {
00344         REPORT
00345             CornerClear();
00346         Real sum1=0.0; Real sum2=0.0; Real* s=store; int i=nrows; int l=lower;
00347         while (i--)
00348         { int j = l; while (j--) sum2 += *s++; sum1 += *s++; }
00349         ((GeneralMatrix&)*this).tDelete(); return sum1 + 2.0 * sum2;
00350     }
00351 
00352 #ifdef use_namespace
00353 }
00354 #endif
Generated on Fri Aug 23 01:37:01 2002 for VirtualU by
1.2.17