1 // Public domain code from Yi-Kuo Yu & Stephen Altschul, NCBI
2
lubksb(double ** a,int n,int * indx,double b[])3 void lubksb(double **a, int n, int *indx, double b[]) {
4 int i, ii = 0, ip, j;
5 double sum;
6
7 for (i = 1; i <= n; i++) {
8 ip = indx[i];
9 sum = b[ip];
10 b[ip] = b[i];
11 if (ii) {
12 for (j = ii; j <= i - 1; j++) {
13 sum -= a[i][j] * b[j];
14 }
15 } else if (sum) {
16 ii = i;
17 }
18 b[i] = sum;
19 }
20 for (i = n; i >= 1; i--) {
21 sum = b[i];
22 for (j = i + 1; j <= n; j++) {
23 sum -= a[i][j] * b[j];
24 }
25 b[i] = sum / a[i][i];
26 }
27 }
28