#############################################################
# Building 'tdgesv' and sample output on the Crays
#############################################################

#############################################################
# First 'rcp' source from SGIs to Cray
#############################################################

einstein% pwd; ls
/usr2/people/phy329/linsys/ex1
Makefile    tdgesv1.f

einstein% rcp * phas761@charon:~/phy329/linsys/ex1

#############################################################
# Now login to the Cray and translate the source.  Note 
# that the script handles the change in name of the 
# real*8 LAPACK routine from 'dgesv' to 'sgesv'.
#############################################################

einstein% rlogin phas761@charon

charon% cd phy329/linsys/ex1; ls
Makefile    tdgesv1.f

charon% f77transcray tdgesv1.f
f77transcray: Translating tdgesv1.f

charon% cat tdgesv1.f
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     Version 1: Uses fixed-size 2-d arrays.
c-----------------------------------------------------------
      program         tdgesv1

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

c-----------------------------------------------------------
c     Storage for arrays 
c-----------------------------------------------------------
      real            a(maxn,maxn),   asave(maxn,maxn),
     &                b(maxn),        x(maxn)
      integer         ipiv(maxn)
      integer         d1a,            d2a       

      integer         i,              nrhs,      
     &                n,              info

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

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

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

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

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.  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 all contiguous in memory.
c-----------------------------------------------------------
      n    = 3
      nrhs = 1

      call sgesv( 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 sgesv() is invalid'
      else 
c-----------------------------------------------------------
c        Matrix is singular.
c-----------------------------------------------------------
         write(0,*) 'tdgesv1: sgesv() detected singular ',
     &              'matrix'
      end if
      

      stop

      end

#############################################################
# The 'blas' library is built-in on the Cray, so we do not 
# have to specify '-lblas' during the load phase.  The 
# Makefile works properly provided LIBBLAS is either 
# unset, or set to nothing (see ~/.cshrc on this account).
#############################################################

charon% make
        cf77 -g -c tdgesv1.f
        cf77 -g -L/hpcf/u0/ph/as/phas761/lib tdgesv1.o \
             -lp329f -llapack  -o tdgesv1

charon% tdgsev1


 5.426364412431667,  -0.325775376817397,  -0.4083508069894641

 STOP executed at line 89 in Fortran routine 'TDGESV1'
 CPU: 0.005s,  Wallclock: 0.024s,  4.5% of 5-CPU Machine
 Memory HWM: 212775, Stack HWM: 16509, Stack segment expansions: 0