C------------------------------------------------------------------------------
C       Example program to show the use of the "zismp" routine for
c       solving Hermitian and symmetric complex systems using the 
c       CSR-LT/CSC-UT format.
C------------------------------------------------------------------------------
C.. This program can be obtained from:
C
C   http://www.cs.umn.edu/~agupta/wsmp
C
C   (C) IBM Corporation, 2012.
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 wssmp_ex1
        implicit none

        integer n, ldb, ldx, ldy, nrhs
        integer iparm(64)
        integer ia(10)
        integer ja(22)
        double precision dparm(64)
        double complex avals(22)
        double complex b(9), x(9), y(9), z(9)
        double complex derr
        double complex rmisc, cvgh      ! just placeholders in this program

        integer i
        double precision waltime, wsmprtc

C.. Fill all arrays containing matrix data.

        data n /9/, ldb /9/, ldx /9/, ldy /9/, nrhs /1/

        data ia /1,2,3,6,8,10,13,17,20,23/

        data ja 
     1        /1,          
     2	             2,    
     3         1,    2,    3,
     4                     3,    4,
     5               2,                5,
     6                           4,    5,    6, 
     7                     3,    4,          6,    7,      
     8         1,                4,                      8,      
     9               2,                5,                      9/
        data avals 
     1       /(14.d0,0.d0),
     2        (15.d0,0.d0), 
     3        (-1.d0,1.d0), (-1.d0,2.d0), (16.d0,0.d0),
     4        (-2.d0,3.d0), (17.d0,0.d0), 
     5        (-3.d0,2.d0), (13.d0,0.d0),
     6        (-1.d0,2.d0), (2.d0,-1.d0), (18.d0,0.d0),
     7        (2.5d0,2.d0), (-2.d0,1.d0), (1.8d0,0.d0), (16.d0,0.d0),
     8        (2.d0,-2.d0), (-1.d0,1.d0), (71.d0,0.d0),
     9        (-2.d0,2.d0), (1.d0,-2.d0), (15.d0,0.d0)/

	data z /(1.1d0,1.9d0), (1.2d0,1.8d0), (1.3d0,1.7d0), 
     1          (1.4d0,1.6d0), (1.5d0,1.5d0), (1.6d0,1.4d0),
     2          (1.7d0,1.3d0), (1.8d0,1.2d0), (1.9d0,1.1d0)/

C.. Executable part of the program starts here.

        waltime = wsmprtc()

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

        call wsetmaxthrds(2)	! Use 2 threads
        iparm(1) = 0
        iparm(2) = 0
        iparm(3) = 0
        call zismp (n, ia, ja, avals, b, ldb, x, ldx, nrhs,
     +              rmisc, cvgh, iparm, dparm)

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

C.. Modify some defaults.

	iparm(4) = 4	! CSR-LT/CSC-UT format
	iparm(7) = 3	! Hermitian matrix

C.. Analyze structure.

        waltime = wsmprtc()
        iparm(2) = 1
        iparm(3) = 1
        call zismp (n, ia, ja, avals, b, ldb, x, ldx, nrhs,
     +              rmisc, cvgh, iparm, dparm)

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

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

C.. Analyze Values.

        waltime = wsmprtc()
        iparm(2) = 2
        iparm(3) = 2
        call zismp (n, ia, ja, avals, b, ldb, x, ldx, nrhs,
     +              rmisc, cvgh, iparm, dparm)

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

C.. Matrix-vector multiplication (compute b = Ax).

        do i = 1, n
          x(i) = z(i)
        end do

        waltime = wsmprtc()
        iparm(2) = 5
        iparm(3) = 5
        call zismp (n, ia, ja, avals, b, ldb, x, ldx, nrhs,
     +              rmisc, cvgh, iparm, dparm)

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

C.. Generate Preconditioner.

        waltime = wsmprtc()
        iparm(2) = 3
        iparm(3) = 3
        call zismp (n, ia, ja, avals, b, ldb, x, ldx, nrhs,
     +              rmisc, cvgh, iparm, dparm)

        waltime = wsmprtc() - waltime
        print *,'Preconditioner generated in time - ',waltime
        if (iparm(64) .ne. 0) then
          print *,'The following ERROR was detected: ',iparm(64)
          stop
        end if

C.. Iterative Solution (solve for y in Ay = b).

        waltime = wsmprtc()
        iparm(2) = 4
        iparm(3) = 4
        call zismp (n, ia, ja, avals, b, ldb, y, ldy, nrhs,
     +              rmisc, cvgh, iparm, dparm)

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

        print *,'Solution computed in',iparm(26),'iterations:'
        rmisc = (0.d0, 0.d0)
        derr = (0.d0, 0.d0)
        do i = 1, n
          print *,i,' : ',y(i)
          derr = derr + (y(i) - x(i)) * (y(i) - x(i))
	  rmisc = rmisc + b(i) * b(i)
        end do
        print *,'Relative error = ',derr / rmisc

C.. Now solve the same sysmtem assuming symmetric non-Hermitian matrix

C.. Modify some defaults.

        iparm(7) = 2	! Symmetric non-Hermitian matrix
        iparm(15) = 2	! SSOR preconditioning

C.. Analyze structure.

        waltime = wsmprtc()
        iparm(2) = 1
        iparm(3) = 1
        call zismp (n, ia, ja, avals, b, ldb, x, ldx, nrhs,
     +              rmisc, cvgh, iparm, dparm)

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

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

C.. Analyze Values.

        waltime = wsmprtc()
        iparm(2) = 2
        iparm(3) = 2
        call zismp (n, ia, ja, avals, b, ldb, x, ldx, nrhs,
     +              rmisc, cvgh, iparm, dparm)

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

C.. Matrix-vector multiplication (compute b = Ax).

        do i = 1, n
          x(i) = z(i)
        end do

        waltime = wsmprtc()
        iparm(2) = 5
        iparm(3) = 5
        call zismp (n, ia, ja, avals, b, ldb, x, ldx, nrhs,
     +              rmisc, cvgh, iparm, dparm)

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

C.. Generate Preconditioner.

        waltime = wsmprtc()
        iparm(2) = 3
        iparm(3) = 3
        call zismp (n, ia, ja, avals, b, ldb, x, ldx, nrhs,
     +              rmisc, cvgh, iparm, dparm)

        waltime = wsmprtc() - waltime
        print *,'Preconditioner generated in time - ',waltime
        if (iparm(64) .ne. 0) then
          print *,'The following ERROR was detected: ',iparm(64)
          stop
        end if

C.. Iterative Solution (solve for y in Ay = b).

        waltime = wsmprtc()
        iparm(2) = 4
        iparm(3) = 4
        call zismp (n, ia, ja, avals, b, ldb, y, ldy, nrhs,
     +              rmisc, cvgh, iparm, dparm)

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

        print *,'Solution computed in',iparm(26),'iterations:'
        rmisc = 0.d0
        derr = 0.d0
        do i = 1, n
          print *,i,' : ',y(i)
          derr = derr + (y(i) - x(i)) * (y(i) - x(i))
          rmisc = rmisc + b(i) * b(i)
        end do
        print *,'Relative error = ',derr / rmisc

        stop
        end

