c======================================================
      program dpa1
c DOUBLE PRECISION VERSION
c======================================================
c This program solves a system of linear equations by using 
c Gaussian Elimination with full, partial, and no pivoting.
c Results from the corresponding IMSL function are also obtained
      implicit real*8(a-h,o-z)
      parameter(nvar=100)
      dimension a(nvar,nvar),b(nvar),x(nvar)

      call read_data(a,b,nvar,nsys)

      write (*,*) 'Results:'
      write (*,*)
      call no_piv_ge(a,b,nsys,nvar,x)
      write (*,*) 'G.E. routine without pivoting'
      do i=1,nsys
        write (*,'(a,i3,a,f16.8)') 'x(',i,') = ',x(i)
      end do
      call part_piv_ge(a,b,nsys,nvar,x,singul)
      write (*,*)
      write (*,*) 'ge routine with partial pivoting'
      do i=1,nsys
        write (*,'(a,i3,a,f16.8)') 'x(',i,') = ',x(i)
      end do
      call full_piv_ge(a,b,nsys,nvar,x,singul)
      write (*,*)
      write (*,*) 'G.E. routine with pivoting'
      do i=1,nsys
        write (*,'(a,i3,a,f16.8)') 'x(',i,') = ',x(i)
      end do
!!    initialize x
      do i=1,nvar
        x(i) = 0.0
      end do
      CALL DLSARG (nsys, A, nvar, B, 1, x)
      write (*,*)
      write (*,*) 'results from IMSL'
      do i=1,nsys
        write (*,'(a,i3,a,f16.8)') 'x(',i,') = ',x(i)
      end do

      stop
      end
      

c======================================================
      subroutine read_data(a,b,n,nsys)
c======================================================
      implicit real*8(a-h,o-z)
      dimension a(n,n),b(n)
      character*20 fname, form*50

	write(*,*)
        write(*,*) 'Enter configuration file name: '
        read (*,'(a)') fname
	open(unit=10,file=fname,status='unknown')
	read(10,*) nsys
        if(nsys.eq.5) then 
           form='(5f10.5)'
        else
           form='(6f10.5)'
        end if
	read(10,*) ((a(i,j), j=1,nsys), i=1,nsys)
	read(10,*) (b(i), i=1,nsys)
	write(*,*)
	write(*,'(a,i3,a,i3)') ' System size : ',nsys,'  x',nsys
	write(*,*)
	write(*,*) 'xmulticient matrix : '
	write(20,fmt=form) ((a(i,j), j=1,nsys), i=1,nsys)
	write(*,fmt=form) ((a(i,j), j=1,nsys), i=1,nsys)
	write(*,*)
	write(*,*) 'constant vector: '
	write(21,fmt=form) (b(i), i=1,nsys)
	write(*,fmt=form) (b(i), i=1,nsys)
	write(*,*)

      return
      end
c======================================================
      subroutine no_piv_ge(aa,bb,n,np,x)
c======================================================
! Perform Gaussian elimination with no pivoting to solve aa*x = bb
! Matrix aa is physically np by np but only n by n is used (n <= np)
      implicit real*8(a-h,o-z)
      parameter( nmax = 100 )
      dimension aa(np,np),bb(np),x(np)
      dimension a(nmax,nmax), b(nmax)

!! initialize x
      do i=1,np
        x(i) = 0.0
      end do
!! Make copies of a(,) and b()
      do i=1,np
        b(i) = bb(i)
        do j=1,np
          a(i,j) = aa(i,j)
        end do
      end do
!! Forward elimination !!
      do k=1,n-1
       do i=k+1,n
        xmult = a(i,k)/a(k,k)
        do j=k+1,n
          a(i,j) = a(i,j) - xmult*a(k,j)
        end do
        a(i,k) = xmult
        b(i) = b(i) - xmult*b(k)
       end do
      end do
