'*********************************************************
'*  Solving a linear system AX = B by LU decomposition   *
'*                                                       *
'*                    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 (test_lu.dat):                             *
'*                                                       *
'*  4                                                    *
'*  8  2    3  12     25.0                               *
'*  2  4    7   0.25  13.25                              *
'*  3  7    3   5     18.0                               *
'* 12  0.25 5   2     19.25                              *
'*                                                       *
'* Output file (test_lu.lst):                            *
'*                                                       *
'* ---------------------------------------------------   *
'*  LINEAR SYSTEM TO BE SOLVED:                          *
'* ---------------------------------------------------   *
'*  N=4                                                  *
'*                                                       *
'*  8.000000  2.000000  3.000000  12.00000  25.00000     *
'*  2.000000  4.000000  7.000000  0.250000  13.25000     *
'*  3.000000  7.000000  3.000000  5.000000  18.00000     *
'*  12.00000  0.250000  5.000000  2.000000  19.25000     *
'*                                                       *
'*  System solution:                                     *
'*                                                       *
'*  X1=  1.000000                                        *
'*  X2=  1.000000                                        *
'*  X3=  1.000000                                        *
'*  X4=  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.                        * 
'*                                                       *
'*********************************************************
 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

  n2 = (N + 1) * (N + 1)
  DIM A(n2), B(N + 1), temp(N + 2), INDX(N + 1)

'------------------------------------------------------------------
'NOTA: Matrix A of size N x N is 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, "      LINEAR SYSTEM TO BE SOLVED:"
  PRINT #2, "--------------------------------------------------"
  PRINT #2,
  PRINT #2, "  N="; N
  PRINT #2,

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

  GOSUB 1000  'LUDCMP(A,n,INDX,D,rc);

  IF rc = 0 THEN GOSUB 2000   'LUBKSB(A,n,INDX,B);

  IF rc = 1 THEN
    PRINT #2, " The system matrix is singular !"
  ELSE
    PRINT #2,
    PRINT #2, "  System solution:"
    PRINT #2,
    FOR I = 1 TO N
      PRINT #2, USING "  X## = ##.######"; I; B(I)
      PRINT #2,
    NEXT I
  END IF
  PRINT #2, "--------------------------------------------------"
  PRINT #2, " End of file "; ou$
  CLOSE #2
  PRINT
  PRINT " Results in file "; ou$
  PRINT

END 'of main program

  '***************************************************************
  '* 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 'SUBROUTINE LUDCMP
  NMAX = 100: TINY = 1.5E-16
  DIM VV(N)
  D = 1: CODE = 0: N1 = N + 1

  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  'SUBROUTINE LUBKSB
  II = 0: N1 = N + 1
  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 lu.bas