ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pzmatgen.f
Go to the documentation of this file.
1  SUBROUTINE pzmatgen( ICTXT, AFORM, DIAG, M, N, MB, NB, A, LDA,
2  $ IAROW, IACOL, ISEED, IROFF, IRNUM, ICOFF,
3  $ ICNUM, MYROW, MYCOL, NPROW, NPCOL )
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 *
579  END
max
#define max(A, B)
Definition: pcgemr.c:180
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
pzmatgen
subroutine pzmatgen(ICTXT, AFORM, DIAG, M, N, MB, NB, A, LDA, IAROW, IACOL, ISEED, IROFF, IRNUM, ICOFF, ICNUM, MYROW, MYCOL, NPROW, NPCOL)
Definition: pzmatgen.f:4
pxerbla
subroutine pxerbla(ICTXT, SRNAME, INFO)
Definition: pxerbla.f:2
setran
subroutine setran(IRAN, IA, IC)
Definition: pmatgeninc.f:142