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