C------------------------------------------------------------------------------
C       Example parallel program to show the use of the "pzgsmp" routine.
C------------------------------------------------------------------------------
C.. This program can be obtained from:
C
C   http://www.cs.umn.edu/~agupta/wsmp
C
C   (C) IBM Corporation, 2004, 2013.
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 pzgsmp_ex2
        implicit none
        include 'mpif.h'
        integer mperr, myid

        call mpi_init_(mperr)
        if (mperr .ne. 0) then
          print *,'MPI initialization error'
          stop
        end if
        call mpi_comm_rank_(MPI_COMM_WORLD,myid,mperr)

C..  This is a small matrix, so use only 1 thread per process.
c    NOTE:  For real problems, use wsetmaxthrds only as needed.

	call wsetmaxthrds(1)

C..  This program must be called on *at least* 3 nodes.

        if (myid .eq. 0) then
          call master
        else
          call slave (myid)
        end if

        call mpi_finalize_(mperr)
        stop
        end
C------------------------------------------------------------------------------
        subroutine master
        implicit none

        integer n, ldb, nrhs
        integer iparm(64)
        integer ia(10)
        integer ja(32)
        double precision dparm(64)
        complex*16 avals(32)
        complex*16 b(9)
        complex*16 rmisc          ! 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/

        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,-7.d0),(-5.d0,2.5d0),(-1.d0,0.5d0),(-6.d0,3.0d0),
     2 (14.d0,-7.d0),(-1.d0,0.5d0),(-3.d0,1.5d0),(-1.d0,0.5d0),
     3 (-1.d0,0.5d0),(16.d0,-8.d0),(-2.d0,1.0d0),
     4 (14.d0,-7.d0),(-3.d0,1.5d0),
     5 (14.d0,-7.d0),(-1.d0,0.5d0),(-1.d0,0.5d0),
     6 (-2.d0,1.0d0),(-1.d0,0.5d0),(16.d0,-8.d0),(-2.d0,1.0d0),
     +                                                  (-4.d0,2.0d0),
     7 (-1.d0,0.5d0),(16.d0,-8.d0),
     8 (-3.d0,1.5d0),(-4.d0,2.0d0),(-3.d0,1.5d0),(-4.d0,2.0d0),
     +                                                  (71.d0,-35.5d0),
     9 (-1.d0,0.5d0),(-2.d0,1.0d0),(-4.d0,2.0d0),(16.d0,-8.d0)/

C.. Executable part of the program starts here.

        waltime = wsmprtc()

C.. Fill 'iparm' and 'dparm' arrays with default values.

C.. As an alternative to this step, the values in 'iparm' and 'dparm' can be
c   filled with values suitable for the application either manually or by 
c   using a "data" statement according to the description in the User's guide.

        iparm(1) = 0
        iparm(2) = 0
        iparm(3) = 0
        call pzgsmp(n, ia, ja, avals, b, ldb, nrhs, rmisc, iparm, dparm)      

        print *,'Initialization complete in time - ',wsmprtc()-waltime
        if (iparm(64) .ne. 0) then
          print *,'The following ERROR was detected: ',iparm(64)
          return
        end if

C.. This is a tiny problem, so use small granularity.
c   NOTE:  Do not reset dparm(25) to a small value for real problems.

        dparm(25) = 1.d0 

C.. Analysis.

        waltime = wsmprtc()
        iparm(2) = 1
        iparm(3) = 1
        call pzgsmp(n, ia, ja, avals, b, ldb, nrhs, rmisc, iparm, dparm)

        print *,'Analysis complete in time - ',wsmprtc()-waltime
        if (iparm(64) .ne. 0) then
          print *,'The following ERROR was detected: ',iparm(64)
          return
        end if
        print *,'Number of expected nonzeros in factors = 1000 X ',
     +		 iparm(24)
        print *,'Number of expected FLOPS in factorization = ',dparm(24)

C.. LU Factorization.

        waltime = wsmprtc()
        iparm(2) = 2
        iparm(3) = 2
        call pzgsmp(n, ia, ja, avals, b, ldb, nrhs, rmisc, iparm, dparm)

        waltime = wsmprtc() - waltime
        print *,'LU complete in time - ',waltime
        print *,'Actual number of nonzeros in factors = 1000 X ',
     +		 iparm(23)
        print *,'Actual number of FLOPS in factorization = ',dparm(23)
        print *,'Factorization MegaFlops = ',
     +           (dparm(23) * 1.d-6) / waltime
        if (iparm(64) .ne. 0) then
          print *,'The following ERROR was detected: ',iparm(64)
	  return
	end if

C.. Back substitution.

        do i = 1, n
          b(i) = (1.d0,-0.5d0)
        end do
        waltime = wsmprtc()
        iparm(2) = 3
        iparm(3) = 3
        call pzgsmp(n, ia, ja, avals, b, ldb, nrhs, rmisc, iparm, dparm)

        print *,'Back substitution done in time - ',wsmprtc()-waltime
        if (iparm(64) .ne. 0) then
          print *,'The following ERROR was detected: ',iparm(64)
          return
        end if

C.. Iterative refinement.

        waltime = wsmprtc()
        iparm(2) = 4
        iparm(3) = 4
        call pzgsmp(n, ia, ja, avals, b, ldb, nrhs, rmisc, iparm, dparm)

        print *,'Iterative ref. done in time - ',wsmprtc()-waltime      
        print *,'Norm of residual - ',dparm(7)
        if (iparm(64) .ne. 0) then
          print *,'The following ERROR was detected: ',iparm(64)
          return
        end if
        print *,'The solution of the system is as follows:'
        do i = 1, n
          print *,i,' : ',b(i)
        end do

        return
        end
C------------------------------------------------------------------------------
        subroutine slave (myid)
        implicit none
        integer myid

        integer n, ldb, nrhs
        integer iparm(64)
        integer ia, ja 
        double precision dparm(64)
        complex*16 avals, b       ! just placeholders in this routine
        complex*16 rmisc          ! just a placeholder in this program

        data n /0/, ldb /1/, nrhs /1/, ia /0/

C.. Fill 'iparm' and 'dparm' arrays with default values.

        iparm(1) = 0
        iparm(2) = 0
        iparm(3) = 0
        call pzgsmp(n, ia, ja, avals, b, ldb, nrhs, rmisc, iparm, dparm)

        if (iparm(64) .ne. 0) return

C.. This is a tiny problem, so use small granularity.
c   NOTE:  Do not reset dparm(25) to a small value for real problems.

        dparm(25) = 1.d0 

C.. Analysis.

        iparm(2) = 1
        iparm(3) = 1
        call pzgsmp(n, ia, ja, avals, b, ldb, nrhs, rmisc, iparm, dparm)

        if (iparm(64) .ne. 0) return

C.. LU Factorization.

        iparm(2) = 2
        iparm(3) = 2
        call pzgsmp(n, ia, ja, avals, b, ldb, nrhs, rmisc, iparm, dparm)

        if (iparm(64) .ne. 0) return

C.. Back substitution.

        iparm(2) = 3
        iparm(3) = 3
        call pzgsmp(n, ia, ja, avals, b, ldb, nrhs, rmisc, iparm, dparm)

        if (iparm(64) .ne. 0) return

C.. Iterative refinement.

        iparm(2) = 4
        iparm(3) = 4
        call pzgsmp(n, ia, ja, avals, b, ldb, nrhs, rmisc, iparm, dparm)

        if (iparm(64) .ne. 0) return

        return
        end
C------------------------------------------------------------------------------


