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

◆ pdmatgen()

subroutine pdmatgen ( integer  ictxt,
character*1  aform,
character*1  diag,
integer  m,
integer  n,
integer  mb,
integer  nb,
double precision, 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 pdmatgen.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 DOUBLE PRECISION A( LDA, * )
18* ..
19*
20* Purpose
21* =======
22*
23* PDMATGEN : Parallel Real 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) DOUBLE PRECISION, pointer into the local
59* memory to an array of dimension ( LDA, * ) containing the
60* local 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
119 parameter( one = 1.0d+0, two = 2.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* ..
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 DOUBLE PRECISION PDRAND
145 EXTERNAL iceil, numroc, lsame, pdrand
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, 'PDMATGEN', 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*pdrand(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*pdrand(0)
268 30 CONTINUE
269 a(ik,jk) = one - two*pdrand(0)
270 DO 40 j = 1, mb-i
271 IF( ik+j .GT. irnum ) GO TO 60
272 a(ik+j,jk) = one - two*pdrand(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*pdrand(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*pdrand(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*pdrand(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 PDMATGEN
546*
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: