#! /bin/sh
# This is a shell archive, meaning:
# 1. Remove everything above the #! /bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create the files:
#	Doc/
#	Fortran/
#	Fortran/Sp/
#	Fortran/Sp/Drivers/
#	Fortran/Sp/Drivers/Makefile
#	Fortran/Sp/Drivers/driver.f
#	Fortran/Sp/Drivers/res
#	Fortran/Sp/Src/
#	Fortran/Sp/Src/src.f
# This archive created: Thu Dec 15 17:15:36 2005
export PATH; PATH=/bin:$PATH
if test ! -d 'Doc'
then
	mkdir 'Doc'
fi
cd 'Doc'
cd ..
if test ! -d 'Fortran'
then
	mkdir 'Fortran'
fi
cd 'Fortran'
if test ! -d 'Sp'
then
	mkdir 'Sp'
fi
cd 'Sp'
if test ! -d 'Drivers'
then
	mkdir 'Drivers'
fi
cd 'Drivers'
if test -f 'Makefile'
then
	echo shar: will not over-write existing file "'Makefile'"
else
cat << "SHAR_EOF" > 'Makefile'
all: Res

src.o: src.f
	$(F77) $(F77OPTS) -c src.f

driver.o: driver.f
	$(F77) $(F77OPTS) -c driver.f

DRIVERS= driver1
RESULTS= Res

Objs1= driver.o src.o
driver1: $(Objs1)
	$(F77) $(F77OPTS) -o driver1 $(Objs1) $(SRCLIBS)
Res: driver1 
	./driver1 >Res

diffres:Res res
	echo "Differences in results from driver"
	$(DIFF) Res res

clean: 
	rm -rf *.o $(DRIVERS) $(CLEANUP) $(RESULTS)
SHAR_EOF
fi # end of overwriting check
if test -f 'driver.f'
then
	echo shar: will not over-write existing file "'driver.f'"
else
cat << "SHAR_EOF" > 'driver.f'
      program main

c*******************************************************************************
c
cc TOMS423_PRB tests the TOMS423 routines DECOMP and SOLVE.
c
c  Discussion:
c
c    Solve A*x = b where A is a given matrix, and B a right hand side.
c
      implicit none

      integer n
      integer ndim

      parameter ( n = 3 )
      parameter ( ndim = 4 )

      real a(ndim,ndim)
      real b(ndim)
      real determ
      integer i
      integer ip(ndim)
      integer j

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TOMS423_PRB'
      write ( *, '(a)' ) '  DECOMP factors a general matrix;'
      write ( *, '(a)' ) '  SOLVE solves a factored linear system;'
      write ( *, '(a)' ) ' '
      write ( *, '(a,i6)' ) '  The matrix dimension NDIM = ', ndim
      write ( *, '(a,i6)' ) '  The number of equations is N = ', n
c
c  Set the values of the matrix A.
c
      a(1,1) = 1.0E+00
      a(1,2) = 2.0E+00
      a(1,3) = 3.0E+00

      a(2,1) = 4.0E+00
      a(2,2) = 5.0E+00
      a(2,3) = 6.0E+00

      a(3,1) = 7.0E+00
      a(3,2) = 8.0E+00
      a(3,3) = 0.0E+00

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  The matrix A:'
      write ( *, '(a)' ) ' '

      do i = 1, n
        write ( *, '(2x,3g14.6)' ) ( a(i,j), j = 1, n )
      end do
c
c  Set the values of the right hand side vector B.
c
      b(1) = 14.0E+00
      b(2) = 32.0E+00
      b(3) = 23.0E+00

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  The right hand side B is '
      write ( *, '(a)' ) ' '
      do i = 1, n
        write ( *, '(2x,g14.6)' ) b(i)
      end do
c
c  Factor the matrix.
c
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  Factor the matrix'

      call decomp ( n, ndim, a, ip )
c
c  Get the determinant.
c
      determ = ip(n)
      do i = 1, n
        determ = determ * a(i,i)
      end do

      write ( *, '(a)' ) ' '
      write ( *, '(a,g14.6)' ) '  Matrix determinant = ', determ

      if ( ip(n) .eq. 0 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) '  DECOMP reports the matrix is singular.'
        stop
      end if
c
c  If no error occurred, now use DGESL to solve the system.
c
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  Solve the linear system.'

      call solve ( n, ndim, a, b, ip )
c
c  Print the results.
c
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  SOLVE returns the solution:'
      write ( *, '(a)' ) '  (Should be (1,2,3))'
      write ( *, '(a)' ) ' '

      do i = 1, n
        write ( *, '(2x,g14.6)' ) b(i)
      end do

      stop
      end
SHAR_EOF
fi # end of overwriting check
if test -f 'res'
then
	echo shar: will not over-write existing file "'res'"
else
cat << "SHAR_EOF" > 'res'
 
