6169
6170
6171
6172
6173
6174
6175
6176 CHARACTER*1 TRANS, UPLO
6177 INTEGER IA, IB, IC, ICTXT, INFO, JA, JB, JC, K, N
6178 DOUBLE PRECISION ERR
6179 COMPLEX*16 ALPHA, BETA
6180
6181
6182 INTEGER DESCA( * ), DESCB( * ), DESCC( * )
6183 DOUBLE PRECISION G( * )
6184 COMPLEX*16 A( * ), B( * ), C( * ), CT( * ),
6185 $ PC( * )
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
6360 INTEGER BLOCK_CYCLIC_2D_INB, CSRC_, CTXT_, DLEN_,
6361 $ DTYPE_, IMB_, INB_, LLD_, MB_, M_, NB_, N_,
6362 $ RSRC_
6363 parameter( block_cyclic_2d_inb = 2, dlen_ = 11,
6364 $ dtype_ = 1, ctxt_ = 2, m_ = 3, n_ = 4,
6365 $ imb_ = 5, inb_ = 6, mb_ = 7, nb_ = 8,
6366 $ rsrc_ = 9, csrc_ = 10, lld_ = 11 )
6367 DOUBLE PRECISION RZERO, RONE
6368 parameter( rzero = 0.0d+0, rone = 1.0d+0 )
6369 COMPLEX*16 ZERO
6370 parameter( zero = ( 0.0d+0, 0.0d+0 ) )
6371
6372
6373 LOGICAL COLREP, HTRAN, NOTRAN, ROWREP, TRAN, UPPER
6374 INTEGER I, IBB, IBEG, ICCOL, ICROW, ICURROW, IEND, IIC,
6375 $ IN, IOFFAK, IOFFAN, IOFFBK, IOFFBN, IOFFC, J,
6376 $ JJC, KK, LDA, LDB, LDC, LDPC, MYCOL, MYROW,
6377 $ NPCOL, NPROW
6378 DOUBLE PRECISION EPS, ERRI
6379 COMPLEX*16 Z
6380
6381
6382 EXTERNAL blacs_gridinfo, dgamx2d, igsum2d,
pb_infog2l
6383
6384
6385 LOGICAL LSAME
6386 DOUBLE PRECISION PDLAMCH
6388
6389
6390 INTRINSIC abs, dble, dconjg, dimag,
max,
min, mod, sqrt
6391
6392
6393 DOUBLE PRECISION ABS1
6394 abs1( z ) = abs( dble( z ) ) + abs( dimag( z ) )
6395
6396
6397
6398 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
6399
6401
6402 upper =
lsame( uplo,
'U' )
6403 htran =
lsame( trans,
'H' )
6404 notran =
lsame( trans,
'N' )
6405 tran =
lsame( trans,
'T' )
6406
6407 lda =
max( 1, desca( m_ ) )
6408 ldb =
max( 1, descb( m_ ) )
6409 ldc =
max( 1, descc( m_ ) )
6410
6411
6412
6413
6414
6415 DO 140 j = 1, n
6416
6417 IF( upper ) THEN
6418 ibeg = 1
6419 iend = j
6420 ELSE
6421 ibeg = j
6422 iend = n
6423 END IF
6424
6425 DO 10 i = 1, n
6426 ct( i ) = zero
6427 g( i ) = rzero
6428 10 CONTINUE
6429
6430 IF( notran ) THEN
6431 DO 30 kk = 1, k
6432 ioffak = ia + j - 1 + ( ja + kk - 2 ) * lda
6433 ioffbk = ib + j - 1 + ( jb + kk - 2 ) * ldb
6434 DO 20 i = ibeg, iend
6435 ioffan = ia + i - 1 + ( ja + kk - 2 ) * lda
6436 ioffbn = ib + i - 1 + ( jb + kk - 2 ) * ldb
6437 ct( i ) = ct( i ) + alpha * (
6438 $ a( ioffan ) * b( ioffbk ) +
6439 $ b( ioffbn ) * a( ioffak ) )
6440 g( i ) = g( i ) + abs( alpha ) * (
6441 $ abs1( a( ioffan ) ) * abs1( b( ioffbk ) ) +
6442 $ abs1( b( ioffbn ) ) * abs1( a( ioffak ) ) )
6443 20 CONTINUE
6444 30 CONTINUE
6445 ELSE IF( tran ) THEN
6446 DO 50 kk = 1, k
6447 ioffak = ia + kk - 1 + ( ja + j - 2 ) * lda
6448 ioffbk = ib + kk - 1 + ( jb + j - 2 ) * ldb
6449 DO 40 i = ibeg, iend
6450 ioffan = ia + kk - 1 + ( ja + i - 2 ) * lda
6451 ioffbn = ib + kk - 1 + ( jb + i - 2 ) * ldb
6452 ct( i ) = ct( i ) + alpha * (
6453 $ a( ioffan ) * b( ioffbk ) +
6454 $ b( ioffbn ) * a( ioffak ) )
6455 g( i ) = g( i ) + abs( alpha ) * (
6456 $ abs1( a( ioffan ) ) * abs1( b( ioffbk ) ) +
6457 $ abs1( b( ioffbn ) ) * abs1( a( ioffak ) ) )
6458 40 CONTINUE
6459 50 CONTINUE
6460 ELSE IF( htran ) THEN
6461 DO 70 kk = 1, k
6462 ioffak = ia + j - 1 + ( ja + kk - 2 ) * lda
6463 ioffbk = ib + j - 1 + ( jb + kk - 2 ) * ldb
6464 DO 60 i = ibeg, iend
6465 ioffan = ia + i - 1 + ( ja + kk - 2 ) * lda
6466 ioffbn = ib + i - 1 + ( jb + kk - 2 ) * ldb
6467 ct( i ) = ct( i ) +
6468 $ alpha * a( ioffan ) * dconjg( b( ioffbk ) ) +
6469 $ b( ioffbn ) * dconjg( alpha * a( ioffak ) )
6470 g( i ) = g( i ) + abs1( alpha ) * (
6471 $ abs1( a( ioffan ) ) * abs1( b( ioffbk ) ) +
6472 $ abs1( b( ioffbn ) ) * abs1( a( ioffak ) ) )
6473 60 CONTINUE
6474 70 CONTINUE
6475 ELSE
6476 DO 90 kk = 1, k
6477 ioffak = ia + kk - 1 + ( ja + j - 2 ) * lda
6478 ioffbk = ib + kk - 1 + ( jb + j - 2 ) * ldb
6479 DO 80 i = ibeg, iend
6480 ioffan = ia + kk - 1 + ( ja + i - 2 ) * lda
6481 ioffbn = ib + kk - 1 + ( jb + i - 2 ) * ldb
6482 ct( i ) = ct( i ) +
6483 $ alpha * dconjg( a( ioffan ) ) * b( ioffbk ) +
6484 $ dconjg( alpha * b( ioffbn ) ) * a( ioffak )
6485 g( i ) = g( i ) + abs1( alpha ) * (
6486 $ abs1( dconjg( a( ioffan ) ) * b( ioffbk ) ) +
6487 $ abs1( dconjg( b( ioffbn ) ) * a( ioffak ) ) )
6488 80 CONTINUE
6489 90 CONTINUE
6490 END IF
6491
6492 ioffc = ic + ibeg - 1 + ( jc + j - 2 ) * ldc
6493
6494 DO 100 i = ibeg, iend
6495 ct( i ) = ct( i ) + beta * c( ioffc )
6496 g( i ) = g( i ) + abs1( beta )*abs1( c( ioffc ) )
6497 c( ioffc ) = ct( i )
6498 ioffc = ioffc + 1
6499 100 CONTINUE
6500
6501
6502
6503 err = rzero
6504 info = 0
6505 ldpc = descc( lld_ )
6506 ioffc = ic + ( jc + j - 2 ) * ldc
6507 CALL pb_infog2l( ic, jc+j-1, descc, nprow, npcol, myrow, mycol,
6508 $ iic, jjc, icrow, iccol )
6509 icurrow = icrow
6510 rowrep = ( icrow.EQ.-1 )
6511 colrep = ( iccol.EQ.-1 )
6512
6513 IF( mycol.EQ.iccol .OR. colrep ) THEN
6514
6515 ibb = descc( imb_ ) - ic + 1
6516 IF( ibb.LE.0 )
6517 $ ibb = ( ( -ibb ) / descc( mb_ ) + 1 )*descc( mb_ ) + ibb
6519 in = ic + ibb - 1
6520
6521 DO 110 i = ic, in
6522
6523 IF( myrow.EQ.icurrow .OR. rowrep ) THEN
6524 erri = abs( pc( iic+(jjc-1)*ldpc ) -
6525 $ c( ioffc ) ) / eps
6526 IF( g( i-ic+1 ).NE.rzero )
6527 $ erri = erri / g( i-ic+1 )
6528 err =
max( err, erri )
6529 IF( err*sqrt( eps ).GE.rone )
6530 $ info = 1
6531 iic = iic + 1
6532 END IF
6533
6534 ioffc = ioffc + 1
6535
6536 110 CONTINUE
6537
6538 icurrow = mod( icurrow+1, nprow )
6539
6540 DO 130 i = in+1, ic+n-1, descc( mb_ )
6541 ibb =
min( ic+n-i, descc( mb_ ) )
6542
6543 DO 120 kk = 0, ibb-1
6544
6545 IF( myrow.EQ.icurrow .OR. rowrep ) THEN
6546 erri = abs( pc( iic+(jjc-1)*ldpc ) -
6547 $ c( ioffc ) )/eps
6548 IF( g( i+kk-ic+1 ).NE.rzero )
6549 $ erri = erri / g( i+kk-ic+1 )
6550 err =
max( err, erri )
6551 IF( err*sqrt( eps ).GE.rone )
6552 $ info = 1
6553 iic = iic + 1
6554 END IF
6555
6556 ioffc = ioffc + 1
6557
6558 120 CONTINUE
6559
6560 icurrow = mod( icurrow+1, nprow )
6561
6562 130 CONTINUE
6563
6564 END IF
6565
6566
6567
6568 CALL igsum2d( ictxt, 'All', ' ', 1, 1, info, 1, -1, mycol )
6569 CALL dgamx2d( ictxt, 'All', ' ', 1, 1, err, 1, i, j, -1, -1,
6570 $ mycol )
6571 IF( info.NE.0 )
6572 $ GO TO 150
6573
6574 140 CONTINUE
6575
6576 150 CONTINUE
6577
6578 RETURN
6579
6580
6581
subroutine pb_infog2l(i, j, desc, nprow, npcol, myrow, mycol, ii, jj, prow, pcol)
double precision function pdlamch(ictxt, cmach)