ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
psmatgen.f
Go to the documentation of this file.
1  SUBROUTINE psmatgen( 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  REAL A( LDA, * )
18 * ..
19 *
20 * Purpose
21 * =======
22 *
23 * PSMATGEN : Parallel Real Single 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) REAL, 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  REAL ONE, TWO
119  PARAMETER ( ONE = 1.0e+0, two = 2.0e+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 * ..
128 * .. Local Arrays ..
129  INTEGER IADD(2), IA1(2), IA2(2), IA3(2), IA4(2),
130  $ IA5(2), IB1(2), IB2(2), IB3(2), IC1(2), IC2(2),
131  $ ic3(2), ic4(2), ic5(2), iran1(2), iran2(2),
132  $ iran3(2), iran4(2), itmp1(2), itmp2(2),
133  $ itmp3(2), jseed(2), mult(2)
134 * ..
135 * .. External Subroutines ..
136  EXTERNAL jumpit, pxerbla, setran, xjumpm
137 * ..
138 * .. Intrinsic Functions ..
139  INTRINSIC abs, max, mod
140 * ..
141 * .. External Functions ..
142  LOGICAL LSAME
143  INTEGER ICEIL, NUMROC
144  REAL PSRAND
145  EXTERNAL iceil, numroc, lsame, psrand
146 * ..
147 * .. Executable Statements ..
148 *
149 * Test the input arguments
150 *
151  mp = numroc( m, mb, myrow, iarow, nprow )
152  nq = numroc( n, nb, mycol, iacol, npcol )
153  symm = lsame( aform, 'S' )
154  herm = lsame( aform, 'H' )
155  tran = lsame( aform, 'T' )
156 *
157  info = 0
158  IF( .NOT.lsame( diag, 'D' ) .AND.
159  $ .NOT.lsame( diag, 'N' ) ) THEN
160  info = 3
161  ELSE IF( symm.OR.herm ) THEN
162  IF( m.NE.n ) THEN
163  info = 5
164  ELSE IF( mb.NE.nb ) THEN
165  info = 7
166  END IF
167  ELSE IF( m.LT.0 ) THEN
168  info = 4
169  ELSE IF( n.LT.0 ) THEN
170  info = 5
171  ELSE IF( mb.LT.1 ) THEN
172  info = 6
173  ELSE IF( nb.LT.1 ) THEN
174  info = 7
175  ELSE IF( lda.LT.0 ) THEN
176  info = 9
177  ELSE IF( ( iarow.LT.0 ).OR.( iarow.GE.nprow ) ) THEN
178  info = 10
179  ELSE IF( ( iacol.LT.0 ).OR.( iacol.GE.npcol ) ) THEN
180  info = 11
181  ELSE IF( mod(iroff,mb).GT.0 ) THEN
182  info = 13
183  ELSE IF( irnum.GT.(mp-iroff) ) THEN
184  info = 14
185  ELSE IF( mod(icoff,nb).GT.0 ) THEN
186  info = 15
187  ELSE IF( icnum.GT.(nq-icoff) ) THEN
188  info = 16
189  ELSE IF( ( myrow.LT.0 ).OR.( myrow.GE.nprow ) ) THEN
190  info = 17
191  ELSE IF( ( mycol.LT.0 ).OR.( mycol.GE.npcol ) ) THEN
192  info = 18
193  END IF
194  IF( info.NE.0 ) THEN
195  CALL pxerbla( ictxt, 'PSMATGEN', info )
196  RETURN
197  END IF
198 *
199  mrrow = mod( nprow+myrow-iarow, nprow )
200  mrcol = mod( npcol+mycol-iacol, npcol )
201  npmb = nprow * mb
202  nqnb = npcol * nb
203  moff = iroff / mb
204  noff = icoff / nb
205  mend = iceil(irnum, mb) + moff
206  nend = iceil(icnum, nb) + noff
207 *
208  mult(1) = mult0
209  mult(2) = mult1
210  iadd(1) = iadd0
211  iadd(2) = iadd1
212  jseed(1) = iseed
213  jseed(2) = 0
214 *
215 * Symmetric or Hermitian matrix will be generated.
216 *
217  IF( symm.OR.herm ) THEN
218 *
219 * First, generate the lower triangular part (with diagonal block)
220 *
221  jump1 = 1
222  jump2 = npmb
223  jump3 = m
224  jump4 = nqnb
225  jump5 = nb
226  jump6 = mrcol
227  jump7 = mb*mrrow
228 *
229  CALL xjumpm( jump1, mult, iadd, jseed, iran1, ia1, ic1 )
230  CALL xjumpm( jump2, mult, iadd, iran1, itmp1, ia2, ic2 )
231  CALL xjumpm( jump3, mult, iadd, iran1, itmp1, ia3, ic3 )
232  CALL xjumpm( jump4, ia3, ic3, iran1, itmp1, ia4, ic4 )
233  CALL xjumpm( jump5, ia3, ic3, iran1, itmp1, ia5, ic5 )
234  CALL xjumpm( jump6, ia5, ic5, iran1, itmp3, itmp1, itmp2 )
235  CALL xjumpm( jump7, mult, iadd, itmp3, iran1, itmp1, itmp2 )
236  CALL xjumpm( noff, ia4, ic4, iran1, itmp1, itmp2, itmp3 )
237  CALL xjumpm( moff, ia2, ic2, itmp1, iran1, itmp2, itmp3 )
238  CALL setran( iran1, ia1, ic1 )
239 *
240  DO 10 i = 1, 2
241  ib1(i) = iran1(i)
242  ib2(i) = iran1(i)
243  ib3(i) = iran1(i)
244  10 CONTINUE
245 *
246  jk = 1
247  DO 80 ic = noff+1, nend
248  ioffc = ((ic-1)*npcol+mrcol) * nb
249  DO 70 i = 1, nb
250  IF( jk .GT. icnum ) GO TO 90
251 *
252  ik = 1
253  DO 50 ir = moff+1, mend
254  ioffr = ((ir-1)*nprow+mrrow) * mb
255 *
256  IF( ioffr .GT. ioffc ) THEN
257  DO 20 j = 1, mb
258  IF( ik .GT. irnum ) GO TO 60
259  a(ik,jk) = one - two*psrand(0)
260  ik = ik + 1
261  20 CONTINUE
262 *
263  ELSE IF( ioffc .EQ. ioffr ) THEN
264  ik = ik + i - 1
265  IF( ik .GT. irnum ) GO TO 60
266  DO 30 j = 1, i-1
267  a(ik,jk) = one - two*psrand(0)
268  30 CONTINUE
269  a(ik,jk) = one - two*psrand(0)
270  DO 40 j = 1, mb-i
271  IF( ik+j .GT. irnum ) GO TO 60
272  a(ik+j,jk) = one - two*psrand(0)
273  a(ik,jk+j) = a(ik+j,jk)
274  40 CONTINUE
275  ik = ik + mb - i + 1
276  ELSE
277  ik = ik + mb
278  END IF
279 *
280  CALL jumpit( ia2, ic2, ib1, iran2 )
281  ib1(1) = iran2(1)
282  ib1(2) = iran2(2)
283  50 CONTINUE
284 *
285  60 CONTINUE
286  jk = jk + 1
287  CALL jumpit( ia3, ic3, ib2, iran3 )
288  ib1(1) = iran3(1)
289  ib1(2) = iran3(2)
290  ib2(1) = iran3(1)
291  ib2(2) = iran3(2)
292  70 CONTINUE
293 *
294  CALL jumpit( ia4, ic4, ib3, iran4 )
295  ib1(1) = iran4(1)
296  ib1(2) = iran4(2)
297  ib2(1) = iran4(1)
298  ib2(2) = iran4(2)
299  ib3(1) = iran4(1)
300  ib3(2) = iran4(2)
301  80 CONTINUE
302 *
303 * Next, generate the upper triangular part.
304 *
305  90 CONTINUE
306  mult(1) = mult0
307  mult(2) = mult1
308  iadd(1) = iadd0
309  iadd(2) = iadd1
310  jseed(1) = iseed
311  jseed(2) = 0
312 *
313  jump1 = 1
314  jump2 = nqnb
315  jump3 = n
316  jump4 = npmb
317  jump5 = mb
318  jump6 = mrrow
319  jump7 = nb*mrcol
320 *
321  CALL xjumpm( jump1, mult, iadd, jseed, iran1, ia1, ic1 )
322  CALL xjumpm( jump2, mult, iadd, iran1, itmp1, ia2, ic2 )
323  CALL xjumpm( jump3, mult, iadd, iran1, itmp1, ia3, ic3 )
324  CALL xjumpm( jump4, ia3, ic3, iran1, itmp1, ia4, ic4 )
325  CALL xjumpm( jump5, ia3, ic3, iran1, itmp1, ia5, ic5 )
326  CALL xjumpm( jump6, ia5, ic5, iran1, itmp3, itmp1, itmp2 )
327  CALL xjumpm( jump7, mult, iadd, itmp3, iran1, itmp1, itmp2 )
328  CALL xjumpm( moff, ia4, ic4, iran1, itmp1, itmp2, itmp3 )
329  CALL xjumpm( noff, ia2, ic2, itmp1, iran1, itmp2, itmp3 )
330  CALL setran( iran1, ia1, ic1 )
331 *
332  DO 100 i = 1, 2
333  ib1(i) = iran1(i)
334  ib2(i) = iran1(i)
335  ib3(i) = iran1(i)
336  100 CONTINUE
337 *
338  ik = 1
339  DO 150 ir = moff+1, mend
340  ioffr = ((ir-1)*nprow+mrrow) * mb
341  DO 140 j = 1, mb
342  IF( ik .GT. irnum ) GO TO 160
343  jk = 1
344  DO 120 ic = noff+1, nend
345  ioffc = ((ic-1)*npcol+mrcol) * nb
346  IF( ioffc .GT. ioffr ) THEN
347  DO 110 i = 1, nb
348  IF( jk .GT. icnum ) GO TO 130
349  a(ik,jk) = one - two*psrand(0)
350  jk = jk + 1
351  110 CONTINUE
352  ELSE
353  jk = jk + nb
354  END IF
355  CALL jumpit( ia2, ic2, ib1, iran2 )
356  ib1(1) = iran2(1)
357  ib1(2) = iran2(2)
358  120 CONTINUE
359 *
360  130 CONTINUE
361  ik = ik + 1
362  CALL jumpit( ia3, ic3, ib2, iran3 )
363  ib1(1) = iran3(1)
364  ib1(2) = iran3(2)
365  ib2(1) = iran3(1)
366  ib2(2) = iran3(2)
367  140 CONTINUE
368 *
369  CALL jumpit( ia4, ic4, ib3, iran4 )
370  ib1(1) = iran4(1)
371  ib1(2) = iran4(2)
372  ib2(1) = iran4(1)
373  ib2(2) = iran4(2)
374  ib3(1) = iran4(1)
375  ib3(2) = iran4(2)
376  150 CONTINUE
377  160 CONTINUE
378 *
379 * (Conjugate) Transposed matrix A will be generated.
380 *
381  ELSE IF( tran .OR. lsame( aform, 'C' ) ) THEN
382 *
383  jump1 = 1
384  jump2 = nqnb
385  jump3 = n
386  jump4 = npmb
387  jump5 = mb
388  jump6 = mrrow
389  jump7 = nb*mrcol
390 *
391  CALL xjumpm( jump1, mult, iadd, jseed, iran1, ia1, ic1 )
392  CALL xjumpm( jump2, mult, iadd, iran1, itmp1, ia2, ic2 )
393  CALL xjumpm( jump3, mult, iadd, iran1, itmp1, ia3, ic3 )
394  CALL xjumpm( jump4, ia3, ic3, iran1, itmp1, ia4, ic4 )
395  CALL xjumpm( jump5, ia3, ic3, iran1, itmp1, ia5, ic5 )
396  CALL xjumpm( jump6, ia5, ic5, iran1, itmp3, itmp1, itmp2 )
397  CALL xjumpm( jump7, mult, iadd, itmp3, iran1, itmp1, itmp2 )
398  CALL xjumpm( moff, ia4, ic4, iran1, itmp1, itmp2, itmp3 )
399  CALL xjumpm( noff, ia2, ic2, itmp1, iran1, itmp2, itmp3 )
400  CALL setran( iran1, ia1, ic1 )
401 *
402  DO 170 i = 1, 2
403  ib1(i) = iran1(i)
404  ib2(i) = iran1(i)
405  ib3(i) = iran1(i)
406  170 CONTINUE
407 *
408  ik = 1
409  DO 220 ir = moff+1, mend
410  ioffr = ((ir-1)*nprow+mrrow) * mb
411  DO 210 j = 1, mb
412  IF( ik .GT. irnum ) GO TO 230
413  jk = 1
414  DO 190 ic = noff+1, nend
415  ioffc = ((ic-1)*npcol+mrcol) * nb
416  DO 180 i = 1, nb
417  IF( jk .GT. icnum ) GO TO 200
418  a(ik,jk) = one - two*psrand(0)
419  jk = jk + 1
420  180 CONTINUE
421  CALL jumpit( ia2, ic2, ib1, iran2 )
422  ib1(1) = iran2(1)
423  ib1(2) = iran2(2)
424  190 CONTINUE
425 *
426  200 CONTINUE
427  ik = ik + 1
428  CALL jumpit( ia3, ic3, ib2, iran3 )
429  ib1(1) = iran3(1)
430  ib1(2) = iran3(2)
431  ib2(1) = iran3(1)
432  ib2(2) = iran3(2)
433  210 CONTINUE
434 *
435  CALL jumpit( ia4, ic4, ib3, iran4 )
436  ib1(1) = iran4(1)
437  ib1(2) = iran4(2)
438  ib2(1) = iran4(1)
439  ib2(2) = iran4(2)
440  ib3(1) = iran4(1)
441  ib3(2) = iran4(2)
442  220 CONTINUE
443  230 CONTINUE
444 *
445 * A random matrix is generated.
446 *
447  ELSE
448 *
449  jump1 = 1
450  jump2 = npmb
451  jump3 = m
452  jump4 = nqnb
453  jump5 = nb
454  jump6 = mrcol
455  jump7 = mb*mrrow
456 *
457  CALL xjumpm( jump1, mult, iadd, jseed, iran1, ia1, ic1 )
458  CALL xjumpm( jump2, mult, iadd, iran1, itmp1, ia2, ic2 )
459  CALL xjumpm( jump3, mult, iadd, iran1, itmp1, ia3, ic3 )
460  CALL xjumpm( jump4, ia3, ic3, iran1, itmp1, ia4, ic4 )
461  CALL xjumpm( jump5, ia3, ic3, iran1, itmp1, ia5, ic5 )
462  CALL xjumpm( jump6, ia5, ic5, iran1, itmp3, itmp1, itmp2 )
463  CALL xjumpm( jump7, mult, iadd, itmp3, iran1, itmp1, itmp2 )
464  CALL xjumpm( noff, ia4, ic4, iran1, itmp1, itmp2, itmp3 )
465  CALL xjumpm( moff, ia2, ic2, itmp1, iran1, itmp2, itmp3 )
466  CALL setran( iran1, ia1, ic1 )
467 *
468  DO 240 i = 1, 2
469  ib1(i) = iran1(i)
470  ib2(i) = iran1(i)
471  ib3(i) = iran1(i)
472  240 CONTINUE
473 *
474  jk = 1
475  DO 290 ic = noff+1, nend
476  ioffc = ((ic-1)*npcol+mrcol) * nb
477  DO 280 i = 1, nb
478  IF( jk .GT. icnum ) GO TO 300
479  ik = 1
480  DO 260 ir = moff+1, mend
481  ioffr = ((ir-1)*nprow+mrrow) * mb
482  DO 250 j = 1, mb
483  IF( ik .GT. irnum ) GO TO 270
484  a(ik,jk) = one - two*psrand(0)
485  ik = ik + 1
486  250 CONTINUE
487  CALL jumpit( ia2, ic2, ib1, iran2 )
488  ib1(1) = iran2(1)
489  ib1(2) = iran2(2)
490  260 CONTINUE
491 *
492  270 CONTINUE
493  jk = jk + 1
494  CALL jumpit( ia3, ic3, ib2, iran3 )
495  ib1(1) = iran3(1)
496  ib1(2) = iran3(2)
497  ib2(1) = iran3(1)
498  ib2(2) = iran3(2)
499  280 CONTINUE
500 *
501  CALL jumpit( ia4, ic4, ib3, iran4 )
502  ib1(1) = iran4(1)
503  ib1(2) = iran4(2)
504  ib2(1) = iran4(1)
505  ib2(2) = iran4(2)
506  ib3(1) = iran4(1)
507  ib3(2) = iran4(2)
508  290 CONTINUE
509  300 CONTINUE
510  END IF
511 *
512 * Diagonally dominant matrix will be generated.
513 *
514  IF( lsame( diag, 'D' ) ) THEN
515  IF( mb.NE.nb ) THEN
516  WRITE(*,*) 'Diagonally dominant matrices with rowNB not'//
517  $ ' equal colNB is not supported!'
518  RETURN
519  END IF
520 *
521  maxmn = max(m, n)
522  jk = 1
523  DO 340 ic = noff+1, nend
524  ioffc = ((ic-1)*npcol+mrcol) * nb
525  ik = 1
526  DO 320 ir = moff+1, mend
527  ioffr = ((ir-1)*nprow+mrrow) * mb
528  IF( ioffc.EQ.ioffr ) THEN
529  DO 310 j = 0, mb-1
530  IF( ik .GT. irnum ) GO TO 330
531  a(ik,jk+j) = abs(a(ik,jk+j)) + maxmn
532  ik = ik + 1
533  310 CONTINUE
534  ELSE
535  ik = ik + mb
536  END IF
537  320 CONTINUE
538  330 CONTINUE
539  jk = jk + nb
540  340 CONTINUE
541  END IF
542 *
543  RETURN
544 *
545 * End of PSMATGEN
546 *
547  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
psmatgen
subroutine psmatgen(ICTXT, AFORM, DIAG, M, N, MB, NB, A, LDA, IAROW, IACOL, ISEED, IROFF, IRNUM, ICOFF, ICNUM, MYROW, MYCOL, NPROW, NPCOL)
Definition: psmatgen.f:4
pxerbla
subroutine pxerbla(ICTXT, SRNAME, INFO)
Definition: pxerbla.f:2
setran
subroutine setran(IRAN, IA, IC)
Definition: pmatgeninc.f:142