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 90 sparse linear algebraic equation subroutines and their utilities to solve this example.
The following is the general sparse matrix A:
* *
| 2.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 |
| 0.0 2.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 |
| 0.0 1.0 2.0 0.0 0.0 0.0 0.0 0.0 0.0 |
| 1.0 0.0 0.0 2.0 -1.0 0.0 0.0 0.0 0.0 |
| 0.0 0.0 0.0 1.0 2.0 -1.0 0.0 0.0 0.0 |
| 0.0 0.0 0.0 0.0 1.0 2.0 -1.0 0.0 0.0 |
| 0.0 0.0 0.0 0.0 0.0 1.0 2.0 -1.0 0.0 |
| 0.0 0.0 0.0 0.0 0.0 0.0 1.0 2.0 -1.0 |
| 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 2.0 |
* *
The following is the dense vector b, containing the right-hand side:
* *
| 2.0 |
| 1.0 |
| 3.0 |
| 2.0 |
| 2.0 |
| 2.0 |
| 2.0 |
| 2.0 |
| 3.0 |
* *
The following is the dense vector x, containing the initial guess to the solution:
* *
| 0.0 |
| 0.0 |
| 0.0 |
| 0.0 |
| 0.0 |
| 0.0 |
| 0.0 |
| 0.0 |
| 0.0 |
* *
Global vector x:
B,D 0
* *
| 1.0 |
0 | 1.0 |
| 1.0 |
| --- |
| 1.0 |
1 | 1.0 |
| 1.0 |
| --- |
| 1.0 |
2 | 1.0 |
| 1.0 |
* *
The following is the 3 × 1 process grid:
B,D | 0 -----| ------- 0 | P00 -----| ------- 1 | P10 -----| ------- 2 | P20
Local vector x:
p,q | 0
------|------
| 1.0
0 | 1.0
| 1.0
------|------
| 1.0
1 | 1.0
| 1.0
------|------
| 1.0
2 | 1.0
| 1.0
ITER = 4
ERR = 0.4071D-15
The value of info is 0 on all processes.
This application program illustrates how to use the Fortran 90 sparse linear algebraic equation subroutines and their utilities.
@process init(f90ptr)
!
! This program illustrates how to use the PESSL F90 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 EXAMPLE90
! Interface module required to use the PESSL F90 Sparse Iterative Solver
USE F90SPARSE
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
! Parameters
CHARACTER, PARAMETER :: ORDER='R'
INTEGER, PARAMETER :: IZERO=0, IONE=1
! Sparse Matrices
TYPE(D_SPMAT) :: A, BLCK
! Preconditioner Data Structure
TYPE(D_PRECN) :: PRC
! Dense Vectors
REAL(KIND(1.D0)), POINTER :: B(:), X(:)
! Communications data structure
TYPE(DESC_TYPE) :: DESC_A
! BLACS parameters
INTEGER :: NPROW, NPCOL, ICTXT, IAM, NP, MYROW, MYCOL
! Solver parameters
INTEGER :: ITER, ITMAX, IERR, ITRACE,
& IPREC, METHD, ISTOPC, IPARM(20)
REAL(KIND(1.D0)) :: ERR, EPS, RPARM(20)
! Other variables
CHARACTER*5 :: AFMT, ATYPE
INTEGER :: IRCODE, IRCODE1, IRCODE2, IRCODE3
INTEGER :: I,J
INTEGER :: N,NNZERO
INTEGER, POINTER :: PV(:)
INTEGER :: LPROCS, NROW, NCOL
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(22),B_GLOBAL(9),XINIT_GLOBAL(9)
INTEGER :: JA(22),IA(10)
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)
!
! Initialize the global problem size
!
N = SIZE(IA)-1
!
! Guess for the local number of nonzeros
!
NNZERO = SIZE(A_GLOBAL)
!
! Allocate and initialize some elements of the sparse matrix A
! its descriptor vector, DESC_A, the rhs vector B, and the
! solution vector X.
!
CALL PADALL(N,PART_BLOCK,DESC_A,ICTXT)
CALL PSPALL(A,DESC_A,NNZ=NNZERO)
CALL PGEALL(B,DESC_A)
CALL PGEALL(X,DESC_A)
!
! Allocate an integer work area to be used as an argument for
! the PART_BLOCK PARTS subroutine
!
NROW = N
NCOL = NROW
LPROCS = MAX(NPROW, NROW + NCOL)
ALLOCATE(PV(LPROCS), STAT = IRCODE)
IF (IRCODE /= 0) THEN
WRITE(6,*) 'PV Allocation failed'
CALL BLACS_ABORT(ICTXT,-1)
STOP
ENDIF
! SETUP BLCK
BLCK%M = 1
BLCK%N = NCOL
BLCK%FIDA = 'CSR'
ALLOCATE(BLCK%AS(BLCK%N),STAT=IRCODE1)
ALLOCATE(BLCK%IA1(BLCK%N),STAT=IRCODE2)
ALLOCATE(BLCK%IA2(BLCK%M+1),STAT=IRCODE3)
IRCODE = IRCODE1 + IRCODE2 + IRCODE3
IF (IRCODE /= 0) THEN
WRITE(6,*) 'Error allocating BLCK'
CALL BLACS_ABORT(ICTXT,-1)
STOP
ENDIF
!
! 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 PSPINS as many times as necessary
! to insert the local rows it owns.
!
! Each process will call PGEINS as many times as necessary
! to insert the local elements it owns.
!
DO GLOBAL_INDX = 1, NROW
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
BLCK%IA2(1) = 1
BLCK%IA2(2) = 1
DO J = IA(GLOBAL_INDX), IA(GLOBAL_INDX+1)-1
BLCK%AS(BLCK%IA2(2)) = A_GLOBAL(J)
BLCK%IA1(BLCK%IA2(2)) = JA(J)
BLCK%IA2(2) =BLCK%IA2(2) + 1
ENDDO
CALL PSPINS(A,GLOBAL_INDX,1,BLCK,DESC_A)
CALL PGEINS(B,B_GLOBAL(GLOBAL_INDX:GLOBAL_INDX),
& DESC_A,GLOBAL_INDX)
CALL PGEINS(X,XINIT_GLOBAL(GLOBAL_INDX:GLOBAL_INDX),
& DESC_A,GLOBAL_INDX)
ENDIF
END DO
END DO
! Assemble A and DESC_A
AFMT = 'DEF'
ATYPE = 'GEN'
CALL PSPASB(A,DESC_A,MTYPE=ATYPE,
& STOR=AFMT,DUPFLAG=2,INFO=IERR)
IF (IERR /= 0) THEN
IF (IAM.EQ.0) THEN
WRITE(6,*) 'Error in assembly :',IERR
CALL BLACS_ABORT(ICTXT,-1)
STOP
END IF
END IF
! Assemble B and X
CALL PGEASB(B,DESC_A)
CALL PGEASB(X,DESC_A)
!
! Deallocate BLCK
!
IF (ASSOCIATED(BLCK%AS)) DEALLOCATE(BLCK%AS)
IF (ASSOCIATED(BLCK%IA1)) DEALLOCATE(BLCK%IA1)
IF (ASSOCIATED(BLCK%IA2)) DEALLOCATE(BLCK%IA2)
!
! Deallocate Work vector
!
IF (ASSOCIATED(PV)) DEALLOCATE(PV)
!
! Preconditioning
!
! We are using ILU for the preconditioner; PESSL
! will allocate PRC.
!
IPREC = 2
CALL PSPGPR(IPREC,A,PRC,DESC_A,INFO=IERR)
IF (IERR /= 0) THEN
IF (IAM.EQ.0) THEN
WRITE(6,*) 'Error in preconditioner :',IERR
CALL BLACS_ABORT(ICTXT,-1)
STOP
END IF
END IF
!
! Iterative Solver - use the BICGSTAB method
!
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 PSPGIS(A,B,X,PRC,DESC_A,IPARM=IPARM,RPARM=RPARM,
& INFO=IERR)
IF (IERR /= 0) THEN
IF (IAM.EQ.0) THEN
WRITE(6,*) 'Error in solver :',IERR
CALL BLACS_ABORT(ICTXT,-1)
STOP
END IF
END IF
ITER = IPARM(5)
ERR = RPARM(2)
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, NROW
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 vectors, the sparse matrix, and
! the preconditioner data structure.
! Finally, deallocate the descriptor vector
!
CALL PGEFREE(B, DESC_A)
CALL PGEFREE(X, DESC_A)
CALL PSPFREE(A, DESC_A)
CALL PSPFREE(PRC, DESC_A)
CALL PADFREE(DESC_A)
!
! Terminate the process grid and the BLACS
!
CALL BLACS_GRIDEXIT(ICTXT)
CALL BLACS_EXIT(0)
END PROGRAM EXAMPLE90