55 lines
1.2 KiB
C++
55 lines
1.2 KiB
C++
#include "ABACUS.h"
|
|
using namespace std;
|
|
|
|
void ABACUS::ludcmp_CX (SQMat_CX& a, Vect_INT& indx, DP& d)
|
|
{
|
|
const complex<DP> TINY = 1.0e-200;
|
|
int i, j, k;
|
|
int imax = 0;
|
|
complex<DP> big, dum, sum, temp;
|
|
|
|
int n = a.size();
|
|
if (n == 0) ABACUSerror("Zero-sized array in ludcmp_CX");
|
|
Vect_CX vv(n);
|
|
d = 1.0;
|
|
for (i = 0; i < n; i++) {
|
|
big = 0.0;
|
|
for (j = 0; j < n; j++)
|
|
if ((abs(temp = a[i][j])) > abs(big)) big = temp;
|
|
if (big == 0.0) throw Divide_by_zero();
|
|
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 ((abs(dum = vv[i]*sum)) >= abs(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;
|
|
}
|
|
}
|
|
}
|