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