23 #define REPORT { static ExeCounter ExeCount(__LINE__,8); ++ExeCount; } 39 Tracer tr(
"Crout(ludcmp)" ); sing =
false;
42 Real big = fabs(*akk);
int mu = 0;
Real* ai = akk;
int k;
44 for (k = 1; k < nrows_val; k++)
46 ai += nrows_val;
const Real trybig = fabs(*ai);
47 if (big < trybig) { big = trybig; mu = k; }
51 if (nrows_val)
for (k = 0;;)
71 Real* a1 = store + nrows_val * k;
72 Real* a2 = store + nrows_val * mu;
d = !
d;
74 while (j--) {
const Real temp = *a1; *a1++ = *a2; *a2++ = temp; }
77 Real diag = *akk; big = 0; mu = k + 1;
80 ai = akk;
int i = nrows_val - k - 1;
83 ai += nrows_val;
Real* al = ai;
84 Real mult = *al / diag; *al = mult;
85 int l = nrows_val - k - 1;
Real* aj = akk;
90 *(++al) -= (mult * *(++aj));
91 const Real trybig = fabs(*al);
92 if (big < trybig) { big = trybig; mu = nrows_val - i - 1; }
93 while (l--) *(++al) -= (mult * *(++aj));
98 if (++k == nrows_val)
break;
112 Tracer tr(
"Crout(lubksb)");
114 int i, j, ii = nrows_val;
118 for (i = 0; i < nrows_val; i++)
120 int ip = indx[i];
Real temp = B[ip]; B[ip] = B[i]; B[i] = temp;
121 if (temp != 0.0) { ii = i;
break; }
129 bi = B + ii; ai = store + ii + i * nrows_val;
132 int ip = indx[i];
Real sum = B[ip]; B[ip] = B[i];
133 Real* aij = ai;
Real* bj = bi; j = i - ii;
134 while (j--) sum -= *aij++ * *bj++;
136 if (++i == nrows_val)
break;
141 ai = store + nrows_val * nrows_val;
143 for (i = nrows_val - 1; i >= mini; i--)
145 Real* bj = B+i; ai -= nrows_val;
Real* ajx = ai+i;
147 j = nrows_val - i;
while(--j) sum -= *(++ajx) * *(++bj);
160 while (i--) sum +=
square(*s++);
168 while (i--) sum += fabs(*s++);
175 Real sm = 0.0;
int i = storage;
Real* s = store;
176 while (i--) sm += *s++;
222 Real maxval = 0.0;
int l = storage;
Real* s = store;
223 while (l--) {
Real a = fabs(*s++);
if (maxval < a) maxval =
a; }
231 Real maxval = 0.0;
int l = storage;
Real* s = store;
int li = storage;
233 {
Real a = fabs(*s++);
if (maxval <= a) { maxval =
a; li = l; } }
242 int l = storage - 1;
Real* s = store;
Real minval = fabs(*s++);
243 while (l--) {
Real a = fabs(*s++);
if (minval > a) minval =
a; }
251 int l = storage - 1;
Real* s = store;
Real minval = fabs(*s++);
int li = l;
253 {
Real a = fabs(*s++);
if (minval >= a) { minval =
a; li = l; } }
262 int l = storage - 1;
Real* s = store;
Real maxval = *s++;
263 while (l--) {
Real a = *s++;
if (maxval < a) maxval =
a; }
271 int l = storage - 1;
Real* s = store;
Real maxval = *s++;
int li = l;
272 while (l--) {
Real a = *s++;
if (maxval <= a) { maxval =
a; li = l; } }
281 int l = storage - 1;
Real* s = store;
Real minval = *s++;
282 while (l--) {
Real a = *s++;
if (minval > a) minval =
a; }
290 int l = storage - 1;
Real* s = store;
Real minval = *s++;
int li = l;
291 while (l--) {
Real a = *s++;
if (minval >= a) { minval =
a; li = l; } }
300 Real maxval = 0.0;
int nr = Nrows();
302 for (
int r = 1; r <= nr; r++)
305 if (c > 0) { i = r; j = c; }
317 for (
int r = 1; r <= nr; r++)
320 if (c > 0) { i = r; j = c; }
332 for (
int r = 1; r <= nr; r++)
334 int c; maxval = mr.
Maximum1(maxval, c);
335 if (c > 0) { i = r; j = c; }
347 for (
int r = 1; r <= nr; r++)
349 int c; minval = mr.
Minimum1(minval, c);
350 if (c > 0) { i = r; j = c; }
360 i = k / Ncols(); j = k - i * Ncols(); i++; j++;
368 i = k / Ncols(); j = k - i * Ncols(); i++; j++;
376 i = k / Ncols(); j = k - i * Ncols(); i++; j++;
384 i = k / Ncols(); j = k - i * Ncols(); i++; j++;
391 Real sum1 = 0.0;
Real sum2 = 0.0;
Real* s = store;
int nr = nrows_val;
392 for (
int i = 0; i<nr; i++)
395 while (j--) sum2 +=
square(*s++);
404 Real sum1 = 0.0;
Real sum2 = 0.0;
Real* s = store;
int nr = nrows_val;
405 for (
int i = 0; i<nr; i++)
408 while (j--) sum2 += fabs(*s++);
420 Real sum1 = 0.0;
Real sum2 = 0.0;
Real* s = store;
int nr = nrows_val;
421 for (
int i = 0; i<nr; i++)
424 while (j--) sum2 += *s++;
432 Real sum = *store * *store * nrows_val;
540 while (n--) sum += *a++ * *b++;
548 int i = nrows_val;
int d = i+1;
552 if (i)
for (;;) { sum += *s;
if (!(--i))
break; s +=
d; }
559 int i = nrows_val;
Real sum = 0.0;
Real* s = store;
560 while (i--) sum += *s++;
567 int i = nrows_val;
Real sum = 0.0;
Real* s = store;
int j = 2;
569 if (i)
for (;;) { sum += *s;
if (!(--i))
break; s += j++; }
576 int i = nrows_val;
Real sum = 0.0;
Real* s = store;
int j = 2;
578 if (i)
for (;;) { sum += *s;
if (!(--i))
break; s += j++; }
585 int i = nrows_val;
Real sum = 0.0;
Real* s = store;
586 while (i) { sum += *s; s += i--; }
593 int i = nrows_val;
int w = lower_val+upper_val+1;
596 if (i)
for (;;) { sum += *s;
if (!(--i))
break; s += w; }
603 int i = nrows_val;
int w = lower_val+1;
606 if (i)
for (;;) { sum += *s;
if (!(--i))
break; s += w; }
627 if (x > 0.0) { log_val += log(x); }
628 else if (x < 0.0) { log_val += log(-x); sign_val = -sign_val; }
637 if ( (k & 1) == 0 ) sign_val = 1;
643 Tracer et(
"LogAndSign::value");
646 return sign_val * exp(log_val);
651 if (f == 0.0) { log_val = 0.0; sign_val = 0;
return; }
652 else if (f < 0.0) { sign_val = -1; f = -f; }
661 while (i--) sum *= *s++;
670 if (i)
for(;;) { sum *= *s;
if (!(--i))
break; s += j++; }
678 while (i) { sum *= *s; s += i--; }
686 if (i > 0) { sum = *store; sum.
PowEq(i); }
699 Tracer tr(
"log_determinant");
707 if (sing)
return 0.0;
730 gm = ( ((
BaseMatrix&)bm).Evaluate() )->MakeSolver();
731 if (gm==&bm) {
REPORT gm = gm->Image(); }
733 else {
REPORT gm->Protect(); }
740 int nr = gm->
nrows();
742 if (gm->size() == 0) {
REPORT ssq = 0.0; }
746 for (
int i = 1; i <= nr; ++i)
751 while (s--) sum +=
square(*in++);
764 int nr = gm->
nrows();
int nc = gm->ncols();
769 for (
int i = 1; i <= nr; ++i)
773 while (s--) *out++ +=
square(*in++);
785 int nr = gm->
nrows();
787 if (gm->size() == 0) {
REPORT sum_vec = 0.0; }
791 for (
int i = 1; i <= nr; ++i)
796 while (s--) sum += *in++;
809 int nr = gm->
nrows();
int nc = gm->ncols();
814 for (
int i = 1; i <= nr; ++i)
818 while (s--) *out++ += *in++;
virtual Real minimum_absolute_value() const
Real minimum2(int &i, int &j) const
LogAndSign log_determinant() const
Real norm_Frobenius() const
Real minimum2(int &i, int &j) const
Miscellaneous exception (details in character string).
Real maximum_absolute_value2(int &i, int &j) const
virtual Real sum_square() const
ReturnMatrix for_return() const
virtual Real maximum() const
LogAndSign log_determinant() const
static void NullMatrixError(const GeneralMatrix *gm)
void operator*=(Real)
multiply by a real
Real minimum_absolute_value1(int &i) const
virtual Real maximum2(int &i, int &j) const
void pow_eq(int k)
raise to power of k
virtual Real trace() const
virtual Real minimum2(int &i, int &j) const
ReturnMatrix sum_rows() const
Real minimum_absolute_value2(int &i, int &j) const
GeneralMatrix * Evaluate(MatrixType mt=MatrixTypeUnSp)
virtual Real minimum_absolute_value2(int &i, int &j) const
ReturnMatrix sum_columns() const
virtual Real maximum_absolute_value1(int &i) const
Real maximum_absolute_value() const
virtual Real minimum_absolute_value1(int &i) const
A matrix is not square exception.
virtual Real minimum() const
LinearEquationSolver(const BaseMatrix &bm)
Real minimum_absolute_value() const
Real sum_absolute_value() const
virtual Real minimum1(int &i) const
LogAndSign log_determinant() const
Real maximum1(int &i) const
Real value() const
the value
Real sum_absolute_value() const
Real maximum2(int &i, int &j) const
void lubksb(Real *, int=0)
Real Maximum1(Real r, int &i)
virtual Real sum_absolute_value() const
Real maximum_absolute_value1(int &i) const
virtual Real maximum_absolute_value() const
Singular matrix exception.
Real maximum2(int &i, int &j) const
Real dotproduct(const Matrix &A, const Matrix &B)
Real maximum_absolute_value2(int &i, int &j) const
Base of the matrix classes.
Real sum_square(const BaseMatrix &B)
LogAndSign log_determinant() const
ReturnMatrix sum_square_rows() const
LogAndSign log_determinant() const
The usual rectangular matrix.
Real MaximumAbsoluteValue1(Real r, int &i)
Incompatible dimensions exception.
virtual Real maximum_absolute_value2(int &i, int &j) const
FloatVector FloatVector * a
LogAndSign log_determinant() const
Real minimum_absolute_value2(int &i, int &j) const
Real trace(const BaseMatrix &B)
ReturnMatrix sum_square_columns() const
Real sum(const BaseMatrix &B)
Real MinimumAbsoluteValue1(Real r, int &i)
The classes for matrices that can contain data are derived from this.
Real sum_absolute_value() const
Real minimum1(int &i) const
Real Minimum1(Real r, int &i)
virtual Real maximum1(int &i) const
virtual LogAndSign log_determinant() const