1 SUBROUTINE pdmatgen2( ICTXT, AFORM, DIAG, M, N, MB, NB, A, LDA,
2 $ IAROW, IACOL, ISEED, IROFF, IRNUM, ICOFF,
3 $ ICNUM, MYROW, MYCOL, NPROW, NPCOL )
11 CHARACTER*1 AFORM, DIAG
12 INTEGER IACOL, IAROW, ICNUM, ICOFF, ICTXT, IRNUM,
13 $ iroff, iseed, lda, m, mb, mycol, myrow, n,
17 DOUBLE PRECISION A( LDA, * )
116 INTEGER MULT0, MULT1, IADD0, IADD1
117 PARAMETER ( MULT0=20077, mult1=16838, iadd0=12345,
119 DOUBLE PRECISION ONE, TWO, ZERO
120 PARAMETER ( ONE = 1.0d+0, two = 2.0d+0, zero = 0.0d+0 )
123 LOGICAL SYMM, HERM, TRAN, UPPR, RANDOM
124 INTEGER I, IC, IK, INFO, IOFFC, IOFFR, IR, J, JK,
125 $ jump1, jump2, jump3, jump4, jump5, jump6,
126 $ jump7, maxmn, mend, moff, mp, mrcol, mrrow,
127 $ nend, noff, npmb, nq, nqnb
130 INTEGER IADD(2), IA1(2), IA2(2), IA3(2), IA4(2),
131 $ IA5(2), IB1(2), IB2(2), IB3(2), IC1(2), IC2(2),
132 $ ic3(2), ic4(2), ic5(2), iran1(2), iran2(2),
133 $ iran3(2), iran4(2), itmp1(2), itmp2(2),
134 $ itmp3(2), jseed(2), mult(2)
140 INTRINSIC abs,
max, mod
144 INTEGER ICEIL, NUMROC
145 DOUBLE PRECISION PDRAND
146 EXTERNAL iceil, numroc, lsame, pdrand
152 mp = numroc( m, mb, myrow, iarow, nprow )
153 nq = numroc( n, nb, mycol, iacol, npcol )
154 symm = lsame( aform,
'S' )
155 uppr = lsame( aform,
'U' )
156 herm = lsame( aform,
'H' )
157 tran = lsame( aform,
'T' )
158 random = lsame( aform,
'R' )
161 IF( .NOT.( uppr.OR.symm.OR.herm.OR.tran.OR.random ) .AND.
162 $ .NOT.lsame( aform,
'C' ) .AND.
163 $ .NOT.lsame( aform,
'N' ) )
THEN
165 ELSE IF( .NOT.lsame( diag,
'D' ) .AND.
166 $ .NOT.lsame( diag,
'N' ) )
THEN
168 ELSE IF( uppr.OR.symm.OR.herm )
THEN
171 ELSE IF( mb.NE.nb )
THEN
174 ELSE IF( m.LT.0 )
THEN
176 ELSE IF( n.LT.0 )
THEN
178 ELSE IF( mb.LT.1 )
THEN
180 ELSE IF( nb.LT.1 )
THEN
182 ELSE IF( lda.LT.0 )
THEN
184 ELSE IF( ( iarow.LT.0 ).OR.( iarow.GE.nprow ) )
THEN
186 ELSE IF( ( iacol.LT.0 ).OR.( iacol.GE.npcol ) )
THEN
188 ELSE IF( mod(iroff,mb).GT.0 )
THEN
190 ELSE IF( irnum.GT.(mp-iroff) )
THEN
192 ELSE IF( mod(icoff,nb).GT.0 )
THEN
194 ELSE IF( icnum.GT.(nq-icoff) )
THEN
196 ELSE IF( ( myrow.LT.0 ).OR.( myrow.GE.nprow ) )
THEN
198 ELSE IF( ( mycol.LT.0 ).OR.( mycol.GE.npcol ) )
THEN
202 CALL pxerbla( ictxt,
'PDMATGEN2', info )
205 mrrow = mod( nprow+myrow-iarow, nprow )
206 mrcol = mod( npcol+mycol-iacol, npcol )
211 mend = iceil(irnum, mb) + moff
212 nend = iceil(icnum, nb) + noff
223 IF( symm.OR.herm )
THEN
235 CALL xjumpm( jump1, mult, iadd, jseed, iran1, ia1, ic1 )
236 CALL xjumpm( jump2, mult, iadd, iran1, itmp1, ia2, ic2 )
237 CALL xjumpm( jump3, mult, iadd, iran1, itmp1, ia3, ic3 )
238 CALL xjumpm( jump4, ia3, ic3, iran1, itmp1, ia4, ic4 )
239 CALL xjumpm( jump5, ia3, ic3, iran1, itmp1, ia5, ic5 )
240 CALL xjumpm( jump6, ia5, ic5, iran1, itmp3, itmp1, itmp2 )
241 CALL xjumpm( jump7, mult, iadd, itmp3, iran1, itmp1, itmp2 )
242 CALL xjumpm( noff, ia4, ic4, iran1, itmp1, itmp2, itmp3 )
243 CALL xjumpm( moff, ia2, ic2, itmp1, iran1, itmp2, itmp3 )
244 CALL setran( iran1, ia1, ic1 )
253 DO 80 ic = noff+1, nend
254 ioffc = ((ic-1)*npcol+mrcol) * nb
256 IF( jk .GT. icnum )
GO TO 90
259 DO 50 ir = moff+1, mend
260 ioffr = ((ir-1)*nprow+mrrow) * mb
262 IF( ioffr .GT. ioffc )
THEN
264 IF( ik .GT. irnum )
GO TO 60
265 a(ik,jk) = one - two*pdrand(0)
269 ELSE IF( ioffc .EQ. ioffr )
THEN
271 IF( ik .GT. irnum )
GO TO 60
273 a(ik,jk) = one - two*pdrand(0)
275 a(ik,jk) = one - two*pdrand(0)
277 IF( ik+j .GT. irnum )
GO TO 60
278 a(ik+j,jk) = one - two*pdrand(0)
279 a(ik,jk+j) = a(ik+j,jk)
286 CALL jumpit( ia2, ic2, ib1, iran2 )
293 CALL jumpit( ia3, ic3, ib2, iran3 )
300 CALL jumpit( ia4, ic4, ib3, iran4 )
327 CALL xjumpm( jump1, mult, iadd, jseed, iran1, ia1, ic1 )
328 CALL xjumpm( jump2, mult, iadd, iran1, itmp1, ia2, ic2 )
329 CALL xjumpm( jump3, mult, iadd, iran1, itmp1, ia3, ic3 )
330 CALL xjumpm( jump4, ia3, ic3, iran1, itmp1, ia4, ic4 )
331 CALL xjumpm( jump5, ia3, ic3, iran1, itmp1, ia5, ic5 )
332 CALL xjumpm( jump6, ia5, ic5, iran1, itmp3, itmp1, itmp2 )
333 CALL xjumpm( jump7, mult, iadd, itmp3, iran1, itmp1, itmp2 )
334 CALL xjumpm( moff, ia4, ic4, iran1, itmp1, itmp2, itmp3 )
335 CALL xjumpm( noff, ia2, ic2, itmp1, iran1, itmp2, itmp3 )
336 CALL setran( iran1, ia1, ic1 )
345 DO 150 ir = moff+1, mend
346 ioffr = ((ir-1)*nprow+mrrow) * mb
348 IF( ik .GT. irnum )
GO TO 160
350 DO 120 ic = noff+1, nend
351 ioffc = ((ic-1)*npcol+mrcol) * nb
352 IF( ioffc .GT. ioffr )
THEN
354 IF( jk .GT. icnum )
GO TO 130
355 a(ik,jk) = one - two*pdrand(0)
361 CALL jumpit( ia2, ic2, ib1, iran2 )
368 CALL jumpit( ia3, ic3, ib2, iran3 )
375 CALL jumpit( ia4, ic4, ib3, iran4 )
387 ELSE IF ( uppr )
THEN
396 CALL xjumpm( jump1, mult, iadd, jseed, iran1, ia1, ic1 )
397 CALL xjumpm( jump2, mult, iadd, iran1, itmp1, ia2, ic2 )
398 CALL xjumpm( jump3, mult, iadd, iran1, itmp1, ia3, ic3 )
399 CALL xjumpm( jump4, ia3, ic3, iran1, itmp1, ia4, ic4 )
400 CALL xjumpm( jump5, ia3, ic3, iran1, itmp1, ia5, ic5 )
401 CALL xjumpm( jump6, ia5, ic5, iran1, itmp3, itmp1, itmp2 )
402 CALL xjumpm( jump7, mult, iadd, itmp3, iran1, itmp1, itmp2 )
403 CALL xjumpm( noff, ia4, ic4, iran1, itmp1, itmp2, itmp3 )
404 CALL xjumpm( moff, ia2, ic2, itmp1, iran1, itmp2, itmp3 )
405 CALL setran( iran1, ia1, ic1 )
414 DO 8000 ic = noff+1, nend
415 ioffc = ((ic-1)*npcol+mrcol) * nb
417 IF( jk .GT. icnum )
GO TO 8000
420 DO 5000 ir = moff+1, mend
421 ioffr = ((ir-1)*nprow+mrrow) * mb
423 IF( ioffc .EQ. ioffr )
THEN
425 IF( ik .GT. irnum )
GO TO 6000
427 a(ik,jk) = one - two*pdrand(0)
429 a(ik,jk) = one - two*pdrand(0)
431 IF( ik+j .GT. irnum )
GO TO 6000
432 a(ik,jk+j) = one - two*pdrand(0)
439 CALL jumpit( ia2, ic2, ib1, iran2 )
446 CALL jumpit( ia3, ic3, ib2, iran3 )
453 CALL jumpit( ia4, ic4, ib3, iran4 )
476 CALL xjumpm( jump1, mult, iadd, jseed, iran1, ia1, ic1 )
477 CALL xjumpm( jump2, mult, iadd, iran1, itmp1, ia2, ic2 )
478 CALL xjumpm( jump3, mult, iadd, iran1, itmp1, ia3, ic3 )
479 CALL xjumpm( jump4, ia3, ic3, iran1, itmp1, ia4, ic4 )
480 CALL xjumpm( jump5, ia3, ic3, iran1, itmp1, ia5, ic5 )
481 CALL xjumpm( jump6, ia5, ic5, iran1, itmp3, itmp1, itmp2 )
482 CALL xjumpm( jump7, mult, iadd, itmp3, iran1, itmp1, itmp2 )
483 CALL xjumpm( moff, ia4, ic4, iran1, itmp1, itmp2, itmp3 )
484 CALL xjumpm( noff, ia2, ic2, itmp1, iran1, itmp2, itmp3 )
485 CALL setran( iran1, ia1, ic1 )
494 DO 1500 ir = moff+1, mend
495 ioffr = ((ir-1)*nprow+mrrow) * mb
497 IF( ik .GT. irnum )
GO TO 1600
499 DO 1200 ic = noff+1, nend
500 ioffc = ((ic-1)*npcol+mrcol) * nb
501 IF( ioffc .GT. ioffr )
THEN
503 IF( jk .GT. icnum )
GO TO 1300
504 a(ik,jk) = one - two*pdrand(0)
510 CALL jumpit( ia2, ic2, ib1, iran2 )
517 CALL jumpit( ia3, ic3, ib2, iran3 )
524 CALL jumpit( ia4, ic4, ib3, iran4 )
536 ELSE IF( tran .OR. lsame( aform,
'C' ) )
THEN
546 CALL xjumpm( jump1, mult, iadd, jseed, iran1, ia1, ic1 )
547 CALL xjumpm( jump2, mult, iadd, iran1, itmp1, ia2, ic2 )
548 CALL xjumpm( jump3, mult, iadd, iran1, itmp1, ia3, ic3 )
549 CALL xjumpm( jump4, ia3, ic3, iran1, itmp1, ia4, ic4 )
550 CALL xjumpm( jump5, ia3, ic3, iran1, itmp1, ia5, ic5 )
551 CALL xjumpm( jump6, ia5, ic5, iran1, itmp3, itmp1, itmp2 )
552 CALL xjumpm( jump7, mult, iadd, itmp3, iran1, itmp1, itmp2 )
553 CALL xjumpm( moff, ia4, ic4, iran1, itmp1, itmp2, itmp3 )
554 CALL xjumpm( noff, ia2, ic2, itmp1, iran1, itmp2, itmp3 )
555 CALL setran( iran1, ia1, ic1 )
564 DO 220 ir = moff+1, mend
565 ioffr = ((ir-1)*nprow+mrrow) * mb
567 IF( ik .GT. irnum )
GO TO 230
569 DO 190 ic = noff+1, nend
570 ioffc = ((ic-1)*npcol+mrcol) * nb
572 IF( jk .GT. icnum )
GO TO 200
573 a(ik,jk) = one - two*pdrand(0)
576 CALL jumpit( ia2, ic2, ib1, iran2 )
583 CALL jumpit( ia3, ic3, ib2, iran3 )
590 CALL jumpit( ia4, ic4, ib3, iran4 )
602 ELSEIF( random )
THEN
612 CALL xjumpm( jump1, mult, iadd, jseed, iran1, ia1, ic1 )
613 CALL xjumpm( jump2, mult, iadd, iran1, itmp1, ia2, ic2 )
614 CALL xjumpm( jump3, mult, iadd, iran1, itmp1, ia3, ic3 )
615 CALL xjumpm( jump4, ia3, ic3, iran1, itmp1, ia4, ic4 )
616 CALL xjumpm( jump5, ia3, ic3, iran1, itmp1, ia5, ic5 )
617 CALL xjumpm( jump6, ia5, ic5, iran1, itmp3, itmp1, itmp2 )
618 CALL xjumpm( jump7, mult, iadd, itmp3, iran1, itmp1, itmp2 )
619 CALL xjumpm( noff, ia4, ic4, iran1, itmp1, itmp2, itmp3 )
620 CALL xjumpm( moff, ia2, ic2, itmp1, iran1, itmp2, itmp3 )
621 CALL setran( iran1, ia1, ic1 )
630 DO 290 ic = noff+1, nend
631 ioffc = ((ic-1)*npcol+mrcol) * nb
633 IF( jk .GT. icnum )
GO TO 300
635 DO 260 ir = moff+1, mend
636 ioffr = ((ir-1)*nprow+mrrow) * mb
638 IF( ik .GT. irnum )
GO TO 270
639 a(ik,jk) = one - two*pdrand(0)
642 CALL jumpit( ia2, ic2, ib1, iran2 )
649 CALL jumpit( ia3, ic3, ib2, iran3 )
656 CALL jumpit( ia4, ic4, ib3, iran4 )
669 IF( lsame( diag,
'D' ) )
THEN
671 WRITE(*,*)
'Diagonally dominant matrices with rowNB not'//
672 $
' equal colNB is not supported!'
678 DO 340 ic = noff+1, nend
679 ioffc = ((ic-1)*npcol+mrcol) * nb
681 DO 320 ir = moff+1, mend
682 ioffr = ((ir-1)*nprow+mrrow) * mb
683 IF( ioffc.EQ.ioffr )
THEN
685 IF( ik .GT. irnum )
GO TO 330
686 a(ik,jk+j) = abs(a(ik,jk+j)) + maxmn