!! Backsubstitution !!
      x(n) = b(n)/a(n,n)
      do i=n-1,1,-1
        sum = b(i)
        do j=i+1,n
          sum = sum - a(i,j)*x(j)
        end do
        x(i) = sum/a(i,i)
      end do
      return
      end

c======================================================
      subroutine part_piv_ge(aa,bb,n,np,x)
c======================================================
! Perform Gaussian elimination with partial pivoting to solve aa*x =bb
! Matrix aa is physically np by np but only n by n is used (n <= np)
      implicit real*8(a-h,o-z)
      parameter( nmax = 100 )
      dimension aa(np,np),bb(np),x(np)
      dimension a(nmax,nmax), b(nmax)
      dimension idx(nmax), scale(nmax)

!! initialize x
      do i=1,np
        x(i) = 0.0
      end do
!! Make copies of a(,) and b()
      do i=1,np
        b(i) = bb(i)
        do j=1,np
          a(i,j) = aa(i,j)
        end do
      end do
!!!!! Forward elimination !!!!!
      do k=1,n-1
        big=abs(a(i,i))
        l=i
        i1=i+1
        do jj=i1,n
           if(abs(a(jj,i)).gt.big) then
              big=abs(a(jj,i))
              l=jj
           end if
        end do
        if(l.ne.i) then
          do jj=1,n
            a(l,jj)=a(i,jj)
          end do
          b(l)=b(i)
        end if
        do i=k+1,n
          xmult = a(i,k)/a(k,k)
          do j=k+1,n
            a(i,j) = a(i,j) - xmult*a(k,j)
          end do
          a(i,k) = xmult
          b(i) = b(i) - xmult*b(k)
        end do
      end do
!!!! Back substitution !!!
      x(n) = b(n)/a(n,n)
      do i=n-1,1,-1
        sum = b(i)
        do j=i+1,n
          sum = sum - a(i,j)*x(j)
        end do
        x(i) = sum/a(i,i)
      end do
          
      return
      end
c======================================================
      subroutine full_piv_ge(aa,bb,n,np,x)
c======================================================
! Perform Gaussian elimination with full pivoting to solve aa*x = bb
! Matrix aa is physically np by np but only n by n is used (n <= np)
      implicit real*8(a-h,o-z)
      parameter( nmax = 100 )
      dimension aa(np,np),bb(np),x(np)
      dimension a(nmax,nmax), b(nmax)
      dimension idx(nmax), scale(nmax)

!! initialize x
      do i=1,np
        x(i) = 0.0
      end do
!! Make copies of a(,) and b()
      do i=1,np
        b(i) = bb(i)
        do j=1,np
          a(i,j) = aa(i,j)
        end do
      end do
!!!!! Forward elimination !!!!!
      do i=1,n
        idx(i) = i
        scalemax = 0.
        do j=1,n
          scalemax = amax1(scalemax,abs(a(i,j)))
        end do
        scale(i) = scalemax
      end do
      do k=1,n-1
        ratiomax = 0.
        do i=k,n
          ratio = abs(a(idx(i),k))/scale(idx(i))
          if( ratio .gt. ratiomax ) then
            j=i
            ratiomax = ratio
          end if
        end do
        idxk = idx(j)
        idx(j) = idx(k)
        idx(k) = idxk
        do i=k+1,n
          xmult = a(idx(i),k)/a(idxk,k)
          do j=k+1,n
            a(idx(i),j) = a(idx(i),j) - xmult*a(idxk,j)
          end do
          a(idx(i),k) = xmult
          b(idx(i)) = b(idx(i)) - a(idx(i),k)*b(idxk)
        end do
      end do
!!!! Back substitution !!!
      x(n) = b(idx(n))/a(idx(n),n)
      do i=n-1,1,-1
        sum = b(idx(i))
        do j=i+1,n
          sum = sum - a(idx(i),j)*x(j)
        end do
        x(i) = sum/a(idx(i),i)
      end do
          
      return
      end