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

Go to the documentation of this file.
00001 //Owner: Fred
00002 //$$ newmat8.cpp         Advanced LU transform, scalar functions
00003 
00004 // Copyright (C) 1991,2,3,4: R B Davies
00005 
00006 #define WANT_MATH
00007 
00008 #include "include.h"
00009 
00010 #include "newmatap.h"
00011 #include "precisio.h"
00012 
00013 #ifdef use_namespace
00014 namespace NEWMAT {
00015 #endif
00016 
00017 #ifdef DO_REPORT
00018 #define REPORT { static ExeCounter ExeCount(__LINE__,8); ++ExeCount; }
00019 #else
00020 #define REPORT {}
00021 #endif
00022 
00023     /************************** LU transformation ****************************/
00024 
00025     void CroutMatrix::ludcmp() {
00026         // LU decomposition from Golub & Van Loan, algorithm 3.4.1, (the "outer
00027         // product" version).
00028         // This replaces the code derived from Numerical Recipes in C in previous
00029         // versions of newmat and being row oriented runs much faster with large
00030         // matrices.
00031         REPORT
00032             Tracer trace( "Crout(ludcmp)" ); sing = false;
00033         Real* akk = store;                            // runs down diagonal
00034 
00035         Real big = fabs(*akk); int mu = 0; Real* ai = akk; int k;
00036 
00037         for (k = 1; k < nrows; k++) {
00038             ai += nrows; const Real trybig = fabs(*ai);
00039             if (big < trybig) { big = trybig; mu = k; }
00040         }
00041 
00042         for (k = 0; k < nrows; k++) {
00043             /*
00044               int mu1;
00045               {
00046               Real big = fabs(*akk); mu1 = k; Real* ai = akk; int i;
00047 
00048               for (i = k+1; i < nrows; i++)
00049               {
00050               ai += nrows; const Real trybig = fabs(*ai);
00051               if (big < trybig) { big = trybig; mu1 = i; }
00052               }
00053               }
00054               if (mu1 != mu) cout << k << " " << mu << " " << mu1 << endl;
00055             */
00056 
00057             indx[k] = mu;
00058 
00059             if (mu != k) {                              //row swap
00060                 Real* a1 = store + nrows * k; Real* a2 = store + nrows * mu; d = !d;
00061                 int j = nrows;
00062                 while (j--) { const Real temp = *a1; *a1++ = *a2; *a2++ = temp; }
00063             }
00064 
00065             Real diag = *akk; big = 0; mu = k + 1;
00066             if (diag != 0) {
00067                 ai = akk; int i = nrows - k - 1;
00068                 while (i--) {
00069                     ai += nrows; Real* al = ai; Real mult = *al / diag; *al = mult;
00070                     int l = nrows - k - 1; Real* aj = akk;
00071                     // work out the next pivot as part of this loop
00072                     // this saves a column operation
00073                     if (l-- != 0) {
00074                         *(++al) -= (mult * *(++aj));
00075                         const Real trybig = fabs(*al);
00076                         if (big < trybig) { big = trybig; mu = nrows - i - 1; }
00077                         while (l--) *(++al) -= (mult * *(++aj));
00078                     }
00079                 }
00080             }
00081             else sing = true;
00082             akk += nrows + 1;
00083         }
00084     }
00085 
00086     void CroutMatrix::lubksb(Real* B, int mini) {
00087         REPORT
00088             // this has been adapted from Numerical Recipes in C. The code has been
00089             // substantially streamlined, so I don't think much of the original
00090             // copyright remains. However there is not much opportunity for
00091             // variation in the code, so it is still similar to the NR code.
00092             // I follow the NR code in skipping over initial zeros in the B vector.
00093 
00094             Tracer trace("Crout(lubksb)");
00095         if (sing) Throw(SingularException(*this));
00096         int i, j, ii;
00097 
00098         // scan for first non-zero in B
00099         for (i = 0; i < nrows; i++) {
00100             int ip = indx[i]; Real temp = B[ip]; B[ip] = B[i]; B[i] = temp;
00101             if (temp != 0.0) { ii = i; break; }
00102         }
00103 
00104         Real* bi = B + ii; Real* ai = store + ii + (ii + 1) * nrows;
00105 
00106         for (i = ii + 1; i < nrows; i++) {
00107             int ip = indx[i]; Real sum = B[ip]; B[ip] = B[i];
00108             Real* aij = ai; Real* bj = bi; j = i - ii;
00109             while (j--) sum -= *aij++ * *bj++;
00110             B[i] = sum; ai += nrows;
00111         }
00112 
00113         ai = store + nrows * nrows;
00114 
00115         for (i = nrows - 1; i >= mini; i--) {
00116             Real* bj = B+i; ai -= nrows; Real* ajx = ai+i;
00117             Real sum = *bj; Real diag = *ajx;
00118             j = nrows - i; while(--j) sum -= *(++ajx) * *(++bj);
00119             B[i] = sum / diag;
00120         }
00121     }
00122 
00123     /****************************** scalar functions ****************************/
00124 
00125     inline Real square(Real x) { return x*x; }
00126 
00127     Real GeneralMatrix::SumSquare() const {
00128         REPORT
00129             Real sum = 0.0; int i = storage; Real* s = store;
00130             while (i--) sum += square(*s++);
00131             ((GeneralMatrix&)*this).tDelete(); return sum;
00132     }
00133 
00134     Real GeneralMatrix::SumAbsoluteValue() const {
00135         REPORT
00136             Real sum = 0.0; int i = storage; Real* s = store;
00137             while (i--) sum += fabs(*s++);
00138             ((GeneralMatrix&)*this).tDelete(); return sum;
00139     }
00140 
00141     Real GeneralMatrix::Sum() const {
00142         REPORT
00143             Real sum = 0.0; int i = storage; Real* s = store;
00144             while (i--) sum += *s++;
00145             ((GeneralMatrix&)*this).tDelete(); return sum;
00146     }
00147 
00148     Real GeneralMatrix::MaximumAbsoluteValue() const {
00149         REPORT
00150             Real maxval = 0.0; int i = storage; Real* s = store;
00151             while (i--) { Real a = fabs(*s++); if (maxval < a) maxval = a; }
00152             ((GeneralMatrix&)*this).tDelete(); return maxval;
00153     }
00154 
00155     Real GeneralMatrix::MaximumValue() const {
00156         REPORT
00157             Real maxval = 0.0; int i = storage; Real* s = store;
00158             while (i--) { Real a = *s++; if (maxval < a) maxval = a; }
00159             ((GeneralMatrix&)*this).tDelete(); return maxval;
00160     }
00161 
00162     Real GeneralMatrix::MinimumValue() const {
00163         REPORT
00164             Real minval = 0.0; int i = storage; Real* s = store;
00165             while (i--) { Real a = *s++; if (minval > a) minval = a; }
00166             ((GeneralMatrix&)*this).tDelete(); return minval;
00167     }
00168 
00169     Real SymmetricMatrix::SumSquare() const {
00170         REPORT
00171             Real sum1 = 0.0; Real sum2 = 0.0; Real* s = store; int nr = nrows;
00172             for (int i = 0; i<nr; i++) {
00173                 int j = i;
00174                 while (j--) sum2 += square(*s++);
00175                 sum1 += square(*s++);
00176             }
00177             ((GeneralMatrix&)*this).tDelete(); return sum1 + 2.0 * sum2;
00178     }
00179 
00180     Real SymmetricMatrix::SumAbsoluteValue() const {
00181         REPORT
00182             Real sum1 = 0.0; Real sum2 = 0.0; Real* s = store; int nr = nrows;
00183             for (int i = 0; i<nr; i++) {
00184                 int j = i;
00185                 while (j--) sum2 += fabs(*s++);
00186                 sum1 += fabs(*s++);
00187             }
00188             ((GeneralMatrix&)*this).tDelete(); return sum1 + 2.0 * sum2;
00189     }
00190 
00191     Real SymmetricMatrix::Sum() const {
00192         REPORT
00193             Real sum1 = 0.0; Real sum2 = 0.0; Real* s = store; int nr = nrows;
00194             for (int i = 0; i<nr; i++) {
00195                 int j = i;
00196                 while (j--) sum2 += *s++;
00197                 sum1 += *s++;
00198             }
00199             ((GeneralMatrix&)*this).tDelete(); return sum1 + 2.0 * sum2;
00200     }
00201 
00202     Real BaseMatrix::SumSquare() const {
00203         REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
00204         Real s = gm->SumSquare(); return s;
00205     }
00206 
00207     Real BaseMatrix::SumAbsoluteValue() const {
00208         REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
00209         Real s = gm->SumAbsoluteValue(); return s;
00210     }
00211 
00212     Real BaseMatrix::Sum() const {
00213         REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
00214         Real s = gm->Sum(); return s;
00215     }
00216 
00217     Real BaseMatrix::MaximumAbsoluteValue() const {
00218         REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
00219         Real s = gm->MaximumAbsoluteValue(); return s;
00220     }
00221 
00222     Real BaseMatrix::MaximumValue() const {
00223         REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
00224         Real s = gm->MaximumValue(); return s;
00225     }
00226     Real BaseMatrix::MinimumValue() const {
00227         REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
00228         Real s = gm->MinimumValue(); return s;
00229     }
00230 
00231     Real Matrix::Trace() const {
00232         REPORT
00233             Tracer trace("Trace");
00234         int i = nrows; int d = i+1;
00235         if (i != ncols) Throw(NotSquareException(*this));
00236         Real sum = 0.0; Real* s = store;
00237         while (i--) { sum += *s; s += d; }
00238         ((GeneralMatrix&)*this).tDelete(); return sum;
00239     }
00240 
00241     Real DiagonalMatrix::Trace() const {
00242         REPORT
00243             int i = nrows; Real sum = 0.0; Real* s = store;
00244             while (i--) sum += *s++;
00245             ((GeneralMatrix&)*this).tDelete(); return sum;
00246     }
00247 
00248     Real SymmetricMatrix::Trace() const {
00249         REPORT
00250             int i = nrows; Real sum = 0.0; Real* s = store; int j = 2;
00251             while (i--) { sum += *s; s += j++; }
00252             ((GeneralMatrix&)*this).tDelete(); return sum;
00253     }
00254 
00255     Real LowerTriangularMatrix::Trace() const {
00256         REPORT
00257             int i = nrows; Real sum = 0.0; Real* s = store; int j = 2;
00258             while (i--) { sum += *s; s += j++; }
00259             ((GeneralMatrix&)*this).tDelete(); return sum;
00260     }
00261 
00262     Real UpperTriangularMatrix::Trace() const {
00263         REPORT
00264             int i = nrows; Real sum = 0.0; Real* s = store;
00265             while (i) { sum += *s; s += i--; }
00266             ((GeneralMatrix&)*this).tDelete(); return sum;
00267     }
00268 
00269     Real BandMatrix::Trace() const {
00270         REPORT
00271             int i = nrows; int w = lower+upper+1;
00272         Real sum = 0.0; Real* s = store+lower;
00273         while (i--) { sum += *s; s += w; }
00274         ((GeneralMatrix&)*this).tDelete(); return sum;
00275     }
00276 
00277     Real SymmetricBandMatrix::Trace() const {
00278         REPORT
00279             int i = nrows; int w = lower+1;
00280         Real sum = 0.0; Real* s = store+lower;
00281         while (i--) { sum += *s; s += w; }
00282         ((GeneralMatrix&)*this).tDelete(); return sum;
00283     }
00284 
00285     Real BaseMatrix::Trace() const {
00286         REPORT
00287             GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate(MatrixType::Dg);
00288         Real sum = gm->Trace(); return sum;
00289     }
00290 
00291     void LogAndSign::operator*=(Real x) {
00292         if (x > 0.0) { log_value += log(x); }
00293         else if (x < 0.0) { log_value += log(-x); sign = -sign; }
00294         else sign = 0;
00295     }
00296 
00297     Real LogAndSign::Value() const {
00298         Tracer et("LogAndSign::Value");
00299         if (log_value >= FloatingPointPrecision::LnMaximum())
00300             Throw(OverflowException("Overflow in exponential"));
00301         return sign * exp(log_value);
00302     }
00303 
00304     LogAndSign::LogAndSign(Real f) {
00305         if (f == 0.0) { log_value = 0.0; sign = 0; return; }
00306         else if (f < 0.0) { sign = -1; f = -f; }
00307         else sign = 1;
00308         log_value = log(f);
00309     }
00310 
00311     LogAndSign DiagonalMatrix::LogDeterminant() const {
00312         REPORT
00313             int i = nrows; LogAndSign sum; Real* s = store;
00314             while (i--) sum *= *s++;
00315             ((GeneralMatrix&)*this).tDelete(); return sum;
00316     }
00317 
00318     LogAndSign LowerTriangularMatrix::LogDeterminant() const {
00319         REPORT
00320             int i = nrows; LogAndSign sum; Real* s = store; int j = 2;
00321             while (i--) { sum *= *s; s += j++; }
00322             ((GeneralMatrix&)*this).tDelete(); return sum;
00323     }
00324 
00325     LogAndSign UpperTriangularMatrix::LogDeterminant() const {
00326         REPORT
00327             int i = nrows; LogAndSign sum; Real* s = store;
00328             while (i) { sum *= *s; s += i--; }
00329             ((GeneralMatrix&)*this).tDelete(); return sum;
00330     }
00331 
00332     LogAndSign BaseMatrix::LogDeterminant() const {
00333         REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
00334         LogAndSign sum = gm->LogDeterminant(); return sum;
00335     }
00336 
00337     LogAndSign GeneralMatrix::LogDeterminant() const {
00338         REPORT
00339             Tracer tr("Determinant");
00340         if (nrows != ncols) Throw(NotSquareException(*this));
00341         CroutMatrix C(*this); return C.LogDeterminant();
00342     }
00343 
00344     LogAndSign CroutMatrix::LogDeterminant() const {
00345         REPORT
00346             if (sing) return 0.0;
00347         int i = nrows; int dd = i+1; LogAndSign sum; Real* s = store;
00348         while (i--) { sum *= *s; s += dd; }
00349         if (!d) sum.ChangeSign(); return sum;
00350 
00351     }
00352 
00353     LinearEquationSolver::LinearEquationSolver(const BaseMatrix& bm)
00354         : gm( ( ((BaseMatrix&)bm).Evaluate() )->MakeSolver() ) {
00355         if (gm==&bm) { REPORT  gm = gm->Image(); }
00356         // want a copy if  *gm is actually bm
00357         else { REPORT  gm->Protect(); }
00358     }
00359 
00360 #ifdef use_namespace
00361 }
00362 #endif
Generated on Fri Aug 23 01:37:10 2002 for VirtualU by
1.2.17