'*********************************************************
'* Inversion of a real square matrix by LU decomposition *
'* with dynamic allocations                              *
'*                                                       *
'*                    Basic Version By J-P Moreau, Paris *
'*                            (www.jpmoreau.fr)          *
'* ----------------------------------------------------- *
'* Reference:                                            *
'*                                                       *
'*  "Numerical Recipes By W.H. Press, B. P. Flannery,    *
'*   S.A. Teukolsky and W.T. Vetterling, Cambridge       *
'*   University Press, 1986" [BIBLI 08].                 *
'* ----------------------------------------------------- *
'* SAMPLE RUN:                                           *
'*                                                       *
'* Input file (inv_lu.dat):                              *
'*                                                       *
'*  4                                                    *
'*  8  2    3  12                                        *
'*  2  4    7   0.25                                     *
'*  3  7    3   5                                        *
'* 12  0.25 5   2                                        *
'*                                                       *
'* Output file (inv_lu.lst):                             *
'*                                                       *
'* ----------------------------------------------------- *
'*  INVERSION OF A REAL SQUARE MATRIX:                   *
'* ----------------------------------------------------- *
'*  N=4                                                  *
'*                                                       *
'*  8.000000  2.000000  3.000000  12.00000               *
'*  2.000000  4.000000  7.000000  0.250000               *
'*  3.000000  7.000000  3.000000  5.000000               *
'*  12.00000  0.250000  5.000000  2.000000               *
'*                                                       *
'*  Inverted matrix Y:                                   *
'*                                                       *
'* -0.040155 -0.085788  0.056557  0.110262               *
'* -0.085788 -0.062018  0.202201  0.016978               *
'*  0.056557  0.202201 -0.130316 -0.038826               *
'*  0.110262  0.016978 -0.038826 -0.066631               *
'*                                                       *
'*  Verification A*Y = I:                                *
'*                                                       *
'*  1.000000 -0.000000  0.000000  0.000000               *
'* -0.000000  1.000000 -0.000000 -0.000000               *
'* -0.000000 -0.000000  1.000000  0.000000               *
'*  0.000000  0.000000 -0.000000  1.000000               *
'* ----------------------------------------------------- *
'* Subroutines used:                                     *
'* 1000 LUDCMP: LU decomposition of a real square matrix *
'* 2000 LUBKSB: Solve linear system AX = B, A being the  *
'*              output of LUDCMP.                        * 
'*********************************************************
'Program Inv_matrix
 DEFINT I-N
 DEFDBL A-H, O-Z

  CLS : PRINT
  INPUT " Data file name (without .dat): ", nom$
  in$ = nom$ + ".dat"
  ou$ = nom$ + ".lst"

  OPEN in$ FOR INPUT AS #1
  OPEN ou$ FOR OUTPUT AS #2

  INPUT #1, N     'size of linear system

  n1 = N + 1
  n2 = n1 * n1
  DIM A(n2), A1(n2), Y(n2), B(n1), temp(n1), INDX(n1)

'------------------------------------------------------------------
'NOTA: Matrices A, A1 and Y of size N x N are stored in a vector,
'      line by line. The element (i,j) of the matrix is at location
'      i(n+1) + j in the vector. The access is faster.
'------------------------------------------------------------------

  PRINT #2, "--------------------------------------------------"
  PRINT #2, "      INVERSION OF A REAL SQUARE MATRIX:"
  PRINT #2, "--------------------------------------------------"
  PRINT #2,
  PRINT #2, "  N="; N
  PRINT #2,

  FOR I = 1 TO N
    FOR J = 1 TO N
      INPUT #1, temp(J)  'read a line
      A(I * n1 + J) = temp(J)
      A1(I * n1 + J) = temp(J)   'copy of A
      Y(I * n1 + J) = 0#
      PRINT #2, USING "###.######"; temp(J);
    NEXT J
    Y(I * n1 + I) = 1#
    PRINT #2,
  NEXT I
  CLOSE #1

'call the LU decomposition routine only once
  GOSUB 1000  'LUDCMP(A,n,INDX,D,rc);

'call solver if previous return code is ok
'to obtain inverse of A one column at a time
  IF rc = 0 THEN
    FOR JJ = 1 TO N   'Caution, J is used in 2000
      FOR I = 1 TO N
	B(I) = Y(I * n1 + JJ)
      NEXT I
      GOSUB 2000  'LUBKSB(A,n,INDX,B);
      FOR I = 1 TO N
	Y(I * n1 + JJ) = B(I)
      NEXT I
    NEXT JJ
  END IF
' inverse of matrix A is now in matrix Y
' the original matrix A is destroyed

  IF rc = 1 THEN
    PRINT #2, " The matrix is singular !"
  ELSE
    PRINT #2,
    PRINT #2, "  Inverted matrix:"
    PRINT #2,
    FOR I = 1 TO N
      FOR J = 1 TO N
	PRINT #2, USING "###.######"; Y(I * n1 + J);
      NEXT J
      PRINT #2,
    NEXT I
  END IF

