6168
 6169
 6170
 6171
 6172
 6173
 6174
 6175      CHARACTER*1        TRANS, UPLO
 6176      INTEGER            IA, IB, IC, ICTXT, INFO, JA, JB, JC, K, N
 6177      REAL               ERR
 6178      COMPLEX            ALPHA, BETA
 6179
 6180
 6181      INTEGER            DESCA( * ), DESCB( * ), DESCC( * )
 6182      REAL               G( * )
 6183      COMPLEX            A( * ), B( * ), C( * ), CT( * ),
 6184     $                   PC( * )
 6185
 6186
 6187
 6188
 6189
 6190
 6191
 6192
 6193
 6194
 6195
 6196
 6197
 6198
 6199
 6200
 6201
 6202
 6203
 6204
 6205
 6206
 6207
 6208
 6209
 6210
 6211
 6212
 6213
 6214
 6215
 6216
 6217
 6218
 6219
 6220
 6221
 6222
 6223
 6224
 6225
 6226
 6227
 6228
 6229
 6230
 6231
 6232
 6233
 6234
 6235
 6236
 6237
 6238
 6239
 6240
 6241
 6242
 6243
 6244
 6245
 6246
 6247
 6248
 6249
 6250
 6251
 6252
 6253
 6254
 6255
 6256
 6257
 6258
 6259
 6260
 6261
 6262
 6263
 6264
 6265
 6266
 6267
 6268
 6269
 6270
 6271
 6272
 6273
 6274
 6275
 6276
 6277
 6278
 6279
 6280
 6281
 6282
 6283
 6284
 6285
 6286
 6287
 6288
 6289
 6290
 6291
 6292
 6293
 6294
 6295
 6296
 6297
 6298
 6299
 6300
 6301
 6302
 6303
 6304
 6305
 6306
 6307
 6308
 6309
 6310
 6311
 6312
 6313
 6314
 6315
 6316
 6317
 6318
 6319
 6320
 6321
 6322
 6323
 6324
 6325
 6326
 6327
 6328
 6329
 6330
 6331
 6332
 6333
 6334
 6335
 6336
 6337
 6338
 6339
 6340
 6341
 6342
 6343
 6344
 6345
 6346
 6347
 6348
 6349
 6350
 6351
 6352
 6353
 6354
 6355
 6356
 6357
 6358
 6359      INTEGER            BLOCK_CYCLIC_2D_INB, CSRC_, CTXT_, DLEN_,
 6360     $                   DTYPE_, IMB_, INB_, LLD_, MB_, M_, NB_, N_,
 6361     $                   RSRC_
 6362      parameter( block_cyclic_2d_inb = 2, dlen_ = 11,
 6363     $                   dtype_ = 1, ctxt_ = 2, m_ = 3, n_ = 4,
 6364     $                   imb_ = 5, inb_ = 6, mb_ = 7, nb_ = 8,
 6365     $                   rsrc_ = 9, csrc_ = 10, lld_ = 11 )
 6366      REAL               RZERO, RONE
 6367      parameter( rzero = 0.0e+0, rone = 1.0e+0 )
 6368      COMPLEX            ZERO
 6369      parameter( zero = ( 0.0e+0, 0.0e+0 ) )
 6370
 6371
 6372      LOGICAL            COLREP, HTRAN, NOTRAN, ROWREP, TRAN, UPPER
 6373      INTEGER            I, IBB, IBEG, ICCOL, ICROW, ICURROW, IEND, IIC,
 6374     $                   IN, IOFFAK, IOFFAN, IOFFBK, IOFFBN, IOFFC, J,
 6375     $                   JJC, KK, LDA, LDB, LDC, LDPC, MYCOL, MYROW,
 6376     $                   NPCOL, NPROW
 6377      REAL               EPS, ERRI
 6378      COMPLEX            Z
 6379
 6380
 6381      EXTERNAL           blacs_gridinfo, igsum2d, 
pb_infog2l, sgamx2d
 
 6382
 6383
 6384      LOGICAL            LSAME
 6385      REAL               PSLAMCH
 6387
 6388
 6389      INTRINSIC          abs, aimag, conjg, 
max, 
min, mod, real, sqrt
 
 6390
 6391
 6392      REAL               ABS1
 6393      abs1( z ) = abs( real( z ) ) + abs( aimag( z ) )
 6394
 6395
 6396
 6397      CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
 6398
 6400
 6401      upper = 
lsame( uplo, 
'U' )
 
 6402      htran = 
