Solution of linear algebraic equations

PROBLEM
We wish to find a solution to the matrix equation A.*X.=B. where operator * denotes matrix multiplications, A. is a square matrix, B. is a known right-hand side vector, and X. is an unknown vector.

ALGORITHM
The LUBKSB algorithm (for forward substitution and backsubstitution) solves the set of N linear equations A.*X.=B., here A. is output from LUDCMP routine. For given matrix A. and vector B. solves the following statements:

 call LUDCMP N; call LUBKSB N

the matrix equation A.*X.=B. and the solution is saved in B. vector.

IMPLEMENTATION
Unit: nonrecursive internal subroutine

Global variables: input N by N matrix A. - LU decomposition, input vector Indx. - permutation vector; A. and Indx. are determined by the LUDCMP routine

Parameters: a positive integer N

Result: the B. vector includes the solution of matrix equation given above

Note:
The A. matrix and Indx. vector are not modified. We can use them for successive calls with different B. vectors.

 LUBKSB: procedure expose A. B. Indx. parse arg N L = 0 do I = 1 to N   P = Indx.I; Sum = B.P; B.P = B.I   if L <> 0     then do       do J = L to I - 1         Sum = Sum - A.I.J * B.J       end     end     else if Sum <> 0 then L = I   B.I = Sum end do I = N to 1 by -1   Sum = B.I   do J = I + 1 to N     Sum = Sum - A.I.J * B.J   end   B.I = Sum / A.I.I end return

EXAMPLE

 N=3 A.1.1 = 2; A.1.2 = -1; A.1.3 = -1; B.1 = 4 A.2.1 = 3; A.2.2 =  4; A.2.3 = -2; B.2 = 11 A.3.1 = 3; A.3.2 = -2; A.3.3 =  4; B.3 = 11 call LUDCMP N; call LUBKSB N S = "" do J = 1 to N; S = S B.J; end say S

displays 3 1 1

CONNECTIONS
Linear Algebraic Equations
Determinant of a matrix
Inverse of a matrix
LU decomposition

Literature
Press W.H., Teukolsky S.A., Vetterling W.T., Flannery B.P. Numerical Recipes in C : the art of scientific computing
- 2nd ed. University Press, Cambridge, 1992