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)