SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ pzmatgen()

subroutine pzmatgen ( integer  ictxt,
character*1  aform,
character*1  diag,
integer  m,
integer  n,
integer  mb,
integer  nb,
complex*16, dimension( lda, * )  a,
integer  lda,
integer  iarow,
integer  iacol,
integer  iseed,
integer  iroff,
integer  irnum,
integer  icoff,
integer  icnum,
integer  myrow,
integer  mycol,
integer  nprow,
integer  npcol 
)

Definition at line 1 of file pzmatgen.f.

4*
5* -- ScaLAPACK routine (version 1.7) --
6* University of Tennessee, Knoxville, Oak Ridge National Laboratory,
7* and University of California, Berkeley.
8* May 1, 1997
9*
10* .. Scalar Arguments ..
11 CHARACTER*1 AFORM, DIAG
12 INTEGER IACOL, IAROW, ICNUM, ICOFF, ICTXT, IRNUM,
13 $ IROFF, ISEED, LDA, M, MB, MYCOL, MYROW, N,
14 $ NB, NPCOL, NPROW
15* ..
16* .. Array Arguments ..
17 COMPLEX*16 A( LDA, * )
18* ..
19*
20* Purpose
21* =======
22*
23* PZMATGEN : Parallel Complex Double precision MATrix GENerator.
24* Generate (or regenerate) a distributed matrix A (or sub-matrix of A).
25*
26* Arguments
27* =========
28*
29* ICTXT (global input) INTEGER
30* The BLACS context handle, indicating the global context of
31* the operation. The context itself is global.
32*
33* AFORM (global input) CHARACTER*1
34* if AFORM = 'S' : A is returned is a symmetric matrix.
35* if AFORM = 'H' : A is returned is a Hermitian matrix.
36* if AFORM = 'T' : A is overwritten with the transpose of
37* what would normally be generated.
38* if AFORM = 'C' : A is overwritten with the conjugate trans-
39* pose of what would normally be generated.
40* otherwise a random matrix is generated.
41*
42* DIAG (global input) CHARACTER*1
43* if DIAG = 'D' : A is diagonally dominant.
44*
45* M (global input) INTEGER
46* The number of rows in the generated distributed matrix.
47*
48* N (global input) INTEGER
49* The number of columns in the generated distributed
50* matrix.
51*
52* MB (global input) INTEGER
53* The row blocking factor of the distributed matrix A.
54*
55* NB (global input) INTEGER
56* The column blocking factor of the distributed matrix A.
57*
58* A (local output) COMPLEX*16, pointer into the local memory
59* to an array of dimension ( LDA, * ) containing the local
60* pieces of the distributed matrix.
61*
62* LDA (local input) INTEGER
63* The leading dimension of the array containing the local
64* pieces of the distributed matrix A.
65*
66* IAROW (global input) INTEGER
67* The row processor coordinate which holds the first block
68* of the distributed matrix A.
69*
70* IACOL (global input) INTEGER
71* The column processor coordinate which holds the first
72* block of the distributed matrix A.
73*
74* ISEED (global input) INTEGER
75* The seed number to generate the distributed matrix A.
76*
77* IROFF (local input) INTEGER
78* The number of local rows of A that have already been
79* generated. It should be a multiple of MB.
80*
81* IRNUM (local input) INTEGER
82* The number of local rows to be generated.
83*
84* ICOFF (local input) INTEGER
85* The number of local columns of A that have already been
86* generated. It should be a multiple of NB.
87*
88* ICNUM (local input) INTEGER
89* The number of local columns to be generated.
90*
91* MYROW (local input) INTEGER
92* The row process coordinate of the calling process.
93*
94* MYCOL (local input) INTEGER
95* The column process coordinate of the calling process.
96*
97* NPROW (global input) INTEGER
98* The number of process rows in the grid.
99*
100* NPCOL (global input) INTEGER
101* The number of process columns in the grid.
102*
103* Notes
104* =====
105*
106* The code is originally developed by David Walker, ORNL,
107* and modified by Jaeyoung Choi, ORNL.
108*
109* Reference: G. Fox et al.
110* Section 12.3 of "Solving problems on concurrent processors Vol. I"
111*
112* =====================================================================
113*
114* .. Parameters ..
115 INTEGER MULT0, MULT1, IADD0, IADD1
116 parameter( mult0=20077, mult1=16838, iadd0=12345,
117 $ iadd1=0 )
118 DOUBLE PRECISION ONE, TWO, ZERO
119 parameter( one = 1.0d+0, two = 2.0d+0, zero = 0.0d+0 )
120* ..
121* .. Local Scalars ..
122 LOGICAL SYMM, HERM, TRAN
123 INTEGER I, IC, IK, INFO, IOFFC, IOFFR, IR, J, JK,
124 $ JUMP1, JUMP2, JUMP3, JUMP4, JUMP5, JUMP6,
125 $ JUMP7, MAXMN, MEND, MOFF, MP, MRCOL, MRROW,
126 $ NEND, NOFF, NPMB, NQ, NQNB
127 DOUBLE PRECISION DUMMY
128* ..
129* .. Local Arrays ..
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)
135* ..
136* .. External Subroutines ..
137 EXTERNAL jumpit, pxerbla, setran, xjumpm
138* ..
139* .. Intrinsic Functions ..
140 INTRINSIC abs, dble, dcmplx, dconjg, dimag, max, mod
141* ..
142* .. External Functions ..
143 LOGICAL LSAME
144 INTEGER ICEIL, NUMROC
145 DOUBLE PRECISION PDRAND
146 EXTERNAL iceil, numroc, lsame, pdrand
147* ..
148* .. Executable Statements ..
149*
150* Test the input arguments
151*
152 mp = numroc( m, mb, myrow, iarow, nprow )
153 nq = numroc( n, nb, mycol, iacol, npcol )
154 symm = lsame( aform, 'S' )
155 herm = lsame( aform, 'H' )
156 tran = lsame( aform, 'T' )
157*
158 info = 0
159 IF( .NOT.lsame( diag, 'D' ) .AND.
160 $ .NOT.lsame( diag, 'N' ) ) THEN
161 info = 3
162 ELSE IF( symm.OR.herm ) THEN
163 IF( m.NE.n ) THEN
164 info = 5
165 ELSE IF( mb.NE.nb ) THEN
166 info = 7
167 END IF
168 ELSE IF( m.LT.0 ) THEN
169 info = 4
170 ELSE IF( n.LT.0 ) THEN
171 info = 5
172 ELSE IF( mb.LT.1 ) THEN
173 info = 6
174 ELSE IF( nb.LT.1 ) THEN
175 info = 7
176 ELSE IF( lda.LT.0 ) THEN
177 info = 9
178 ELSE IF( ( iarow.LT.0 ).OR.( iarow.GE.nprow ) ) THEN
179 info = 10
180 ELSE IF( ( iacol.LT.0 ).OR.( iacol.GE.npcol ) ) THEN
181 info = 11
182 ELSE IF( mod(iroff,mb).GT.0 ) THEN
183 info = 13
184 ELSE IF( irnum.GT.(mp-iroff) ) THEN
185 info = 14
186 ELSE IF( mod(icoff,nb).GT.0 ) THEN
187 info = 15
188 ELSE IF( icnum.GT.(nq-icoff) ) THEN
189 info = 16
190 ELSE IF( ( myrow.LT.0 ).OR.( myrow.GE.nprow ) ) THEN
191 info = 17
192 ELSE IF( ( mycol.LT.0 ).OR.( mycol.GE.npcol ) ) THEN
193 info = 18
194 END IF
195 IF( info.NE.0 ) THEN
196 CALL pxerbla( ictxt, 'PZMATGEN', info )
197 RETURN
198 END IF
199*
200 mrrow = mod( nprow+myrow-iarow, nprow )
201 mrcol = mod( npcol+mycol-iacol, npcol )
202 npmb = nprow * mb
203 nqnb = npcol * nb
204 moff = iroff / mb
205 noff = icoff / nb
206 mend = iceil(irnum, mb) + moff
207 nend = iceil(icnum, nb) + noff
208*
209 mult(1) = mult0
210 mult(2) = mult1
211 iadd(1) = iadd0
212 iadd(2) = iadd1
213 jseed(1) = iseed
214 jseed(2) = 0
215*
216* Symmetric or Hermitian matrix will be generated.
217*
218 IF( symm.OR.herm ) THEN
219*
220* First, generate the lower triangular part (with diagonal block)
221*
222 jump1 = 1
223 jump2 = 2*npmb
224 jump3 = 2*m
225 jump4 = nqnb
226 jump5 = nb
227 jump6 = mrcol
228 jump7 = 2*mb*mrrow
229*
230 CALL xjumpm( jump1, mult, iadd, jseed, iran1, ia1, ic1 )
231 CALL xjumpm( jump2, mult, iadd, iran1, itmp1, ia2, ic2 )
232 CALL xjumpm( jump3, mult, iadd, iran1, itmp1, ia3, ic3 )
233 CALL xjumpm( jump4, ia3, ic3, iran1, itmp1, ia4, ic4 )
234 CALL xjumpm( jump5, ia3, ic3, iran1, itmp1, ia5, ic5 )
235 CALL xjumpm( jump6, ia5, ic5, iran1, itmp3, itmp1, itmp2 )
236 CALL xjumpm( jump7, mult, iadd, itmp3, iran1, itmp1, itmp2 )
237 CALL xjumpm( noff, ia4, ic4, iran1, itmp1, itmp2, itmp3 )
238 CALL xjumpm( moff, ia2, ic2, itmp1, iran1, itmp2, itmp3 )
239 CALL setran( iran1, ia1, ic1 )
240*
241 DO 10 i = 1, 2
242 ib1(i) = iran1(i)
243 ib2(i) = iran1(i)
244 ib3(i) = iran1(i)
245 10 CONTINUE
246*
247 jk = 1
248 DO 80 ic = noff+1, nend
249 ioffc = ((ic-1)*npcol+mrcol) * nb
250 DO 70 i = 1, nb
251 IF( jk .GT. icnum ) GO TO 90
252*
253 ik = 1
254 DO 50 ir = moff+1, mend
255 ioffr = ((ir-1)*nprow+mrrow) * mb
256*
257 IF( ioffr .GT. ioffc ) THEN
258 DO 20 j = 1, mb
259 IF( ik .GT. irnum ) GO TO 60
260 a(ik,jk) = dcmplx( one - two*pdrand(0),
261 $ one - two*pdrand(0) )
262 ik = ik + 1
263 20 CONTINUE
264*
265 ELSE IF( ioffc .EQ. ioffr ) THEN
266 ik = ik + i - 1
267 IF( ik .GT. irnum ) GO TO 60
268 DO 30 j = 1, i-1
269 a(ik,jk) = dcmplx( pdrand(0), pdrand(0) )
270 30 CONTINUE
271 IF( symm ) THEN
272 a(ik,jk) = dcmplx( one - two*pdrand(0),
273 $ one - two*pdrand(0) )
274 ELSE
275 a(ik,jk) = dcmplx( one - two*pdrand(0), zero )
276 dummy = pdrand(0)
277 END IF
278 DO 40 j = 1, mb-i
279 IF( ik+j .GT. irnum ) GO TO 60
280 a(ik+j,jk) = dcmplx( one - two*pdrand(0),
281 $ one - two*pdrand(0) )
282 IF( herm ) THEN
283 a(ik,jk+j) = dconjg( a(ik+j,jk) )
284 ELSE
285 a(ik,jk+j) = a(ik+j,jk)
286 END IF
287 40 CONTINUE
288 ik = ik + mb - i + 1
289 ELSE
290 ik = ik + mb
291 END IF
292*
293 CALL jumpit( ia2, ic2, ib1, iran2 )
294 ib1(1) = iran2(1)
295 ib1(2) = iran2(2)
296 50 CONTINUE
297*
298 60 CONTINUE
299 jk = jk + 1
300 CALL jumpit( ia3, ic3, ib2, iran3 )
301 ib1(1) = iran3(1)
302 ib1(2) = iran3(2)
303 ib2(1) = iran3(1)
304 ib2(2) = iran3(2)
305 70 CONTINUE
306*
307 CALL jumpit( ia4, ic4, ib3, iran4 )
308 ib1(1) = iran4(1)
309 ib1(2) = iran4(2)
310 ib2(1) = iran4(1)
311 ib2(2) = iran4(2)
312 ib3(1) = iran4(1)
313 ib3(2) = iran4(2)
314 80 CONTINUE
315*
316* Next, generate the upper triangular part.
317*
318 90 CONTINUE
319 mult(1) = mult0
320 mult(2) = mult1
321 iadd(1) = iadd0
322 iadd(2) = iadd1
323 jseed(1) = iseed
324 jseed(2) = 0
325*
326 jump1 = 1
327 jump2 = 2*nqnb
328 jump3 = 2*n
329 jump4 = npmb
330 jump5 = mb
331 jump6 = mrrow
332 jump7 = 2*nb*mrcol
333*
334 CALL xjumpm( jump1, mult, iadd, jseed, iran1, ia1, ic1 )
335 CALL xjumpm( jump2, mult, iadd, iran1, itmp1, ia2, ic2 )
336 CALL xjumpm( jump3, mult, iadd, iran1, itmp1, ia3, ic3 )
337 CALL xjumpm( jump4, ia3, ic3, iran1, itmp1, ia4, ic4 )
338 CALL xjumpm( jump5, ia3, ic3, iran1, itmp1, ia5, ic5 )
339 CALL xjumpm( jump6, ia5, ic5, iran1, itmp3, itmp1, itmp2 )
340 CALL xjumpm( jump7, mult, iadd, itmp3, iran1, itmp1, itmp2 )
341 CALL xjumpm( moff, ia4, ic4, iran1, itmp1, itmp2, itmp3 )
342 CALL xjumpm( noff, ia2, ic2, itmp1, iran1, itmp2, itmp3 )
343 CALL setran( iran1, ia1, ic1 )
344*
345 DO 100 i = 1, 2
346 ib1(i) = iran1(i)
347 ib2(i) = iran1(i)
348 ib3(i) = iran1(i)
349 100 CONTINUE
350*
351 ik = 1
352 DO 150 ir = moff+1, mend
353 ioffr = ((ir-1)*nprow+mrrow) * mb
354 DO 140 j = 1, mb
355 IF( ik .GT. irnum ) GO TO 160
356 jk = 1
357 DO 120 ic = noff+1, nend
358 ioffc = ((ic-1)*npcol+mrcol) * nb
359 IF( ioffc .GT. ioffr ) THEN
360 DO 110 i = 1, nb
361 IF( jk .GT. icnum ) GO TO 130
362 IF( symm ) THEN
363 a(ik,jk) = dcmplx( one - two*pdrand(0),
364 $ one - two*pdrand(0) )
365 ELSE
366 a(ik,jk) = dcmplx( one - two*pdrand(0),
367 $ two*pdrand(0) - one )
368 END IF
369 jk = jk + 1
370 110 CONTINUE
371 ELSE
372 jk = jk + nb
373 END IF
374 CALL jumpit( ia2, ic2, ib1, iran2 )
375 ib1(1) = iran2(1)
376 ib1(2) = iran2(2)
377 120 CONTINUE
378*
379 130 CONTINUE
380 ik = ik + 1
381 CALL jumpit( ia3, ic3, ib2, iran3 )
382 ib1(1) = iran3(1)
383 ib1(2) = iran3(2)
384 ib2(1) = iran3(1)
385 ib2(2) = iran3(2)
386 140 CONTINUE
387*
388 CALL jumpit( ia4, ic4, ib3, iran4 )
389 ib1(1) = iran4(1)
390 ib1(2) = iran4(2)
391 ib2(1) = iran4(1)
392 ib2(2) = iran4(2)
393 ib3(1) = iran4(1)
394 ib3(2) = iran4(2)
395 150 CONTINUE
396 160 CONTINUE
397*
398* (Conjugate) Transposed matrix A will be generated.
399*
400 ELSE IF( tran .OR. lsame( aform, 'C' ) ) THEN
401*
402 jump1 = 1
403 jump2 = 2*nqnb
404 jump3 = 2*n
405 jump4 = npmb
406 jump5 = mb
407 jump6 = mrrow
408 jump7 = 2*nb*mrcol
409*
410 CALL xjumpm( jump1, mult, iadd, jseed, iran1, ia1, ic1 )
411 CALL xjumpm( jump2, mult, iadd, iran1, itmp1, ia2, ic2 )
412 CALL xjumpm( jump3, mult, iadd, iran1, itmp1, ia3, ic3 )
413 CALL xjumpm( jump4, ia3, ic3, iran1, itmp1, ia4, ic4 )
414 CALL xjumpm( jump5, ia3, ic3, iran1, itmp1, ia5, ic5 )
415 CALL xjumpm( jump6, ia5, ic5, iran1, itmp3, itmp1, itmp2 )
416 CALL xjumpm( jump7, mult, iadd, itmp3, iran1, itmp1, itmp2 )
417 CALL xjumpm( moff, ia4, ic4, iran1, itmp1, itmp2, itmp3 )
418 CALL xjumpm( noff, ia2, ic2, itmp1, iran1, itmp2, itmp3 )
419 CALL setran( iran1, ia1, ic1 )
420*
421 DO 170 i = 1, 2
422 ib1(i) = iran1(i)
423 ib2(i) = iran1(i)
424 ib3(i) = iran1(i)
425 170 CONTINUE
426*
427 ik = 1
428 DO 220 ir = moff+1, mend
429 ioffr = ((ir-1)*nprow+mrrow) * mb
430 DO 210 j = 1, mb
431 IF( ik .GT. irnum ) GO TO 230
432 jk = 1
433 DO 190 ic = noff+1, nend
434 ioffc = ((ic-1)*npcol+mrcol) * nb
435 DO 180 i = 1, nb
436 IF( jk .GT. icnum ) GO TO 200
437 IF( tran ) THEN
438 a(ik,jk) = dcmplx( one - two*pdrand(0),
439 $ one - two*pdrand(0) )
440 ELSE
441 a(ik,jk) = dcmplx( one - two*pdrand(0),
442 $ two*pdrand(0) - one )
443 END IF
444 jk = jk + 1
445 180 CONTINUE
446 CALL jumpit( ia2, ic2, ib1, iran2 )
447 ib1(1) = iran2(1)
448 ib1(2) = iran2(2)
449 190 CONTINUE
450*
451 200 CONTINUE
452 ik = ik + 1
453 CALL jumpit( ia3, ic3, ib2, iran3 )
454 ib1(1) = iran3(1)
455 ib1(2) = iran3(2)
456 ib2(1) = iran3(1)
457 ib2(2) = iran3(2)
458 210 CONTINUE
459*
460 CALL jumpit( ia4, ic4, ib3, iran4 )
461 ib1(1) = iran4(1)
462 ib1(2) = iran4(2)
463 ib2(1) = iran4(1)
464 ib2(2) = iran4(2)
465 ib3(1) = iran4(1)
466 ib3(2) = iran4(2)
467 220 CONTINUE
468 230 CONTINUE
469*
470* A random matrix is generated.
471*
472 ELSE
473*
474 jump1 = 1
475 jump2 = 2*npmb
476 jump3 = 2*m
477 jump4 = nqnb
478 jump5 = nb
479 jump6 = mrcol
480 jump7 = 2*mb*mrrow
481*
482 CALL xjumpm( jump1, mult, iadd, jseed, iran1, ia1, ic1 )
483 CALL xjumpm( jump2, mult, iadd, iran1, itmp1, ia2, ic2 )
484 CALL xjumpm( jump3, mult, iadd, iran1, itmp1, ia3, ic3 )
485 CALL xjumpm( jump4, ia3, ic3, iran1, itmp1, ia4, ic4 )
486 CALL xjumpm( jump5, ia3, ic3, iran1, itmp1, ia5, ic5 )
487 CALL xjumpm( jump6, ia5, ic5, iran1, itmp3, itmp1, itmp2 )
488 CALL xjumpm( jump7, mult, iadd, itmp3, iran1, itmp1, itmp2 )
489 CALL xjumpm( noff, ia4, ic4, iran1, itmp1, itmp2, itmp3 )
490 CALL xjumpm( moff, ia2, ic2, itmp1, iran1, itmp2, itmp3 )
491 CALL setran( iran1, ia1, ic1 )
492*
493 DO 240 i = 1, 2
494 ib1(i) = iran1(i)
495 ib2(i) = iran1(i)
496 ib3(i) = iran1(i)
497 240 CONTINUE
498*
499 jk = 1
500 DO 290 ic = noff+1, nend
501 ioffc = ((ic-1)*npcol+mrcol) * nb
502 DO 280 i = 1, nb
503 IF( jk .GT. icnum ) GO TO 300
504 ik = 1
505 DO 260 ir = moff+1, mend
506 ioffr = ((ir-1)*nprow+mrrow) * mb
507 DO 250 j = 1, mb
508 IF( ik .GT. irnum ) GO TO 270
509 a(ik,jk) = dcmplx( one - two*pdrand(0),
510 $ one - two*pdrand(0) )
511 ik = ik + 1
512 250 CONTINUE
513 CALL jumpit( ia2, ic2, ib1, iran2 )
514 ib1(1) = iran2(1)
515 ib1(2) = iran2(2)
516 260 CONTINUE
517*
518 270 CONTINUE
519 jk = jk + 1
520 CALL jumpit( ia3, ic3, ib2, iran3 )
521 ib1(1) = iran3(1)
522 ib1(2) = iran3(2)
523 ib2(1) = iran3(1)
524 ib2(2) = iran3(2)
525 280 CONTINUE
526*
527 CALL jumpit( ia4, ic4, ib3, iran4 )
528 ib1(1) = iran4(1)
529 ib1(2) = iran4(2)
530 ib2(1) = iran4(1)
531 ib2(2) = iran4(2)
532 ib3(1) = iran4(1)
533 ib3(2) = iran4(2)
534 290 CONTINUE
535 300 CONTINUE
536 END IF
537*
538* Diagonally dominant matrix will be generated.
539*
540 IF( lsame( diag, 'D' ) ) THEN
541 IF( mb.NE.nb ) THEN
542 WRITE(*,*) 'Diagonally dominant matrices with rowNB not'//
543 $ ' equal colNB is not supported!'
544 RETURN
545 END IF
546*
547 maxmn = max(m, n)
548 jk = 1
549 DO 340 ic = noff+1, nend
550 ioffc = ((ic-1)*npcol+mrcol) * nb
551 ik = 1
552 DO 320 ir = moff+1, mend
553 ioffr = ((ir-1)*nprow+mrrow) * mb
554 IF( ioffc.EQ.ioffr ) THEN
555 DO 310 j = 0, mb-1
556 IF( ik .GT. irnum ) GO TO 330
557 IF( herm ) THEN
558 a(ik,jk+j) = dcmplx(
559 $ abs(dble(a(ik,jk+j)))+2*maxmn, zero )
560 ELSE
561 a(ik,jk+j)= dcmplx( abs(dble(a(ik,jk+j)))+maxmn,
562 $ abs(dimag(a(ik,jk+j)))+ maxmn )
563 END IF
564 ik = ik + 1
565 310 CONTINUE
566 ELSE
567 ik = ik + mb
568 END IF
569 320 CONTINUE
570 330 CONTINUE
571 jk = jk + nb
572 340 CONTINUE
573 END IF
574*
575 RETURN
576*
577* End of PZMATGEN
578*
pdrand
double precision function pdrand(idumm)
Definition pmatgeninc.f:272
jumpit
subroutine jumpit(mult, iadd, irann, iranm)
Definition pmatgeninc.f:183
xjumpm
subroutine xjumpm(jumpm, mult, iadd, irann, iranm, iam, icm)
Definition pmatgeninc.f:85
setran
subroutine setran(iran, ia, ic)
Definition pmatgeninc.f:142
iceil
integer function iceil(inum, idenom)
Definition iceil.f:2
numroc
integer function numroc(n, nb, iproc, isrcproc, nprocs)
Definition numroc.f:2
max
#define max(A, B)
Definition pcgemr.c:180
pxerbla
subroutine pxerbla(ictxt, srname, info)
Definition pxerbla.f:2
lsame
logical function lsame(ca, cb)
Definition tools.f:1724
Here is the call graph for this function:
Here is the caller graph for this function: