Implementation

The subroutine SOLVE is designed to follow the steps of the Algorithm 1

      subroutine solve(a,maxrow,neq,flag,b)
      implicit none
!
      integer :: maxrow,neq,flag
      double precision a(maxrow,*),b(*)
!-----------------------------------------------------------------------
!
!  SOLVE solves the linear system A*x=b using the factorization
!  obtained from FACTOR.  Do not use SOLVE if a zero pivot has
!  been detected in FACTOR.
!
!  Input variables:
!     A      = an array containing the coefficient matrix A.
!              On output it has the effect of the elimination process.
!     MAXROW = maximum number of equations allowed; the declared row 
!              dimension of A
!     NEQ    = actual number of equations to be solved;  NEQ cannot
!              exceed MAXROW.
!     B      = right hand side vector b.
!
!  Output variables:
!     A      = it has the effect of the elimination process.
!
!     B      = solution vector x.
!
!     FLAG   = an integer variable that reports whether or not the
!              matrix A has a zero pivot. 
!            - A value of FLAG = 0 means all pivots were nonzero;
!            - if positive, the first zero pivot occurred at equation
!              FLAG and the decomposition could not be completed. 
!            - If FLAG = -1 then there is an input error (NEQ or
!              MAXROW not positive or NEQ > MAXROW).
!             
!
!  Declare local variables and initialize:
      double precision :: t
      integer :: i,j,k,m
      double precision :: zero,one
      data zero/0.d0/,one/1.d0/
!
      if ((neq .le. 0) .or. (maxrow .le. 0) .or. (neq .gt. maxrow)) then
         flag = -1
         return
      endif

      flag = 0
      if (neq .eq. 1) then
!
!        neq = 1 is a special case.
!
         if (a(1,1) .eq. zero) then
            flag = 1
            return
         else
            b(1) = b(1)/a(1,1)
         endif
      else
!
!     Gaussian elimination
!
      do k = 1,neq-1
!
!        Check for a nonzero pivot; if pivot is zero return,
!
         if (a(k,k) .eq. zero) then
            flag = k
            return
         endif             
!
!     Forward elimination.
!
         do i = k+1,neq
            t = a(i,k)/a(k,k)
            a(i,k) = -t
            if (t .ne. zero) then
               do j = k+1,neq
                  a(i,j) = a(i,j)-t*a(k,j)
               end do
               b(i) = b(i) -t*b(k)
            endif
         end do
      end do
!
      if (a(neq,neq) .eq. zero) then
         flag = neq
         return
      endif
!
!
!        Back substitution.
!
         do i = neq,1,-1

         .
         .  you must complete this portion
         .

         end do
      endif

      return
      end

      double precision :: a(10,10),b(10),temp(10)

      data zero/0.d0/,one/1.d0/

• Driver Program