lin1.f

      program linsolve
      real      a(6,6),at(6,6),b(6),bt(6),bc(6)
      integer ipiv(6)
c
c    Driver program to test linear equation solvers
c    Sets up a 6x6 system of linear equations
c    This is more specific than might be desired.  A better form
c    would declare these arrays to be allocatable.
c
c
c    a  -  the array (2-dimensional)of linear equation coefficients
c          the first index gives equation number
c          the second index gives the unknown number
c
c    at -  array to keep a copy of a
c
c    b  -  array (1-dimensional) containing the right hand sides
c          of the linear equations
c
c    bt -  array to keep a copy of bt, in the end it contains the solutions
c          to the equations.
c
c    bc -  an array that should match the original right hand side array
c          generated by multiplying the solutions times the appropriate
c          coefficients
c
c    ipiv- special "pivoting" array used to provide an improved solution
c          in the Gauss Elimination
c
c
c    Set up the A matrix and Right Hand Side Generation
c
c    We start by loading a desired solution in the bt array
c    and creating a matrix "a" with random element values
c
      do 20 i=1,6
         bt(i)=float(i)
        do 20 j=1,6
c         	a(i,j)=float(i+j-(i*j)/7*4+1)
      	   a(i,j)=rand()
 20   continue
c
c     Generate the RHS
c
      do 30 i=1,6
      	b(i)=0.
      	do 30 j=1,6
 30   	   b(i)=b(i)+a(i,j)*bt(j)
      write(6,*) 'Right Hand Side of the linear equation Ax=b'
      write(6,2000)(b(i),i=1,6)
      write(6,*) ' Matrix A'
      write(6,2000)((a(i,j),j=1,6),i=1,6)
 2000 format(1x,6f10.4)
c
c     LINPACK or LAPACK ROUTINES overwrite the A matrix and b array
c     If you need them again make a copy
c
      do 40 i=1,6
      	bt(i)=b(i)
      	do 40 j=1,6
   40	   at(i,j)=a(i,j)	
c
c    Insert calls to LINPACK or LAPACK routines here
c
      call sgefa(at,6,6,ipiv,info)
c
c     Check status by writing the info variable from the LU routine
c
      write(6,*) ' INFO = ',info
      call sgesl(at,6,6,ipiv,bt,0)
c
c    Write out the solution   (this assumes you used bt in the last LINPACK
c    call.
c
      write(6,*) ' Solution of the Linear System:'
      write(6,2000) (bt(i),i=1,6)
c
c     Check Solution by Multipling Ax
c
      do 50 i=1,6
         bc(i)=0.
         do 50 j=1,6
            bc(i)=bc(i)+a(i,j)*bt(j)
   50 continue
      write(6,*)' Ax = '
      write(6,2000) (bc(i),i=1,6)
      stop
      end