ccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c  MAIN
c
c  N: number of rows and/or columns in A,x,b matrices
c  A(,): A matrix
c  x():  solution matrix
c  b():  b matrix
c  i,j:  iterator vars
c  
ccccccccccccccccccccccccccccccccccccccccccccccccccccc

      implicit none

      INTEGER N,i,j
      REAL A(10,10),x(10),b(10)

      N=10

cccc define A and b matrices
      CALL AandB(A,b,N)

cccc solving Ax=b
      CALL GAUSS(A,x,b,N)

cccc writing output to file
      OPEN(UNIT=90,file='hw2_data',STATUS='OLD')
      WRITE(90,*) 'The solution to Ax=b with'

      WRITE(90,*) ' '
      WRITE(90,*) ' '

      WRITE(90,*)'MATRIX A = '
      DO 91, i=1,N
         WRITE(90,*) ( A(i,j),j=1,N)
 91   CONTINUE

      WRITE(90,*) ' '
      WRITE(90,*) ' '

      WRITE(90,*) 'MATRIX b = ' 
      DO 92, i=1,N
         WRITE(90,*) b(i)
 92   CONTINUE

      WRITE(90,*) ' '
      WRITE(90,*) ' '

      DO 94, i=1,N
         WRITE(90,*) 'x(',i, ' ) = ',  x(i)
 94   CONTINUE

      CLOSE(90)

      END


cccccccccccccccccccccccccccccccccccccccccc
c
c  SUBROUTINE GAUSS: GAUSSIAN ELIMINATION
c
c  A(N,N),x(N),b(N): Ax=b arrays
c  Ap(N,N+1):  Augmented Matrix, A and b
c  row, column:  row and column numbers
c  norm:  normalizing value
c  i,j:  iterators 
c  test,epsil: checks for singular matrices 
c  
cccccccccccccccccccccccccccccccccccccccccc

      SUBROUTINE GAUSS(A,x,b,N)

      implicit none

      INTEGER N,row,i,j,column,bcol
      REAL    Ap(N,N+1),A(N,N),x(N),b(N),norm,epsil,test

      bcol = N+1
      epsil = 1.E-4
      test=1.
      
cccc reassign A to Ap
      DO 10 i=1, N
	 DO 12 j=1, N
	    Ap(i,j) = A(i,j)
 12      CONTINUE
 10   CONTINUE

cccc reassign b to Ap
      DO 16 i=1, N
	 Ap(i,bcol) = b(i)
 16   CONTINUE


      DO 20 row=1, N

	 CALL PIVOT(Ap,N,row,bcol)
       
cccc check for singularity  
         test = abs( Ap(row,row) )
         IF(test.LT.epsil) THEN
            PRINT*,'Matrix may be singular... program terminated'
            STOP
         END IF

cccc normalize all elements
	 column = row

	 DO 30 i = row+1, N
	    norm = Ap(i,column)/Ap(row,column)
            Ap(i,column) = 0

            DO 40 j=column+1,bcol
               Ap(i,j) = Ap(i,j) - norm*Ap(row,j)
 40         CONTINUE

 30      CONTINUE

 20   CONTINUE

cccc back substitution
      CALL BACKSUB(Ap,N,x,bcol)

      RETURN
      END	 



ccccccccccccccccccccccccccccccccccccccccc
c
c  SUBROUTINE PIVOT: from GAUSS
c 
c  row0:  initial row
c  bcol:  last column of Ap (b vector)
c  val:   value of element
c  maxval:  maximum value in column
c  maxrow:  row with maxval
c  tempvar: temporary variable
c
ccccccccccccccccccccccccccccccccccccccccc

      SUBROUTINE PIVOT(Ap,N,row0,bcol)

      implicit none

      INTEGER N,row0,bcol,col,i,j,maxrow
      REAL    Ap(N,bcol),maxval,val,tempvar

      col = row0
      maxval = abs( Ap(row0,col) )
      maxrow = row0

      DO 100 i=row0+1, N
	 val = abs( Ap(i,col) )
         IF(val.GT.maxval) THEN
            maxval = val
	    maxrow = i
	 END IF
 100  CONTINUE

      IF(maxrow.NE.row0) THEN
         DO 110 j=col, bcol
	    tempvar = Ap(row0,j)
            Ap(row0,j) = Ap(maxrow,j)
	    Ap(maxrow,j) = tempvar
 110     CONTINUE
      END IF

      RETURN
      END



