'*****************************************************
'* Calculate the determinant of a real square matrix *
'* A(n,n) by Gauss method with full pivoting.        *
'* ------------------------------------------------- *
'* Ref.: "Algebre - Algorithmes et programmes en     *
'*        Pascal By Jean-Louis Jardrin, Dunod -      *
'*        Bordas Paris, 1988 p. 76-79" [BIBLI 10].   *
'* ------------------------------------------------- *
'* SAMPLE RUN:                                       *
'* (Calculate the determinant of matrix:             *
'*            10 18  1 14 22                         *
'*             4 12 25  8 16                         *
'*            23  6 19  2 15                         *
'*            17  5 13 21  9                         *
'*            11 24  7 20  3  )                      *
'*                                                   *
'* Input size of square real matrix: 5               *
'* Input desired precision: 1e-10                    *
'*                                                   *
'* Line 1                                            *
'* Element 1: 10                                     *
'* Element 2: 18                                     *
'* Element 3: 1                                      *
'* Element 4: 14                                     *
'* Element 5: 22                                     *
'*                                                   *
'* Line 2                                            *
'* Element 1: 4                                      *
'* Element 2: 12                                     *
'* Element 3: 25                                     *
'* Element 4:  8                                     *
'* Element 5: 16                                     *
'*                                                   *
'* Line 3                                            *
'* Element 1: 23                                     *
'* Element 2:  6                                     *
'* Element 3: 19                                     *
'* Element 4:  2                                     *
'* Element 5: 15                                     *
'*                                                   *
'* Line 4                                            *
'* Element 1: 17                                     *
'* Element 2:  5                                     *
'* Element 3: 13                                     *
'* Element 4: 21                                     *
'* Element 5:  9                                     *
'*                                                   *
'* Line 5                                            *
'* Element 1: 11                                     *
'* Element 2: 24                                     *
'* Element 3:  7                                     *
'* Element 4: 20                                     *
'* Element 5:  3                                     *
'*                                                   *
'* Determinant = -4680000                            *
'*                                                   *
'*                                                   *
'*       Basic version by Jean-Pierre Moreau, Paris. *
'*                    (www.jpmoreau.fr)              *
'*****************************************************
'PROGRAM Calculate_deter
defint i-n
defdbl a-h,o-z

MAXSIZE=25

'integer n                   size of matrix A
dim A(MAXSIZE,MAXSIZE)      'pointer to input matrix
'real d0                     determinant of matrix A
'real eps                    desired precision

  cls
  print
  print " Calculate the determinant of a real square matrix"
  print "        by Gauss method with full pivoting."
  print
  input " Input size of square real matrix: ", n
' read precision
  input " Input desired precision: ", eps
' read matrix A(n,n)
  print
  for i=1 to n
    print " Line "; i
    for j=1 to n
      print " Element "; j; ": ";
      input "", A(i,j)
    next j
    print
  next i

  gosub 2000  'd0=DMGT(eps,n,A)

  print
  print  " Determinant = "; d0
  print
  print

END


'The subroutine TSRGT applies to input real square matrix A(n,n) the upper
'triangularization algorithm of Gauss method with full pivoting and keeps
'trace of successive transformations done in integer vectors KP and LP.
'-----------------------------------------------------------------------------
'  Input parameters:
'  eps        precision (real*8)
'  n          size of A matrix (integer)
'  A          pointer to input real square matrix (real*8)
'  Output parameters:
'  it         flag=1 if A matrix ok, =0 if A matrix is singular (integer)
'  C          pointer to table storing main diagonal elements and supra-
'             diagonal elements of upper triangular matrix and the multi-
'             plying coefficients used during triangularization process
'  KP         table storing informations concerning the column exchanges
'             during process (integer)
'  LP         table storing informations concerning the line exchanges
'             during process (integer)
'-----------------------------------------------------------------------------
'The table C is first initialized to A matrix, then receives at each step k
'of the triangularization process, usefull elements of A matrix at step k for
'k=1,2,...n.
'The variables po(real*8), lo and ko(integer) store respectively pivot at step k,
'its line number and its column number.
'------------------------------------------------------------------------------
1000 'Subroutine TSRGT(eps, n, A, it, C, Kp, Lp)
dim C(MAXSIZE,MAXSIZE), Kp(MAXSIZE),Lp(MAXSIZE)
' set C to A
  for i=1 to n
    for j=1 to n
      C(i,j)=A(i,j)
    next j
  next i	
  it=1: k=1
1010 po=C(k,k): lo=k: ko=k
    for i=k to n
      for j=k to n
        if abs(C(i,j))>abs(po) then
          po=C(i,j): lo=i: ko=j
        end if
      next j
    next i
    Lp(k)=lo: Kp(k)=ko
    if abs(po)<eps then
      it=0
    else
      if lo <> k then
        for j=k to n
          t0=C(k,j): C(k,j)=C(lo,j): C(lo,j)=t0
        next j
      end if
      if ko <> k then
        for i=1 to n
          t0=C(i,k): C(i,k)=C(i,ko): C(i,ko)=t0
        next i
      end if
      for i=k+1 to n
        C(i,k)=C(i,k)/po
        for j=k+1 to n
          C(i,j)=C(i,j)-C(i,k)*C(k,j)
        next j
      next i
      k=k+1
    end if
  if it=1 and k<n then goto 1010
  if it=1 and abs(C(n,n))<eps then it=0
return 'TSRGT

'The function DMGT returns the determinant of a real square matrix
'A(n,n) by Gauss method with full pivoting.
'----------------------------------------------------------------------------
'  Input parameters:
'  eps        precision (real*8)
'  n          size of A matrix (integer)
'  A          pointer to input real square matrix
'  Output parameters:
'  None
'-----------------------------------------------------------------------------
'The procedure TSRGT is used to reduce A matrix to an upper triangular matrix.
'Output variables are it(integer), C(n,n), Kp(n) and Lp(n).
'If it=0, matrix A is singular, if it=1, matrix A is regular. Table C contains
'at location i,j (j>=i) the corresponding element of the upper triangular matrix.
'Tables Lp and Kp contain informations relative to exchanges of line or column
'that occured during the process. For instance, the element number k of Lp is
'an integer <> k if an exchange of line has been made at step k (k=1,2,...,n).
'The number of exchanges of lines and columns is stored in integer L. the
'determinant of A matrix is stored in d0 (real*8).
'-----------------------------------------------------------------------------
2000 'd0=DMGT(eps, n, A)

  gosub 1000  'call TSRGT(eps,n,A,it,C,Kp,Lp)
  if it=0 then
    d0=0#    'matrix singular, det=0
  else       'matrix regular, det<>0
    d0=1#
    for k=1 to n
	  d0=d0*C(k,k)
    next k
    l=0
    for k=1 to n-1
      if Lp(k) <> k then l=l+1
      if Kp(k) <> k then l=l+1
    next k
	if l MOD 2 <> 0 then d0=-d0  'l is odd
  end if
return

'End of file deter.bas