lsame( trans, 
'H' )
 
 6403      notran = 
lsame( trans, 
'N' )
 
 6404      tran = 
lsame( trans, 
'T' )
 
 6405
 6406      lda = 
max( 1, desca( m_ ) )
 
 6407      ldb = 
max( 1, descb( m_ ) )
 
 6408      ldc = 
max( 1, descc( m_ ) )
 
 6409
 6410
 6411
 6412
 6413
 6414      DO 140 j = 1, n
 6415
 6416         IF( upper ) THEN
 6417            ibeg = 1
 6418            iend = j
 6419         ELSE
 6420            ibeg = j
 6421            iend = n
 6422         END IF
 6423
 6424         DO 10 i = 1, n
 6425            ct( i ) = zero
 6426            g( i )  = rzero
 6427   10    CONTINUE
 6428
 6429         IF( notran ) THEN
 6430            DO 30 kk = 1, k
 6431               ioffak = ia + j - 1 + ( ja + kk - 2 ) * lda
 6432               ioffbk = ib + j - 1 + ( jb + kk - 2 ) * ldb
 6433               DO 20 i = ibeg, iend
 6434                  ioffan = ia + i - 1 + ( ja + kk - 2 ) * lda
 6435                  ioffbn = ib + i - 1 + ( jb + kk - 2 ) * ldb
 6436                  ct( i ) = ct( i ) + alpha * (
 6437     $                      a( ioffan ) * b( ioffbk ) +
 6438     $                      b( ioffbn ) * a( ioffak ) )
 6439                  g( i ) = g( i ) + abs( alpha ) * (
 6440     $                     abs1( a( ioffan ) ) * abs1( b( ioffbk ) ) +
 6441     $                     abs1( b( ioffbn ) ) * abs1( a( ioffak ) ) )
 6442   20          CONTINUE
 6443   30       CONTINUE
 6444         ELSE IF( tran ) THEN
 6445            DO 50 kk = 1, k
 6446               ioffak = ia + kk - 1 + ( ja + j - 2 ) * lda
 6447               ioffbk = ib + kk - 1 + ( jb + j - 2 ) * ldb
 6448               DO 40 i = ibeg, iend
 6449                  ioffan = ia + kk - 1 + ( ja + i - 2 ) * lda
 6450                  ioffbn = ib + kk - 1 + ( jb + i - 2 ) * ldb
 6451                  ct( i ) = ct( i ) + alpha * (
 6452     $                      a( ioffan ) * b( ioffbk ) +
 6453     $                      b( ioffbn ) * a( ioffak ) )
 6454                  g( i ) = g( i ) + abs( alpha ) * (
 6455     $                     abs1( a( ioffan ) ) * abs1( b( ioffbk ) ) +
 6456     $                     abs1( b( ioffbn ) ) * abs1( a( ioffak ) ) )
 6457   40          CONTINUE
 6458   50       CONTINUE
 6459         ELSE IF( htran ) THEN
 6460            DO 70 kk = 1, k
 6461               ioffak = ia + j - 1 + ( ja + kk - 2 ) * lda
 6462               ioffbk = ib + j - 1 + ( jb + kk - 2 ) * ldb
 6463               DO 60 i = ibeg, iend
 6464                  ioffan = ia + i - 1 + ( ja + kk - 2 ) * lda
 6465                  ioffbn = ib + i - 1 + ( jb + kk - 2 ) * ldb
 6466                  ct( i ) = ct( i ) +
 6467     $                alpha * a( ioffan ) * conjg( b( ioffbk ) ) +
 6468     $                b( ioffbn ) * conjg( alpha * a( ioffak ) )
 6469                  g( i ) = g( i ) + abs1( alpha ) * (
 6470     $                     abs1( a( ioffan ) ) * abs1( b( ioffbk ) ) +
 6471     $                     abs1( b( ioffbn ) ) * abs1( a( ioffak ) ) )
 6472   60          CONTINUE
 6473   70       CONTINUE
 6474         ELSE
 6475            DO 90 kk = 1, k
 6476               ioffak = ia + kk - 1 + ( ja + j - 2 ) * lda
 6477               ioffbk = ib + kk - 1 + ( jb + j - 2 ) * ldb
 6478               DO 80 i = ibeg, iend
 6479                  ioffan = ia + kk - 1 + ( ja + i - 2 ) * lda
 6480                  ioffbn = ib + kk - 1 + ( jb + i - 2 ) * ldb
 6481                  ct( i ) = ct( i ) +
 6482     $                   alpha * conjg( a( ioffan ) ) * b( ioffbk ) +
 6483     $                   conjg( alpha * b( ioffbn ) ) * a( ioffak )
 6484                  g( i ) = g( i ) + abs1( alpha ) * (
 6485     $                   abs1( conjg( a( ioffan ) ) * b( ioffbk ) ) +
 6486     $                   abs1( conjg( b( ioffbn ) ) * a( ioffak ) ) )
 6487   80          CONTINUE
 6488   90       CONTINUE
 6489         END IF
 6490
 6491         ioffc = ic + ibeg - 1 + ( jc + j - 2 ) * ldc
 6492
 6493         DO 100 i = ibeg, iend
 6494            ct( i ) = ct( i ) + beta * c( ioffc )
 6495            g( i ) = g( i ) + abs1( beta )*abs1( c( ioffc ) )
 6496            c( ioffc ) = ct( i )
 6497            ioffc = ioffc + 1
 6498  100    CONTINUE
 6499
 6500
 6501
 6502         err  = rzero
 6503         info = 0
 6504         ldpc = descc( lld_ )
 6505         ioffc = ic + ( jc + j - 2 ) * ldc
 6506         CALL pb_infog2l( ic, jc+j-1, descc, nprow, npcol, myrow, mycol,
 
 6507     $                    iic, jjc, icrow, iccol )
 6508         icurrow = icrow
 6509         rowrep  = ( icrow.EQ.-1 )
 6510         colrep  = ( iccol.EQ.-1 )
 6511
 6512         IF( mycol.EQ.iccol .OR. colrep ) THEN
 6513
 6514            ibb = descc( imb_ ) - ic + 1
 6515            IF( ibb.LE.0 )
 6516     $         ibb = ( ( -ibb ) / descc( mb_ ) + 1 )*descc( mb_ ) + ibb
 6518            in = ic + ibb - 1
 6519
 6520            DO 110 i = ic, in
 6521
 6522               IF( myrow.EQ.icurrow .OR. rowrep ) THEN
 6523                  erri = abs( pc( iic+(jjc-1)*ldpc ) -
 6524     $                        c( ioffc ) ) / eps
 6525                  IF( g( i-ic+1 ).NE.rzero )
 6526     $               erri = erri / g( i-ic+1 )
 6527                  err = 