cccccccccccccccccccccccccccccccccccccccc
c
c  SUBROUTINE BACKSUB: from GAUSS
c
c  all variables defined above
c
cccccccccccccccccccccccccccccccccccccccc

      SUBROUTINE BACKSUB(Ap,N,x,bcol)

      implicit none

      INTEGER N,bcol,i,j
      REAL    Ap(N,bcol),x(N),temp

      x(N) = Ap(N,bcol) / Ap(N,N)

      DO 200 i=N-1, 1, -1

	 temp = Ap(i,bcol)
         DO 210 j=i+1, N
            temp = temp - x(j)*Ap(i,j)
 210     CONTINUE

         x(i) = temp / Ap(i,i)

 200  CONTINUE

      RETURN
      END 

ccccccccccccccccccccccccccccccccccccccc
c
c  SUBROUTINE AandB: Sets A and b matrices
c
ccccccccccccccccccccccccccccccccccccccc
      SUBROUTINE AandB(A,b,N)
      
      implicit none
      
      INTEGER N
      REAL A(N,N), b(N)

      A(1,1) = 1.
      A(1,2) = 0.
      A(1,3) = 1.
      A(1,4) = 0.
      A(1,5) = 1.
      A(1,6) = 0.
      A(1,7) = 1.
      A(1,8) = 0.
      A(1,9) = 1.
      A(1,10) =0.
      
      A(2,1) = 0.
      A(2,2) = 1.
      A(2,3) = 0.		        
      A(2,4) = 1.
      A(2,5) = 0.
      A(2,6) = 1.
      A(2,7) = 0.
      A(2,8) = 1.
      A(2,9) = 0.
      A(2,10) =1.

      A(3,1) = 1.
      A(3,2) = 2.
      A(3,3) = 3.
      A(3,4) = 3.
      A(3,5) = 1.
      A(3,6) = 1.
      A(3,7) = 1.
      A(3,8) = 0.
      A(3,9) = 0.
      A(3,10) =0.

      A(4,1) = 0.
      A(4,2) = 3.
      A(4,3) = 0.
      A(4,4) = 0.
      A(4,5) = 0.
      A(4,6) = 3.
      A(4,7) = 3.
      A(4,8) = 3.
      A(4,9) = 3.
      A(4,10) =4.

      A(5,1) = 0.
      A(5,2) = 1.
      A(5,3) = 0.
      A(5,4) = 0.
      A(5,5) = 1.
      A(5,6) = 0.
      A(5,7) = 0.
      A(5,8) = 1.
      A(5,9) = 0.
      A(5,10) =0.
      
      A(6,1) = 0.
      A(6,2) = 0.
      A(6,3) = 1.
      A(6,4) = 0.
      A(6,5) = 1.
      A(6,6) = 0.
      A(6,7) = 0.
      A(6,8) = 0.
      A(6,9) = 0.
      A(6,10) =1.
      
      A(7,1) = 4.
      A(7,2) = 5.
      A(7,3) = 0.
      A(7,4) = 0.
      A(7,5) = 7.
      A(7,6) = 6.
      A(7,7) = 0.
      A(7,8) = 0.
      A(7,9) = 6.
      A(7,10) =6.
      
      A(8,1) = 4.
      A(8,2) = 0.
      A(8,3) = 0.
      A(8,4) = 0.
      A(8,5) = 4.
      A(8,6) = 0.
      A(8,7) = 0.
      A(8,8) = 4.
      A(8,9) = 0.
      A(8,10) =0.
 
      A(9,1) = 0.
      A(9,2) = 4.
      A(9,3) = 0.
      A(9,4) = 0.
      A(9,5) = 0.
      A(9,6) = 4.
      A(9,7) = 0
      A(9,8) = 0.
      A(9,9) = 0.
      A(9,10) =4.
      
      A(10,1) = 0.
      A(10,2) = 0.
      A(10,3) = 0.
      A(10,4) = 0.
      A(10,5) = 0.
      A(10,6) = 0.
      A(10,7) = 0.
      A(10,8) = 1.
      A(10,9) = 2.
      A(10,10) =3.

      b(1) = 1.
      b(2) = 0.
      b(3) = 0.
      b(4) = 0.
      b(5) = 0.
      b(6) = 0.
      b(7) = 0.
      b(8) = 0.
      b(9) = 0.
      b(10) =0.
      
      RETURN
      END
								  

      
