PROGRAM PDPOSVXEXAMPLE
*
* Example Program solving Ax=b via ScaLAPACK routine PDPOSV
*
* The global symmetric positive definite matrix is
* __ __
* | 9 0 0 1 1 |
* | 0 9 0 2 -2 |
* A = | 0 0 9 1 -1 |
* | 1 2 1 9 0 |
* | 1 -2 -1 0 9 |
* -- --
* The righthand side is
*
* __ __
* | 11 |
* | 9 |
* b = | 9 |
* | 13 |
* | 7 |
* -- --
* In this program, only upper triangular part of this matrix is
* stored in 2 by 2 process grid with 2 by 2 block.
* Each process stores elements as follows
*
* --- P0 --- --- P1 ---
* A = 9 0 1 b = 11 A = 0 1
* 0 9 -2 9 0 2
* 1 -2 9 7 -1 9
*
* --- P3 --- --- P4 ---
* A = 1 0 -1 b = 9 A = 9 1
* 0 2 0 13 1 9
*
* Further Details
* ===============
*
* Contributed by Keita Teranishi, University of Tennessee, 1996.
*
* =====================================================================
*
* .. Parameters ..
INTEGER DLEN_, IA, JA, IB, JB, N, MB, NB, RSRC, CSRC,
$ MXLLDA, MXLLDB, NRHS, NBRHS, NOUT, LWORK,
$ MXLOCC, MXRHSC
PARAMETER ( DLEN_ = 9, IA = 1, JA = 1, IB = 1, JB = 1,
$ N = 5, MB = 2, NB = 2, RSRC = 0, CSRC = 0,
$ MXLLDA = 3, MXLLDB = 3, NRHS = 1, NBRHS = 1,
$ NOUT = 6, LWORK = 20, MXLOCC = 3, 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_ )
DOUBLE PRECISION A( MXLLDA, MXLOCC ), A0( MXLLDA, MXLOCC ),
$ B( MXLLDB, MXRHSC ), B0( MXLLDB, MXRHSC ),
$ WORK( LWORK )
* ..
* .. External Functions ..
DOUBLE PRECISION PDLAMCH, PDLANGE, PDLANSY
EXTERNAL PDLAMCH, PDLANGE, PDLANSY
* ..
* .. External Subroutines ..
EXTERNAL BLACS_EXIT, BLACS_GRIDEXIT, BLACS_GRIDINFO,
$ DESCINIT, MATINIT, PDLACPY, PDLAPRNT, PDPOSV,
$ PDSYMM, SL_INIT
* ..
* .. Intrinsic Functions ..
INTRINSIC DBLE
* ..
* .. Data statements ..
DATA NPROW / 2 / , NPCOL / 2 /
* ..
* .. 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, N, 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 )
*
IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
WRITE( NOUT, FMT = 9992 )
END IF
*
* ... Print out the elements of b ...
*
CALL PDLAPRNT( N, NRHS, B, IB, JB, DESCB, RSRC, CSRC, 'B', NOUT,
$ WORK )
*
* CALL THE SCALAPACK ROUTINE
* Solve the linear system A * X = B
*
CALL PDPOSV( 'U', N, NRHS, A, IA, JA, DESCA, B, IB, JB, DESCB,
$ INFO )
*
IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
WRITE( NOUT, FMT = 9999 )
WRITE( NOUT, FMT = 9998 )N, N, NB
WRITE( NOUT, FMT = 9997 )NPROW*NPCOL, NPROW, NPCOL
WRITE( NOUT, FMT = 9996 )INFO
END IF
IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
WRITE( NOUT, FMT = 9991 )
END IF
*
* ... Print out the elements of x ...
*
CALL PDLAPRNT( N, NRHS, B, IB, JB, DESCB, RSRC, CSRC, 'X', NOUT,
$ WORK )
*
* Compute residual ||A * X - B|| / ( ||X|| * ||A|| * eps * N )
*
EPS = PDLAMCH( ICTXT, 'Epsilon' )
ANORM = PDLANSY( 'I', 'U', N, A0, 1, 1, DESCA, WORK )
BNORM = PDLANGE( 'I', N, NRHS, B, 1, 1, DESCB, WORK )
CALL PDSYMM( 'L', 'U', N, NRHS, 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*DBLE( N ) )
*
IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
IF( RESID.LT.10.0D+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 PDPOSV = ', 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 )
9992 FORMAT( / 'The elements of b ' )
9991 FORMAT( / 'The elements of x' )
STOP
END
*
SUBROUTINE MATINIT( AA, DESCA, B, DESCB )
*
* 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
* ..
* .. External Subroutines ..
EXTERNAL BLACS_GRIDINFO
* ..
* .. Executable Statements ..
*
ICTXT = DESCA( CTXT_ )
CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )
*
*
MXLLDA = DESCA( LLD_ )
*
IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
AA( 1 ) = 9.0D0
AA( 2 ) = 0.0D0
AA( 3 ) = 1.0D0
AA( 1+MXLLDA ) = 0.0D0
AA( 2+MXLLDA ) = 9.0D0
AA( 3+MXLLDA ) = -2.0D0
AA( 1+2*MXLLDA ) = 1.0D0
AA( 2+2*MXLLDA ) = -2.0D0
AA( 3+2*MXLLDA ) = 9.0D0
B( 1 ) = 11.0D0
B( 2 ) = 9.0D0
B( 3 ) = 7.0D0
ELSE IF( MYROW.EQ.0 .AND. MYCOL.EQ.1 ) THEN
AA( 1 ) = 0.0D0
AA( 2 ) = 0.0D0
AA( 3 ) = -1.0D0
AA( 1+MXLLDA ) = 1.0D0
AA( 2+MXLLDA ) = 2.0D0
AA( 3+MXLLDA ) = 0.0D0
ELSE IF( MYROW.EQ.1 .AND. MYCOL.EQ.0 ) THEN
AA( 1 ) = 0.0D0
AA( 2 ) = 1.0D0
AA( 1+MXLLDA ) = 0.0D0
AA( 2+MXLLDA ) = 2.0D0
AA( 1+2*MXLLDA ) = -1.0D0
AA( 2+2*MXLLDA ) = 0.0D0
B( 1 ) = 9.0D0
B( 2 ) = 13.0D0
ELSE IF( MYROW.EQ.1 .AND. MYCOL.EQ.1 ) THEN
AA( 1 ) = 9.0D0
AA( 2 ) = 1.0D0
AA( 1+MXLLDA ) = 1.0D0
AA( 2+MXLLDA ) = 9.0D0
END IF
RETURN
END