TOMS423_PRB
  DECOMP factors a general matrix;
  SOLVE solves a factored linear system;
 
  The matrix dimension NDIM =      4
  The number of equations is N =      3
 
  The matrix A:
 
     1.00000       2.00000       3.00000    
     4.00000       5.00000       6.00000    
     7.00000       8.00000       0.00000    
 
  The right hand side B is 
 
     14.0000    
     32.0000    
     23.0000    
 
  Factor the matrix
 
  Matrix determinant =    27.0000    
 
  Solve the linear system.
 
  SOLVE returns the solution:
  (Should be (1,2,3))
 
     1.00000    
     2.00000    
     3.00000    
SHAR_EOF
fi # end of overwriting check
cd ..
if test ! -d 'Src'
then
	mkdir 'Src'
fi
cd 'Src'
if test -f 'src.f'
then
	echo shar: will not over-write existing file "'src.f'"
else
cat << "SHAR_EOF" > 'src.f'
      SUBROUTINE DECOMP(N, NDIM, A, IP)
      INTEGER N, NDIM, K, KP1, M, I, J
      REAL A(NDIM,NDIM),T
      INTEGER IP(NDIM)
C
C  MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION.
C  INPUT..
C    N = ORDER OF MATRIX.
C    NDIM = DECLARED DIMENSION OF ARRAY A.
C    A = MATRIX TO BE TRIANGULARIZED.
C  OUTPUT..
C    A(I,J), I.LE.J = UPPER TRIANGULAR FACTOR U.
C    A(I,J), I.GT.J = MULTIPLIERS = LOWER TRIANGULAR
C                     FACTOR, I-L.
C    IP(K), K.LT.N = INDEX OF K-TH PIVOT ROW.
C    IP(N) = (-1)**(NUMBER OF INTERCHANGES) OR 0.
C  USE 'SOLVE' TO OBTRAIN SOLUTION OF LINEAR SYSTEM.
C  DETERM(A) = IP(N)*A(1,1)*A(2,2)*...*A(N,N).
C  IF IP(N) = 0, A IS SINGULAR, SOLVE WILL DIVIDE BY ZERO.
C  INTERCHANGES FINISHED IN U, ONLY PARTLY IN L.
C
      IP(N) = 1
      DO 6 K = 1,N
        IF (K.EQ.N) GO TO 5
        KP1 = K+1
        M = K
        DO 1 I = KP1,N
          IF(ABS(A(I,K)).GT.ABS(A(M,K)))M=I
1       CONTINUE
        IP(K) = M
        IF(M.NE.K) IP(N) = -IP(N)
        T = A(M,K)
        A(M,K) = A(K,K)
        A(K,K) = T
        IF(T.EQ.0.) GO TO 5
        DO 2 I = KP1,N
          A(I,K) = -A(I,K)/T
2       CONTINUE
        DO 4 J = KP1,N
          T = A(M,J)
          A(M,J) = A(K,J)
          A(K,J) = T
          IF(T.EQ.0.) GO TO 4
          DO 3 I = KP1,N
            A(I,J) = A(I,J) + A(I,K)*T
3         CONTINUE
4       CONTINUE
5       IF(A(K,K).EQ.0.) IP(N) = 0
6     CONTINUE
      RETURN
      END

      SUBROUTINE SOLVE(N, NDIM, A, B, IP)
      INTEGER N, NDIM, NM1, K, KP1, M, I, KB, KM1
      REAL A(NDIM,NDIM),B(NDIM),T
      INTEGER IP(NDIM)
C
C  SOLUTION OF LINEAR SYSTEM, A*X = B.
C  INPUT..
C    N = ORDER 0F MATRIX.
C    NDIM = DECLARED DIMENSION OF ARRAY A.
C    A = TRIANGULARIZED MATRIX OBTAINED FROM 'DECOMP'.
C    B = RIGHT HAND SIDE VECTOR.
C    IP = PIVOT VECTOR OBTAINED FROM 'DECOMP'.
C  DO NOT USE IF DECOMP HAS SET IP(N) = 0.
C  OUTPUT..
C    B = SOLUTION VECTOR, X.
C
      IF(N.EQ.1) GO TO 9
      NM1 = N-1
      DO 75 K = 1,NM1
        KP1 = K+1
        M = IP(K)
        T = B(M)
        B(M) = B(K)
        B(K) = T
        DO 7 I = KP1, N
          B(I) = B(I) + A(I,K)*T
7       CONTINUE
75    CONTINUE
      DO 85 KB = 1,NM1
        KM1 = N-KB
        K = KM1+1
        B(K) = B(K)/A(K,K)
        T = -B(K)
        DO 8 I = 1,KM1
          B(I) = B(I) + A(I,K)*T
8       CONTINUE
85    CONTINUE
9     B(1) = B(1)/A(1,1)
      RETURN
      END
SHAR_EOF
fi # end of overwriting check
cd ..
cd ..
cd ..
#       End of shell archive
exit 0