Tekkotsu Homepage
Demos
Overview
Downloads
Dev. Resources
Reference
Credits

newmat8.cpp

Go to the documentation of this file.
00001 //$$ newmat8.cpp         Advanced LU transform, scalar functions
00002 
00003 // Copyright (C) 1991,2,3,4,8: R B Davies
00004 
00005 #ifndef WANT_MATH
00006 #define WANT_MATH
00007 #endif
00008 
00009 #include "include.h"
00010 
00011 #include "newmat.h"
00012 #include "newmatrc.h"
00013 #include "precisio.h"
00014 
00015 #ifdef use_namespace
00016 namespace NEWMAT {
00017 #endif
00018 
00019 
00020 #ifdef DO_REPORT
00021 #define REPORT { static ExeCounter ExeCount(__LINE__,8); ++ExeCount; }
00022 #else
00023 #define REPORT {}
00024 #endif
00025 
00026 
00027 /************************** LU transformation ****************************/
00028 
00029 void CroutMatrix::ludcmp()
00030 // LU decomposition from Golub & Van Loan, algorithm 3.4.1, (the "outer
00031 // product" version).
00032 // This replaces the code derived from Numerical Recipes in C in previous
00033 // versions of newmat and being row oriented runs much faster with large
00034 // matrices.
00035 {
00036    REPORT
00037    Tracer trace( "Crout(ludcmp)" ); sing = false;
00038    Real* akk = store;                    // runs down diagonal
00039 
00040    Real big = fabs(*akk); int mu = 0; Real* ai = akk; int k;
00041 
00042    for (k = 1; k < nrows_value; k++)
00043    {
00044       ai += nrows_value; const Real trybig = fabs(*ai);
00045       if (big < trybig) { big = trybig; mu = k; }
00046    }
00047 
00048 
00049    if (nrows_value) for (k = 0;;)
00050    {
00051       /*
00052       int mu1;
00053       {
00054          Real big = fabs(*akk); mu1 = k; Real* ai = akk; int i;
00055 
00056          for (i = k+1; i < nrows_value; i++)
00057          {
00058             ai += nrows_value; const Real trybig = fabs(*ai);
00059             if (big < trybig) { big = trybig; mu1 = i; }
00060          }
00061       }
00062       if (mu1 != mu) cout << k << " " << mu << " " << mu1 << endl;
00063       */
00064 
00065       indx[k] = mu;
00066 
00067       if (mu != k)                       //row swap
00068       {
00069          Real* a1 = store + nrows_value * k;
00070          Real* a2 = store + nrows_value * mu; d = !d;
00071          int j = nrows_value;
00072          while (j--) { const Real temp = *a1; *a1++ = *a2; *a2++ = temp; }
00073       }
00074 
00075       Real diag = *akk; big = 0; mu = k + 1;
00076       if (diag != 0)
00077       {
00078          ai = akk; int i = nrows_value - k - 1;
00079          while (i--)
00080          {
00081             ai += nrows_value; Real* al = ai;
00082             Real mult = *al / diag; *al = mult;
00083             int l = nrows_value - k - 1; Real* aj = akk;
00084             // work out the next pivot as part of this loop
00085             // this saves a column operation
00086             if (l-- != 0)
00087             {
00088                *(++al) -= (mult * *(++aj));
00089                const Real trybig = fabs(*al);
00090                if (big < trybig) { big = trybig; mu = nrows_value - i - 1; }
00091                while (l--) *(++al) -= (mult * *(++aj));
00092             }
00093          }
00094       }
00095       else sing = true;
00096       if (++k == nrows_value) break;          // so next line won't overflow
00097       akk += nrows_value + 1;
00098    }
00099 }
00100 
00101 void CroutMatrix::lubksb(Real* B, int mini)
00102 {
00103    REPORT
00104    // this has been adapted from Numerical Recipes in C. The code has been
00105    // substantially streamlined, so I do not think much of the original
00106    // copyright remains. However there is not much opportunity for
00107    // variation in the code, so it is still similar to the NR code.
00108    // I follow the NR code in skipping over initial zeros in the B vector.
00109 
00110    Tracer trace("Crout(lubksb)");
00111    if (sing) Throw(SingularException(*this));
00112    int i, j, ii = nrows_value;       // ii initialised : B might be all zeros
00113 
00114 
00115    // scan for first non-zero in B
00116    for (i = 0; i < nrows_value; i++)
00117    {
00118       int ip = indx[i]; Real temp = B[ip]; B[ip] = B[i]; B[i] = temp;
00119       if (temp != 0.0) { ii = i; break; }
00120    }
00121 
00122    Real* bi; Real* ai;
00123    i = ii + 1;
00124 
00125    if (i < nrows_value)
00126    {
00127       bi = B + ii; ai = store + ii + i * nrows_value;
00128       for (;;)
00129       {
00130          int ip = indx[i]; Real sum = B[ip]; B[ip] = B[i];
00131          Real* aij = ai; Real* bj = bi; j = i - ii;
00132          while (j--) sum -= *aij++ * *bj++;
00133          B[i] = sum;
00134          if (++i == nrows_value) break;
00135          ai += nrows_value;
00136       }
00137    }
00138 
00139    ai = store + nrows_value * nrows_value;
00140 
00141    for (i = nrows_value - 1; i >= mini; i--)
00142    {
00143       Real* bj = B+i; ai -= nrows_value; Real* ajx = ai+i;
00144       Real sum = *bj; Real diag = *ajx;
00145       j = nrows_value - i; while(--j) sum -= *(++ajx) * *(++bj);
00146       B[i] = sum / diag;
00147    }
00148 }
00149 
00150 /****************************** scalar functions ****************************/
00151 
00152 inline Real square(Real x) { return x*x; }
00153 
00154 Real GeneralMatrix::SumSquare() const
00155 {
00156    REPORT
00157    Real sum = 0.0; int i = storage; Real* s = store;
00158    while (i--) sum += square(*s++);
00159    (const_cast<GeneralMatrix&>(*this)).tDelete(); return sum;
00160 }
00161 
00162 Real GeneralMatrix::SumAbsoluteValue() const
00163 {
00164    REPORT
00165    Real sum = 0.0; int i = storage; Real* s = store;
00166    while (i--) sum += fabs(*s++);
00167    (const_cast<GeneralMatrix&>(*this)).tDelete(); return sum;
00168 }
00169 
00170 Real GeneralMatrix::Sum() const
00171 {
00172    REPORT
00173    Real sum = 0.0; int i = storage; Real* s = store;
00174    while (i--) sum += *s++;
00175    (const_cast<GeneralMatrix&>(*this)).tDelete(); return sum;
00176 }
00177 
00178 // maxima and minima
00179 
00180 // There are three sets of routines
00181 // MaximumAbsoluteValue, MinimumAbsoluteValue, Maximum, Minimum
00182 // ... these find just the maxima and minima
00183 // MaximumAbsoluteValue1, MinimumAbsoluteValue1, Maximum1, Minimum1
00184 // ... these find the maxima and minima and their locations in a
00185 //     one dimensional object
00186 // MaximumAbsoluteValue2, MinimumAbsoluteValue2, Maximum2, Minimum2
00187 // ... these find the maxima and minima and their locations in a
00188 //     two dimensional object
00189 
00190 // If the matrix has no values throw an exception
00191 
00192 // If we do not want the location find the maximum or minimum on the
00193 // array stored by GeneralMatrix
00194 // This won't work for BandMatrices. We call ClearCorner for
00195 // MaximumAbsoluteValue but for the others use the AbsoluteMinimumValue2
00196 // version and discard the location.
00197 
00198 // For one dimensional objects, when we want the location of the
00199 // maximum or minimum, work with the array stored by GeneralMatrix
00200 
00201 // For two dimensional objects where we want the location of the maximum or
00202 // minimum proceed as follows:
00203 
00204 // For rectangular matrices use the array stored by GeneralMatrix and
00205 // deduce the location from the location in the GeneralMatrix
00206 
00207 // For other two dimensional matrices use the Matrix Row routine to find the
00208 // maximum or minimum for each row.
00209 
00210 static void NullMatrixError(const GeneralMatrix* gm)
00211 {
00212    (const_cast<GeneralMatrix&>(*gm)).tDelete();
00213    Throw(ProgramException("Maximum or minimum of null matrix"));
00214 }
00215 
00216 Real GeneralMatrix::MaximumAbsoluteValue() const
00217 {
00218    REPORT
00219    if (storage == 0) NullMatrixError(this);
00220    Real maxval = 0.0; int l = storage; Real* s = store;
00221    while (l--) { Real a = fabs(*s++); if (maxval < a) maxval = a; }
00222    (const_cast<GeneralMatrix&>(*this)).tDelete(); return maxval;
00223 }
00224 
00225 Real GeneralMatrix::MaximumAbsoluteValue1(int& i) const
00226 {
00227    REPORT
00228    if (storage == 0) NullMatrixError(this);
00229    Real maxval = 0.0; int l = storage; Real* s = store; int li = storage;
00230    while (l--)
00231       { Real a = fabs(*s++); if (maxval <= a) { maxval = a; li = l; }  }
00232    i = storage - li;
00233    (const_cast<GeneralMatrix&>(*this)).tDelete(); return maxval;
00234 }
00235 
00236 Real GeneralMatrix::MinimumAbsoluteValue() const
00237 {
00238    REPORT
00239    if (storage == 0) NullMatrixError(this);
00240    int l = storage - 1; Real* s = store; Real minval = fabs(*s++);
00241    while (l--) { Real a = fabs(*s++); if (minval > a) minval = a; }
00242    (const_cast<GeneralMatrix&>(*this)).tDelete(); return minval;
00243 }
00244 
00245 Real GeneralMatrix::MinimumAbsoluteValue1(int& i) const
00246 {
00247    REPORT
00248    if (storage == 0) NullMatrixError(this);
00249    int l = storage - 1; Real* s = store; Real minval = fabs(*s++); int li = l;
00250    while (l--)
00251       { Real a = fabs(*s++); if (minval >= a) { minval = a; li = l; }  }
00252    i = storage - li;
00253    (const_cast<GeneralMatrix&>(*this)).tDelete(); return minval;
00254 }
00255 
00256 Real GeneralMatrix::Maximum() const
00257 {
00258    REPORT
00259    if (storage == 0) NullMatrixError(this);
00260    int l = storage - 1; Real* s = store; Real maxval = *s++;
00261    while (l--) { Real a = *s++; if (maxval < a) maxval = a; }
00262    (const_cast<GeneralMatrix&>(*this)).tDelete(); return maxval;
00263 }
00264 
00265 Real GeneralMatrix::Maximum1(int& i) const
00266 {
00267    REPORT
00268    if (storage == 0) NullMatrixError(this);
00269    int l = storage - 1; Real* s = store; Real maxval = *s++; int li = l;
00270    while (l--) { Real a = *s++; if (maxval <= a) { maxval = a; li = l; } }
00271    i = storage - li;
00272    (const_cast<GeneralMatrix&>(*this)).tDelete(); return maxval;
00273 }
00274 
00275 Real GeneralMatrix::Minimum() const
00276 {
00277    REPORT
00278    if (storage == 0) NullMatrixError(this);
00279    int l = storage - 1; Real* s = store; Real minval = *s++;
00280    while (l--) { Real a = *s++; if (minval > a) minval = a; }
00281    (const_cast<GeneralMatrix&>(*this)).tDelete(); return minval;
00282 }
00283 
00284 Real GeneralMatrix::Minimum1(int& i) const
00285 {
00286    REPORT
00287    if (storage == 0) NullMatrixError(this);
00288    int l = storage - 1; Real* s = store; Real minval = *s++; int li = l;
00289    while (l--) { Real a = *s++; if (minval >= a) { minval = a; li = l; } }
00290    i = storage - li;
00291    (const_cast<GeneralMatrix&>(*this)).tDelete(); return minval;
00292 }
00293 
00294 Real GeneralMatrix::MaximumAbsoluteValue2(int& i, int& j) const
00295 {
00296    REPORT
00297    if (storage == 0) NullMatrixError(this);
00298    Real maxval = 0.0; int nr = Nrows();
00299    MatrixRow mr(const_cast<GeneralMatrix*>(this), LoadOnEntry+DirectPart);
00300    for (int r = 1; r <= nr; r++)
00301    {
00302       int c; maxval = mr.MaximumAbsoluteValue1(maxval, c);
00303       if (c > 0) { i = r; j = c; }
00304       mr.Next();
00305    }
00306    (const_cast<GeneralMatrix&>(*this)).tDelete(); return maxval;
00307 }
00308 
00309 Real GeneralMatrix::MinimumAbsoluteValue2(int& i, int& j) const
00310 {
00311    REPORT
00312    if (storage == 0)  NullMatrixError(this);
00313    Real minval = FloatingPointPrecision::Maximum(); int nr = Nrows();
00314    MatrixRow mr(const_cast<GeneralMatrix*>(this), LoadOnEntry+DirectPart);
00315    for (int r = 1; r <= nr; r++)
00316    {
00317       int c; minval = mr.MinimumAbsoluteValue1(minval, c);
00318       if (c > 0) { i = r; j = c; }
00319       mr.Next();
00320    }
00321    (const_cast<GeneralMatrix&>(*this)).tDelete(); return minval;
00322 }
00323 
00324 Real GeneralMatrix::Maximum2(int& i, int& j) const
00325 {
00326    REPORT
00327    if (storage == 0) NullMatrixError(this);
00328    Real maxval = -FloatingPointPrecision::Maximum(); int nr = Nrows();
00329    MatrixRow mr(const_cast<GeneralMatrix*>(this), LoadOnEntry+DirectPart);
00330    for (int r = 1; r <= nr; r++)
00331    {
00332       int c; maxval = mr.Maximum1(maxval, c);
00333       if (c > 0) { i = r; j = c; }
00334       mr.Next();
00335    }
00336    (const_cast<GeneralMatrix&>(*this)).tDelete(); return maxval;
00337 }
00338 
00339 Real GeneralMatrix::Minimum2(int& i, int& j) const
00340 {
00341    REPORT
00342    if (storage == 0) NullMatrixError(this);
00343    Real minval = FloatingPointPrecision::Maximum(); int nr = Nrows();
00344    MatrixRow mr(const_cast<GeneralMatrix*>(this), LoadOnEntry+DirectPart);
00345    for (int r = 1; r <= nr; r++)
00346    {
00347       int c; minval = mr.Minimum1(minval, c);
00348       if (c > 0) { i = r; j = c; }
00349       mr.Next();
00350    }
00351    (const_cast<GeneralMatrix&>(*this)).tDelete(); return minval;
00352 }
00353 
00354 Real Matrix::MaximumAbsoluteValue2(int& i, int& j) const
00355 {
00356    REPORT
00357    int k; Real m = GeneralMatrix::MaximumAbsoluteValue1(k); k--;
00358    i = k / Ncols(); j = k - i * Ncols(); i++; j++;
00359    return m;
00360 }
00361 
00362 Real Matrix::MinimumAbsoluteValue2(int& i, int& j) const
00363 {
00364    REPORT
00365    int k; Real m = GeneralMatrix::MinimumAbsoluteValue1(k); k--;
00366    i = k / Ncols(); j = k - i * Ncols(); i++; j++;
00367    return m;
00368 }
00369 
00370 Real Matrix::Maximum2(int& i, int& j) const
00371 {
00372    REPORT
00373    int k; Real m = GeneralMatrix::Maximum1(k); k--;
00374    i = k / Ncols(); j = k - i * Ncols(); i++; j++;
00375    return m;
00376 }
00377 
00378 Real Matrix::Minimum2(int& i, int& j) const
00379 {
00380    REPORT
00381    int k; Real m = GeneralMatrix::Minimum1(k); k--;
00382    i = k / Ncols(); j = k - i * Ncols(); i++; j++;
00383    return m;
00384 }
00385 
00386 Real SymmetricMatrix::SumSquare() const
00387 {
00388    REPORT
00389    Real sum1 = 0.0; Real sum2 = 0.0; Real* s = store; int nr = nrows_value;
00390    for (int i = 0; i<nr; i++)
00391    {
00392       int j = i;
00393       while (j--) sum2 += square(*s++);
00394       sum1 += square(*s++);
00395    }
00396    (const_cast<SymmetricMatrix&>(*this)).tDelete(); return sum1 + 2.0 * sum2;
00397 }
00398 
00399 Real SymmetricMatrix::SumAbsoluteValue() const
00400 {
00401    REPORT
00402    Real sum1 = 0.0; Real sum2 = 0.0; Real* s = store; int nr = nrows_value;
00403    for (int i = 0; i<nr; i++)
00404    {
00405       int j = i;
00406       while (j--) sum2 += fabs(*s++);
00407       sum1 += fabs(*s++);
00408    }
00409    (const_cast<SymmetricMatrix&>(*this)).tDelete(); return sum1 + 2.0 * sum2;
00410 }
00411 
00412 Real IdentityMatrix::SumAbsoluteValue() const
00413    { REPORT  return fabs(Trace()); }    // no need to do tDelete?
00414 
00415 Real SymmetricMatrix::Sum() const
00416 {
00417    REPORT
00418    Real sum1 = 0.0; Real sum2 = 0.0; Real* s = store; int nr = nrows_value;
00419    for (int i = 0; i<nr; i++)
00420    {
00421       int j = i;
00422       while (j--) sum2 += *s++;
00423       sum1 += *s++;
00424    }
00425    (const_cast<SymmetricMatrix&>(*this)).tDelete(); return sum1 + 2.0 * sum2;
00426 }
00427 
00428 Real IdentityMatrix::SumSquare() const
00429 {
00430    Real sum = *store * *store * nrows_value;
00431    (const_cast<IdentityMatrix&>(*this)).tDelete(); return sum;
00432 }
00433 
00434 
00435 Real BaseMatrix::SumSquare() const
00436 {
00437    REPORT GeneralMatrix* gm = (const_cast<BaseMatrix&>(*this)).Evaluate();
00438    Real s = gm->SumSquare(); return s;
00439 }
00440 
00441 Real BaseMatrix::NormFrobenius() const
00442    { REPORT  return sqrt(SumSquare()); }
00443 
00444 Real BaseMatrix::SumAbsoluteValue() const
00445 {
00446    REPORT GeneralMatrix* gm = (const_cast<BaseMatrix&>(*this)).Evaluate();
00447    Real s = gm->SumAbsoluteValue(); return s;
00448 }
00449 
00450 Real BaseMatrix::Sum() const
00451 {
00452    REPORT GeneralMatrix* gm = (const_cast<BaseMatrix&>(*this)).Evaluate();
00453    Real s = gm->Sum(); return s;
00454 }
00455 
00456 Real BaseMatrix::MaximumAbsoluteValue() const
00457 {
00458    REPORT GeneralMatrix* gm = (const_cast<BaseMatrix&>(*this)).Evaluate();
00459    Real s = gm->MaximumAbsoluteValue(); return s;
00460 }
00461 
00462 Real BaseMatrix::MaximumAbsoluteValue1(int& i) const
00463 {
00464    REPORT GeneralMatrix* gm = (const_cast<BaseMatrix&>(*this)).Evaluate();
00465    Real s = gm->MaximumAbsoluteValue1(i); return s;
00466 }
00467 
00468 Real BaseMatrix::MaximumAbsoluteValue2(int& i, int& j) const
00469 {
00470    REPORT GeneralMatrix* gm = (const_cast<BaseMatrix&>(*this)).Evaluate();
00471    Real s = gm->MaximumAbsoluteValue2(i, j); return s;
00472 }
00473 
00474 Real BaseMatrix::MinimumAbsoluteValue() const
00475 {
00476    REPORT GeneralMatrix* gm = (const_cast<BaseMatrix&>(*this)).Evaluate();
00477    Real s = gm->MinimumAbsoluteValue(); return s;
00478 }
00479 
00480 Real BaseMatrix::MinimumAbsoluteValue1(int& i) const
00481 {
00482    REPORT GeneralMatrix* gm = (const_cast<BaseMatrix&>(*this)).Evaluate();
00483    Real s = gm->MinimumAbsoluteValue1(i); return s;
00484 }
00485 
00486 Real BaseMatrix::MinimumAbsoluteValue2(int& i, int& j) const
00487 {
00488    REPORT GeneralMatrix* gm = (const_cast<BaseMatrix&>(*this)).Evaluate();
00489    Real s = gm->MinimumAbsoluteValue2(i, j); return s;
00490 }
00491 
00492 Real BaseMatrix::Maximum() const
00493 {
00494    REPORT GeneralMatrix* gm = (const_cast<BaseMatrix&>(*this)).Evaluate();
00495    Real s = gm->Maximum(); return s;
00496 }
00497 
00498 Real BaseMatrix::Maximum1(int& i) const
00499 {
00500    REPORT GeneralMatrix* gm = (const_cast<BaseMatrix&>(*this)).Evaluate();
00501    Real s = gm->Maximum1(i); return s;
00502 }
00503 
00504 Real BaseMatrix::Maximum2(int& i, int& j) const
00505 {
00506    REPORT GeneralMatrix* gm = (const_cast<BaseMatrix&>(*this)).Evaluate();
00507    Real s = gm->Maximum2(i, j); return s;
00508 }
00509 
00510 Real BaseMatrix::Minimum() const
00511 {
00512    REPORT GeneralMatrix* gm = (const_cast<BaseMatrix&>(*this)).Evaluate();
00513    Real s = gm->Minimum(); return s;
00514 }
00515 
00516 Real BaseMatrix::Minimum1(int& i) const
00517 {
00518    REPORT GeneralMatrix* gm = (const_cast<BaseMatrix&>(*this)).Evaluate();
00519    Real s = gm->Minimum1(i); return s;
00520 }
00521 
00522 Real BaseMatrix::Minimum2(int& i, int& j) const
00523 {
00524    REPORT GeneralMatrix* gm = (const_cast<BaseMatrix&>(*this)).Evaluate();
00525    Real s = gm->Minimum2(i, j); return s;
00526 }
00527 
00528 Real DotProduct(const Matrix& A, const Matrix& B)
00529 {
00530    REPORT
00531    int n = A.storage;
00532    if (n != B.storage) Throw(IncompatibleDimensionsException(A,B));
00533    Real sum = 0.0; Real* a = A.store; Real* b = B.store;
00534    while (n--) sum += *a++ * *b++;
00535    return sum;
00536 }
00537 
00538 Real Matrix::Trace() const
00539 {
00540    REPORT
00541    Tracer trace("Trace");
00542    int i = nrows_value; int d = i+1;
00543    if (i != ncols_value) Throw(NotSquareException(*this));
00544    Real sum = 0.0; Real* s = store;
00545 //   while (i--) { sum += *s; s += d; }
00546    if (i) for (;;) { sum += *s; if (!(--i)) break; s += d; }
00547    (const_cast<Matrix&>(*this)).tDelete(); return sum;
00548 }
00549 
00550 Real DiagonalMatrix::Trace() const
00551 {
00552    REPORT
00553    int i = nrows_value; Real sum = 0.0; Real* s = store;
00554    while (i--) sum += *s++;
00555    (const_cast<DiagonalMatrix&>(*this)).tDelete(); return sum;
00556 }
00557 
00558 Real SymmetricMatrix::Trace() const
00559 {
00560    REPORT
00561    int i = nrows_value; Real sum = 0.0; Real* s = store; int j = 2;
00562    // while (i--) { sum += *s; s += j++; }
00563    if (i) for (;;) { sum += *s; if (!(--i)) break; s += j++; }
00564    (const_cast<SymmetricMatrix&>(*this)).tDelete(); return sum;
00565 }
00566 
00567 Real LowerTriangularMatrix::Trace() const
00568 {
00569    REPORT
00570    int i = nrows_value; Real sum = 0.0; Real* s = store; int j = 2;
00571    // while (i--) { sum += *s; s += j++; }
00572    if (i) for (;;) { sum += *s; if (!(--i)) break; s += j++; }
00573    (const_cast<LowerTriangularMatrix&>(*this)).tDelete(); return sum;
00574 }
00575 
00576 Real UpperTriangularMatrix::Trace() const
00577 {
00578    REPORT
00579    int i = nrows_value; Real sum = 0.0; Real* s = store;
00580    while (i) { sum += *s; s += i--; }             // won t cause a problem
00581    (const_cast<UpperTriangularMatrix&>(*this)).tDelete(); return sum;
00582 }
00583 
00584 Real BandMatrix::Trace() const
00585 {
00586    REPORT
00587    int i = nrows_value; int w = lower+upper+1;
00588    Real sum = 0.0; Real* s = store+lower;
00589    // while (i--) { sum += *s; s += w; }
00590    if (i) for (;;) { sum += *s; if (!(--i)) break; s += w; }
00591    (const_cast<BandMatrix&>(*this)).tDelete(); return sum;
00592 }
00593 
00594 Real SymmetricBandMatrix::Trace() const
00595 {
00596    REPORT
00597    int i = nrows_value; int w = lower+1;
00598    Real sum = 0.0; Real* s = store+lower;
00599    // while (i--) { sum += *s; s += w; }
00600    if (i) for (;;) { sum += *s; if (!(--i)) break; s += w; }
00601    (const_cast<SymmetricBandMatrix&>(*this)).tDelete(); return sum;
00602 }
00603 
00604 Real IdentityMatrix::Trace() const
00605 {
00606    Real sum = *store * nrows_value;
00607    (const_cast<IdentityMatrix&>(*this)).tDelete(); return sum;
00608 }
00609 
00610 
00611 Real BaseMatrix::Trace() const
00612 {
00613    REPORT
00614    MatrixType Diag = MatrixType::Dg; Diag.SetDataLossOK();
00615    GeneralMatrix* gm = (const_cast<BaseMatrix&>(*this)).Evaluate(Diag);
00616    Real sum = gm->Trace(); return sum;
00617 }
00618 
00619 void LogAndSign::operator*=(Real x)
00620 {
00621    if (x > 0.0) { log_value += log(x); }
00622    else if (x < 0.0) { log_value += log(-x); sign = -sign; }
00623    else sign = 0;
00624 }
00625 
00626 void LogAndSign::PowEq(int k)
00627 {
00628    if (sign)
00629    {
00630       log_value *= k;
00631       if ( (k & 1) == 0 ) sign = 1;
00632    }
00633 }
00634 
00635 Real LogAndSign::Value() const
00636 {
00637    Tracer et("LogAndSign::Value");
00638    if (log_value >= FloatingPointPrecision::LnMaximum())
00639       Throw(OverflowException("Overflow in exponential"));
00640    return sign * exp(log_value);
00641 }
00642 
00643 LogAndSign::LogAndSign(Real f)
00644 {
00645    if (f == 0.0) { log_value = 0.0; sign = 0; return; }
00646    else if (f < 0.0) { sign = -1; f = -f; }
00647    else sign = 1;
00648    log_value = log(f);
00649 }
00650 
00651 LogAndSign DiagonalMatrix::LogDeterminant() const
00652 {
00653    REPORT
00654    int i = nrows_value; LogAndSign sum; Real* s = store;
00655    while (i--) sum *= *s++;
00656    (const_cast<DiagonalMatrix&>(*this)).tDelete(); return sum;
00657 }
00658 
00659 LogAndSign LowerTriangularMatrix::LogDeterminant() const
00660 {
00661    REPORT
00662    int i = nrows_value; LogAndSign sum; Real* s = store; int j = 2;
00663    // while (i--) { sum *= *s; s += j++; }
00664    if (i) for(;;) { sum *= *s; if (!(--i)) break; s += j++; }
00665    (const_cast<LowerTriangularMatrix&>(*this)).tDelete(); return sum;
00666 }
00667 
00668 LogAndSign UpperTriangularMatrix::LogDeterminant() const
00669 {
00670    REPORT
00671    int i = nrows_value; LogAndSign sum; Real* s = store;
00672    while (i) { sum *= *s; s += i--; }
00673    (const_cast<UpperTriangularMatrix&>(*this)).tDelete(); return sum;
00674 }
00675 
00676 LogAndSign IdentityMatrix::LogDeterminant() const
00677 {
00678    REPORT
00679    int i = nrows_value; LogAndSign sum;
00680    if (i > 0) { sum = *store; sum.PowEq(i); }
00681    (const_cast<IdentityMatrix&>(*this)).tDelete(); return sum;
00682 }
00683 
00684 LogAndSign BaseMatrix::LogDeterminant() const
00685 {
00686    REPORT GeneralMatrix* gm = (const_cast<BaseMatrix&>(*this)).Evaluate();
00687    LogAndSign sum = gm->LogDeterminant(); return sum;
00688 }
00689 
00690 LogAndSign GeneralMatrix::LogDeterminant() const
00691 {
00692    REPORT
00693    Tracer tr("LogDeterminant");
00694    if (nrows_value != ncols_value) Throw(NotSquareException(*this));
00695    CroutMatrix C(*this); return C.LogDeterminant();
00696 }
00697 
00698 LogAndSign CroutMatrix::LogDeterminant() const
00699 {
00700    REPORT
00701    if (sing) return 0.0;
00702    int i = nrows_value; int dd = i+1; LogAndSign sum; Real* s = store;
00703    if (i) for(;;)
00704    {
00705       sum *= *s;
00706       if (!(--i)) break;
00707       s += dd;
00708    }
00709    if (!d) sum.ChangeSign(); return sum;
00710 
00711 }
00712 
00713 Real BaseMatrix::Determinant() const
00714 {
00715    REPORT
00716    Tracer tr("Determinant");
00717    REPORT GeneralMatrix* gm = (const_cast<BaseMatrix&>(*this)).Evaluate();
00718    LogAndSign ld = gm->LogDeterminant();
00719    return ld.Value();
00720 }
00721 
00722 
00723 
00724 
00725 
00726 LinearEquationSolver::LinearEquationSolver(const BaseMatrix& bm)
00727 : gm( ( (const_cast<BaseMatrix&>(bm)).Evaluate() )->MakeSolver() )
00728 {
00729    if (gm==&bm) { REPORT  gm = gm->Image(); }
00730    // want a copy if  *gm is actually bm
00731    else { REPORT  gm->Protect(); }
00732 }
00733 
00734 
00735 #ifdef use_namespace
00736 }
00737 #endif
00738 

newmat11b
Generated Mon May 9 04:54:18 2016 by Doxygen 1.6.3