c===========================================================
c     Test program for LAPACK "driver" routine 'dgesv' 
c     which computes the solution of a real system
c     of linear equations:  A x = b
c
c     This version uses fixed-size 2-d arrays (size fixed at 
c     some maximum value commensurate with needs and/or 
c     available memory), illustrating another commonly used 
c     Fortran technique to implement run-time dimensioning, 
c     PARTICULARLY FOR RANK-2 ARRAYS.
c     
c     This time the rules are as follows:  All subroutines and
c     functions which manipulate the array must be passed:
c
c     (1)  The array itself.
c     (2)  The "true" or "physical" dimensions;
c          i.e. the dimensions in MAIN (*).
c     (3)  The "run-time" or "logical" dimensions (*).
c 
c     (*)  More precisely, due to the nature of the FORTRAN 
c          subscripting computation, the leading d-1 dimensions 
c          must be passed for a rank-d array.  In particular, 
c          for rank-2 array (matrices), THE leading physical
c          dimension (often denoted 'LDA' in LAPACK code), and 
c          THE leading logical dimension (often denoted 'N')
c          must BOTH be passed.
c
c     Passing the physical dimensions ensures that the 
c     linearization/subscripting calculation is identical 
c     in all program units INCUDING MAIN---so that, e.g.,
c     one can safely and conveniently use a(i,j) etc. in 
c     MAIN.
c
c     Passing the logical dimensions allows us to write 
c     routines which function for a general case (here,
c     typically for N x N matrices).
c
c     Passing BOTH sets of dimensions is slightly cumbersome, 
c     but is the price we pay in this case for convenience 
c     and generality.  
c===========================================================
      program         tdgesv

      implicit        none
c-----------------------------------------------------------
c     Maximum size of linear system.
c-----------------------------------------------------------
      integer         maxn
      parameter     ( maxn = 100 )

c-----------------------------------------------------------
c     Storage for arrays.
c-----------------------------------------------------------
      real*8          a(maxn,maxn),   
     &                b(maxn)
      integer         ipiv(maxn)

      integer         i,              nrhs,      
     &                n,              info

c-----------------------------------------------------------
c     Set up sample 3 x 3 system ...
c-----------------------------------------------------------
      a(1,1) =  1.23d0
      a(1,2) =  0.24d0
      a(1,3) = -0.45d0

      a(2,1) = -0.43d0
      a(2,2) =  2.45d0
      a(2,3) =  0.78d0

      a(3,1) =  0.51d0
      a(3,2) = -0.68d0
      a(3,3) =  3.23d0

      b(1)   =  6.78d0
      b(2)   = -3.45d0
      b(3)   =  1.67d0

c-----------------------------------------------------------
c     ... and solve it.  Note that 'dgsev' is general 
c     enough to solve A x_i = b_i for multiple right-hand-
c     sides b_i.  Here we have only one right-hand-side.
c     Also note that the procedure performs the LU 
c     decomposition in place, thus destroying the 
c     input-matrix, it also overwrites the right-hand-side(s)
c     with the solution(s).  Finally, observe that we 
c     pass the "leading dimension" (maxn) of both 'a' and 
c     'b' to the routine.  Again, this allows us to load array
c     elements in the main program as we have just done,
c     without running into troubles due to the fact that 
c     these elements ARE NOT, in general all contiguous in 
c     memory.  This certainly includes the current 3 x 3 case.
c-----------------------------------------------------------
      n    = 3
      nrhs = 1

      call dgesv( n, nrhs, a, maxn, ipiv, b, maxn, info )

      if(      info .eq. 0 ) then
c-----------------------------------------------------------
c        Solution successful, write soln to stdout.       
c        Note the use of "implied-do-loop" to write a 
c        sequence of elements: the enclosing parenthesis 
c        around the  "loop" are required.
c-----------------------------------------------------------
         write(*,*) ( b(i) , i = 1 , n )
      else if( info .lt. 0 ) then
c-----------------------------------------------------------
c        Bad argument detected.
c-----------------------------------------------------------
         write(0,*) 'tdgesv1: Argument ', abs(info), 
     &              ' to dgesv() is invalid'
      else 
c-----------------------------------------------------------
c        Matrix is singular.
c-----------------------------------------------------------
         write(0,*) 'tdgesv1: dgesv() detected singular ',
     &              'matrix'
      end if
      

      stop

      end