max( err, erri )
 
 6528                  IF( err*sqrt( eps ).GE.rone )
 6529     $               info = 1
 6530                  iic = iic + 1
 6531               END IF
 6532
 6533               ioffc = ioffc + 1
 6534
 6535  110       CONTINUE
 6536
 6537            icurrow = mod( icurrow+1, nprow )
 6538
 6539            DO 130 i = in+1, ic+n-1, descc( mb_ )
 6540               ibb = 
min( ic+n-i, descc( mb_ ) )
 
 6541
 6542               DO 120 kk = 0, ibb-1
 6543
 6544                  IF( myrow.EQ.icurrow .OR. rowrep ) THEN
 6545                     erri = abs( pc( iic+(jjc-1)*ldpc ) -
 6546     $                           c( ioffc ) )/eps
 6547                     IF( g( i+kk-ic+1 ).NE.rzero )
 6548     $                  erri = erri / g( i+kk-ic+1 )
 6549                     err = 
max( err, erri )
 
 6550                     IF( err*sqrt( eps ).GE.rone )
 6551     $                  info = 1
 6552                     iic = iic + 1
 6553                  END IF
 6554
 6555                  ioffc = ioffc + 1
 6556
 6557  120          CONTINUE
 6558
 6559               icurrow = mod( icurrow+1, nprow )
 6560
 6561  130       CONTINUE
 6562
 6563         END IF
 6564
 6565
 6566
 6567         CALL igsum2d( ictxt, 'All', ' ', 1, 1, info, 1, -1, mycol )
 6568         CALL sgamx2d( ictxt, 'All', ' ', 1, 1, err, 1, i, j, -1, -1,
 6569     $                 mycol )
 6570         IF( info.NE.0 )
 6571     $      GO TO 150
 6572
 6573  140 CONTINUE
 6574
 6575  150 CONTINUE
 6576
 6577      RETURN
 6578
 6579
 6580
subroutine pb_infog2l(i, j, desc, nprow, npcol, myrow, mycol, ii, jj, prow, pcol)
 
real function pslamch(ictxt, cmach)