ScaLAPACK  2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
pdscaex.f
Go to the documentation of this file.
00001       PROGRAM PDSCAEX
00002 *
00003 *  -- ScaLAPACK example code --
00004 *     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
00005 *     and University of California, Berkeley.
00006 *
00007 *     Written by Antoine Petitet, August 1995 (petitet@cs.utk.edu)
00008 *
00009 *     This program solves a linear system by calling the ScaLAPACK 
00010 *     routine PDGESV. The input matrix and right-and-sides are
00011 *     read from a file. The solution is written to a file.
00012 *
00013 *     .. Parameters ..
00014       INTEGER            DBLESZ, INTGSZ, MEMSIZ, TOTMEM
00015       PARAMETER          ( DBLESZ = 8, INTGSZ = 4, TOTMEM = 2000000,
00016      $                     MEMSIZ = TOTMEM / DBLESZ )
00017       INTEGER            BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DT_,
00018      $                   LLD_, MB_, M_, NB_, N_, RSRC_
00019       PARAMETER          ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DT_ = 1,
00020      $                     CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6,
00021      $                     RSRC_ = 7, CSRC_ = 8, LLD_ = 9 )
00022       DOUBLE PRECISION   ONE
00023       PARAMETER          ( ONE = 1.0D+0 )
00024 *     ..
00025 *     .. Local Scalars ..
00026       CHARACTER*80       OUTFILE
00027       INTEGER            IAM, ICTXT, INFO, IPA, IPACPY, IPB, IPPIV, IPX,
00028      $                   IPW, LIPIV, MYCOL, MYROW, N, NB, NOUT, NPCOL,
00029      $                   NPROCS, NPROW, NP, NQ, NQRHS, NRHS, WORKSIZ
00030       DOUBLE PRECISION   ANORM, BNORM, EPS, XNORM, RESID
00031 *     ..
00032 *     .. Local Arrays ..
00033       INTEGER            DESCA( DLEN_ ), DESCB( DLEN_ ), DESCX( DLEN_ )
00034       DOUBLE PRECISION   MEM( MEMSIZ )
00035 *     ..
00036 *     .. External Subroutines ..
00037       EXTERNAL           BLACS_EXIT, BLACS_GET, BLACS_GRIDEXIT,
00038      $                   BLACS_GRIDINFO, BLACS_GRIDINIT, BLACS_PINFO,
00039      $                   DESCINIT, IGSUM2D, PDSCAEXINFO, PDGESV,
00040      $                   PDGEMM, PDLACPY, PDLAPRNT, PDLAREAD, PDLAWRITE
00041 *     ..
00042 *     .. External Functions ..
00043       INTEGER            ICEIL, NUMROC
00044       DOUBLE PRECISION   PDLAMCH, PDLANGE
00045       EXTERNAL           ICEIL, NUMROC, PDLAMCH, PDLANGE
00046 *     ..
00047 *     .. Intrinsic Functions ..
00048       INTRINSIC          DBLE, MAX
00049 *     ..
00050 *     .. Executable Statements ..
00051 *
00052 *     Get starting information
00053 *
00054       CALL BLACS_PINFO( IAM, NPROCS )
00055       CALL PDSCAEXINFO( OUTFILE, NOUT, N, NRHS, NB, NPROW, NPCOL, MEM,
00056      $                  IAM, NPROCS )
00057 *
00058 *     Define process grid
00059 *
00060       CALL BLACS_GET( -1, 0, ICTXT )
00061       CALL BLACS_GRIDINIT( ICTXT, 'Row-major', NPROW, NPCOL )
00062       CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )
00063 *
00064 *     Go to bottom of process grid loop if this case doesn't use my
00065 *     process
00066 *
00067       IF( MYROW.GE.NPROW .OR. MYCOL.GE.NPCOL )
00068      $   GO TO 20
00069 *
00070       NP    = NUMROC( N, NB, MYROW, 0, NPROW )
00071       NQ    = NUMROC( N, NB, MYCOL, 0, NPCOL )
00072       NQRHS = NUMROC( NRHS, NB, MYCOL, 0, NPCOL )
00073 *
00074 *     Initialize the array descriptor for the matrix A and B
00075 *
00076       CALL DESCINIT( DESCA, N, N, NB, NB, 0, 0, ICTXT, MAX( 1, NP ),
00077      $               INFO )
00078       CALL DESCINIT( DESCB, N, NRHS, NB, NB, 0, 0, ICTXT, MAX( 1, NP ),
00079      $               INFO )
00080       CALL DESCINIT( DESCX, N, NRHS, NB, NB, 0, 0, ICTXT, MAX( 1, NP ),
00081      $               INFO )
00082 *
00083 *     Assign pointers into MEM for SCALAPACK arrays, A is
00084 *     allocated starting at position MEM( 1 )
00085 *
00086       IPA = 1
00087       IPACPY = IPA + DESCA( LLD_ )*NQ
00088       IPB = IPACPY + DESCA( LLD_ )*NQ
00089       IPX = IPB + DESCB( LLD_ )*NQRHS
00090       IPPIV = IPX + DESCB( LLD_ )*NQRHS
00091       LIPIV = ICEIL( INTGSZ*( NP+NB ), DBLESZ )
00092       IPW = IPPIV + MAX( NP, LIPIV )
00093 *
00094       WORKSIZ = NB
00095 *
00096 *     Check for adequate memory for problem size
00097 *
00098       INFO = 0
00099       IF( IPW+WORKSIZ.GT.MEMSIZ ) THEN
00100          IF( IAM.EQ.0 )
00101      $      WRITE( NOUT, FMT = 9998 ) 'test', ( IPW+WORKSIZ )*DBLESZ
00102          INFO = 1
00103       END IF
00104 *
00105 *     Check all processes for an error
00106 *
00107       CALL IGSUM2D( ICTXT, 'All', ' ', 1, 1, INFO, 1, -1, 0 )
00108       IF( INFO.GT.0 ) THEN
00109          IF( IAM.EQ.0 )
00110      $      WRITE( NOUT, FMT = 9999 ) 'MEMORY'
00111          GO TO 10
00112       END IF
00113 *
00114 *     Read from file and distribute matrices A and B
00115 *
00116       CALL PDLAREAD( 'DSCAEXMAT.dat', MEM( IPA ), DESCA, 0, 0,
00117      $               MEM( IPW ) )
00118       CALL PDLAREAD( 'DSCAEXRHS.dat', MEM( IPB ), DESCB, 0, 0,
00119      $               MEM( IPW ) )
00120 *
00121 *     Make a copy of A and the rhs for checking purposes
00122 *
00123       CALL PDLACPY( 'All', N, N, MEM( IPA ), 1, 1, DESCA,
00124      $              MEM( IPACPY ), 1, 1, DESCA )
00125       CALL PDLACPY( 'All', N, NRHS, MEM( IPB ), 1, 1, DESCB,
00126      $              MEM( IPX ), 1, 1, DESCX )
00127 *
00128 **********************************************************************
00129 *     Call ScaLAPACK PDGESV routine
00130 **********************************************************************
00131 *
00132       IF( IAM.EQ.0 ) THEN
00133          WRITE( NOUT, FMT = * )
00134          WRITE( NOUT, FMT = * )
00135      $         '***********************************************'
00136          WRITE( NOUT, FMT = * )
00137      $         'Example of ScaLAPACK routine call: (PDGESV)'
00138          WRITE( NOUT, FMT = * )
00139      $         '***********************************************'
00140          WRITE( NOUT, FMT = * )
00141          WRITE( NOUT, FMT = * ) 'A * X = B, Matrix A:'
00142          WRITE( NOUT, FMT = * )
00143       END IF
00144       CALL PDLAPRNT( N, N, MEM( IPA ), 1, 1, DESCA, 0, 0,
00145      $               'A', NOUT, MEM( IPW ) )
00146       IF( IAM.EQ.0 ) THEN
00147          WRITE( NOUT, FMT = * )
00148          WRITE( NOUT, FMT = * ) 'Matrix B:'
00149          WRITE( NOUT, FMT = * )
00150       END IF
00151       CALL PDLAPRNT( N, NRHS, MEM( IPB ), 1, 1, DESCB, 0, 0,
00152      $               'B', NOUT, MEM( IPW ) )
00153 *
00154       CALL PDGESV( N, NRHS, MEM( IPA ), 1, 1, DESCA, MEM( IPPIV ),
00155      $             MEM( IPB ), 1, 1, DESCB, INFO )
00156 *
00157       IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
00158          WRITE( NOUT, FMT = * )
00159          WRITE( NOUT, FMT = * ) 'INFO code returned by PDGESV = ', INFO
00160          WRITE( NOUT, FMT = * )
00161          WRITE( NOUT, FMT = * ) 'Matrix X = A^{-1} * B'
00162          WRITE( NOUT, FMT = * )
00163       END IF
00164       CALL PDLAPRNT( N, NRHS, MEM( IPB ), 1, 1, DESCB, 0, 0, 'X', NOUT,
00165      $               MEM( IPW ) )
00166       CALL PDLAWRITE( 'DSCAEXSOL.dat', N, NRHS, MEM( IPB ), 1, 1, DESCB,
00167      $                0, 0, MEM( IPW ) )
00168 *
00169 *     Compute residual ||A * X  - B|| / ( ||X|| * ||A|| * eps * N )
00170 *
00171       EPS = PDLAMCH( ICTXT, 'Epsilon' )
00172       ANORM = PDLANGE( 'I', N, N, MEM( IPA ), 1, 1, DESCA, MEM( IPW ) )
00173       BNORM = PDLANGE( 'I', N, NRHS, MEM( IPB ), 1, 1, DESCB,
00174      $                 MEM( IPW ) )
00175       CALL PDGEMM( 'No transpose', 'No transpose', N, NRHS, N, ONE,
00176      $             MEM( IPACPY ), 1, 1, DESCA, MEM( IPB ), 1, 1, DESCB,
00177      $             -ONE, MEM( IPX ), 1, 1, DESCX )
00178       XNORM = PDLANGE( 'I', N, NRHS, MEM( IPX ), 1, 1, DESCX,
00179      $                 MEM( IPW ) )
00180       RESID = XNORM / ( ANORM * BNORM * EPS * DBLE( N ) )
00181 *
00182       IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
00183          WRITE( NOUT, FMT = * )
00184          WRITE( NOUT, FMT = * )
00185      $     '||A * X  - B|| / ( ||X|| * ||A|| * eps * N ) = ', RESID
00186          WRITE( NOUT, FMT = * )
00187          IF( RESID.LT.10.0D+0 ) THEN
00188             WRITE( NOUT, FMT = * ) 'The answer is correct.'
00189          ELSE
00190             WRITE( NOUT, FMT = * ) 'The answer is suspicious.'
00191          END IF
00192       END IF
00193 *
00194    10 CONTINUE
00195 *
00196       CALL BLACS_GRIDEXIT( ICTXT )
00197 *
00198    20 CONTINUE
00199 *
00200 *     Print ending messages and close output file
00201 *
00202       IF( IAM.EQ.0 ) THEN
00203          WRITE( NOUT, FMT = * )
00204          WRITE( NOUT, FMT = * )
00205          WRITE( NOUT, FMT = 9997 )
00206          WRITE( NOUT, FMT = * )
00207          IF( NOUT.NE.6 .AND. NOUT.NE.0 )
00208      $      CLOSE ( NOUT )
00209       END IF
00210 *
00211       CALL BLACS_EXIT( 0 )
00212 *
00213  9999 FORMAT( 'Bad ', A6, ' parameters: going on to next test case.' )
00214  9998 FORMAT( 'Unable to perform ', A, ': need TOTMEM of at least',
00215      $        I11 )
00216  9997 FORMAT( 'END OF TESTS.' )
00217 *
00218       STOP
00219 *
00220 *     End of PDSCAEX
00221 *
00222       END