IBM Books

Parallel Engineering and Scientific Subroutine Library for AIX Version 2 Release 3: Guide and Reference

Example--Using the Fortran 77 Sparse Subroutines

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.

Application Program

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


[ Top of Page | Previous Page | Next Page | Table of Contents | Index ]