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


kni
Author(s): Martin Günther
autogenerated on Thu Aug 27 2015 13:40:07