C------------------------------------------------------------------------------
C       Example program to show the use of the simple interface for
c       solving general sparse systems in the WSMP library.
C------------------------------------------------------------------------------
C.. This program can be obtained from:
C
C   http://www.cs.umn.edu/~agupta/wsmp
C
C   (C) IBM Corporation, 2000.
C
c       Acceptance and use of this program constitutes the user's understanding
c       that he/she will have no recourse to IBM for any actual or consequential
c       damages, including, but not limited to, lost profits or savings, arising
c       out of the use or inability to use this program.
C------------------------------------------------------------------------------
        program wgsmp_ex2
        implicit none

        integer n, ldb, nrhs, niter, nnz, info
        integer ia(10)
        integer ja(32)
        double precision avals(32)
        double precision b(9)
        double precision berr, thresh, opc
        double precision darray          ! just a placeholder in this program

        integer i
        double precision waltime, wsmprtc

C.. Fill all arrays containing matrix data.

        data n /9/, ldb /9/, nrhs /1/, niter /3/
        data thresh /5.d-2/
        data ia /1,5,9,12,14,17,22,24,29,33/

        data ja
     1        /1,          3,                      7,    8,
     2               2,    3,          5,                      9,
     3         1,          3,                      7,            
     4                           4,                      8,
     5                                 5,    6,                9,
     6               2,          4,          6,    7,    8,      
     7                           4,                7,      
     8         1,          3,          5,          7,    8,      
     9               2,    3,                            8,    9/
        data avals
     1      /14.d0,      -5.d0,                  -1.d0,-6.d0,
     2             14.d0,-1.d0,      -3.d0,                  -1.d0,
     3       -1.d0,      16.d0,                  -2.d0,
     4                         14.d0,                  -3.d0,
     5                               14.d0,-1.d0,            -1.d0,
     6             -2.d0,      -1.d0,      16.d0,-2.d0,-4.d0,
     7                         -1.d0,            16.d0,
     8       -3.d0,      -4.d0,      -3.d0,      -4.d0,71.d0,
     9             -1.d0,-2.d0,                        -4.d0,16.d0/

C.. Analysis.

      waltime = wsmprtc()
      call wgralz (n, ia, ja, avals, nnz, opc, info)      

      print *,'Analysis complete in time - ',wsmprtc()-waltime

      if (info .ne. 0) then
        print *,'The following ERROR was detected: ',info
        stop
      end if

      print *,'Number of nonzeros in LU factors = ',nnz
      print *,'Number of FLOPS in factorization = ',opc

C.. Factorization.

      waltime = wsmprtc()
      call wgrluf (n, ia, ja, avals, thresh, info)

      waltime = wsmprtc() - waltime
      print *,'Factorization complete in time - ',waltime
      if (info .ne. 0) then
        print *,'The following ERROR was detected: ',info
        stop
      end if
      print *,'Factorization MegaFlops = ',(opc*1.d-6)/waltime

C.. Back substitution with iterative refinement.

      do i = 1, n
        b(i) = 1.d0
      end do

      waltime = wsmprtc()
      call wgrslv (n, ia, ja, avals, b, ldb, nrhs, niter, berr, info)

      print *,'Back substitution complete in time - ',wsmprtc()-waltime
      if (info .ne. 0) then
        print *,'The following ERROR was detected: ',info
        stop
      end if

      print *,'Maximum relative error = ',berr
      print *,'The solution of the system is as follows:'
      do i = 1, n
        print *,i,' : ',b(i)
      end do

      stop
      end