|
ScaLAPACK
2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
|
|
ScaLAPACK
2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
|
00001 PROGRAM PCSCAEX 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 PCGESV. 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 CPLXSZ, INTGSZ, MEMSIZ, TOTMEM 00015 PARAMETER ( CPLXSZ = 8, INTGSZ = 4, TOTMEM = 2000000, 00016 $ MEMSIZ = TOTMEM / CPLXSZ ) 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 COMPLEX ONE 00023 PARAMETER ( ONE = (1.0D+0,0.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 REAL ANORM, BNORM, EPS, XNORM, RESID 00031 * .. 00032 * .. Local Arrays .. 00033 INTEGER DESCA( DLEN_ ), DESCB( DLEN_ ), DESCX( DLEN_ ) 00034 COMPLEX MEM( MEMSIZ ) 00035 * .. 00036 * .. External Subroutines .. 00037 EXTERNAL BLACS_EXIT, BLACS_GET, BLACS_GRIDEXIT, 00038 $ BLACS_GRIDINFO, BLACS_GRIDINIT, BLACS_PINFO, 00039 $ DESCINIT, IGSUM2D, PCGESV, 00040 $ PCGEMM, PCLACPY, PCLAPRNT, PCLAREAD, PCLAWRITE 00041 * .. 00042 * .. External Functions .. 00043 INTEGER ICEIL, NUMROC 00044 REAL PSLAMCH, PCLANGE 00045 EXTERNAL ICEIL, NUMROC, PSLAMCH, PCLANGE 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 ), CPLXSZ ) 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 )*CPLXSZ 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 PCLAREAD( 'CSCAEXMAT.dat', MEM( IPA ), DESCA, 0, 0, 00117 $ MEM( IPW ) ) 00118 CALL PCLAREAD( 'CSCAEXRHS.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 PCLACPY( 'All', N, N, MEM( IPA ), 1, 1, DESCA, 00124 $ MEM( IPACPY ), 1, 1, DESCA ) 00125 CALL PCLACPY( 'All', N, NRHS, MEM( IPB ), 1, 1, DESCB, 00126 $ MEM( IPX ), 1, 1, DESCX ) 00127 * 00128 ********************************************************************** 00129 * Call ScaLAPACK PCGESV 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: (PCGESV)' 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 PCLAPRNT( 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 PCLAPRNT( N, NRHS, MEM( IPB ), 1, 1, DESCB, 0, 0, 00152 $ 'B', NOUT, MEM( IPW ) ) 00153 * 00154 CALL PCGESV( 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 PCGESV = ', INFO 00160 WRITE( NOUT, FMT = * ) 00161 WRITE( NOUT, FMT = * ) 'Matrix X = A^{-1} * B' 00162 WRITE( NOUT, FMT = * ) 00163 END IF 00164 CALL PCLAPRNT( N, NRHS, MEM( IPB ), 1, 1, DESCB, 0, 0, 'X', NOUT, 00165 $ MEM( IPW ) ) 00166 CALL PCLAWRITE( 'CSCAEXSOL.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 = PSLAMCH( ICTXT, 'Epsilon' ) 00172 ANORM = PCLANGE( 'I', N, N, MEM( IPA ), 1, 1, DESCA, MEM( IPW ) ) 00173 BNORM = PCLANGE( 'I', N, NRHS, MEM( IPB ), 1, 1, DESCB, 00174 $ MEM( IPW ) ) 00175 CALL PCGEMM( '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 = PCLANGE( '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 PCSCAEX 00221 * 00222 END
