This example finds the solution to the linear system Ax = b. It also contains an application program that shows how you can use the Fortran 77 sparse linear algebraic equation subroutines and their utilities to solve the problem shown in Example--Using the Fortran 90 Sparse Subroutines.
This application program illustrates how to use the Fortran 77 sparse linear algebraic equation subroutines and their utilities.
!
! This program illustrates how to use the PESSL F77 Sparse Iterative
! Solver and its supporting utility subroutines. A very simple problem
! (DSRIS Example 1 from the ESSL Guide and Reference) using an
! HPF BLOCK data distribution is solved.
!
PROGRAM EXAMPLE77
IMPLICIT NONE
! Interface definition for the PARTS subroutine PART_BLOCK
INTERFACE PART_BLOCK
SUBROUTINE PART_BLOCK(GLOBAL_INDX,N,NP,PV,NV)
IMPLICIT NONE
INTEGER, INTENT(IN) :: GLOBAL_INDX, N, NP
INTEGER, INTENT(OUT) :: NV
INTEGER, INTENT(OUT) :: PV(*)
END SUBROUTINE PART_BLOCK
END INTERFACE
! External declaration for the PARTS subroutine PART_BLOCK
EXTERNAL PART_BLOCK
! Parameters
CHARACTER*1 ORDER
CHARACTER*5 STOR
CHARACTER*5 MTYPE
INTEGER*4 IZERO, IONE, DUPFLAG, N, NNZ
PARAMETER (ORDER='R')
PARAMETER (STOR='DEF')
PARAMETER (MTYPE='GEN')
PARAMETER (IZERO=0)
PARAMETER (IONE=1)
PARAMETER (N=9)
PARAMETER (NNZ=22)
PARAMETER (DUPFLAG=2)
! Descriptor Vector
INTEGER*4, ALLOCATABLE :: DESC_A(:)
! Sparse Matrices and related information
REAL*8 AS(NNZ)
INTEGER*4 IA1(NNZ+N), IA2(NNZ+N)
INTEGER*4 INFOA(30)
REAL*8 BS(NNZ)
INTEGER*4 IB1(N+1), IB2(NNZ)
INTEGER*4 INFOB(30)
! Preconditioner Data Structure
REAL*8 PRCS(2*NNZ+2*N+41)
! Dense Vectors
REAL*8 B(N), X(N)
! BLACS parameters
INTEGER*4 NPROW, NPCOL, ICTXT, IAM, NP, MYROW, MYCOL
! Solver parameters
INTEGER*4 ITER, ITMAX, INFO, ITRACE,
& IPREC, METHD, ISTOPC, IPARM(20)
REAL*8 ERR, EPS, RPARM(20)
! We will not have duplicates so PV used by the PARTS subroutine
! PART_BLOCK only needs to be of length 1.
INTEGER PV(1)
! Other variables
INTEGER IERR
INTEGER NB, LDB, LDBG
INTEGER NX, LDX, LDXG
INTEGER NRHS
INTEGER I,J
INTEGER GLOBAL_INDX, NV_COUNT
INTEGER GLOBAL_INDX_OWNER, NV
INTEGER LOCAL_INDX
!
! Global Problem
! DSRIS Example 1 from the ESSL Guide and Reference
!
REAL*8 A_GLOBAL(NNZ),B_GLOBAL(N),XINIT_GLOBAL(N)
INTEGER JA(NNZ),IA(N+1)
DATA A_GLOBAL /2.D0,2.D0,-1.D0,1.D0,2.D0,1.D0,2.D0,-1.D0,
$ 1.D0,2.D0,-1.D0,1.D0,2.D0,-1.D0,1.D0,2.D0,
$ -1.D0,1.D0,2.D0,-1.D0,1.D0,2.D0/
DATA JA /1,2,3,2,3,1,4,5,4,5,6,5,6,7,6,7,8,
$ 7,8,9,8,9/
DATA IA /1,2,4,6,9,12,15,18,21,23/
DATA B_GLOBAL /2.D0,1.D0,3.D0,2.D0,2.D0,2.D0,2.D0,2.D0,
$ 3.D0/
DATA XINIT_GLOBAL /0.D0,0.D0,0.D0,0.D0,0.D0,0.D0,0.D0,0.D0,
$ 0.D0/
! Initialize BLACS
! Define a NP x 1 Process Grid
CALL BLACS_PINFO(IAM, NP)
CALL BLACS_GET(IZERO, IZERO, ICTXT)
CALL BLACS_GRIDINIT(ICTXT, ORDER, NP, IONE)
CALL BLACS_GRIDINFO(ICTXT, NPROW, NPCOL, MYROW, MYCOL)
!
! Allocate the descriptor vector
!
ALLOCATE(DESC_A(30 + 3*NP + 4*N + 3),STAT=IERR)
IF (IERR .NE. 0) THEN
WRITE(6,*) 'Error allocating DESC_A :',IERR
CALL BLACS_ABORT(ICTXT,-1)
STOP
END IF
! Initialize some elements of the sparse matrix A
! and its descriptor vector, DESC_A
!
DESC_A(11) = SIZE(DESC_A)
CALL PADINIT(N,PART_BLOCK,DESC_A,ICTXT)
INFOA(1) = SIZE(AS)
INFOA(2) = SIZE(IA1)
INFOA(3) = SIZE(IA2)
CALL PDSPINIT(AS,IA1,IA2,INFOA,DESC_A)
!
! In this simple example, all processes have a copy of
! the global sparse matrix, A, the global rhs vector B,
! and the global initial guess vector, X
!
! Each process will call PDSPINS as many times as necessary
! to insert the local rows it owns.
!
! Each process will call PDGEINS as many times as necessary
! to insert the local elements it owns.
!
NB = 1
LDB = SIZE(B,1)
LDBG = SIZE(B_GLOBAL,1)
NX = 1
LDX = SIZE(X,1)
LDXG = SIZE(XINIT_GLOBAL,1)
DO GLOBAL_INDX = 1, N
CALL PART_BLOCK(GLOBAL_INDX,N,NP,PV,NV)
!
! In this simple example, NV will always be 1
! since there will not be duplicate coefficients
!
DO NV_COUNT = 1, NV
GLOBAL_INDX_OWNER = PV(NV_COUNT)
IF (GLOBAL_INDX_OWNER == MYROW) THEN
IB2(1) = 1
IB2(2) = 1
DO J = IA(GLOBAL_INDX), IA(GLOBAL_INDX+1)-1
BS(IB2(2)) = A_GLOBAL(J)
IB1(IB2(2)) = JA(J)
IB2(2) = IB2(2) + 1
ENDDO
INFOB(1) = IB2(2) - 1
INFOB(2) = IB2(2) - 1
INFOB(3) = 2
INFOB(4) = 1
INFOB(5) = 1
INFOB(6) = 1
INFOB(7) = N
CALL PDSPINS(AS,IA1,IA2,INFOA,DESC_A,GLOBAL_INDX, 1,
& BS,IB1,IB2,INFOB)
CALL PDGEINS(NB,B,LDB,GLOBAL_INDX,1,1,1,
& B_GLOBAL(GLOBAL_INDX),LDBG,DESC_A)
CALL PDGEINS(NX,X,LDX,GLOBAL_INDX,1,1,1,
& XINIT_GLOBAL(GLOBAL_INDX),LDXG,DESC_A)
ENDIF
END DO
END DO
! Assemble A and DESC_A
CALL PDSPASB(AS,IA1,IA2,INFOA,DESC_A,
& MTYPE,STOR,DUPFLAG,INFO)
IF (INFO .NE. 0) THEN
IF (IAM.EQ.0) THEN
WRITE(6,*) 'Error in assembly :',INFO
CALL BLACS_ABORT(ICTXT,-1)
STOP
END IF
END IF
! Assemble B and X
CALL PDGEASB(NB,B,LDB,DESC_A)
CALL PDGEASB(NX,X,LDX,DESC_A)
!
! Preconditioning
!
! We are using ILU for the preconditioner
!
IPREC = 2
CALL PDSPGPR(IPREC,AS,IA1,IA2,INFOA,
& PRCS,SIZE(PRCS),DESC_A,INFO)
IF (INFO .NE. 0) THEN
IF (IAM.EQ.0) THEN
WRITE(6,*) 'Error in preconditioner :',INFO
CALL BLACS_ABORT(ICTXT,-1)
STOP
END IF
END IF
!
! Iterative Solver - use the BICGSTAB method
!
NRHS = 1
ITMAX = 1000
EPS = 1.D-8
METHD = 1
ISTOPC = 1
ITRACE = 0
IPARM = 0
IPARM(1) = METHD
IPARM(2) = ISTOPC
IPARM(3) = ITMAX
IPARM(4) = ITRACE
RPARM = 0.0D0
RPARM(1) = EPS
CALL PDSPGIS(AS,IA1,IA2,INFOA,NRHS,B,LDB,X,LDX,PRCS,DESC_A,
& IPARM,RPARM,INFO)
IF (INFO .NE. 0) THEN
IF (IAM.EQ.0) THEN
WRITE(6,*) 'Error in solver :',INFO
CALL BLACS_ABORT(ICTXT,-1)
STOP
END IF
END IF
ERR = RPARM(2)
ITER = IPARM(5)
IF (IAM.EQ.0) THEN
WRITE(6,*) 'Number of iterations : ',ITER
WRITE(6,*) 'Error on exit : ',ERR
END IF
!
! Each process prints their local piece of the solution vector
!
IF (IAM.EQ.0) THEN
Write(6,*) 'Solution Vector X'
END IF
LOCAL_INDX = 1
Do GLOBAL_INDX = 1, N
CALL PART_BLOCK(GLOBAL_INDX,N,NP,PV,NV)
!
! In this simple example, NV will always be 1
! since there will not be duplicate coefficients
!
DO NV_COUNT = 1, NV
GLOBAL_INDX_OWNER = PV(NV_COUNT)
IF (GLOBAL_INDX_OWNER == MYROW) THEN
Write(6,*) GLOBAL_INDX, X(LOCAL_INDX)
LOCAL_INDX = LOCAL_INDX +1
ENDIF
END DO
END DO
!
! Deallocate the descriptor vector
!
DEALLOCATE(DESC_A, STAT=IERR)
IF (IERR .NE. 0) THEN
WRITE(6,*) 'Error deallocating DESC_A :',IERR
CALL BLACS_ABORT(ICTXT,-1)
STOP
END IF
!
! Terminate the process grid and the BLACS
!
CALL BLACS_GRIDEXIT(ICTXT)
CALL BLACS_EXIT(0)
END PROGRAM EXAMPLE77