'*****************************************************
'* Calculate the determinant of a real square matrix *
'* by the LU decomposition method                    *
'* ------------------------------------------------- *
'* SAMPLE RUN:                                       *
'* (Calculate the determinant of real square matrix: *
'*        | 1  0.5  0  0    0  0    0  0    0 |      *
'*        |-1  0.5  1  0    0  0    0  0    0 |      *
'*        | 0 -1    0  2.5  0  0    0  0    0 |      *
'*        | 0  0   -1 -1.5  4  0    0  0    0 |      *
'*        | 0  0    0 -1   -3  5.5  0  0    0 |      *
'*        | 0  0    0  0   -1 -4.5  7  0    0 |      *
'*        | 0  0    0  0    0 -1   -6  8.5  0 |      *
'*        | 0  0    0  0    0  0   -1 -7.5 10 |      *
'*        | 0  0    0  0    0  0    0 -1   -9 |  )   *
'*                                                   *
'*                                                   *
'*  Determinant =  1                                 *
'*                                                   *
'*                                                   *
'*               Basic version by J-P Moreau, Paris. *
'*                        (www.jpmoreau.fr)          *
'*****************************************************
' Exact value is: 1
'-----------------------------------------------------
DEFINT I-N
DEFDBL A-H, O-Z

CLS : PRINT
n = 9

DIM A(n, n)
DIM INDX(n)

'matrix A(i,j) column by column
DATA  1,-1,0,0,0,0,0,0,0
DATA  0.5,0.5,-1,0,0,0,0,0,0
DATA  0,1,0,-1,0,0,0,0,0
DATA  0,0,2.5,-1.5,-1,0,0,0,0
DATA  0,0,0,4,-3,-1,0,0,0
DATA  0,0,0,0,5.5,-4.5,-1,0,0
DATA  0,0,0,0,0,7,-6,-1,0
DATA  0,0,0,0,0,0,8.5,-7.5,-1
DATA  0,0,0,0,0,0,0,10,-9

'read matrix A(i,j)
FOR j = 1 TO n
  FOR i = 1 TO n
    READ A(i, j)
  NEXT i
NEXT j

GOSUB 1000  'call LUDCMP(A,N,INDX,ID,CODE)

FOR j = 1 TO n
  D = D * A(j, j)
NEXT j

PRINT
PRINT " Determinant = "; D
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, j)) > AMAX THEN AMAX = ABS(A(i, 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, j)
      FOR K = 1 TO i - 1
        SUM = SUM - A(i, K) * A(K, j)
      NEXT K
      A(i, j) = SUM
    NEXT i  'i loop
    AMAX = 0#
    FOR i = j TO n
      SUM = A(i, j)
      FOR K = 1 TO j - 1
        SUM = SUM - A(i, K) * A(K, j)
      NEXT K
      A(i, 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, K)
        A(IMAX, K) = A(j, K)
        A(j, K) = DUM
      NEXT K  'k loop
      D = -D
      VV(IMAX) = VV(j)
    END IF

    INDX(j) = IMAX
    IF ABS(A(j, j)) < TINY THEN
      A(j, j) = TINY
    END IF

    IF j <> n THEN
      DUM = 1# / A(j, j)
      FOR i = j + 1 TO n
        A(i, j) = A(i, j) * DUM
      NEXT i
    END IF
  NEXT j  'j loop

RETURN  'subroutine LUDCMP

'End of file deter1.bas