C------------------------------------------------------------------------------
C       Example program to show the use of the "wismp" routine for
c       solving a symmetric positive definite system using the default
c       preconditioner and two RHS vectors.
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(49)
        double precision dparm(64)
        double precision avals(49)
        double precision b(9,2), x(9,2), y(9,2)
        double precision derr
        double precision 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 /2/

        data ia /1,5,9,15,19,23,29,35,44,50/

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

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 wismp (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(7) = 3	! symmetric positive definite matrix

C.. Analyze structure.

        waltime = wsmprtc()
        iparm(2) = 1
        iparm(3) = 1
        call wismp (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 wismp (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,1) = 1.d0 + 1.d-1 * dfloat(i-(i/10)*10)
        end do
        do i = 1, n
	  x(i,2) = x(n-i+1,1)
	end do

        waltime = wsmprtc()
        iparm(2) = 5
        iparm(3) = 5
        call wismp (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 wismp (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 wismp (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,1),y(i,2)
          derr = derr + (y(i,1) - x(i,1)) * (y(i,1) - x(i,1))
          derr = derr + (y(i,2) - x(i,2)) * (y(i,2) - x(i,2))
	  rmisc = rmisc + b(i,1) * b(i,1)
	  rmisc = rmisc + b(i,2) * b(i,2)
        end do
        print *,'Relative error = ',derr / rmisc

        stop
        end

