scalapack/example1.f

   1       PROGRAM EXAMPLE1
   2 *
   3 *     Example Program solving Ax=b via ScaLAPACK routine PDGESV
   4 *
   5 *     .. Parameters ..
   6       INTEGER            DLEN_, IA, JA, IB, JB, M, N, MB, NB, RSRC,
   7      $                   CSRC, MXLLDA, MXLLDB, NRHS, NBRHS, NOUT,
   8      $                   MXLOCR, MXLOCC, MXRHSC
   9       PARAMETER          ( DLEN_ = 9, IA = 1, JA = 1, IB = 1, JB = 1,
  10      $                   M = 9, N = 9, MB = 2, NB = 2, RSRC = 0,
  11      $                   CSRC = 0, MXLLDA = 5, MXLLDB = 5, NRHS = 1,
  12      $                   NBRHS = 1, NOUT = 6, MXLOCR = 5, MXLOCC = 4,
  13      $                   MXRHSC = 1 )
  14       DOUBLE PRECISION   ONE
  15       PARAMETER          ( ONE = 1.0D+0 )
  16 *     ..
  17 *     .. Local Scalars ..
  18       INTEGER            ICTXT, INFO, MYCOL, MYROW, NPCOL, NPROW
  19       DOUBLE PRECISION   ANORM, BNORM, EPS, RESID, XNORM
  20 *     ..
  21 *     .. Local Arrays ..
  22       INTEGER            DESCA( DLEN_ ), DESCB( DLEN_ ),
  23      $                   IPIV( MXLOCR+NB )
  24       DOUBLE PRECISION   A( MXLLDA, MXLOCC ), A0( MXLLDA, MXLOCC ),
  25      $                   B( MXLLDB, MXRHSC ), B0( MXLLDB, MXRHSC ),
  26      $                   WORK( MXLOCR )
  27 *     ..
  28 *     .. External Functions ..
  29       DOUBLE PRECISION   PDLAMCH, PDLANGE
  30       EXTERNAL           PDLAMCH, PDLANGE
  31 *     ..
  32 *     .. External Subroutines ..
  33       EXTERNAL           BLACS_EXIT, BLACS_GRIDEXIT, BLACS_GRIDINFO,
  34      $                   DESCINIT, MATINIT, PDGEMM, PDGESV, PDLACPY,
  35      $                   SL_INIT
  36 *     ..
  37 *     .. Intrinsic Functions ..
  38       INTRINSIC          DBLE
  39 *     ..
  40 *     .. Data statements ..
  41       DATA               NPROW / 2 / , NPCOL / 3 /
  42 *     ..
  43 *     .. Executable Statements ..
  44 *
  45 *     INITIALIZE THE PROCESS GRID
  46 *
  47       CALL SL_INIT( ICTXT, NPROW, NPCOL )
  48       CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )
  49 *
  50 *     If I'm not in the process grid, go to the end of the program
  51 *
  52       IF( MYROW.EQ.-1 )
  53      $   GO TO 10
  54 *
  55 *     DISTRIBUTE THE MATRIX ON THE PROCESS GRID
  56 *     Initialize the array descriptors for the matrices A and B
  57 *
  58       CALL DESCINIT( DESCA, M, N, MB, NB, RSRC, CSRC, ICTXT, MXLLDA,
  59      $               INFO )
  60       CALL DESCINIT( DESCB, N, NRHS, NB, NBRHS, RSRC, CSRC, ICTXT,
  61      $               MXLLDB, INFO )
  62 *
  63 *     Generate matrices A and B and distribute to the process grid
  64 *
  65       CALL MATINIT( A, DESCA, B, DESCB )
  66 *
  67 *     Make a copy of A and B for checking purposes
  68 *
  69       CALL PDLACPY( 'All', N, N, A, 1, 1, DESCA, A0, 1, 1, DESCA )
  70       CALL PDLACPY( 'All', N, NRHS, B, 1, 1, DESCB, B0, 1, 1, DESCB )
  71 *
  72 *     CALL THE SCALAPACK ROUTINE
  73 *     Solve the linear system A * X = B
  74 *
  75       CALL PDGESV( N, NRHS, A, IA, JA, DESCA, IPIV, B, IB, JB, DESCB,
  76      $             INFO )
  77 *
  78       IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
  79          WRITE( NOUT, FMT = 9999 )
  80          WRITE( NOUT, FMT = 9998 )M, N, NB
  81          WRITE( NOUT, FMT = 9997 )NPROW*NPCOL, NPROW, NPCOL
  82          WRITE( NOUT, FMT = 9996 )INFO
  83       END IF
  84 *
  85 *     Compute residual ||A * X  - B|| / ( ||X|| * ||A|| * eps * N )
  86 *
  87       EPS = PDLAMCH( ICTXT, 'Epsilon' )
  88       ANORM = PDLANGE( 'I', N, N, A, 1, 1, DESCA, WORK )
  89       BNORM = PDLANGE( 'I', N, NRHS, B, 1, 1, DESCB, WORK )
  90       CALL PDGEMM( 'N', 'N', N, NRHS, N, ONE, A0, 1, 1, DESCA, B, 1, 1,
  91      $             DESCB, -ONE, B0, 1, 1, DESCB )
  92       XNORM = PDLANGE( 'I', N, NRHS, B0, 1, 1, DESCB, WORK )
  93       RESID = XNORM / ( ANORM*BNORM*EPS*DBLE( N ) )
  94 *
  95       IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
  96          IF( RESID.LT.10.0D+0 ) THEN
  97             WRITE( NOUT, FMT = 9995 )
  98             WRITE( NOUT, FMT = 9993 )RESID
  99          ELSE
 100             WRITE( NOUT, FMT = 9994 )
 101             WRITE( NOUT, FMT = 9993 )RESID
 102          END IF
 103       END IF
 104 *
 105 *     RELEASE THE PROCESS GRID
 106 *     Free the BLACS context
 107 *
 108       CALL BLACS_GRIDEXIT( ICTXT )
 109    10 CONTINUE
 110 *
 111 *     Exit the BLACS
 112 *
 113       CALL BLACS_EXIT( 0 )
 114 *
 115  9999 FORMAT( / 'ScaLAPACK Example Program #1 -- May 1, 1997' )
 116  9998 FORMAT( / 'Solving Ax=b where A is a ', I3, ' by ', I3,
 117      $      ' matrix with a block size of ', I3 )
 118  9997 FORMAT( 'Running on ', I3, ' processes, where the process grid',
 119      $      ' is ', I3, ' by ', I3 )
 120  9996 FORMAT( / 'INFO code returned by PDGESV = ', I3 )
 121  9995 FORMAT( /
 122      $   'According to the normalized residual the solution is correct.'
 123      $       )
 124  9994 FORMAT( /
 125      $ 'According to the normalized residual the solution is incorrect.'
 126      $       )
 127  9993 FORMAT( / '||A*x - b|| / ( ||x||*||A||*eps*N ) = ', 1P, E16.8 )
 128       STOP
 129       END
 130       SUBROUTINE MATINIT( AA, DESCA, B, DESCB )
 131 *
 132 *     MATINIT generates and distributes matrices A and B (depicted in
 133 *     Figures 2.5 and 2.6) to a 2 x 3 process grid
 134 *
 135 *     .. Array Arguments ..
 136       INTEGER            DESCA( * ), DESCB( * )
 137       DOUBLE PRECISION   AA( * ), B( * )
 138 *     ..
 139 *     .. Parameters ..
 140       INTEGER            CTXT_, LLD_
 141       PARAMETER          ( CTXT_ = 2, LLD_ = 9 )
 142 *     ..
 143 *     .. Local Scalars ..
 144       INTEGER            ICTXT, MXLLDA, MYCOL, MYROW, NPCOL, NPROW
 145       DOUBLE PRECISION   A, C, K, L, P, S
 146 *     ..
 147 *     .. External Subroutines ..
 148       EXTERNAL           BLACS_GRIDINFO
 149 *     ..
 150 *     .. Executable Statements ..
 151 *
 152       ICTXT = DESCA( CTXT_ )
 153       CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )
 154 *
 155       S = 19.0D0
 156       C = 3.0D0
 157       A = 1.0D0
 158       L = 12.0D0
 159       P = 16.0D0
 160       K = 11.0D0
 161 *
 162       MXLLDA = DESCA( LLD_ )
 163 *
 164       IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
 165          AA( 1 ) = S
 166          AA( 2 ) = -S
 167          AA( 3 ) = -S
 168          AA( 4 ) = -S
 169          AA( 5 ) = -S
 170          AA( 1+MXLLDA ) = C
 171          AA( 2+MXLLDA ) = C
 172          AA( 3+MXLLDA ) = -C
 173          AA( 4+MXLLDA ) = -C
 174          AA( 5+MXLLDA ) = -C
 175          AA( 1+2*MXLLDA ) = A
 176          AA( 2+2*MXLLDA ) = A
 177          AA( 3+2*MXLLDA ) = A
 178          AA( 4+2*MXLLDA ) = A
 179          AA( 5+2*MXLLDA ) = -A
 180          AA( 1+3*MXLLDA ) = C
 181          AA( 2+3*MXLLDA ) = C
 182          AA( 3+3*MXLLDA ) = C
 183          AA( 4+3*MXLLDA ) = C
 184          AA( 5+3*MXLLDA ) = -C
 185          B( 1 ) = 0.0D0
 186          B( 2 ) = 0.0D0
 187          B( 3 ) = 0.0D0
 188          B( 4 ) = 0.0D0
 189          B( 5 ) = 0.0D0
 190       ELSE IF( MYROW.EQ.0 .AND. MYCOL.EQ.1 ) THEN
 191          AA( 1 ) = A
 192          AA( 2 ) = A
 193          AA( 3 ) = -A
 194          AA( 4 ) = -A
 195          AA( 5 ) = -A
 196          AA( 1+MXLLDA ) = L
 197          AA( 2+MXLLDA ) = L
 198          AA( 3+MXLLDA ) = -L
 199          AA( 4+MXLLDA ) = -L
 200          AA( 5+MXLLDA ) = -L
 201          AA( 1+2*MXLLDA ) = K
 202          AA( 2+2*MXLLDA ) = K
 203          AA( 3+2*MXLLDA ) = K
 204          AA( 4+2*MXLLDA ) = K
 205          AA( 5+2*MXLLDA ) = K
 206       ELSE IF( MYROW.EQ.0 .AND. MYCOL.EQ.2 ) THEN
 207          AA( 1 ) = A
 208          AA( 2 ) = A
 209          AA( 3 ) = A
 210          AA( 4 ) = -A
 211          AA( 5 ) = -A
 212          AA( 1+MXLLDA ) = P
 213          AA( 2+MXLLDA ) = P
 214          AA( 3+MXLLDA ) = P
 215          AA( 4+MXLLDA ) = P
 216          AA( 5+MXLLDA ) = -P
 217       ELSE IF( MYROW.EQ.1 .AND. MYCOL.EQ.0 ) THEN
 218          AA( 1 ) = -S
 219          AA( 2 ) = -S
 220          AA( 3 ) = -S
 221          AA( 4 ) = -S
 222          AA( 1+MXLLDA ) = -C
 223          AA( 2+MXLLDA ) = -C
 224          AA( 3+MXLLDA ) = -C
 225          AA( 4+MXLLDA ) = C
 226          AA( 1+2*MXLLDA ) = A
 227          AA( 2+2*MXLLDA ) = A
 228          AA( 3+2*MXLLDA ) = A
 229          AA( 4+2*MXLLDA ) = -A
 230          AA( 1+3*MXLLDA ) = C
 231          AA( 2+3*MXLLDA ) = C
 232          AA( 3+3*MXLLDA ) = C
 233          AA( 4+3*MXLLDA ) = C
 234          B( 1 ) = 1.0D0
 235          B( 2 ) = 0.0D0
 236          B( 3 ) = 0.0D0
 237          B( 4 ) = 0.0D0
 238       ELSE IF( MYROW.EQ.1 .AND. MYCOL.EQ.1 ) THEN
 239          AA( 1 ) = A
 240          AA( 2 ) = -A
 241          AA( 3 ) = -A
 242          AA( 4 ) = -A
 243          AA( 1+MXLLDA ) = L
 244          AA( 2+MXLLDA ) = L
 245          AA( 3+MXLLDA ) = -L
 246          AA( 4+MXLLDA ) = -L
 247          AA( 1+2*MXLLDA ) = K
 248          AA( 2+2*MXLLDA ) = K
 249          AA( 3+2*MXLLDA ) = K
 250          AA( 4+2*MXLLDA ) = K
 251       ELSE IF( MYROW.EQ.1 .AND. MYCOL.EQ.2 ) THEN
 252          AA( 1 ) = A
 253          AA( 2 ) = A
 254          AA( 3 ) = -A
 255          AA( 4 ) = -A
 256          AA( 1+MXLLDA ) = P
 257          AA( 2+MXLLDA ) = P
 258          AA( 3+MXLLDA ) = -P
 259          AA( 4+MXLLDA ) = -P
 260       END IF
 261       RETURN
 262       END
 263       SUBROUTINE SL_INIT( ICTXT, NPROW, NPCOL )
 264 *
 265 *     .. Scalar Arguments ..
 266       INTEGER            ICTXT, NPCOL, NPROW
 267 *     ..
 268 *
 269 *  Purpose
 270 *  =======
 271 *
 272 *  SL_INIT initializes an NPROW x NPCOL process grid using a row-major
 273 *  ordering  of  the  processes. This routine retrieves a default system
 274 *  context  which  will  include all available processes. In addition it
 275 *  spawns the processes if needed.
 276 *
 277 *  Arguments
 278 *  =========
 279 *
 280 *  ICTXT   (global output) INTEGER
 281 *          ICTXT specifies the BLACS context handle identifying the
 282 *          created process grid.  The context itself is global.
 283 *
 284 *  NPROW   (global input) INTEGER
 285 *          NPROW specifies the number of process rows in the grid
 286 *          to be created.
 287 *
 288 *  NPCOL   (global input) INTEGER
 289 *          NPCOL specifies the number of process columns in the grid
 290 *          to be created.
 291 *
 292 *  =====================================================================
 293 *
 294 *     .. Local Scalars ..
 295       INTEGER            IAM, NPROCS
 296 *     ..
 297 *     .. External Subroutines ..
 298       EXTERNAL           BLACS_GET, BLACS_GRIDINIT, BLACS_PINFO,
 299      $                   BLACS_SETUP
 300 *     ..
 301 *     .. Executable Statements ..
 302 *
 303 *     Get starting information
 304 *
 305       CALL BLACS_PINFO( IAM, NPROCS )
 306 *
 307 *     If machine needs additional set up, do it now
 308 *
 309       IF( NPROCS.LT.1 ) THEN
 310          IF( IAM.EQ.0 )
 311      $      NPROCS = NPROW*NPCOL
 312          CALL BLACS_SETUP( IAM, NPROCS )
 313       END IF
 314 *
 315 *     Define process grid
 316 *
 317       CALL BLACS_GET( -1, 0, ICTXT )
 318       CALL BLACS_GRIDINIT( ICTXT, 'Row-major', NPROW, NPCOL )
 319 *
 320       RETURN
 321 *
 322 *     End of SL_INIT
 323 *
 324       END