diff options
Diffstat (limited to 'pkg/utilities/nttools/stxtools/lubksb.f')
| -rw-r--r-- | pkg/utilities/nttools/stxtools/lubksb.f | 50 |
1 files changed, 50 insertions, 0 deletions
diff --git a/pkg/utilities/nttools/stxtools/lubksb.f b/pkg/utilities/nttools/stxtools/lubksb.f new file mode 100644 index 00000000..429cf2e4 --- /dev/null +++ b/pkg/utilities/nttools/stxtools/lubksb.f @@ -0,0 +1,50 @@ + SUBROUTINE LUBKSB (A, N, NP, INDX, B) + +C Solves a matrix equation AX = B. Before using this routine you must +C call ludcmp to decompose the matrix A (in-place) into lower and upper +C triangular portions. The vector B is input to this routine, and the +C answer X is returned in B. +C +C real a(np,np) i: matrix returned by ludcmp +C int n i: logical size of a is n x n +C int np i: space allocated for a is np x np +C int indx(n) i: index returned by ludcmp +C real b(n) io: input b, output x in equation ax = b +C +C 1988 Oct 28 From Numerical Recipes + + INTEGER N, NP + REAL A(NP,NP) + INTEGER INDX(N) + REAL B(N) + + REAL SUM + INTEGER II, LL, I, J + + II = 0 + DO 20 I = 1, N + LL = INDX(I) + SUM = B(LL) + B(LL) = B(I) + IF (II .NE. 0) THEN + DO 10 J = II, I-1 + SUM = SUM - A(I,J) * B(J) + 10 CONTINUE + ELSE IF (SUM .NE. 0.) THEN + II = I + ENDIF + B(I) = SUM + 20 CONTINUE + + DO 40 I = N, 1, -1 + SUM = B(I) + IF (I .LT. N) THEN + DO 30 J = I+1, N + SUM = SUM - A(I,J) * B(J) + 30 CONTINUE + ENDIF + B(I) = SUM / A(I,I) + 40 CONTINUE + + RETURN + END |