OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
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.NOT. IF( LSAME( DIAG, 'd.AND.' )
159.NOT. $ LSAME( DIAG, 'n' ) ) THEN
160 INFO = 3
161.OR. ELSE IF( SYMMHERM ) THEN
162.NE. IF( MN ) THEN
163 INFO = 5
164.NE. ELSE IF( MBNB ) THEN
165 INFO = 7
166 END IF
167.LT. ELSE IF( M0 ) THEN
168 INFO = 4
169.LT. ELSE IF( N0 ) THEN
170 INFO = 5
171.LT. ELSE IF( MB1 ) THEN
172 INFO = 6
173.LT. ELSE IF( NB1 ) THEN
174 INFO = 7
175.LT. ELSE IF( LDA0 ) THEN
176 INFO = 9
177.LT..OR..GE. ELSE IF( ( IAROW0 )( IAROWNPROW ) ) THEN
178 INFO = 10
179.LT..OR..GE. ELSE IF( ( IACOL0 )( IACOLNPCOL ) ) THEN
180 INFO = 11
181.GT. ELSE IF( MOD(IROFF,MB)0 ) THEN
182 INFO = 13
183.GT. ELSE IF( IRNUM(MP-IROFF) ) THEN
184 INFO = 14
185.GT. ELSE IF( MOD(ICOFF,NB)0 ) THEN
186 INFO = 15
187.GT. ELSE IF( ICNUM(NQ-ICOFF) ) THEN
188 INFO = 16
189.LT..OR..GE. ELSE IF( ( MYROW0 )( MYROWNPROW ) ) THEN
190 INFO = 17
191.LT..OR..GE. ELSE IF( ( MYCOL0 )( MYCOLNPCOL ) ) THEN
192 INFO = 18
193 END IF
194.NE. IF( INFO0 ) 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.OR. IF( SYMMHERM ) 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.GT. IF( JK ICNUM ) GO TO 90
251*
252 IK = 1
253 DO 50 IR = MOFF+1, MEND
254 IOFFR = ((IR-1)*NPROW+MRROW) * MB
255*
256.GT. IF( IOFFR IOFFC ) THEN
257 DO 20 J = 1, MB
258.GT. IF( IK IRNUM ) GO TO 60
259 A(IK,JK) = ONE - TWO*PSRAND(0)
260 IK = IK + 1
261 20 CONTINUE
262*
263.EQ. ELSE IF( IOFFC IOFFR ) THEN
264 IK = IK + I - 1
265.GT. IF( IK 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.GT. IF( IK+J 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.GT. IF( IK IRNUM ) GO TO 160
343 JK = 1
344 DO 120 IC = NOFF+1, NEND
345 IOFFC = ((IC-1)*NPCOL+MRCOL) * NB
346.GT. IF( IOFFC IOFFR ) THEN
347 DO 110 I = 1, NB
348.GT. IF( JK 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.OR. ELSE IF( TRAN 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.GT. IF( IK 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.GT. IF( JK 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.GT. IF( JK 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.GT. IF( IK 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.NE. IF( MBNB ) 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
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
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
max
#define max(a, b)
Definition macros.h:21
pxerbla
subroutine pxerbla(contxt, srname, info)
Definition mpi.f:1600