newmat8.cpp
Go to the documentation of this file.00001
00002
00003
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
00028
00029 void CroutMatrix::ludcmp()
00030
00031
00032
00033
00034
00035 {
00036 REPORT
00037 Tracer trace( "Crout(ludcmp)" ); sing = false;
00038 Real* akk = store;
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
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065 indx[k] = mu;
00066
00067 if (mu != k)
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
00085
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;
00097 akk += nrows_value + 1;
00098 }
00099 }
00100
00101 void CroutMatrix::lubksb(Real* B, int mini)
00102 {
00103 REPORT
00104
00105
00106
00107
00108
00109
00110 Tracer trace("Crout(lubksb)");
00111 if (sing) Throw(SingularException(*this));
00112 int i, j, ii = nrows_value;
00113
00114
00115
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
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
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
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()); }
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
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
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
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--; }
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
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
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
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
00731 else { REPORT gm->Protect(); }
00732 }
00733
00734
00735 #ifdef use_namespace
00736 }
00737 #endif
00738