      program main
* This program compares the performance of the three Fortran oriented
* methodes to do LU factorization, JIK-SDOT, JKI_GAXPY, and KJI-SAXPY.
* Programmer A. Gaber Mohamed, Ph.D., Jan. 7, 1991
* NPAC, Syracuse University
      double precision a(:,:),x(:),b(:)
      double precision bb(:,:)
      integer ipvt(:)
      double precision emax
      integer nb
      character*1 trans
      real etime, tt1(2), tt2(2), t1, t2, toh
      real mflops1,mflops2,mflops3
c 
      print *, 'Enter order of progression of matrix dimension :'
      read *, nn
      print *, 'Enter Block Size :'
      read *, nb
      trans = 'n'
      nr = 1
c
      do 10 i=1,3
         t1 = etime(tt1)
        toh = etime(tt2) -t1
  10  continue
      print *, ' Overhead Time = ', toh
      print *, ' Leading Dimension = ', nn+1

       do 200 n = 100,nn,100
           lda = n+1
           ldb = lda
* Allocate arrays
         allocate (a(lda,n), x(n), b(n), bb(ldb,1), ipvt(n))
         write(6,50)lda,n
   50    format(/' lda size of the arrays',i5,' and order is ',i5/)
* JIK-SDOT
         call matgen(a,lda,n,b)
	 bb(1:n,1)=b(1:n)
         t1 = etime(tt1)
         t1 = etime(tt1)
         call jik_block(n,n,a,lda,ipvt,ierr,nb)
         call dgetrs('N',n,nr,a,lda,ipvt,bb,ldb,ierr)
         t2 = etime(tt2) -t1 -toh
         mflops1 =(n*(5.0+n*(4.0*n-3.0)))/(6.0*t2*1.0e6)
         emax = 0.0d0
         do 110 i = 1,n
            emax = dmax1(emax,abs(bb(i,1)))
  110    continue
         write(6,121) t2,mflops1,emax
  121    format(' Using JIK-SDOT and FX/BLAS time = ',
     $   e12.3/'     MFLOPS = ',f9.4,' check = ',e12.3)
* JKI-GAXPY
         call matgen(a,lda,n,b)
	 bb(1:n,1)=b(1:n)
         t1 = etime(tt1)
         t1 = etime(tt1)
         call jki_block(n,n,a,lda,ipvt,ierr,nb)
         call dgetrs('N',n,nr,a,lda,ipvt,bb,ldb,ierr)
         t2 = etime(tt2) -t1 -toh
         mflops2 =(n*(5.0+n*(4.0*n-3.0)))/(6.0*t2*1.0e6)
         emax = 0.0d0
         do 112 i = 1,n
            emax = dmax1(emax,abs(bb(i,1)))
  112    continue
         write(6,122) t2,mflops2,emax
  122    format(' Using JKI-GAXPY and FX/BLAS time = ',
     $   e12.3/'     MFLOPS = ',f9.4,' check = ',e12.3)
* KJI-SAXPY
         call matgen(a,lda,n,b)
	 bb(1:n,1)=b(1:n)
         t1 = etime(tt1)
         t1 = etime(tt1)
         call kji_block(n,n,a,lda,ipvt,ierr,nb)
         call dgetrs('N',n,nr,a,lda,ipvt,bb,ldb,ierr)
         t2 = etime(tt2) -t1 -toh
         mflops3 =(n*(5.0+n*(4.0*n-3.0)))/(6.0*t2*1.0e6)
         emax = 0.0d0
         do 113 i = 1,n
            emax = dmax1(emax,abs(bb(i,1)))
  113    continue
         write(6,123) t2,mflops3,emax
  123    format(' Using KJI-SAXPY and FX/BLAS time = ',
     $   e12.3/'     MFLOPS = ',f9.4,' check = ',e12.3)
         write(1,*) n,mflops1,mflops2,mflops3
* Deallocate arrays
      deallocate (a, x, b, bb, ipvt)
  200 continue
      stop
      end

      subroutine matgen(a,lda,n,b)
      double precision a(lda,n),b(n)
c
      init = 1325
      do 30 j = 1,n
         do 20 i = 1,n
            init = mod(3125*init,65536)
            a(i,j) = (init - 32768.0d0)/16384.0d0
   20    continue
   30 continue
      do 35 i = 1,n
          b(i) = 0.0d0
   35 continue
      do 50 j = 1,n
         do 40 i = 1,n
            b(i) = b(i) + a(i,j)
   40    continue
   50 continue
      return
      end

