1 SUBROUTINE pzlascal( TYPE, M, N, ALPHA, A, IA, JA, DESCA )
136 INTEGER BLOCK_CYCLIC_2D_INB, CSRC_, CTXT_, DLEN_,
137 $ DTYPE_, IMB_, INB_, LLD_, MB_, M_, NB_, N_,
139 parameter( block_cyclic_2d_inb = 2, dlen_ = 11,
140 $ dtype_ = 1, ctxt_ = 2, m_ = 3, n_ = 4,
141 $ imb_ = 5, inb_ = 6, mb_ = 7, nb_ = 8,
142 $ rsrc_ = 9, csrc_ = 10, lld_ = 11 )
146 LOGICAL GODOWN, GOLEFT, LOWER, UPPER
147 INTEGER IACOL, IAROW, ICTXT, IIA, IIMAX, ILOW, IMB1,
148 $ IMBLOC, INB1, INBLOC, IOFFA, IOFFD, ITYPE,
149 $ IUPP, JJA, JJMAX, JOFFA, JOFFD, LCMT, LCMT00,
150 $ LDA, LMBLOC, LNBLOC, LOW, M1, MB, MBLKD, MBLKS,
151 $ MBLOC, MP, MRCOL, MRROW, MYCOL, MYROW, N1, NB,
152 $ NBLKD, NBLKS, NBLOC, NPCOL, NPROW, NQ, PMB,
156 INTEGER DESCA2( DLEN_ )
165 EXTERNAL lsame, pb_numroc
178 ictxt = desca2( ctxt_ )
179 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
183 IF( m.EQ.0 .OR. n.EQ.0 )
186 IF( lsame(
TYPE,
'L' ) ) then
192 ELSE IF( lsame(
TYPE,
'U' ) ) then
198 ELSE IF( lsame(
TYPE,
'H' ) ) then
214 IF( itype.EQ.0 )
THEN
218 CALL pb_infog2l( ia, ja, desca2, nprow, npcol, myrow, mycol,
219 $ iia, jja, iarow, iacol )
220 mp = pb_numroc( m, ia, desca2( imb_ ), desca2( mb_ ), myrow,
221 $ desca2( rsrc_ ), nprow )
222 nq = pb_numroc( n, ja, desca2( inb_ ), desca2( nb_ ), mycol,
223 $ desca2( csrc_ ), npcol )
225 IF( mp.LE.0 .OR. nq.LE.0 )
229 ioffa = iia + ( jja - 1 ) * lda
231 CALL pb_zlascal(
'All', mp, nq, 0, alpha, a( ioffa ), lda )
237 CALL pb_ainfog2l( m, n, ia, ja, desca2, nprow, npcol, myrow,
238 $ mycol, imb1, inb1, mp, nq, iia, jja, iarow,
239 $ iacol, mrrow, mrcol )
241 IF( mp.LE.0 .OR. nq.LE.0 )
251 CALL pb_binfo( ioffd, mp, nq, imb1, inb1, mb, nb, mrrow,
252 $ mrcol, lcmt00, mblks, nblks, imbloc, inbloc,
253 $ lmbloc, lnbloc, ilow, low, iupp, upp )
262 IF( desca2( rsrc_ ).LT.0 )
THEN
267 IF( desca2( csrc_ ).LT.0 )
THEN
276 godown = ( lcmt00.GT.iupp )
277 goleft = ( lcmt00.LT.ilow )
279 IF( .NOT.godown .AND. .NOT.goleft )
THEN
283 goleft = ( ( lcmt00 - ( iupp - upp + pmb ) ).LT.ilow )
286 CALL pb_zlascal( uplo, imbloc, inbloc, lcmt00, alpha,
287 $ a( iia+joffa*lda ), lda )
289 IF( upper .AND. nq.GT.inbloc )
290 $
CALL pb_zlascal(
'All', imbloc, nq-inbloc, 0, alpha,
291 $ a( iia+(joffa+inbloc)*lda ), lda )
295 IF( lower .AND. mp.GT.imbloc )
296 $
CALL pb_zlascal(
'All', mp-imbloc, inbloc, 0, alpha,
297 $ a( iia+imbloc+joffa*lda ), lda )
306 lcmt00 = lcmt00 - ( iupp - upp + pmb )
308 ioffa = ioffa + imbloc
311 IF( mblks.GT.0 .AND. lcmt00.GT.upp )
THEN
312 lcmt00 = lcmt00 - pmb
318 tmp1 =
min( ioffa, iimax ) - iia + 1
319 IF( upper .AND. tmp1.GT.0 )
THEN
321 $ a( iia+joffa*lda ), lda )
335 IF( mblkd.GT.0 .AND. lcmt.GE.ilow )
THEN
338 CALL pb_zlascal( uplo, mbloc, inbloc, lcmt, alpha,
339 $ a( ioffd+1+joffa*lda ), lda )
345 ioffd = ioffd + mbloc
349 tmp1 = m1 - ioffd + iia - 1
350 IF( lower .AND. tmp1.GT.0 )
351 $
CALL pb_zlascal(
'All', tmp1, inbloc, 0, alpha,
352 $ a( ioffd+1+joffa*lda ), lda )
354 tmp1 = ioffa - iia + 1
357 lcmt00 = lcmt00 + low - ilow + qnb
359 joffa = joffa + inbloc
361 IF( upper .AND. tmp1.GT.0 .AND. n1.GT.0 )
363 $ a( iia+joffa*lda ), lda )
368 ELSE IF( goleft )
THEN
370 lcmt00 = lcmt00 + low - ilow + qnb
372 joffa = joffa + inbloc
375 IF( nblks.GT.0 .AND. lcmt00.LT.low )
THEN
376 lcmt00 = lcmt00 + qnb
382 tmp1 =
min( joffa, jjmax ) - jja + 1
383 IF( lower .AND. tmp1.GT.0 )
THEN
385 $ a( iia+(jja-1)*lda ), lda )
399 IF( nblkd.GT.0 .AND. lcmt.LE.iupp )
THEN
402 CALL pb_zlascal( uplo, imbloc, nbloc, lcmt, alpha,
403 $ a( iia+joffd*lda ), lda )
409 joffd = joffd + nbloc
413 tmp1 = n1 - joffd + jja - 1
414 IF( upper .AND. tmp1.GT.0 )
415 $
CALL pb_zlascal(
'All', imbloc, tmp1, 0, alpha,
416 $ a( iia+joffd*lda ), lda )
418 tmp1 = joffa - jja + 1
421 lcmt00 = lcmt00 - ( iupp - upp + pmb )
423 ioffa = ioffa + imbloc
425 IF( lower .AND. m1.GT.0 .AND. tmp1.GT.0 )
427 $ a( ioffa+1+(jja-1)*lda ), lda )
436 IF( nblks.GT.0 )
THEN
440 IF( mblks.GT.0 .AND. lcmt00.GT.upp )
THEN
441 lcmt00 = lcmt00 - pmb
447 tmp1 =
min( ioffa, iimax ) - iia + 1
448 IF( upper .AND. tmp1.GT.0 )
THEN
450 $ a( iia+joffa*lda ), lda )
464 IF( mblkd.GT.0 .AND. lcmt.GE.low )
THEN
467 CALL pb_zlascal( uplo, mbloc, nbloc, lcmt, alpha,
468 $ a( ioffd+1+joffa*lda ), lda )
474 ioffd = ioffd + mbloc
478 tmp1 = m1 - ioffd + iia - 1
479 IF( lower .AND. tmp1.GT.0 )
480 $
CALL pb_zlascal(
'All', tmp1, nbloc, 0, alpha,
481 $ a( ioffd+1+joffa*lda ), lda )
483 tmp1 =
min( ioffa, iimax ) - iia + 1
486 lcmt00 = lcmt00 + qnb
488 joffa = joffa + nbloc
490 IF( upper .AND. tmp1.GT.0 .AND. n1.GT.0 )
492 $ a( iia+joffa*lda ), lda )
508 SUBROUTINE pzlagen( INPLACE, AFORM, DIAG, OFFA, M, N, IA, JA,
509 $ DESCA, IASEED, A, LDA )
518 CHARACTER*1 aform, diag
519 INTEGER ia, iaseed, ja, lda, m, n, offa
523 COMPLEX*16 A( LDA, * )
701 INTEGER BLOCK_CYCLIC_2D_INB, CSRC_, CTXT_, DLEN_,
702 $ dtype_, imb_, inb_, lld_, mb_, m_, nb_, n_,
704 parameter( block_cyclic_2d_inb = 2, dlen_ = 11,
705 $ dtype_ = 1, ctxt_ = 2, m_ = 3, n_ = 4,
706 $ imb_ = 5, inb_ = 6, mb_ = 7, nb_ = 8,
707 $ rsrc_ = 9, csrc_ = 10, lld_ = 11 )
708 INTEGER JMP_1, JMP_COL, JMP_IMBV, JMP_INBV, JMP_LEN,
709 $ jmp_mb, jmp_nb, jmp_npimbloc, jmp_npmb,
710 $ jmp_nqinbloc, jmp_nqnb, jmp_row
711 parameter( jmp_1 = 1, jmp_row = 2, jmp_col = 3,
712 $ jmp_mb = 4, jmp_imbv = 5, jmp_npmb = 6,
713 $ jmp_npimbloc = 7, jmp_nb = 8, jmp_inbv = 9,
714 $ jmp_nqnb = 10, jmp_nqinbloc = 11,
716 DOUBLE PRECISION ZERO
717 parameter( zero = 0.0d+0 )
720 LOGICAL DIAGDO, SYMM, HERM, NOTRAN
721 INTEGER CSRC, I, IACOL, IAROW, ICTXT, IIA, ILOCBLK,
722 $ ilocoff, ilow, imb, imb1, imbloc, imbvir, inb,
723 $ inb1, inbloc, inbvir, info, ioffda, itmp, iupp,
724 $ ivir, jja, jlocblk, jlocoff, jvir, lcmt00,
725 $ lmbloc, lnbloc, low, maxmn, mb, mblks, mp,
726 $ mrcol, mrrow, mycdist, mycol, myrdist, myrow,
727 $ nb, nblks, npcol, nprow, nq, nvir, rsrc, upp
731 INTEGER DESCA2( DLEN_ ), IMULADD( 4, JMP_LEN ),
732 $ iran( 2 ), jmp( jmp_len ), muladd0( 4 )
746 INTRINSIC dble, dcmplx,
max,
min
749 DATA ( muladd0( i ), i = 1, 4 ) / 20077, 16838,
760 ictxt = desca2( ctxt_ )
761 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
766 IF( nprow.EQ.-1 )
THEN
767 info = -( 1000 + ctxt_ )
769 symm = lsame( aform,
'S' )
770 herm = lsame( aform,
'H' )
771 notran = lsame( aform,
'N' )
772 diagdo = lsame( diag,
'D' )
773 IF( .NOT.( symm.OR.herm.OR.notran ) .AND.
774 $ .NOT.( lsame( aform,
'T' ) ) .AND.
775 $ .NOT.( lsame( aform,
'C' ) ) )
THEN
777 ELSE IF( ( .NOT.diagdo ) .AND.
778 $ ( .NOT.lsame( diag,
'N' ) ) )
THEN
781 CALL pb_chkmat( ictxt, m, 5, n, 6, ia, ja, desca2, 10, info )
785 CALL pxerbla( ictxt,
'PZLAGEN', -info )
791 IF( ( m.LE.0 ).OR.( n.LE.0 ) )
800 rsrc = desca2( rsrc_ )
801 csrc = desca2( csrc_ )
805 CALL pb_ainfog2l( m, n, ia, ja, desca2, nprow, npcol, myrow,
806 $ mycol, imb1, inb1, mp, nq, iia, jja, iarow,
807 $ iacol, mrrow, mrcol )
819 ioffda = ja + offa - ia
820 CALL pb_binfo( ioffda, mp, nq, imb1, inb1, mb, nb, mrrow,
821 $ mrcol, lcmt00, mblks, nblks, imbloc, inbloc,
822 $ lmbloc, lnbloc, ilow, low, iupp, upp )
830 itmp =
max( 0, -offa )
833 nvir = desca2( m_ ) + itmp
835 CALL pb_locinfo( ivir, imbvir, mb, myrow, rsrc, nprow, ilocblk,
838 itmp =
max( 0, offa )
841 nvir =
max(
max( nvir, desca2( n_ ) + itmp ),
842 $ desca2( m_ ) + desca2( n_ ) - 1 )
844 CALL pb_locinfo( jvir, inbvir, nb, mycol, csrc, npcol, jlocblk,
847 IF( symm .OR. herm .OR. notran )
THEN
849 CALL pb_initjmp( .true., nvir, imbvir, inbvir, imbloc, inbloc,
850 $ mb, nb, rsrc, csrc, nprow, npcol, 2, jmp )
858 CALL pb_setlocran( iaseed, ilocblk, jlocblk, ilocoff, jlocoff,
859 $ myrdist, mycdist, nprow, npcol, jmp,
862 CALL pb_zlagen(
'Lower', aform, a( iia, jja ), lda, lcmt00,
863 $ iran, mblks, imbloc, mb, lmbloc, nblks, inbloc,
864 $ nb, lnbloc, jmp, imuladd )
868 IF( symm .OR. herm .OR. ( .NOT. notran ) )
THEN
870 CALL pb_initjmp( .false., nvir, imbvir, inbvir, imbloc, inbloc,
871 $ mb, nb, rsrc, csrc, nprow, npcol, 2, jmp )
879 CALL pb_setlocran( iaseed, ilocblk, jlocblk, ilocoff, jlocoff,
880 $ myrdist, mycdist, nprow, npcol, jmp,
883 CALL pb_zlagen(
'Upper', aform, a( iia, jja ), lda, lcmt00,
884 $ iran, mblks, imbloc, mb, lmbloc, nblks, inbloc,
885 $ nb, lnbloc, jmp, imuladd )
891 maxmn =
max( desca2( m_ ), desca2( n_ ) )
893 alpha = dcmplx( dble( 2 * maxmn ), zero )
895 alpha = dcmplx( dble( nvir ), dble( maxmn ) )
898 IF( ioffda.GE.0 )
THEN
900 $ a,
min( ia+ioffda, ia+m-1 ), ja, desca )
903 $ a, ia,
min( ja-ioffda, ja+n-1 ), desca )
913 SUBROUTINE pzladom( INPLACE, N, ALPHA, A, IA, JA, DESCA )
1041 INTEGER BLOCK_CYCLIC_2D_INB, CSRC_, CTXT_, DLEN_,
1042 $ dtype_, imb_, inb_, lld_, mb_, m_, nb_, n_,
1044 parameter( block_cyclic_2d_inb = 2, dlen_ = 11,
1045 $ dtype_ = 1, ctxt_ = 2, m_ = 3, n_ = 4,
1046 $ imb_ = 5, inb_ = 6, mb_ = 7, nb_ = 8,
1047 $ rsrc_ = 9, csrc_ = 10, lld_ = 11 )
1050 LOGICAL GODOWN, GOLEFT
1051 INTEGER I, IACOL, IAROW, ICTXT, IIA, IJOFFA, ILOW,
1052 $ imb1, imbloc, inb1, inbloc, ioffa, ioffd, iupp,
1053 $ jja, joffa, joffd, lcmt, lcmt00, lda, ldap1,
1054 $ lmbloc, lnbloc, low, mb, mblkd, mblks, mbloc,
1055 $ mrcol, mrrow, mycol, myrow, nb, nblkd, nblks,
1056 $ nbloc, np, npcol, nprow, nq, pmb, qnb, upp
1060 INTEGER DESCA2( DLEN_ )
1067 INTRINSIC abs, dble, dcmplx, dimag,
max,
min
1077 ictxt = desca2( ctxt_ )
1078 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
1083 CALL pb_ainfog2l( n, n, ia, ja, desca2, nprow, npcol, myrow,
1084 $ mycol, imb1, inb1, np, nq, iia, jja, iarow,
1085 $ iacol, mrrow, mrcol )
1100 CALL pb_binfo( 0, np, nq, imb1, inb1, mb, nb, mrrow, mrcol,
1101 $ lcmt00, mblks, nblks, imbloc, inbloc, lmbloc,
1102 $ lnbloc, ilow, low, iupp, upp )
1106 lda = desca2( lld_ )
1109 IF( desca2( rsrc_ ).LT.0 )
THEN
1114 IF( desca2( csrc_ ).LT.0 )
THEN
1123 godown = ( lcmt00.GT.iupp )
1124 goleft = ( lcmt00.LT.ilow )
1126 IF( .NOT.godown .AND. .NOT.goleft )
THEN
1130 IF( lcmt00.GE.0 )
THEN
1131 ijoffa = ioffa+lcmt00 + ( joffa - 1 ) * lda
1132 DO 10 i = 1,
min( inbloc,
max( 0, imbloc - lcmt00 ) )
1133 atmp = a( ijoffa + i*ldap1 )
1134 a( ijoffa + i*ldap1 ) = alpha +
1135 $ dcmplx( abs( dble( atmp ) ),
1136 $ abs( dimag( atmp ) ) )
1139 ijoffa = ioffa + ( joffa - lcmt00 - 1 ) * lda
1140 DO 20 i = 1,
min( imbloc,
max( 0, inbloc + lcmt00 ) )
1141 atmp = a( ijoffa + i*ldap1 )
1142 a( ijoffa + i*ldap1 ) = alpha +
1143 $ dcmplx( abs( dble( atmp ) ),
1144 $ abs( dimag( atmp ) ) )
1147 goleft = ( ( lcmt00 - ( iupp - upp + pmb ) ).LT.ilow )
1148 godown = .NOT.goleft
1154 lcmt00 = lcmt00 - ( iupp - upp + pmb )
1156 ioffa = ioffa + imbloc
1159 IF( mblks.GT.0 .AND. lcmt00.GT.upp )
THEN
1160 lcmt00 = lcmt00 - pmb
1172 IF( mblkd.GT.0 .AND. lcmt.GE.ilow )
THEN
1175 IF( lcmt.GE.0 )
THEN
1176 ijoffa = ioffd + lcmt + ( joffa - 1 ) * lda
1177 DO 50 i = 1,
min( inbloc,
max( 0, mbloc - lcmt ) )
1178 atmp = a( ijoffa + i*ldap1 )
1179 a( ijoffa + i*ldap1 ) = alpha +
1180 $ dcmplx( abs( dble( atmp ) ),
1181 $ abs( dimag( atmp ) ) )
1184 ijoffa = ioffd + ( joffa - lcmt - 1 ) * lda
1185 DO 60 i = 1,
min( mbloc,
max( 0, inbloc + lcmt ) )
1186 atmp = a( ijoffa + i*ldap1 )
1187 a( ijoffa + i*ldap1 ) = alpha +
1188 $ dcmplx( abs( dble( atmp ) ),
1189 $ abs( dimag( atmp ) ) )
1197 ioffd = ioffd + mbloc
1201 lcmt00 = lcmt00 + low - ilow + qnb
1203 joffa = joffa + inbloc
1205 ELSE IF( goleft )
THEN
1207 lcmt00 = lcmt00 + low - ilow + qnb
1209 joffa = joffa + inbloc
1212 IF( nblks.GT.0 .AND. lcmt00.LT.low )
THEN
1213 lcmt00 = lcmt00 + qnb
1225 IF( nblkd.GT.0 .AND. lcmt.LE.iupp )
THEN
1228 IF( lcmt.GE.0 )
THEN
1229 ijoffa = ioffa + lcmt + ( joffd - 1 ) * lda
1230 DO 90 i = 1,
min( nbloc,
max( 0, imbloc - lcmt ) )
1231 atmp = a( ijoffa + i*ldap1 )
1232 a( ijoffa + i*ldap1 ) = alpha +
1233 $ dcmplx( abs( dble( atmp ) ),
1234 $ abs( dimag( atmp ) ) )
1237 ijoffa = ioffa + ( joffd - lcmt - 1 ) * lda
1238 DO 100 i = 1,
min( imbloc,
max( 0, nbloc + lcmt ) )
1239 atmp = a( ijoffa + i*ldap1 )
1240 a( ijoffa + i*ldap1 ) = alpha +
1241 $ dcmplx( abs( dble( atmp ) ),
1242 $ abs( dimag( atmp ) ) )
1250 joffd = joffd + nbloc
1254 lcmt00 = lcmt00 - ( iupp - upp + pmb )
1256 ioffa = ioffa + imbloc
1262 IF( nblks.GT.0 )
THEN
1266 IF( mblks.GT.0 .AND. lcmt00.GT.upp )
THEN
1267 lcmt00 = lcmt00 - pmb
1279 IF( mblkd.GT.0 .AND. lcmt.GE.low )
THEN
1282 IF( lcmt.GE.0 )
THEN
1283 ijoffa = ioffd + lcmt + ( joffa - 1 ) * lda
1284 DO 140 i = 1,
min( nbloc,
max( 0, mbloc - lcmt ) )
1285 atmp = a( ijoffa + i*ldap1 )
1286 a( ijoffa + i*ldap1 ) = alpha +
1287 $ dcmplx( abs( dble( atmp ) ),
1288 $ abs( dimag( atmp ) ) )
1291 ijoffa = ioffd + ( joffa - lcmt - 1 ) * lda
1292 DO 150 i = 1,
min( mbloc,
max( 0, nbloc + lcmt ) )
1293 atmp = a( ijoffa + i*ldap1 )
1294 a( ijoffa + i*ldap1 ) = alpha +
1295 $ dcmplx( abs( dble( atmp ) ),
1296 $ abs( dimag( atmp ) ) )
1304 ioffd = ioffd + mbloc
1308 lcmt00 = lcmt00 + qnb
1310 joffa = joffa + nbloc
1320 SUBROUTINE pb_zlascal( UPLO, M, N, IOFFD, ALPHA, A, LDA )
1329 INTEGER IOFFD, LDA, M, N
1333 COMPLEX*16 A( LDA, * )
1423 INTEGER I, J, JTMP, MN
1436 IF( m.LE.0 .OR. n.LE.0 )
1441 IF( lsame( uplo,
'L' ) )
THEN
1445 mn =
max( 0, -ioffd )
1446 DO 20 j = 1,
min( mn, n )
1448 a( i, j ) = alpha * a( i, j )
1451 DO 40 j = mn + 1,
min( m - ioffd, n )
1452 DO 30 i = j + ioffd, m
1453 a( i, j ) = alpha * a( i, j )
1457 ELSE IF( lsame( uplo,
'U' ) )
THEN
1461 mn =
min( m - ioffd, n )
1462 DO 60 j =
max( 0, -ioffd ) + 1, mn
1463 DO 50 i = 1, j + ioffd
1464 a( i, j ) = alpha * a( i, j )
1467 DO 80 j =
max( 0, mn ) + 1, n
1469 a( i, j ) = alpha * a( i, j )
1473 ELSE IF( lsame( uplo,
'D' ) )
THEN
1477 DO 90 j =
max( 0, -ioffd ) + 1,
min( m - ioffd, n )
1479 a( jtmp, j ) = alpha * a( jtmp, j )
1488 a( i, j ) = alpha * a( i, j )
1499 SUBROUTINE pb_zlagen( UPLO, AFORM, A, LDA, LCMT00, IRAN, MBLKS,
1500 $ IMBLOC, MB, LMBLOC, NBLKS, INBLOC, NB,
1501 $ LNBLOC, JMP, IMULADD )
1509 CHARACTER*1 UPLO, AFORM
1510 INTEGER IMBLOC, INBLOC, LCMT00, LDA, LMBLOC, LNBLOC,
1511 $ MB, MBLKS, NB, NBLKS
1514 INTEGER IMULADD( 4, * ), IRAN( * ), JMP( * )
1515 COMPLEX*16 A( LDA, * )
1618 INTEGER JMP_1, JMP_COL, JMP_IMBV, JMP_INBV, JMP_LEN,
1619 $ JMP_MB, JMP_NB, JMP_NPIMBLOC, JMP_NPMB,
1620 $ JMP_NQINBLOC, JMP_NQNB, JMP_ROW
1621 parameter( jmp_1 = 1, jmp_row = 2, jmp_col = 3,
1622 $ jmp_mb = 4, jmp_imbv = 5, jmp_npmb = 6,
1623 $ jmp_npimbloc = 7, jmp_nb = 8, jmp_inbv = 9,
1624 $ jmp_nqnb = 10, jmp_nqinbloc = 11,
1626 DOUBLE PRECISION ZERO
1627 PARAMETER ( ZERO = 0.0d+0 )
1630 INTEGER I, IB, IBLK, II, IK, ITMP, JB, JBLK, JJ, JK,
1631 $ JTMP, LCMTC, LCMTR, LOW, MNB, UPP
1635 INTEGER IB0( 2 ), IB1( 2 ), IB2( 2 ), IB3( 2 )
1642 DOUBLE PRECISION PB_DRAND
1643 EXTERNAL LSAME, PB_DRAND
1646 INTRINSIC dble, dcmplx,
max,
min
1651 ib1( i ) = iran( i )
1652 ib2( i ) = iran( i )
1653 ib3( i ) = iran( i )
1656 IF( lsame( aform,
'N' ) )
THEN
1662 DO 50 jblk = 1, nblks
1664 IF( jblk.EQ.1 )
THEN
1666 ELSE IF( jblk.EQ.nblks )
THEN
1672 DO 40 jk = jj, jj + jb - 1
1676 DO 30 iblk = 1, mblks
1678 IF( iblk.EQ.1 )
THEN
1680 ELSE IF( iblk.EQ.mblks )
THEN
1688 DO 20 ik = ii, ii + ib - 1
1689 a( ik, jk ) = dcmplx( pb_drand( 0 ),
1695 IF( iblk.EQ.1 )
THEN
1699 CALL pb_jumpit( imuladd( 1, jmp_npimbloc ), ib1,
1706 CALL pb_jumpit( imuladd( 1, jmp_npmb ), ib1, ib0 )
1717 CALL pb_jumpit( imuladd( 1, jmp_col ), ib2, ib0 )
1728 IF( jblk.EQ.1 )
THEN
1732 CALL pb_jumpit( imuladd( 1, jmp_nqinbloc ), ib3, ib0 )
1738 CALL pb_jumpit( imuladd( 1, jmp_nqnb ), ib3, ib0 )
1751 ELSE IF( lsame( aform,
'T' ) )
THEN
1758 DO 90 iblk = 1, mblks
1760 IF( iblk.EQ.1 )
THEN
1762 ELSE IF( iblk.EQ.mblks )
THEN
1768 DO 80 ik = ii, ii + ib - 1
1772 DO 70 jblk = 1, nblks
1774 IF( jblk.EQ.1 )
THEN
1776 ELSE IF( jblk.EQ.nblks )
THEN
1784 DO 60 jk = jj, jj + jb - 1
1785 a( ik, jk ) = dcmplx( pb_drand( 0 ),
1791 IF( jblk.EQ.1 )
THEN
1795 CALL pb_jumpit( imuladd( 1, jmp_nqinbloc ), ib1,
1802 CALL pb_jumpit( imuladd( 1, jmp_nqnb ), ib1, ib0 )
1813 CALL pb_jumpit( imuladd( 1, jmp_row ), ib2, ib0 )
1824 IF( iblk.EQ.1 )
THEN
1828 CALL pb_jumpit( imuladd( 1, jmp_npimbloc ), ib3, ib0 )
1834 CALL pb_jumpit( imuladd( 1, jmp_npmb ), ib3, ib0 )
1847 ELSE IF( lsame( aform,
'S' ) )
THEN
1851 IF( lsame( uplo,
'L' ) )
THEN
1858 DO 170 jblk = 1, nblks
1860 IF( jblk.EQ.1 )
THEN
1863 ELSE IF( jblk.EQ.nblks )
THEN
1871 DO 160 jk = jj, jj + jb - 1
1876 DO 150 iblk = 1, mblks
1878 IF( iblk.EQ.1 )
THEN
1881 ELSE IF( iblk.EQ.mblks )
THEN
1891 IF( lcmtr.GT.upp )
THEN
1893 DO 100 ik = ii, ii + ib - 1
1894 dummy = dcmplx( pb_drand( 0 ),
1898 ELSE IF( lcmtr.GE.low )
THEN
1901 mnb =
max( 0, -lcmtr )
1903 IF( jtmp.LE.
min( mnb, jb ) )
THEN
1905 DO 110 ik = ii, ii + ib - 1
1906 a( ik, jk ) = dcmplx( pb_drand( 0 ),
1910 ELSE IF( ( jtmp.GE.( mnb + 1 ) ) .AND.
1911 $ ( jtmp.LE.
min( ib-lcmtr, jb ) ) )
THEN
1913 itmp = ii + jtmp + lcmtr - 1
1915 DO 120 ik = ii, itmp - 1
1916 dummy = dcmplx( pb_drand( 0 ),
1920 DO 130 ik = itmp, ii + ib - 1
1921 a( ik, jk ) = dcmplx( pb_drand( 0 ),
1929 DO 140 ik = ii, ii + ib - 1
1930 a( ik, jk ) = dcmplx( pb_drand( 0 ),
1938 IF( iblk.EQ.1 )
THEN
1942 lcmtr = lcmtr - jmp( jmp_npimbloc )
1943 CALL pb_jumpit( imuladd( 1, jmp_npimbloc ), ib1,
1950 lcmtr = lcmtr - jmp( jmp_npmb )
1951 CALL pb_jumpit( imuladd( 1, jmp_npmb ), ib1,
1963 CALL pb_jumpit( imuladd( 1, jmp_col ), ib2, ib0 )
1974 IF( jblk.EQ.1 )
THEN
1978 lcmtc = lcmtc + jmp( jmp_nqinbloc )
1979 CALL pb_jumpit( imuladd( 1, jmp_nqinbloc ), ib3, ib0 )
1985 lcmtc = lcmtc + jmp( jmp_nqnb )
1986 CALL pb_jumpit( imuladd( 1, jmp_nqnb ), ib3, ib0 )
2006 DO 250 iblk = 1, mblks
2008 IF( iblk.EQ.1 )
THEN
2011 ELSE IF( iblk.EQ.mblks )
THEN
2019 DO 240 ik = ii, ii + ib - 1
2024 DO 230 jblk = 1, nblks
2026 IF( jblk.EQ.1 )
THEN
2029 ELSE IF( jblk.EQ.nblks )
THEN
2039 IF( lcmtc.LT.low )
THEN
2041 DO 180 jk = jj, jj + jb - 1
2042 dummy = dcmplx( pb_drand( 0 ),
2046 ELSE IF( lcmtc.LE.upp )
THEN
2049 mnb =
max( 0, lcmtc )
2051 IF( itmp.LE.
min( mnb, ib ) )
THEN
2053 DO 190 jk = jj, jj + jb - 1
2054 a( ik, jk ) = dcmplx( pb_drand( 0 ),
2058 ELSE IF( ( itmp.GE.( mnb + 1 ) ) .AND.
2059 $ ( itmp.LE.
min( jb+lcmtc, ib ) ) )
THEN
2061 jtmp = jj + itmp - lcmtc - 1
2063 DO 200 jk = jj, jtmp - 1
2064 dummy = dcmplx( pb_drand( 0 ),
2068 DO 210 jk = jtmp, jj + jb - 1
2069 a( ik, jk ) = dcmplx( pb_drand( 0 ),
2077 DO 220 jk = jj, jj + jb - 1
2078 a( ik, jk ) = dcmplx( pb_drand( 0 ),
2086 IF( jblk.EQ.1 )
THEN
2090 lcmtc = lcmtc + jmp( jmp_nqinbloc )
2091 CALL pb_jumpit( imuladd( 1, jmp_nqinbloc ), ib1,
2098 lcmtc = lcmtc + jmp( jmp_nqnb )
2099 CALL pb_jumpit( imuladd( 1, jmp_nqnb ), ib1,
2111 CALL pb_jumpit( imuladd( 1, jmp_row ), ib2, ib0 )
2122 IF( iblk.EQ.1 )
THEN
2126 lcmtr = lcmtr - jmp( jmp_npimbloc )
2127 CALL pb_jumpit( imuladd( 1, jmp_npimbloc ), ib3, ib0 )
2133 lcmtr = lcmtr - jmp( jmp_npmb )
2134 CALL pb_jumpit( imuladd( 1, jmp_npmb ), ib3, ib0 )
2149 ELSE IF( lsame( aform,
'C' ) )
THEN
2156 DO 290 iblk = 1, mblks
2158 IF( iblk.EQ.1 )
THEN
2160 ELSE IF( iblk.EQ.mblks )
THEN
2166 DO 280 ik = ii, ii + ib - 1
2170 DO 270 jblk = 1, nblks
2172 IF( jblk.EQ.1 )
THEN
2174 ELSE IF( jblk.EQ.nblks )
THEN
2182 DO 260 jk = jj, jj + jb - 1
2183 a( ik, jk ) = dcmplx( pb_drand( 0 ),
2189 IF( jblk.EQ.1 )
THEN
2193 CALL pb_jumpit( imuladd( 1, jmp_nqinbloc ), ib1,
2200 CALL pb_jumpit( imuladd( 1, jmp_nqnb ), ib1,
2212 CALL pb_jumpit( imuladd( 1, jmp_row ), ib2, ib0 )
2223 IF( iblk.EQ.1 )
THEN
2227 CALL pb_jumpit( imuladd( 1, jmp_npimbloc ), ib3, ib0 )
2233 CALL pb_jumpit( imuladd( 1, jmp_npmb ), ib3, ib0 )
2246 ELSE IF( lsame( aform,
'H' ) )
THEN
2250 IF( lsame( uplo,
'L' ) )
THEN
2257 DO 370 jblk = 1, nblks
2259 IF( jblk.EQ.1 )
THEN
2262 ELSE IF( jblk.EQ.nblks )
THEN
2270 DO 360 jk = jj, jj + jb - 1
2275 DO 350 iblk = 1, mblks
2277 IF( iblk.EQ.1 )
THEN
2280 ELSE IF( iblk.EQ.mblks )
THEN
2290 IF( lcmtr.GT.upp )
THEN
2292 DO 300 ik = ii, ii + ib - 1
2293 dummy = dcmplx( pb_drand( 0 ),
2297 ELSE IF( lcmtr.GE.low )
THEN
2300 mnb =
max( 0, -lcmtr )
2302 IF( jtmp.LE.
min( mnb, jb ) )
THEN
2304 DO 310 ik = ii, ii + ib - 1
2305 a( ik, jk ) = dcmplx( pb_drand( 0 ),
2309 ELSE IF( ( jtmp.GE.( mnb + 1 ) ) .AND.
2310 $ ( jtmp.LE.
min( ib-lcmtr, jb ) ) )
THEN
2312 itmp = ii + jtmp + lcmtr - 1
2314 DO 320 ik = ii, itmp - 1
2315 dummy = dcmplx( pb_drand( 0 ),
2319 IF( itmp.LE.( ii + ib - 1 ) )
THEN
2320 dummy = dcmplx( pb_drand( 0 ),
2322 a( itmp, jk ) = dcmplx( dble( dummy ),
2326 DO 330 ik = itmp + 1, ii + ib - 1
2327 a( ik, jk ) = dcmplx( pb_drand( 0 ),
2335 DO 340 ik = ii, ii + ib - 1
2336 a( ik, jk ) = dcmplx( pb_drand( 0 ),
2344 IF( iblk.EQ.1 )
THEN
2348 lcmtr = lcmtr - jmp( jmp_npimbloc )
2349 CALL pb_jumpit( imuladd( 1, jmp_npimbloc ), ib1,
2356 lcmtr = lcmtr - jmp( jmp_npmb )
2357 CALL pb_jumpit( imuladd( 1, jmp_npmb ), ib1,
2369 CALL pb_jumpit( imuladd( 1, jmp_col ), ib2, ib0 )
2380 IF( jblk.EQ.1 )
THEN
2384 lcmtc = lcmtc + jmp( jmp_nqinbloc )
2385 CALL pb_jumpit( imuladd( 1, jmp_nqinbloc ), ib3, ib0 )
2391 lcmtc = lcmtc + jmp( jmp_nqnb )
2392 CALL pb_jumpit( imuladd( 1, jmp_nqnb ), ib3, ib0 )
2412 DO 450 iblk = 1, mblks
2414 IF( iblk.EQ.1 )
THEN
2417 ELSE IF( iblk.EQ.mblks )
THEN
2425 DO 440 ik = ii, ii + ib - 1
2430 DO 430 jblk = 1, nblks
2432 IF( jblk.EQ.1 )
THEN
2435 ELSE IF( jblk.EQ.nblks )
THEN
2445 IF( lcmtc.LT.low )
THEN
2447 DO 380 jk = jj, jj + jb - 1
2448 dummy = dcmplx( pb_drand( 0 ),
2452 ELSE IF( lcmtc.LE.upp )
THEN
2455 mnb =
max( 0, lcmtc )
2457 IF( itmp.LE.
min( mnb, ib ) )
THEN
2459 DO 390 jk = jj, jj + jb - 1
2460 a( ik, jk ) = dcmplx( pb_drand( 0 ),
2464 ELSE IF( ( itmp.GE.( mnb + 1 ) ) .AND.
2465 $ ( itmp.LE.
min( jb+lcmtc, ib ) ) )
THEN
2467 jtmp = jj + itmp - lcmtc - 1
2469 DO 400 jk = jj, jtmp - 1
2470 dummy = dcmplx( pb_drand( 0 ),
2474 IF( jtmp.LE.( jj + jb - 1 ) )
THEN
2475 dummy = dcmplx( pb_drand( 0 ),
2477 a( ik, jtmp ) = dcmplx( dble( dummy ),
2481 DO 410 jk = jtmp + 1, jj + jb - 1
2482 a( ik, jk ) = dcmplx( pb_drand( 0 ),
2490 DO 420 jk = jj, jj + jb - 1
2491 a( ik, jk ) = dcmplx( pb_drand( 0 ),
2499 IF( jblk.EQ.1 )
THEN
2503 lcmtc = lcmtc + jmp( jmp_nqinbloc )
2504 CALL pb_jumpit( imuladd( 1, jmp_nqinbloc ), ib1,
2511 lcmtc = lcmtc + jmp( jmp_nqnb )
2512 CALL pb_jumpit( imuladd( 1, jmp_nqnb ), ib1,
2524 CALL pb_jumpit( imuladd( 1, jmp_row ), ib2, ib0 )
2535 IF( iblk.EQ.1 )
THEN
2539 lcmtr = lcmtr - jmp( jmp_npimbloc )
2540 CALL pb_jumpit( imuladd( 1, jmp_npimbloc ), ib3, ib0 )
2546 lcmtr = lcmtr - jmp( jmp_npmb )
2547 CALL pb_jumpit( imuladd( 1, jmp_npmb ), ib3, ib0 )
2569 DOUBLE PRECISION FUNCTION pb_drand( IDUMM )
2615 DOUBLE PRECISION one, two
2616 PARAMETER ( one = 1.0d+0, two = 2.0d+0 )
2631 DOUBLE PRECISION FUNCTION pb_dran( IDUMM )
2676 DOUBLE PRECISION divfac, pow16
2677 PARAMETER ( divfac = 2.147483648d+9,
2678 $ pow16 = 6.5536d+4 )
2690 INTEGER iacs( 4 ), irand( 2 )
2691 common /rancom/ irand, iacs
2698 pb_dran = ( dble( irand( 1 ) ) + pow16 * dble( irand( 2 ) ) ) /
2701 CALL pb_lmul( irand, iacs, j )
2702 CALL pb_ladd( j, iacs( 3 ), irand )