123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354 |
- #include "ABACUS.h"
- using namespace std;
-
- void ABACUS::ludcmp (SQMat_DP& a, Vect_INT& indx, DP& d)
- {
- const DP TINY = 1.0e-200;
- int i, j, k;
- int imax = 0;
- DP big, dum, sum, temp;
-
- int n = a.size();
- if (n == 0) ABACUSerror("Zero-sized array in ludcmp");
- Vect_DP vv(n);
- d = 1.0;
- for (i = 0; i < n; i++) {
- big = 0.0;
- for (j = 0; j < n; j++)
- if ((temp = fabs(a[i][j])) > big) big = temp;
- if (big == 0.0) throw Divide_by_zero(); //ABACUSerror("Singular matrix in routine ludcmp.");
- vv[i] = 1.0/big;
- }
- for (j = 0; j < n; j++) {
- for (i = 0; i < j; i++) {
- sum = a[i][j];
- for (k = 0; k < i; k++) sum -= a[i][k] * a[k][j];
- a[i][j] = sum;
- }
- big = 0.0;
- for (i = j; i < n; i++) {
- sum = a[i][j];
- for (k = 0; k < j; k++) sum -= a[i][k] * a[k][j];
- a[i][j] = sum;
- if ((dum = vv[i]*fabs(sum)) >= big) {
- big = dum;
- imax = i;
- }
- }
- if (j != imax) {
- for (k = 0; k < n; k++) {
- dum = a[imax][k];
- a[imax][k] = a[j][k];
- a[j][k] = dum;
- }
- d = -d;
- vv[imax] = vv[j];
- }
- indx[j] = imax;
- if (a[j][j] == 0.0) a[j][j] = TINY;
- if (j != n-1) {
- dum = 1.0/(a[j][j]);
- for (i = j + 1; i < n; i++) a[i][j] *= dum;
- }
- }
- }
|