!
! Example Program solving Ax=b via ScaLAPACK routine PDGESV
!
! .. Parameters ..
INTEGER DLEN_, IA, JA, IB, JB, M, N, MB, NB, RSRC, &
CSRC, MXLLDA, MXLLDB, NRHS, NBRHS, NOUT, &
MXLOCR, MXLOCC, MXRHSC
PARAMETER ( DLEN_ = 9, IA = 1, JA = 1, IB = 1, JB = 1, &
M = 9, N = 9, MB = 2, NB = 2, RSRC = 0, &
CSRC = 0, MXLLDA = 5, MXLLDB = 5, NRHS = 1, &
NBRHS = 1, NOUT = 6, MXLOCR = 5, MXLOCC = 4, &
MXRHSC = 1 )
DOUBLE PRECISION ONE
PARAMETER ( ONE = 1.0D+0 )
! ..
! .. Local Scalars ..
INTEGER ICTXT, INFO, MYCOL, MYROW, NPCOL, NPROW
DOUBLE PRECISION ANORM, BNORM, EPS, RESID, XNORM
! ..
! .. Local Arrays ..
INTEGER DESCA( DLEN_ ), DESCB( DLEN_ ), &
IPIV( MXLOCR+NB )
DOUBLE PRECISION A( MXLLDA, MXLOCC ), A0( MXLLDA, MXLOCC ), &
B( MXLLDB, MXRHSC ), B0( MXLLDB, MXRHSC ), &
WORK( MXLOCR )
! ..
! .. External Functions ..
DOUBLE PRECISION PDLAMCH, PDLANGE
EXTERNAL PDLAMCH, PDLANGE
! ..
! .. External Subroutines ..
EXTERNAL BLACS_EXIT, BLACS_GRIDEXIT, BLACS_GRIDINFO, &
DESCINIT, MATINIT, PDGEMM, PDGESV, PDLACPY, &
SL_INIT
! ..
! .. Intrinsic Functions ..
INTRINSIC REAL
! ..
! .. Data statements ..
DATA NPROW / 2 / , NPCOL / 3 /
! ..
! .. Executable Statements ..
!
! INITIALIZE THE PROCESS GRID
!
CALL SL_INIT( ICTXT, NPROW, NPCOL )
CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )
!
! If I'm not in the process grid, go to the end of the program
!
IF( MYROW.EQ.-1 ) &
GO TO 10
!
! DISTRIBUTE THE MATRIX ON THE PROCESS GRID
! Initialize the array descriptors for the matrices A and B
!
CALL DESCINIT( DESCA, M, N, MB, NB, RSRC, CSRC, ICTXT, MXLLDA, &
INFO )
CALL DESCINIT( DESCB, N, NRHS, NB, NBRHS, RSRC, CSRC, ICTXT, &
MXLLDB, INFO )
!
! Generate matrices A and B and distribute to the process grid
!
CALL MATINIT( A, DESCA, B, DESCB )
!
! Make a copy of A and B for checking purposes
!
CALL PDLACPY( 'All', N, N, A, 1, 1, DESCA, A0, 1, 1, DESCA )
CALL PDLACPY( 'All', N, NRHS, B, 1, 1, DESCB, B0, 1, 1, DESCB )
!
! CALL THE SCALAPACK ROUTINE
! Solve the linear system A * X = B
!
CALL PDGESV( N, NRHS, A, IA, JA, DESCA, IPIV, B, IB, JB, DESCB, &
INFO )
!
IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
WRITE( NOUT, FMT = 9999 )
WRITE( NOUT, FMT = 9998 )M, N, NB
WRITE( NOUT, FMT = 9997 )NPROW*NPCOL, NPROW, NPCOL
WRITE( NOUT, FMT = 9996 )INFO
END IF
!
! Compute residual ||A * X - B|| / ( ||X|| * ||A|| * eps * N )
!
EPS = PDLAMCH( ICTXT, 'Epsilon' )
ANORM = PDLANGE( 'I', N, N, A, 1, 1, DESCA, WORK )
BNORM = PDLANGE( 'I', N, NRHS, B, 1, 1, DESCB, WORK )
CALL PDGEMM( 'N', 'N', N, NRHS, N, ONE, A0, 1, 1, DESCA, B, 1, 1, &
DESCB, -ONE, B0, 1, 1, DESCB )
XNORM = PDLANGE( 'I', N, NRHS, B0, 1, 1, DESCB, WORK )
RESID = XNORM / ( ANORM*BNORM*EPS*REAL( N ) )
!
IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
IF( RESID.LT.10.0E+0 ) THEN
WRITE( NOUT, FMT = 9995 )
WRITE( NOUT, FMT = 9993 )RESID
ELSE
WRITE( NOUT, FMT = 9994 )
WRITE( NOUT, FMT = 9993 )RESID
END IF
END IF
!
! RELEASE THE PROCESS GRID
! Free the BLACS context
!
CALL BLACS_GRIDEXIT( ICTXT )
10 CONTINUE
!
! Exit the BLACS
!
CALL BLACS_EXIT( 0 )
!
9999 FORMAT( / 'ScaLAPACK Example Program #1 -- May 1, 1997' )
9998 FORMAT( / 'Solving Ax=b where A is a ', I3, ' by ', I3, &
' matrix with a block size of ', I3 )
9997 FORMAT( 'Running on ', I3, ' processes, where the process grid', &
' is ', I3, ' by ', I3 )
9996 FORMAT( / 'INFO code returned by PDGESV = ', I3 )
9995 FORMAT( / &
'According to the normalized residual the solution is correct.' &
)
9994 FORMAT( / &
'According to the normalized residual the solution is incorrect.' &
)
9993 FORMAT( / '||A*x - b|| / ( ||x||*||A||*eps*N ) = ', 1P, E16.8 )
STOP
END
!
! MATINIT generates and distributes matrices A and B (depicted in
! Figures 2.5 and 2.6) to a 2 x 3 process grid
!
! .. Array Arguments ..
INTEGER DESCA( * ), DESCB( * )
DOUBLE PRECISION AA( * ), B( * )
! ..
! .. Parameters ..
INTEGER CTXT_, LLD_
PARAMETER ( CTXT_ = 2, LLD_ = 9 )
! ..
! .. Local Scalars ..
INTEGER ICTXT, MXLLDA, MYCOL, MYROW, NPCOL, NPROW
DOUBLE PRECISION A, C, K, L, P, S
! ..
! .. External Subroutines ..
EXTERNAL BLACS_GRIDINFO
! ..
! .. Executable Statements ..
!
ICTXT = DESCA( CTXT_ )
CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )
!
S = 19.0D0
C = 3.0D0
A = 1.0D0
L = 12.0D0
P = 16.0D0
K = 11.0D0
!
MXLLDA = DESCA( LLD_ )
!
IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
AA( 1 ) = S
AA( 2 ) = -S
AA( 3 ) = -S
AA( 4 ) = -S
AA( 5 ) = -S
AA( 1+MXLLDA ) = C
AA( 2+MXLLDA ) = C
AA( 3+MXLLDA ) = -C
AA( 4+MXLLDA ) = -C
AA( 5+MXLLDA ) = -C
AA( 1+2*MXLLDA ) = A
AA( 2+2*MXLLDA ) = A
AA( 3+2*MXLLDA ) = A
AA( 4+2*MXLLDA ) = A
AA( 5+2*MXLLDA ) = -A
AA( 1+3*MXLLDA ) = C
AA( 2+3*MXLLDA ) = C
AA( 3+3*MXLLDA ) = C
AA( 4+3*MXLLDA ) = C
AA( 5+3*MXLLDA ) = -C
B( 1 ) = 0.0D0
B( 2 ) = 0.0D0
B( 3 ) = 0.0D0
B( 4 ) = 0.0D0
B( 5 ) = 0.0D0
ELSE IF( MYROW.EQ.0 .AND. MYCOL.EQ.1 ) THEN
AA( 1 ) = A
AA( 2 ) = A
AA( 3 ) = -A
AA( 4 ) = -A
AA( 5 ) = -A
AA( 1+MXLLDA ) = L
AA( 2+MXLLDA ) = L
AA( 3+MXLLDA ) = -L
AA( 4+MXLLDA ) = -L
AA( 5+MXLLDA ) = -L
AA( 1+2*MXLLDA ) = K
AA( 2+2*MXLLDA ) = K
AA( 3+2*MXLLDA ) = K
AA( 4+2*MXLLDA ) = K
AA( 5+2*MXLLDA ) = K
ELSE IF( MYROW.EQ.0 .AND. MYCOL.EQ.2 ) THEN
AA( 1 ) = A
AA( 2 ) = A
AA( 3 ) = A
AA( 4 ) = -A
AA( 5 ) = -A
AA( 1+MXLLDA ) = P
AA( 2+MXLLDA ) = P
AA( 3+MXLLDA ) = P
AA( 4+MXLLDA ) = P
AA( 5+MXLLDA ) = -P
ELSE IF( MYROW.EQ.1 .AND. MYCOL.EQ.0 ) THEN
AA( 1 ) = -S
AA( 2 ) = -S
AA( 3 ) = -S
AA( 4 ) = -S
AA( 1+MXLLDA ) = -C
AA( 2+MXLLDA ) = -C
AA( 3+MXLLDA ) = -C
AA( 4+MXLLDA ) = C
AA( 1+2*MXLLDA ) = A
AA( 2+2*MXLLDA ) = A
AA( 3+2*MXLLDA ) = A
AA( 4+2*MXLLDA ) = -A
AA( 1+3*MXLLDA ) = C
AA( 2+3*MXLLDA ) = C
AA( 3+3*MXLLDA ) = C
AA( 4+3*MXLLDA ) = C
B( 1 ) = 1.0D0
B( 2 ) = 0.0D0
B( 3 ) = 0.0D0
B( 4 ) = 0.0D0
ELSE IF( MYROW.EQ.1 .AND. MYCOL.EQ.1 ) THEN
AA( 1 ) = A
AA( 2 ) = -A
AA( 3 ) = -A
AA( 4 ) = -A
AA( 1+MXLLDA ) = L
AA( 2+MXLLDA ) = L
AA( 3+MXLLDA ) = -L
AA( 4+MXLLDA ) = -L
AA( 1+2*MXLLDA ) = K
AA( 2+2*MXLLDA ) = K
AA( 3+2*MXLLDA ) = K
AA( 4+2*MXLLDA ) = K
ELSE IF( MYROW.EQ.1 .AND. MYCOL.EQ.2 ) THEN
AA( 1 ) = A
AA( 2 ) = A
AA( 3 ) = -A
AA( 4 ) = -A
AA( 1+MXLLDA ) = P
AA( 2+MXLLDA ) = P
AA( 3+MXLLDA ) = -P
AA( 4+MXLLDA ) = -P
END IF
RETURN
END
!
! .. Scalar Arguments ..
INTEGER ICTXT, NPCOL, NPROW
! ..
!
! Purpose
! =======
!
! SL_INIT initializes an NPROW x NPCOL process grid using a row-major
! ordering of the processes. This routine retrieves a default system
! context which will include all available processes. In addition it
! spawns the processes if needed.
!
! Arguments
! =========
!
! ICTXT (global output) INTEGER
! ICTXT specifies the BLACS context handle identifying the
! created process grid. The context itself is global.
!
! NPROW (global input) INTEGER
! NPROW specifies the number of process rows in the grid
! to be created.
!
! NPCOL (global input) INTEGER
! NPCOL specifies the number of process columns in the grid
! to be created.
!
! =====================================================================
!
! .. Local Scalars ..
INTEGER IAM, NPROCS
! ..
! .. External Subroutines ..
EXTERNAL BLACS_GET, BLACS_GRIDINIT, BLACS_PINFO, &
BLACS_SETUP
! ..
! .. Executable Statements ..
!
! Get starting information
!
CALL BLACS_PINFO( IAM, NPROCS )
!
! If machine needs additional set up, do it now
!
IF( NPROCS.LT.1 ) THEN
IF( IAM.EQ.0 ) &
NPROCS = NPROW*NPCOL
CALL BLACS_SETUP( IAM, NPROCS )
END IF
!
! Define process grid
!
CALL BLACS_GET( -1, 0, ICTXT )
CALL BLACS_GRIDINIT( ICTXT, 'Row-major', NPROW, NPCOL )
!
RETURN
!
! End of SL_INIT
!
END
HTML generated by f90tohtml