' verify A1 x Y = I (result put in A)
' Multiply A1 by Y and put in A
  FOR I = 1 TO N
    FOR J = 1 TO N
      SUM = 0#
      FOR K = 1 TO N
	SUM = SUM + A1(I * n1 + K) * Y(K * n1 + J)
	A(I * n1 + J) = SUM
      NEXT K
    NEXT J
  NEXT I
' A should now contain identity matrix
  PRINT #2,
  PRINT #2, "  Verification A*Y = I:"
  PRINT #2,
  FOR I = 1 TO N
    FOR J = 1 TO N
      PRINT #2, USING "###.######"; A(I * n1 + J);
    NEXT J
    PRINT #2,
  NEXT I

  PRINT #2, "--------------------------------------------------"
  PRINT #2,
  CLOSE #2
  PRINT
  PRINT " Results in file "; ou$
  PRINT

END

  '***************************************************************
  '* Given an N x N matrix A, this routine replaces it by the LU *
  '* decomposition of a rowwise permutation of itself. A and N   *
  '* are input. INDX is an output vector which records the row   *
  '* permutation effected by the partial pivoting; D is output   *
  '* as -1 or 1, depending on whether the number of row inter-   *
  '* changes was even or odd, respectively. This routine is used *
  '* in combination with LUBKSB to solve linear equations or to  *
  '* invert a matrix. Return code is 1, if matrix is singular.   *
  '***************************************************************
1000 'PROCEDURE LUDCMP;
  NMAX = 100: TINY = 1.5E-16
  DIM VV(N)
  D = 1: CODE = 0:

  FOR I = 1 TO N
    AMAX = 0#
    FOR J = 1 TO N
      IF ABS(A(I * n1 + J)) > AMAX THEN AMAX = ABS(A(I * n1 + J))
    NEXT J
    IF (AMAX < TINY) THEN
      CODE = 1
      RETURN
    END IF
    VV(I) = 1# / AMAX
  NEXT I  'i loop

  FOR J = 1 TO N
    FOR I = 1 TO J - 1
      SUM = A(I * n1 + J)
      FOR K = 1 TO I - 1
	SUM = SUM - A(I * n1 + K) * A(K * n1 + J)
      NEXT K
      A(I * n1 + J) = SUM
    NEXT I  'i loop
    AMAX = 0#
    FOR I = J TO N
      SUM = A(I * n1 + J)
      FOR K = 1 TO J - 1
	SUM = SUM - A(I * n1 + K) * A(K * n1 + J)
      NEXT K
      A(I * n1 + J) = SUM
      DUM = VV(I) * ABS(SUM)
      IF DUM >= AMAX THEN
	IMAX = I
	AMAX = DUM
      END IF
    NEXT I  'i loop

    IF J <> IMAX THEN
      FOR K = 1 TO N
	DUM = A(IMAX * n1 + K)
	A(IMAX * n1 + K) = A(J * n1 + K)
	A(J * n1 + K) = DUM
      NEXT K  'k loop
      D = -D
      VV(IMAX) = VV(J)
    END IF

    INDX(J) = IMAX
    IF ABS(A(J * n1 + J)) < TINY THEN
      A(J * n1 + J) = TINY
    END IF

    IF J <> N THEN
      DUM = 1# / A(J * n1 + J)
      FOR I = J + 1 TO N
	A(I * n1 + J) = A(I * n1 + J) * DUM
      NEXT I
    END IF
  NEXT J  'j loop

RETURN  'subroutine LUDCMP


  '******************************************************************
  '* Solves the set of N linear equations A . X := B.  Here A is    *
  '* input, not as the matrix A but rather as its LU decomposition, *
  '* determined by the routine LUDCMP. INDX is input as the permuta-*
  '* tion vector returned by LUDCMP. B is input as the right-hand   *
  '* side vector B, and returns with the solution vector X. A, N and*
  '* INDX are not modified by this routine and can be used for suc- *
  '* cessive calls with different right-hand sides. This routine is *
  '* also efficient for plain matrix inversion.                     *
  '******************************************************************
2000  'PROCEDURE LUBKSB;
  II = 0
  FOR I = 1 TO N
    LL = INDX(I)
    SUM = B(LL)
    B(LL) = B(I)
    IF II <> 0 THEN
      FOR J = II TO I - 1
	SUM = SUM - A(I * n1 + J) * B(J)
      NEXT J
    ELSEIF SUM <> 0 THEN
      II = I
    END IF
    B(I) = SUM
  NEXT I  'i loop

  FOR I = N TO 1 STEP -1
    SUM = B(I)
    IF I < N THEN
      FOR J = I + 1 TO N
	SUM = SUM - A(I * n1 + J) * B(J)
      NEXT J
    END IF
    B(I) = SUM / A(I * n1 + I)
  NEXT I  'i loop
RETURN

'End of file inv_lu.bas