ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pdblastim.f
Go to the documentation of this file.
1  SUBROUTINE pdlascal( TYPE, M, N, ALPHA, A, IA, JA, DESCA )
2 *
3 * -- PBLAS test routine (version 2.0) --
4 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
5 * and University of California, Berkeley.
6 * April 1, 1998
7 *
8 * .. Scalar Arguments ..
9  CHARACTER*1 TYPE
10  INTEGER IA, JA, M, N
11  DOUBLE PRECISION ALPHA
12 * ..
13 * .. Array Arguments ..
14  INTEGER DESCA( * )
15  DOUBLE PRECISION A( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * PDLASCAL scales the m by n submatrix A(IA:IA+M-1,JA:JA+N-1) denoted
22 * by sub( A ) by the scalar alpha. TYPE specifies if sub( A ) is full,
23 * upper triangular, lower triangular or upper Hessenberg.
24 *
25 * Notes
26 * =====
27 *
28 * A description vector is associated with each 2D block-cyclicly dis-
29 * tributed matrix. This vector stores the information required to
30 * establish the mapping between a matrix entry and its corresponding
31 * process and memory location.
32 *
33 * In the following comments, the character _ should be read as
34 * "of the distributed matrix". Let A be a generic term for any 2D
35 * block cyclicly distributed matrix. Its description vector is DESCA:
36 *
37 * NOTATION STORED IN EXPLANATION
38 * ---------------- --------------- ------------------------------------
39 * DTYPE_A (global) DESCA( DTYPE_ ) The descriptor type.
40 * CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
41 * the NPROW x NPCOL BLACS process grid
42 * A is distributed over. The context
43 * itself is global, but the handle
44 * (the integer value) may vary.
45 * M_A (global) DESCA( M_ ) The number of rows in the distribu-
46 * ted matrix A, M_A >= 0.
47 * N_A (global) DESCA( N_ ) The number of columns in the distri-
48 * buted matrix A, N_A >= 0.
49 * IMB_A (global) DESCA( IMB_ ) The number of rows of the upper left
50 * block of the matrix A, IMB_A > 0.
51 * INB_A (global) DESCA( INB_ ) The number of columns of the upper
52 * left block of the matrix A,
53 * INB_A > 0.
54 * MB_A (global) DESCA( MB_ ) The blocking factor used to distri-
55 * bute the last M_A-IMB_A rows of A,
56 * MB_A > 0.
57 * NB_A (global) DESCA( NB_ ) The blocking factor used to distri-
58 * bute the last N_A-INB_A columns of
59 * A, NB_A > 0.
60 * RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
61 * row of the matrix A is distributed,
62 * NPROW > RSRC_A >= 0.
63 * CSRC_A (global) DESCA( CSRC_ ) The process column over which the
64 * first column of A is distributed.
65 * NPCOL > CSRC_A >= 0.
66 * LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
67 * array storing the local blocks of
68 * the distributed matrix A,
69 * IF( Lc( 1, N_A ) > 0 )
70 * LLD_A >= MAX( 1, Lr( 1, M_A ) )
71 * ELSE
72 * LLD_A >= 1.
73 *
74 * Let K be the number of rows of a matrix A starting at the global in-
75 * dex IA,i.e, A( IA:IA+K-1, : ). Lr( IA, K ) denotes the number of rows
76 * that the process of row coordinate MYROW ( 0 <= MYROW < NPROW ) would
77 * receive if these K rows were distributed over NPROW processes. If K
78 * is the number of columns of a matrix A starting at the global index
79 * JA, i.e, A( :, JA:JA+K-1, : ), Lc( JA, K ) denotes the number of co-
80 * lumns that the process MYCOL ( 0 <= MYCOL < NPCOL ) would receive if
81 * these K columns were distributed over NPCOL processes.
82 *
83 * The values of Lr() and Lc() may be determined via a call to the func-
84 * tion PB_NUMROC:
85 * Lr( IA, K ) = PB_NUMROC( K, IA, IMB_A, MB_A, MYROW, RSRC_A, NPROW )
86 * Lc( JA, K ) = PB_NUMROC( K, JA, INB_A, NB_A, MYCOL, CSRC_A, NPCOL )
87 *
88 * Arguments
89 * =========
90 *
91 * TYPE (global input) CHARACTER*1
92 * On entry, TYPE specifies the type of the input submatrix as
93 * follows:
94 * = 'L' or 'l': sub( A ) is a lower triangular matrix,
95 * = 'U' or 'u': sub( A ) is an upper triangular matrix,
96 * = 'H' or 'h': sub( A ) is an upper Hessenberg matrix,
97 * otherwise sub( A ) is a full matrix.
98 *
99 * M (global input) INTEGER
100 * On entry, M specifies the number of rows of the submatrix
101 * sub( A ). M must be at least zero.
102 *
103 * N (global input) INTEGER
104 * On entry, N specifies the number of columns of the submatrix
105 * sub( A ). N must be at least zero.
106 *
107 * ALPHA (global input) DOUBLE PRECISION
108 * On entry, ALPHA specifies the scalar alpha.
109 *
110 * A (local input/local output) DOUBLE PRECISION array
111 * On entry, A is an array of dimension (LLD_A, Ka), where Ka is
112 * at least Lc( 1, JA+N-1 ). Before entry, this array contains
113 * the local entries of the matrix A.
114 * On exit, the local entries of this array corresponding to the
115 * to the entries of the submatrix sub( A ) are overwritten by
116 * the local entries of the m by n scaled submatrix.
117 *
118 * IA (global input) INTEGER
119 * On entry, IA specifies A's global row index, which points to
120 * the beginning of the submatrix sub( A ).
121 *
122 * JA (global input) INTEGER
123 * On entry, JA specifies A's global column index, which points
124 * to the beginning of the submatrix sub( A ).
125 *
126 * DESCA (global and local input) INTEGER array
127 * On entry, DESCA is an integer array of dimension DLEN_. This
128 * is the array descriptor for the matrix A.
129 *
130 * -- Written on April 1, 1998 by
131 * Antoine Petitet, University of Tennessee, Knoxville 37996, USA.
132 *
133 * =====================================================================
134 *
135 * .. Parameters ..
136  INTEGER BLOCK_CYCLIC_2D_INB, CSRC_, CTXT_, DLEN_,
137  $ DTYPE_, IMB_, INB_, LLD_, MB_, M_, NB_, N_,
138  $ RSRC_
139  parameter( block_cyclic_2d_inb = 2, dlen_ = 11,
140  $ dtype_ = 1, ctxt_ = 2, m_ = 3, n_ = 4,
141  $ imb_ = 5, inb_ = 6, mb_ = 7, nb_ = 8,
142  $ rsrc_ = 9, csrc_ = 10, lld_ = 11 )
143 * ..
144 * .. Local Scalars ..
145  CHARACTER*1 UPLO
146  LOGICAL GODOWN, GOLEFT, LOWER, UPPER
147  INTEGER IACOL, IAROW, ICTXT, IIA, IIMAX, ILOW, IMB1,
148  $ IMBLOC, INB1, INBLOC, IOFFA, IOFFD, ITYPE,
149  $ IUPP, JJA, JJMAX, JOFFA, JOFFD, LCMT, LCMT00,
150  $ LDA, LMBLOC, LNBLOC, LOW, M1, MB, MBLKD, MBLKS,
151  $ MBLOC, MP, MRCOL, MRROW, MYCOL, MYROW, N1, NB,
152  $ NBLKD, NBLKS, NBLOC, NPCOL, NPROW, NQ, PMB,
153  $ QNB, TMP1, UPP
154 * ..
155 * .. Local Arrays ..
156  INTEGER DESCA2( DLEN_ )
157 * ..
158 * .. External Subroutines ..
159  EXTERNAL blacs_gridinfo, pb_ainfog2l, pb_binfo,
160  $ pb_desctrans, pb_dlascal, pb_infog2l
161 * ..
162 * .. External Functions ..
163  LOGICAL LSAME
164  INTEGER PB_NUMROC
165  EXTERNAL lsame, pb_numroc
166 * ..
167 * .. Intrinsic Functions ..
168  INTRINSIC min
169 * ..
170 * .. Executable Statements ..
171 *
172 * Convert descriptor
173 *
174  CALL pb_desctrans( desca, desca2 )
175 *
176 * Get grid parameters
177 *
178  ictxt = desca2( ctxt_ )
179  CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
180 *
181 * Quick return if possible
182 *
183  IF( m.EQ.0 .OR. n.EQ.0 )
184  $ RETURN
185 *
186  IF( lsame( TYPE, 'L' ) ) then
187  itype = 1
188  uplo = TYPE
189  upper = .false.
190  lower = .true.
191  ioffd = 0
192  ELSE IF( lsame( TYPE, 'U' ) ) then
193  itype = 2
194  uplo = TYPE
195  upper = .true.
196  lower = .false.
197  ioffd = 0
198  ELSE IF( lsame( TYPE, 'H' ) ) then
199  itype = 3
200  uplo = 'U'
201  upper = .true.
202  lower = .false.
203  ioffd = 1
204  ELSE
205  itype = 0
206  uplo = 'A'
207  upper = .true.
208  lower = .true.
209  ioffd = 0
210  END IF
211 *
212 * Compute local indexes
213 *
214  IF( itype.EQ.0 ) THEN
215 *
216 * Full matrix
217 *
218  CALL pb_infog2l( ia, ja, desca2, nprow, npcol, myrow, mycol,
219  $ iia, jja, iarow, iacol )
220  mp = pb_numroc( m, ia, desca2( imb_ ), desca2( mb_ ), myrow,
221  $ desca2( rsrc_ ), nprow )
222  nq = pb_numroc( n, ja, desca2( inb_ ), desca2( nb_ ), mycol,
223  $ desca2( csrc_ ), npcol )
224 *
225  IF( mp.LE.0 .OR. nq.LE.0 )
226  $ RETURN
227 *
228  lda = desca2( lld_ )
229  ioffa = iia + ( jja - 1 ) * lda
230 *
231  CALL pb_dlascal( 'All', mp, nq, 0, alpha, a( ioffa ), lda )
232 *
233  ELSE
234 *
235 * Trapezoidal matrix
236 *
237  CALL pb_ainfog2l( m, n, ia, ja, desca2, nprow, npcol, myrow,
238  $ mycol, imb1, inb1, mp, nq, iia, jja, iarow,
239  $ iacol, mrrow, mrcol )
240 *
241  IF( mp.LE.0 .OR. nq.LE.0 )
242  $ RETURN
243 *
244 * Initialize LCMT00, MBLKS, NBLKS, IMBLOC, INBLOC, LMBLOC,
245 * LNBLOC, ILOW, LOW, IUPP, and UPP.
246 *
247  mb = desca2( mb_ )
248  nb = desca2( nb_ )
249  lda = desca2( lld_ )
250 *
251  CALL pb_binfo( ioffd, mp, nq, imb1, inb1, mb, nb, mrrow,
252  $ mrcol, lcmt00, mblks, nblks, imbloc, inbloc,
253  $ lmbloc, lnbloc, ilow, low, iupp, upp )
254 *
255  m1 = mp
256  n1 = nq
257  ioffa = iia - 1
258  joffa = jja - 1
259  iimax = ioffa + mp
260  jjmax = joffa + nq
261 *
262  IF( desca2( rsrc_ ).LT.0 ) THEN
263  pmb = mb
264  ELSE
265  pmb = nprow * mb
266  END IF
267  IF( desca2( csrc_ ).LT.0 ) THEN
268  qnb = nb
269  ELSE
270  qnb = npcol * nb
271  END IF
272 *
273 * Handle the first block of rows or columns separately, and
274 * update LCMT00, MBLKS and NBLKS.
275 *
276  godown = ( lcmt00.GT.iupp )
277  goleft = ( lcmt00.LT.ilow )
278 *
279  IF( .NOT.godown .AND. .NOT.goleft ) THEN
280 *
281 * LCMT00 >= ILOW && LCMT00 <= IUPP
282 *
283  goleft = ( ( lcmt00 - ( iupp - upp + pmb ) ).LT.ilow )
284  godown = .NOT.goleft
285 *
286  CALL pb_dlascal( uplo, imbloc, inbloc, lcmt00, alpha,
287  $ a( iia+joffa*lda ), lda )
288  IF( godown ) THEN
289  IF( upper .AND. nq.GT.inbloc )
290  $ CALL pb_dlascal( 'All', imbloc, nq-inbloc, 0, alpha,
291  $ a( iia+(joffa+inbloc)*lda ), lda )
292  iia = iia + imbloc
293  m1 = m1 - imbloc
294  ELSE
295  IF( lower .AND. mp.GT.imbloc )
296  $ CALL pb_dlascal( 'All', mp-imbloc, inbloc, 0, alpha,
297  $ a( iia+imbloc+joffa*lda ), lda )
298  jja = jja + inbloc
299  n1 = n1 - inbloc
300  END IF
301 *
302  END IF
303 *
304  IF( godown ) THEN
305 *
306  lcmt00 = lcmt00 - ( iupp - upp + pmb )
307  mblks = mblks - 1
308  ioffa = ioffa + imbloc
309 *
310  10 CONTINUE
311  IF( mblks.GT.0 .AND. lcmt00.GT.upp ) THEN
312  lcmt00 = lcmt00 - pmb
313  mblks = mblks - 1
314  ioffa = ioffa + mb
315  GO TO 10
316  END IF
317 *
318  tmp1 = min( ioffa, iimax ) - iia + 1
319  IF( upper .AND. tmp1.GT.0 ) THEN
320  CALL pb_dlascal( 'All', tmp1, n1, 0, alpha,
321  $ a( iia+joffa*lda ), lda )
322  iia = iia + tmp1
323  m1 = m1 - tmp1
324  END IF
325 *
326  IF( mblks.LE.0 )
327  $ RETURN
328 *
329  lcmt = lcmt00
330  mblkd = mblks
331  ioffd = ioffa
332 *
333  mbloc = mb
334  20 CONTINUE
335  IF( mblkd.GT.0 .AND. lcmt.GE.ilow ) THEN
336  IF( mblkd.EQ.1 )
337  $ mbloc = lmbloc
338  CALL pb_dlascal( uplo, mbloc, inbloc, lcmt, alpha,
339  $ a( ioffd+1+joffa*lda ), lda )
340  lcmt00 = lcmt
341  lcmt = lcmt - pmb
342  mblks = mblkd
343  mblkd = mblkd - 1
344  ioffa = ioffd
345  ioffd = ioffd + mbloc
346  GO TO 20
347  END IF
348 *
349  tmp1 = m1 - ioffd + iia - 1
350  IF( lower .AND. tmp1.GT.0 )
351  $ CALL pb_dlascal( 'All', tmp1, inbloc, 0, alpha,
352  $ a( ioffd+1+joffa*lda ), lda )
353 *
354  tmp1 = ioffa - iia + 1
355  m1 = m1 - tmp1
356  n1 = n1 - inbloc
357  lcmt00 = lcmt00 + low - ilow + qnb
358  nblks = nblks - 1
359  joffa = joffa + inbloc
360 *
361  IF( upper .AND. tmp1.GT.0 .AND. n1.GT.0 )
362  $ CALL pb_dlascal( 'All', tmp1, n1, 0, alpha,
363  $ a( iia+joffa*lda ), lda )
364 *
365  iia = ioffa + 1
366  jja = joffa + 1
367 *
368  ELSE IF( goleft ) THEN
369 *
370  lcmt00 = lcmt00 + low - ilow + qnb
371  nblks = nblks - 1
372  joffa = joffa + inbloc
373 *
374  30 CONTINUE
375  IF( nblks.GT.0 .AND. lcmt00.LT.low ) THEN
376  lcmt00 = lcmt00 + qnb
377  nblks = nblks - 1
378  joffa = joffa + nb
379  GO TO 30
380  END IF
381 *
382  tmp1 = min( joffa, jjmax ) - jja + 1
383  IF( lower .AND. tmp1.GT.0 ) THEN
384  CALL pb_dlascal( 'All', m1, tmp1, 0, alpha,
385  $ a( iia+(jja-1)*lda ), lda )
386  jja = jja + tmp1
387  n1 = n1 - tmp1
388  END IF
389 *
390  IF( nblks.LE.0 )
391  $ RETURN
392 *
393  lcmt = lcmt00
394  nblkd = nblks
395  joffd = joffa
396 *
397  nbloc = nb
398  40 CONTINUE
399  IF( nblkd.GT.0 .AND. lcmt.LE.iupp ) THEN
400  IF( nblkd.EQ.1 )
401  $ nbloc = lnbloc
402  CALL pb_dlascal( uplo, imbloc, nbloc, lcmt, alpha,
403  $ a( iia+joffd*lda ), lda )
404  lcmt00 = lcmt
405  lcmt = lcmt + qnb
406  nblks = nblkd
407  nblkd = nblkd - 1
408  joffa = joffd
409  joffd = joffd + nbloc
410  GO TO 40
411  END IF
412 *
413  tmp1 = n1 - joffd + jja - 1
414  IF( upper .AND. tmp1.GT.0 )
415  $ CALL pb_dlascal( 'All', imbloc, tmp1, 0, alpha,
416  $ a( iia+joffd*lda ), lda )
417 *
418  tmp1 = joffa - jja + 1
419  m1 = m1 - imbloc
420  n1 = n1 - tmp1
421  lcmt00 = lcmt00 - ( iupp - upp + pmb )
422  mblks = mblks - 1
423  ioffa = ioffa + imbloc
424 *
425  IF( lower .AND. m1.GT.0 .AND. tmp1.GT.0 )
426  $ CALL pb_dlascal( 'All', m1, tmp1, 0, alpha,
427  $ a( ioffa+1+(jja-1)*lda ), lda )
428 *
429  iia = ioffa + 1
430  jja = joffa + 1
431 *
432  END IF
433 *
434  nbloc = nb
435  50 CONTINUE
436  IF( nblks.GT.0 ) THEN
437  IF( nblks.EQ.1 )
438  $ nbloc = lnbloc
439  60 CONTINUE
440  IF( mblks.GT.0 .AND. lcmt00.GT.upp ) THEN
441  lcmt00 = lcmt00 - pmb
442  mblks = mblks - 1
443  ioffa = ioffa + mb
444  GO TO 60
445  END IF
446 *
447  tmp1 = min( ioffa, iimax ) - iia + 1
448  IF( upper .AND. tmp1.GT.0 ) THEN
449  CALL pb_dlascal( 'All', tmp1, n1, 0, alpha,
450  $ a( iia+joffa*lda ), lda )
451  iia = iia + tmp1
452  m1 = m1 - tmp1
453  END IF
454 *
455  IF( mblks.LE.0 )
456  $ RETURN
457 *
458  lcmt = lcmt00
459  mblkd = mblks
460  ioffd = ioffa
461 *
462  mbloc = mb
463  70 CONTINUE
464  IF( mblkd.GT.0 .AND. lcmt.GE.low ) THEN
465  IF( mblkd.EQ.1 )
466  $ mbloc = lmbloc
467  CALL pb_dlascal( uplo, mbloc, nbloc, lcmt, alpha,
468  $ a( ioffd+1+joffa*lda ), lda )
469  lcmt00 = lcmt
470  lcmt = lcmt - pmb
471  mblks = mblkd
472  mblkd = mblkd - 1
473  ioffa = ioffd
474  ioffd = ioffd + mbloc
475  GO TO 70
476  END IF
477 *
478  tmp1 = m1 - ioffd + iia - 1
479  IF( lower .AND. tmp1.GT.0 )
480  $ CALL pb_dlascal( 'All', tmp1, nbloc, 0, alpha,
481  $ a( ioffd+1+joffa*lda ), lda )
482 *
483  tmp1 = min( ioffa, iimax ) - iia + 1
484  m1 = m1 - tmp1
485  n1 = n1 - nbloc
486  lcmt00 = lcmt00 + qnb
487  nblks = nblks - 1
488  joffa = joffa + nbloc
489 *
490  IF( upper .AND. tmp1.GT.0 .AND. n1.GT.0 )
491  $ CALL pb_dlascal( 'All', tmp1, n1, 0, alpha,
492  $ a( iia+joffa*lda ), lda )
493 *
494  iia = ioffa + 1
495  jja = joffa + 1
496 *
497  GO TO 50
498 *
499  END IF
500 *
501  END IF
502 *
503  RETURN
504 *
505 * End of PDLASCAL
506 *
507  END
508  SUBROUTINE pdlagen( INPLACE, AFORM, DIAG, OFFA, M, N, IA, JA,
509  $ DESCA, IASEED, A, LDA )
510 *
511 * -- PBLAS test routine (version 2.0) --
512 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
513 * and University of California, Berkeley.
514 * April 1, 1998
515 *
516 * .. Scalar Arguments ..
517  LOGICAL inplace
518  CHARACTER*1 aform, diag
519  INTEGER ia, iaseed, ja, lda, m, n, offa
520 * ..
521 * .. Array Arguments ..
522  INTEGER desca( * )
523  DOUBLE PRECISION A( LDA, * )
524 * ..
525 *
526 * Purpose
527 * =======
528 *
529 * PDLAGEN generates (or regenerates) a submatrix sub( A ) denoting
530 * A(IA:IA+M-1,JA:JA+N-1).
531 *
532 * Notes
533 * =====
534 *
535 * A description vector is associated with each 2D block-cyclicly dis-
536 * tributed matrix. This vector stores the information required to
537 * establish the mapping between a matrix entry and its corresponding
538 * process and memory location.
539 *
540 * In the following comments, the character _ should be read as
541 * "of the distributed matrix". Let A be a generic term for any 2D
542 * block cyclicly distributed matrix. Its description vector is DESCA:
543 *
544 * NOTATION STORED IN EXPLANATION
545 * ---------------- --------------- ------------------------------------
546 * DTYPE_A (global) DESCA( DTYPE_ ) The descriptor type.
547 * CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
548 * the NPROW x NPCOL BLACS process grid
549 * A is distributed over. The context
550 * itself is global, but the handle
551 * (the integer value) may vary.
552 * M_A (global) DESCA( M_ ) The number of rows in the distribu-
553 * ted matrix A, M_A >= 0.
554 * N_A (global) DESCA( N_ ) The number of columns in the distri-
555 * buted matrix A, N_A >= 0.
556 * IMB_A (global) DESCA( IMB_ ) The number of rows of the upper left
557 * block of the matrix A, IMB_A > 0.
558 * INB_A (global) DESCA( INB_ ) The number of columns of the upper
559 * left block of the matrix A,
560 * INB_A > 0.
561 * MB_A (global) DESCA( MB_ ) The blocking factor used to distri-
562 * bute the last M_A-IMB_A rows of A,
563 * MB_A > 0.
564 * NB_A (global) DESCA( NB_ ) The blocking factor used to distri-
565 * bute the last N_A-INB_A columns of
566 * A, NB_A > 0.
567 * RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
568 * row of the matrix A is distributed,
569 * NPROW > RSRC_A >= 0.
570 * CSRC_A (global) DESCA( CSRC_ ) The process column over which the
571 * first column of A is distributed.
572 * NPCOL > CSRC_A >= 0.
573 * LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
574 * array storing the local blocks of
575 * the distributed matrix A,
576 * IF( Lc( 1, N_A ) > 0 )
577 * LLD_A >= MAX( 1, Lr( 1, M_A ) )
578 * ELSE
579 * LLD_A >= 1.
580 *
581 * Let K be the number of rows of a matrix A starting at the global in-
582 * dex IA,i.e, A( IA:IA+K-1, : ). Lr( IA, K ) denotes the number of rows
583 * that the process of row coordinate MYROW ( 0 <= MYROW < NPROW ) would
584 * receive if these K rows were distributed over NPROW processes. If K
585 * is the number of columns of a matrix A starting at the global index
586 * JA, i.e, A( :, JA:JA+K-1, : ), Lc( JA, K ) denotes the number of co-
587 * lumns that the process MYCOL ( 0 <= MYCOL < NPCOL ) would receive if
588 * these K columns were distributed over NPCOL processes.
589 *
590 * The values of Lr() and Lc() may be determined via a call to the func-
591 * tion PB_NUMROC:
592 * Lr( IA, K ) = PB_NUMROC( K, IA, IMB_A, MB_A, MYROW, RSRC_A, NPROW )
593 * Lc( JA, K ) = PB_NUMROC( K, JA, INB_A, NB_A, MYCOL, CSRC_A, NPCOL )
594 *
595 * Arguments
596 * =========
597 *
598 * INPLACE (global input) LOGICAL
599 * On entry, INPLACE specifies if the matrix should be generated
600 * in place or not. If INPLACE is .TRUE., the local random array
601 * to be generated will start in memory at the local memory lo-
602 * cation A( 1, 1 ), otherwise it will start at the local posi-
603 * tion induced by IA and JA.
604 *
605 * AFORM (global input) CHARACTER*1
606 * On entry, AFORM specifies the type of submatrix to be genera-
607 * ted as follows:
608 * AFORM = 'S', sub( A ) is a symmetric matrix,
609 * AFORM = 'H', sub( A ) is a Hermitian matrix,
610 * AFORM = 'T', sub( A ) is overrwritten with the transpose
611 * of what would normally be generated,
612 * AFORM = 'C', sub( A ) is overwritten with the conjugate
613 * transpose of what would normally be genera-
614 * ted.
615 * AFORM = 'N', a random submatrix is generated.
616 *
617 * DIAG (global input) CHARACTER*1
618 * On entry, DIAG specifies if the generated submatrix is diago-
619 * nally dominant or not as follows:
620 * DIAG = 'D' : sub( A ) is diagonally dominant,
621 * DIAG = 'N' : sub( A ) is not diagonally dominant.
622 *
623 * OFFA (global input) INTEGER
624 * On entry, OFFA specifies the offdiagonal of the underlying
625 * matrix A(1:DESCA(M_),1:DESCA(N_)) of interest when the subma-
626 * trix is symmetric, Hermitian or diagonally dominant. OFFA = 0
627 * specifies the main diagonal, OFFA > 0 specifies a subdiago-
628 * nal, and OFFA < 0 specifies a superdiagonal (see further de-
629 * tails).
630 *
631 * M (global input) INTEGER
632 * On entry, M specifies the global number of matrix rows of the
633 * submatrix sub( A ) to be generated. M must be at least zero.
634 *
635 * N (global input) INTEGER
636 * On entry, N specifies the global number of matrix columns of
637 * the submatrix sub( A ) to be generated. N must be at least
638 * zero.
639 *
640 * IA (global input) INTEGER
641 * On entry, IA specifies A's global row index, which points to
642 * the beginning of the submatrix sub( A ).
643 *
644 * JA (global input) INTEGER
645 * On entry, JA specifies A's global column index, which points
646 * to the beginning of the submatrix sub( A ).
647 *
648 * DESCA (global and local input) INTEGER array
649 * On entry, DESCA is an integer array of dimension DLEN_. This
650 * is the array descriptor for the matrix A.
651 *
652 * IASEED (global input) INTEGER
653 * On entry, IASEED specifies the seed number to generate the
654 * matrix A. IASEED must be at least zero.
655 *
656 * A (local output) DOUBLE PRECISION array
657 * On entry, A is an array of dimension (LLD_A, Ka), where Ka is
658 * at least Lc( 1, JA+N-1 ). On exit, this array contains the
659 * local entries of the randomly generated submatrix sub( A ).
660 *
661 * LDA (local input) INTEGER
662 * On entry, LDA specifies the local leading dimension of the
663 * array A. When INPLACE is .FALSE., LDA is usually DESCA(LLD_).
664 * This restriction is however not enforced, and this subroutine
665 * requires only that LDA >= MAX( 1, Mp ) where
666 *
667 * Mp = PB_NUMROC( M, IA, IMB_A, MB_A, MYROW, RSRC_A, NPROW ).
668 *
669 * PB_NUMROC is a ScaLAPACK tool function; MYROW, MYCOL, NPROW
670 * and NPCOL can be determined by calling the BLACS subroutine
671 * BLACS_GRIDINFO.
672 *
673 * Further Details
674 * ===============
675 *
676 * OFFD is tied to the matrix described by DESCA, as opposed to the
677 * piece that is currently (re)generated. This is a global information
678 * independent from the distribution parameters. Below are examples of
679 * the meaning of OFFD for a global 7 by 5 matrix:
680 *
681 * ---------------------------------------------------------------------
682 * OFFD | 0 -1 -2 -3 -4 0 -1 -2 -3 -4 0 -1 -2 -3 -4
683 * -------|-------------------------------------------------------------
684 * | | OFFD=-1 | OFFD=0 OFFD=2
685 * | V V
686 * 0 | . d . . . -> d . . . . . . . . .
687 * 1 | . . d . . . d . . . . . . . .
688 * 2 | . . . d . . . d . . -> d . . . .
689 * 3 | . . . . d . . . d . . d . . .
690 * 4 | . . . . . . . . . d . . d . .
691 * 5 | . . . . . . . . . . . . . d .
692 * 6 | . . . . . . . . . . . . . . d
693 * ---------------------------------------------------------------------
694 *
695 * -- Written on April 1, 1998 by
696 * Antoine Petitet, University of Tennessee, Knoxville 37996, USA.
697 *
698 * =====================================================================
699 *
700 * .. Parameters ..
701  INTEGER BLOCK_CYCLIC_2D_INB, CSRC_, CTXT_, DLEN_,
702  $ dtype_, imb_, inb_, lld_, mb_, m_, nb_, n_,
703  $ rsrc_
704  parameter( block_cyclic_2d_inb = 2, dlen_ = 11,
705  $ dtype_ = 1, ctxt_ = 2, m_ = 3, n_ = 4,
706  $ imb_ = 5, inb_ = 6, mb_ = 7, nb_ = 8,
707  $ rsrc_ = 9, csrc_ = 10, lld_ = 11 )
708  INTEGER JMP_1, JMP_COL, JMP_IMBV, JMP_INBV, JMP_LEN,
709  $ jmp_mb, jmp_nb, jmp_npimbloc, jmp_npmb,
710  $ jmp_nqinbloc, jmp_nqnb, jmp_row
711  parameter( jmp_1 = 1, jmp_row = 2, jmp_col = 3,
712  $ jmp_mb = 4, jmp_imbv = 5, jmp_npmb = 6,
713  $ jmp_npimbloc = 7, jmp_nb = 8, jmp_inbv = 9,
714  $ jmp_nqnb = 10, jmp_nqinbloc = 11,
715  $ jmp_len = 11 )
716 * ..
717 * .. Local Scalars ..
718  LOGICAL DIAGDO, SYMM, HERM, NOTRAN
719  INTEGER CSRC, I, IACOL, IAROW, ICTXT, IIA, ILOCBLK,
720  $ ilocoff, ilow, imb, imb1, imbloc, imbvir, inb,
721  $ inb1, inbloc, inbvir, info, ioffda, itmp, iupp,
722  $ ivir, jja, jlocblk, jlocoff, jvir, lcmt00,
723  $ lmbloc, lnbloc, low, maxmn, mb, mblks, mp,
724  $ mrcol, mrrow, mycdist, mycol, myrdist, myrow,
725  $ nb, nblks, npcol, nprow, nq, nvir, rsrc, upp
726  DOUBLE PRECISION ALPHA
727 * ..
728 * .. Local Arrays ..
729  INTEGER DESCA2( DLEN_ ), IMULADD( 4, JMP_LEN ),
730  $ iran( 2 ), jmp( jmp_len ), muladd0( 4 )
731 * ..
732 * .. External Subroutines ..
733  EXTERNAL blacs_gridinfo, pb_ainfog2l, pb_binfo,
734  $ pb_chkmat, pb_desctrans, pb_dlagen, pb_initjmp,
735  $ pb_initmuladd, pb_jump, pb_jumpit, pb_locinfo,
736  $ pb_setlocran, pb_setran, pdladom, pxerbla
737 * ..
738 * .. External Functions ..
739  LOGICAL LSAME
740  EXTERNAL lsame
741 * ..
742 * .. Intrinsic Functions ..
743  INTRINSIC dble, max, min
744 * ..
745 * .. Data Statements ..
746  DATA ( muladd0( i ), i = 1, 4 ) / 20077, 16838,
747  $ 12345, 0 /
748 * ..
749 * .. Executable Statements ..
750 *
751 * Convert descriptor
752 *
753  CALL pb_desctrans( desca, desca2 )
754 *
755 * Test the input arguments
756 *
757  ictxt = desca2( ctxt_ )
758  CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
759 *
760 * Test the input parameters
761 *
762  info = 0
763  IF( nprow.EQ.-1 ) THEN
764  info = -( 1000 + ctxt_ )
765  ELSE
766  symm = lsame( aform, 'S' )
767  herm = lsame( aform, 'H' )
768  notran = lsame( aform, 'N' )
769  diagdo = lsame( diag, 'D' )
770  IF( .NOT.( symm.OR.herm.OR.notran ) .AND.
771  $ .NOT.( lsame( aform, 'T' ) ) .AND.
772  $ .NOT.( lsame( aform, 'C' ) ) ) THEN
773  info = -2
774  ELSE IF( ( .NOT.diagdo ) .AND.
775  $ ( .NOT.lsame( diag, 'N' ) ) ) THEN
776  info = -3
777  END IF
778  CALL pb_chkmat( ictxt, m, 5, n, 6, ia, ja, desca2, 10, info )
779  END IF
780 *
781  IF( info.NE.0 ) THEN
782  CALL pxerbla( ictxt, 'PDLAGEN', -info )
783  RETURN
784  END IF
785 *
786 * Quick return if possible
787 *
788  IF( ( m.LE.0 ).OR.( n.LE.0 ) )
789  $ RETURN
790 *
791 * Start the operations
792 *
793  mb = desca2( mb_ )
794  nb = desca2( nb_ )
795  imb = desca2( imb_ )
796  inb = desca2( inb_ )
797  rsrc = desca2( rsrc_ )
798  csrc = desca2( csrc_ )
799 *
800 * Figure out local information about the distributed matrix operand
801 *
802  CALL pb_ainfog2l( m, n, ia, ja, desca2, nprow, npcol, myrow,
803  $ mycol, imb1, inb1, mp, nq, iia, jja, iarow,
804  $ iacol, mrrow, mrcol )
805 *
806 * Decide where the entries shall be stored in memory
807 *
808  IF( inplace ) THEN
809  iia = 1
810  jja = 1
811  END IF
812 *
813 * Initialize LCMT00, MBLKS, NBLKS, IMBLOC, INBLOC, LMBLOC, LNBLOC,
814 * ILOW, LOW, IUPP, and UPP.
815 *
816  ioffda = ja + offa - ia
817  CALL pb_binfo( ioffda, mp, nq, imb1, inb1, mb, nb, mrrow,
818  $ mrcol, lcmt00, mblks, nblks, imbloc, inbloc,
819  $ lmbloc, lnbloc, ilow, low, iupp, upp )
820 *
821 * Initialize ILOCBLK, ILOCOFF, MYRDIST, JLOCBLK, JLOCOFF, MYCDIST
822 * This values correspond to the square virtual underlying matrix
823 * of size MAX( M_ + MAX( 0, -OFFA ), N_ + MAX( 0, OFFA ) ) used
824 * to set up the random sequence. For practical purposes, the size
825 * of this virtual matrix is upper bounded by M_ + N_ - 1.
826 *
827  itmp = max( 0, -offa )
828  ivir = ia + itmp
829  imbvir = imb + itmp
830  nvir = desca2( m_ ) + itmp
831 *
832  CALL pb_locinfo( ivir, imbvir, mb, myrow, rsrc, nprow, ilocblk,
833  $ ilocoff, myrdist )
834 *
835  itmp = max( 0, offa )
836  jvir = ja + itmp
837  inbvir = inb + itmp
838  nvir = max( max( nvir, desca2( n_ ) + itmp ),
839  $ desca2( m_ ) + desca2( n_ ) - 1 )
840 *
841  CALL pb_locinfo( jvir, inbvir, nb, mycol, csrc, npcol, jlocblk,
842  $ jlocoff, mycdist )
843 *
844  IF( symm .OR. herm .OR. notran ) THEN
845 *
846  CALL pb_initjmp( .true., nvir, imbvir, inbvir, imbloc, inbloc,
847  $ mb, nb, rsrc, csrc, nprow, npcol, 1, jmp )
848 *
849 * Compute constants to jump JMP( * ) numbers in the sequence
850 *
851  CALL pb_initmuladd( muladd0, jmp, imuladd )
852 *
853 * Compute and set the random value corresponding to A( IA, JA )
854 *
855  CALL pb_setlocran( iaseed, ilocblk, jlocblk, ilocoff, jlocoff,
856  $ myrdist, mycdist, nprow, npcol, jmp,
857  $ imuladd, iran )
858 *
859  CALL pb_dlagen( 'Lower', aform, a( iia, jja ), lda, lcmt00,
860  $ iran, mblks, imbloc, mb, lmbloc, nblks, inbloc,
861  $ nb, lnbloc, jmp, imuladd )
862 *
863  END IF
864 *
865  IF( symm .OR. herm .OR. ( .NOT. notran ) ) THEN
866 *
867  CALL pb_initjmp( .false., nvir, imbvir, inbvir, imbloc, inbloc,
868  $ mb, nb, rsrc, csrc, nprow, npcol, 1, jmp )
869 *
870 * Compute constants to jump JMP( * ) numbers in the sequence
871 *
872  CALL pb_initmuladd( muladd0, jmp, imuladd )
873 *
874 * Compute and set the random value corresponding to A( IA, JA )
875 *
876  CALL pb_setlocran( iaseed, ilocblk, jlocblk, ilocoff, jlocoff,
877  $ myrdist, mycdist, nprow, npcol, jmp,
878  $ imuladd, iran )
879 *
880  CALL pb_dlagen( 'Upper', aform, a( iia, jja ), lda, lcmt00,
881  $ iran, mblks, imbloc, mb, lmbloc, nblks, inbloc,
882  $ nb, lnbloc, jmp, imuladd )
883 *
884  END IF
885 *
886  IF( diagdo ) THEN
887 *
888  maxmn = max( desca2( m_ ), desca2( n_ ) )
889  alpha = dble( maxmn )
890 *
891  IF( ioffda.GE.0 ) THEN
892  CALL pdladom( inplace, min( max( 0, m-ioffda ), n ), alpha,
893  $ a, min( ia+ioffda, ia+m-1 ), ja, desca )
894  ELSE
895  CALL pdladom( inplace, min( m, max( 0, n+ioffda ) ), alpha,
896  $ a, ia, min( ja-ioffda, ja+n-1 ), desca )
897  END IF
898 *
899  END IF
900 *
901  RETURN
902 *
903 * End of PDLAGEN
904 *
905  END
906  SUBROUTINE pdladom( INPLACE, N, ALPHA, A, IA, JA, DESCA )
907 *
908 * -- PBLAS test routine (version 2.0) --
909 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
910 * and University of California, Berkeley.
911 * April 1, 1998
912 *
913 * .. Scalar Arguments ..
914  LOGICAL INPLACE
915  INTEGER IA, JA, N
916  DOUBLE PRECISION ALPHA
917 * ..
918 * .. Array Arguments ..
919  INTEGER DESCA( * )
920  DOUBLE PRECISION A( * )
921 * ..
922 *
923 * Purpose
924 * =======
925 *
926 * PDLADOM adds alpha to the diagonal entries of an n by n submatrix
927 * sub( A ) denoting A( IA:IA+N-1, JA:JA+N-1 ).
928 *
929 * Notes
930 * =====
931 *
932 * A description vector is associated with each 2D block-cyclicly dis-
933 * tributed matrix. This vector stores the information required to
934 * establish the mapping between a matrix entry and its corresponding
935 * process and memory location.
936 *
937 * In the following comments, the character _ should be read as
938 * "of the distributed matrix". Let A be a generic term for any 2D
939 * block cyclicly distributed matrix. Its description vector is DESCA:
940 *
941 * NOTATION STORED IN EXPLANATION
942 * ---------------- --------------- ------------------------------------
943 * DTYPE_A (global) DESCA( DTYPE_ ) The descriptor type.
944 * CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
945 * the NPROW x NPCOL BLACS process grid
946 * A is distributed over. The context
947 * itself is global, but the handle
948 * (the integer value) may vary.
949 * M_A (global) DESCA( M_ ) The number of rows in the distribu-
950 * ted matrix A, M_A >= 0.
951 * N_A (global) DESCA( N_ ) The number of columns in the distri-
952 * buted matrix A, N_A >= 0.
953 * IMB_A (global) DESCA( IMB_ ) The number of rows of the upper left
954 * block of the matrix A, IMB_A > 0.
955 * INB_A (global) DESCA( INB_ ) The number of columns of the upper
956 * left block of the matrix A,
957 * INB_A > 0.
958 * MB_A (global) DESCA( MB_ ) The blocking factor used to distri-
959 * bute the last M_A-IMB_A rows of A,
960 * MB_A > 0.
961 * NB_A (global) DESCA( NB_ ) The blocking factor used to distri-
962 * bute the last N_A-INB_A columns of
963 * A, NB_A > 0.
964 * RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
965 * row of the matrix A is distributed,
966 * NPROW > RSRC_A >= 0.
967 * CSRC_A (global) DESCA( CSRC_ ) The process column over which the
968 * first column of A is distributed.
969 * NPCOL > CSRC_A >= 0.
970 * LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
971 * array storing the local blocks of
972 * the distributed matrix A,
973 * IF( Lc( 1, N_A ) > 0 )
974 * LLD_A >= MAX( 1, Lr( 1, M_A ) )
975 * ELSE
976 * LLD_A >= 1.
977 *
978 * Let K be the number of rows of a matrix A starting at the global in-
979 * dex IA,i.e, A( IA:IA+K-1, : ). Lr( IA, K ) denotes the number of rows
980 * that the process of row coordinate MYROW ( 0 <= MYROW < NPROW ) would
981 * receive if these K rows were distributed over NPROW processes. If K
982 * is the number of columns of a matrix A starting at the global index
983 * JA, i.e, A( :, JA:JA+K-1, : ), Lc( JA, K ) denotes the number of co-
984 * lumns that the process MYCOL ( 0 <= MYCOL < NPCOL ) would receive if
985 * these K columns were distributed over NPCOL processes.
986 *
987 * The values of Lr() and Lc() may be determined via a call to the func-
988 * tion PB_NUMROC:
989 * Lr( IA, K ) = PB_NUMROC( K, IA, IMB_A, MB_A, MYROW, RSRC_A, NPROW )
990 * Lc( JA, K ) = PB_NUMROC( K, JA, INB_A, NB_A, MYCOL, CSRC_A, NPCOL )
991 *
992 * Arguments
993 * =========
994 *
995 * INPLACE (global input) LOGICAL
996 * On entry, INPLACE specifies if the matrix should be generated
997 * in place or not. If INPLACE is .TRUE., the local random array
998 * to be generated will start in memory at the local memory lo-
999 * cation A( 1, 1 ), otherwise it will start at the local posi-
1000 * tion induced by IA and JA.
1001 *
1002 * N (global input) INTEGER
1003 * On entry, N specifies the global order of the submatrix
1004 * sub( A ) to be modified. N must be at least zero.
1005 *
1006 * ALPHA (global input) DOUBLE PRECISION
1007 * On entry, ALPHA specifies the scalar alpha.
1008 *
1009 * A (local input/local output) DOUBLE PRECISION array
1010 * On entry, A is an array of dimension (LLD_A, Ka), where Ka is
1011 * at least Lc( 1, JA+N-1 ). Before entry, this array contains
1012 * the local entries of the matrix A. On exit, the local entries
1013 * of this array corresponding to the main diagonal of sub( A )
1014 * have been updated.
1015 *
1016 * IA (global input) INTEGER
1017 * On entry, IA specifies A's global row index, which points to
1018 * the beginning of the submatrix sub( A ).
1019 *
1020 * JA (global input) INTEGER
1021 * On entry, JA specifies A's global column index, which points
1022 * to the beginning of the submatrix sub( A ).
1023 *
1024 * DESCA (global and local input) INTEGER array
1025 * On entry, DESCA is an integer array of dimension DLEN_. This
1026 * is the array descriptor for the matrix A.
1027 *
1028 * -- Written on April 1, 1998 by
1029 * Antoine Petitet, University of Tennessee, Knoxville 37996, USA.
1030 *
1031 * =====================================================================
1032 *
1033 * .. Parameters ..
1034  INTEGER BLOCK_CYCLIC_2D_INB, CSRC_, CTXT_, DLEN_,
1035  $ dtype_, imb_, inb_, lld_, mb_, m_, nb_, n_,
1036  $ rsrc_
1037  parameter( block_cyclic_2d_inb = 2, dlen_ = 11,
1038  $ dtype_ = 1, ctxt_ = 2, m_ = 3, n_ = 4,
1039  $ imb_ = 5, inb_ = 6, mb_ = 7, nb_ = 8,
1040  $ rsrc_ = 9, csrc_ = 10, lld_ = 11 )
1041 * ..
1042 * .. Local Scalars ..
1043  LOGICAL GODOWN, GOLEFT
1044  INTEGER I, IACOL, IAROW, ICTXT, IIA, IJOFFA, ILOW,
1045  $ imb1, imbloc, inb1, inbloc, ioffa, ioffd, iupp,
1046  $ jja, joffa, joffd, lcmt, lcmt00, lda, ldap1,
1047  $ lmbloc, lnbloc, low, mb, mblkd, mblks, mbloc,
1048  $ mrcol, mrrow, mycol, myrow, nb, nblkd, nblks,
1049  $ nbloc, np, npcol, nprow, nq, pmb, qnb, upp
1050  DOUBLE PRECISION ATMP
1051 * ..
1052 * .. Local Scalars ..
1053  INTEGER DESCA2( DLEN_ )
1054 * ..
1055 * .. External Subroutines ..
1056  EXTERNAL blacs_gridinfo, pb_ainfog2l, pb_binfo,
1057  $ pb_desctrans
1058 * ..
1059 * .. Intrinsic Functions ..
1060  INTRINSIC abs, max, min
1061 * ..
1062 * .. Executable Statements ..
1063 *
1064 * Convert descriptor
1065 *
1066  CALL pb_desctrans( desca, desca2 )
1067 *
1068 * Get grid parameters
1069 *
1070  ictxt = desca2( ctxt_ )
1071  CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
1072 *
1073  IF( n.EQ.0 )
1074  $ RETURN
1075 *
1076  CALL pb_ainfog2l( n, n, ia, ja, desca2, nprow, npcol, myrow,
1077  $ mycol, imb1, inb1, np, nq, iia, jja, iarow,
1078  $ iacol, mrrow, mrcol )
1079 *
1080 * Decide where the entries shall be stored in memory
1081 *
1082  IF( inplace ) THEN
1083  iia = 1
1084  jja = 1
1085  END IF
1086 *
1087 * Initialize LCMT00, MBLKS, NBLKS, IMBLOC, INBLOC, LMBLOC, LNBLOC,
1088 * ILOW, LOW, IUPP, and UPP.
1089 *
1090  mb = desca2( mb_ )
1091  nb = desca2( nb_ )
1092 *
1093  CALL pb_binfo( 0, np, nq, imb1, inb1, mb, nb, mrrow, mrcol,
1094  $ lcmt00, mblks, nblks, imbloc, inbloc, lmbloc,
1095  $ lnbloc, ilow, low, iupp, upp )
1096 *
1097  ioffa = iia - 1
1098  joffa = jja - 1
1099  lda = desca2( lld_ )
1100  ldap1 = lda + 1
1101 *
1102  IF( desca2( rsrc_ ).LT.0 ) THEN
1103  pmb = mb
1104  ELSE
1105  pmb = nprow * mb
1106  END IF
1107  IF( desca2( csrc_ ).LT.0 ) THEN
1108  qnb = nb
1109  ELSE
1110  qnb = npcol * nb
1111  END IF
1112 *
1113 * Handle the first block of rows or columns separately, and update
1114 * LCMT00, MBLKS and NBLKS.
1115 *
1116  godown = ( lcmt00.GT.iupp )
1117  goleft = ( lcmt00.LT.ilow )
1118 *
1119  IF( .NOT.godown .AND. .NOT.goleft ) THEN
1120 *
1121 * LCMT00 >= ILOW && LCMT00 <= IUPP
1122 *
1123  IF( lcmt00.GE.0 ) THEN
1124  ijoffa = ioffa+lcmt00 + ( joffa - 1 ) * lda
1125  DO 10 i = 1, min( inbloc, max( 0, imbloc - lcmt00 ) )
1126  atmp = a( ijoffa + i*ldap1 )
1127  a( ijoffa + i*ldap1 ) = abs( atmp ) + alpha
1128  10 CONTINUE
1129  ELSE
1130  ijoffa = ioffa + ( joffa - lcmt00 - 1 ) * lda
1131  DO 20 i = 1, min( imbloc, max( 0, inbloc + lcmt00 ) )
1132  atmp = a( ijoffa + i*ldap1 )
1133  a( ijoffa + i*ldap1 ) = abs( atmp ) + alpha
1134  20 CONTINUE
1135  END IF
1136  goleft = ( ( lcmt00 - ( iupp - upp + pmb ) ).LT.ilow )
1137  godown = .NOT.goleft
1138 *
1139  END IF
1140 *
1141  IF( godown ) THEN
1142 *
1143  lcmt00 = lcmt00 - ( iupp - upp + pmb )
1144  mblks = mblks - 1
1145  ioffa = ioffa + imbloc
1146 *
1147  30 CONTINUE
1148  IF( mblks.GT.0 .AND. lcmt00.GT.upp ) THEN
1149  lcmt00 = lcmt00 - pmb
1150  mblks = mblks - 1
1151  ioffa = ioffa + mb
1152  GO TO 30
1153  END IF
1154 *
1155  lcmt = lcmt00
1156  mblkd = mblks
1157  ioffd = ioffa
1158 *
1159  mbloc = mb
1160  40 CONTINUE
1161  IF( mblkd.GT.0 .AND. lcmt.GE.ilow ) THEN
1162  IF( mblkd.EQ.1 )
1163  $ mbloc = lmbloc
1164  IF( lcmt.GE.0 ) THEN
1165  ijoffa = ioffd + lcmt + ( joffa - 1 ) * lda
1166  DO 50 i = 1, min( inbloc, max( 0, mbloc - lcmt ) )
1167  atmp = a( ijoffa + i*ldap1 )
1168  a( ijoffa + i*ldap1 ) = abs( atmp ) + alpha
1169  50 CONTINUE
1170  ELSE
1171  ijoffa = ioffd + ( joffa - lcmt - 1 ) * lda
1172  DO 60 i = 1, min( mbloc, max( 0, inbloc + lcmt ) )
1173  atmp = a( ijoffa + i*ldap1 )
1174  a( ijoffa + i*ldap1 ) = abs( atmp ) + alpha
1175  60 CONTINUE
1176  END IF
1177  lcmt00 = lcmt
1178  lcmt = lcmt - pmb
1179  mblks = mblkd
1180  mblkd = mblkd - 1
1181  ioffa = ioffd
1182  ioffd = ioffd + mbloc
1183  GO TO 40
1184  END IF
1185 *
1186  lcmt00 = lcmt00 + low - ilow + qnb
1187  nblks = nblks - 1
1188  joffa = joffa + inbloc
1189 *
1190  ELSE IF( goleft ) THEN
1191 *
1192  lcmt00 = lcmt00 + low - ilow + qnb
1193  nblks = nblks - 1
1194  joffa = joffa + inbloc
1195 *
1196  70 CONTINUE
1197  IF( nblks.GT.0 .AND. lcmt00.LT.low ) THEN
1198  lcmt00 = lcmt00 + qnb
1199  nblks = nblks - 1
1200  joffa = joffa + nb
1201  GO TO 70
1202  END IF
1203 *
1204  lcmt = lcmt00
1205  nblkd = nblks
1206  joffd = joffa
1207 *
1208  nbloc = nb
1209  80 CONTINUE
1210  IF( nblkd.GT.0 .AND. lcmt.LE.iupp ) THEN
1211  IF( nblkd.EQ.1 )
1212  $ nbloc = lnbloc
1213  IF( lcmt.GE.0 ) THEN
1214  ijoffa = ioffa + lcmt + ( joffd - 1 ) * lda
1215  DO 90 i = 1, min( nbloc, max( 0, imbloc - lcmt ) )
1216  atmp = a( ijoffa + i*ldap1 )
1217  a( ijoffa + i*ldap1 ) = abs( atmp ) + alpha
1218  90 CONTINUE
1219  ELSE
1220  ijoffa = ioffa + ( joffd - lcmt - 1 ) * lda
1221  DO 100 i = 1, min( imbloc, max( 0, nbloc + lcmt ) )
1222  atmp = a( ijoffa + i*ldap1 )
1223  a( ijoffa + i*ldap1 ) = abs( atmp ) + alpha
1224  100 CONTINUE
1225  END IF
1226  lcmt00 = lcmt
1227  lcmt = lcmt + qnb
1228  nblks = nblkd
1229  nblkd = nblkd - 1
1230  joffa = joffd
1231  joffd = joffd + nbloc
1232  GO TO 80
1233  END IF
1234 *
1235  lcmt00 = lcmt00 - ( iupp - upp + pmb )
1236  mblks = mblks - 1
1237  ioffa = ioffa + imbloc
1238 *
1239  END IF
1240 *
1241  nbloc = nb
1242  110 CONTINUE
1243  IF( nblks.GT.0 ) THEN
1244  IF( nblks.EQ.1 )
1245  $ nbloc = lnbloc
1246  120 CONTINUE
1247  IF( mblks.GT.0 .AND. lcmt00.GT.upp ) THEN
1248  lcmt00 = lcmt00 - pmb
1249  mblks = mblks - 1
1250  ioffa = ioffa + mb
1251  GO TO 120
1252  END IF
1253 *
1254  lcmt = lcmt00
1255  mblkd = mblks
1256  ioffd = ioffa
1257 *
1258  mbloc = mb
1259  130 CONTINUE
1260  IF( mblkd.GT.0 .AND. lcmt.GE.low ) THEN
1261  IF( mblkd.EQ.1 )
1262  $ mbloc = lmbloc
1263  IF( lcmt.GE.0 ) THEN
1264  ijoffa = ioffd + lcmt + ( joffa - 1 ) * lda
1265  DO 140 i = 1, min( nbloc, max( 0, mbloc - lcmt ) )
1266  atmp = a( ijoffa + i*ldap1 )
1267  a( ijoffa + i*ldap1 ) = abs( atmp ) + alpha
1268  140 CONTINUE
1269  ELSE
1270  ijoffa = ioffd + ( joffa - lcmt - 1 ) * lda
1271  DO 150 i = 1, min( mbloc, max( 0, nbloc + lcmt ) )
1272  atmp = a( ijoffa + i*ldap1 )
1273  a( ijoffa + i*ldap1 ) = abs( atmp ) + alpha
1274  150 CONTINUE
1275  END IF
1276  lcmt00 = lcmt
1277  lcmt = lcmt - pmb
1278  mblks = mblkd
1279  mblkd = mblkd - 1
1280  ioffa = ioffd
1281  ioffd = ioffd + mbloc
1282  GO TO 130
1283  END IF
1284 *
1285  lcmt00 = lcmt00 + qnb
1286  nblks = nblks - 1
1287  joffa = joffa + nbloc
1288  GO TO 110
1289 *
1290  END IF
1291 *
1292  RETURN
1293 *
1294 * End of PDLADOM
1295 *
1296  END
1297  SUBROUTINE pb_dlascal( UPLO, M, N, IOFFD, ALPHA, A, LDA )
1298 *
1299 * -- PBLAS test routine (version 2.0) --
1300 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
1301 * and University of California, Berkeley.
1302 * April 1, 1998
1303 *
1304 * .. Scalar Arguments ..
1305  CHARACTER*1 UPLO
1306  INTEGER IOFFD, LDA, M, N
1307  DOUBLE PRECISION ALPHA
1308 * ..
1309 * .. Array Arguments ..
1310  DOUBLE PRECISION A( LDA, * )
1311 * ..
1312 *
1313 * Purpose
1314 * =======
1315 *
1316 * PB_DLASCAL scales a two-dimensional array A by the scalar alpha.
1317 *
1318 * Arguments
1319 * =========
1320 *
1321 * UPLO (input) CHARACTER*1
1322 * On entry, UPLO specifies which trapezoidal part of the ar-
1323 * ray A is to be scaled as follows:
1324 * = 'L' or 'l': the lower trapezoid of A is scaled,
1325 * = 'U' or 'u': the upper trapezoid of A is scaled,
1326 * = 'D' or 'd': diagonal specified by IOFFD is scaled,
1327 * Otherwise: all of the array A is scaled.
1328 *
1329 * M (input) INTEGER
1330 * On entry, M specifies the number of rows of the array A. M
1331 * must be at least zero.
1332 *
1333 * N (input) INTEGER
1334 * On entry, N specifies the number of columns of the array A.
1335 * N must be at least zero.
1336 *
1337 * IOFFD (input) INTEGER
1338 * On entry, IOFFD specifies the position of the offdiagonal de-
1339 * limiting the upper and lower trapezoidal part of A as follows
1340 * (see the notes below):
1341 *
1342 * IOFFD = 0 specifies the main diagonal A( i, i ),
1343 * with i = 1 ... MIN( M, N ),
1344 * IOFFD > 0 specifies the subdiagonal A( i+IOFFD, i ),
1345 * with i = 1 ... MIN( M-IOFFD, N ),
1346 * IOFFD < 0 specifies the superdiagonal A( i, i-IOFFD ),
1347 * with i = 1 ... MIN( M, N+IOFFD ).
1348 *
1349 * ALPHA (input) DOUBLE PRECISION
1350 * On entry, ALPHA specifies the scalar alpha.
1351 *
1352 * A (input/output) DOUBLE PRECISION array
1353 * On entry, A is an array of dimension (LDA,N). Before entry
1354 * with UPLO = 'U' or 'u', the leading m by n part of the array
1355 * A must contain the upper trapezoidal part of the matrix as
1356 * specified by IOFFD to be scaled, and the strictly lower tra-
1357 * pezoidal part of A is not referenced; When UPLO = 'L' or 'l',
1358 * the leading m by n part of the array A must contain the lower
1359 * trapezoidal part of the matrix as specified by IOFFD to be
1360 * scaled, and the strictly upper trapezoidal part of A is not
1361 * referenced. On exit, the entries of the trapezoid part of A
1362 * determined by UPLO and IOFFD are scaled.
1363 *
1364 * LDA (input) INTEGER
1365 * On entry, LDA specifies the leading dimension of the array A.
1366 * LDA must be at least max( 1, M ).
1367 *
1368 * Notes
1369 * =====
1370 * N N
1371 * ---------------------------- -----------
1372 * | d | | |
1373 * M | d 'U' | | 'U' |
1374 * | 'L' 'D' | |d |
1375 * | d | M | d |
1376 * ---------------------------- | 'D' |
1377 * | d |
1378 * IOFFD < 0 | 'L' d |
1379 * | d|
1380 * N | |
1381 * ----------- -----------
1382 * | d 'U'|
1383 * | d | IOFFD > 0
1384 * M | 'D' |
1385 * | d| N
1386 * | 'L' | ----------------------------
1387 * | | | 'U' |
1388 * | | |d |
1389 * | | | 'D' |
1390 * | | | d |
1391 * | | |'L' d |
1392 * ----------- ----------------------------
1393 *
1394 * -- Written on April 1, 1998 by
1395 * Antoine Petitet, University of Tennessee, Knoxville 37996, USA.
1396 *
1397 * =====================================================================
1398 *
1399 * .. Local Scalars ..
1400  INTEGER I, J, JTMP, MN
1401 * ..
1402 * .. External Functions ..
1403  LOGICAL LSAME
1404  EXTERNAL lsame
1405 * ..
1406 * .. Intrinsic Functions ..
1407  INTRINSIC max, min
1408 * ..
1409 * .. Executable Statements ..
1410 *
1411 * Quick return if possible
1412 *
1413  IF( m.LE.0 .OR. n.LE.0 )
1414  $ RETURN
1415 *
1416 * Start the operations
1417 *
1418  IF( lsame( uplo, 'L' ) ) THEN
1419 *
1420 * Scales the lower triangular part of the array by ALPHA.
1421 *
1422  mn = max( 0, -ioffd )
1423  DO 20 j = 1, min( mn, n )
1424  DO 10 i = 1, m
1425  a( i, j ) = alpha * a( i, j )
1426  10 CONTINUE
1427  20 CONTINUE
1428  DO 40 j = mn + 1, min( m - ioffd, n )
1429  DO 30 i = j + ioffd, m
1430  a( i, j ) = alpha * a( i, j )
1431  30 CONTINUE
1432  40 CONTINUE
1433 *
1434  ELSE IF( lsame( uplo, 'U' ) ) THEN
1435 *
1436 * Scales the upper triangular part of the array by ALPHA.
1437 *
1438  mn = min( m - ioffd, n )
1439  DO 60 j = max( 0, -ioffd ) + 1, mn
1440  DO 50 i = 1, j + ioffd
1441  a( i, j ) = alpha * a( i, j )
1442  50 CONTINUE
1443  60 CONTINUE
1444  DO 80 j = max( 0, mn ) + 1, n
1445  DO 70 i = 1, m
1446  a( i, j ) = alpha * a( i, j )
1447  70 CONTINUE
1448  80 CONTINUE
1449 *
1450  ELSE IF( lsame( uplo, 'D' ) ) THEN
1451 *
1452 * Scales the diagonal entries by ALPHA.
1453 *
1454  DO 90 j = max( 0, -ioffd ) + 1, min( m - ioffd, n )
1455  jtmp = j + ioffd
1456  a( jtmp, j ) = alpha * a( jtmp, j )
1457  90 CONTINUE
1458 *
1459  ELSE
1460 *
1461 * Scales the entire array by ALPHA.
1462 *
1463  DO 110 j = 1, n
1464  DO 100 i = 1, m
1465  a( i, j ) = alpha * a( i, j )
1466  100 CONTINUE
1467  110 CONTINUE
1468 *
1469  END IF
1470 *
1471  RETURN
1472 *
1473 * End of PB_DLASCAL
1474 *
1475  END
1476  SUBROUTINE pb_dlagen( UPLO, AFORM, A, LDA, LCMT00, IRAN, MBLKS,
1477  $ IMBLOC, MB, LMBLOC, NBLKS, INBLOC, NB,
1478  $ LNBLOC, JMP, IMULADD )
1479 *
1480 * -- PBLAS test routine (version 2.0) --
1481 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
1482 * and University of California, Berkeley.
1483 * April 1, 1998
1484 *
1485 * .. Scalar Arguments ..
1486  CHARACTER*1 UPLO, AFORM
1487  INTEGER IMBLOC, INBLOC, LCMT00, LDA, LMBLOC, LNBLOC,
1488  $ MB, MBLKS, NB, NBLKS
1489 * ..
1490 * .. Array Arguments ..
1491  INTEGER IMULADD( 4, * ), IRAN( * ), JMP( * )
1492  DOUBLE PRECISION A( LDA, * )
1493 * ..
1494 *
1495 * Purpose
1496 * =======
1497 *
1498 * PB_DLAGEN locally initializes an array A.
1499 *
1500 * Arguments
1501 * =========
1502 *
1503 * UPLO (global input) CHARACTER*1
1504 * On entry, UPLO specifies whether the lower (UPLO='L') trape-
1505 * zoidal part or the upper (UPLO='U') trapezoidal part is to be
1506 * generated when the matrix to be generated is symmetric or
1507 * Hermitian. For all the other values of AFORM, the value of
1508 * this input argument is ignored.
1509 *
1510 * AFORM (global input) CHARACTER*1
1511 * On entry, AFORM specifies the type of submatrix to be genera-
1512 * ted as follows:
1513 * AFORM = 'S', sub( A ) is a symmetric matrix,
1514 * AFORM = 'H', sub( A ) is a Hermitian matrix,
1515 * AFORM = 'T', sub( A ) is overrwritten with the transpose
1516 * of what would normally be generated,
1517 * AFORM = 'C', sub( A ) is overwritten with the conjugate
1518 * transpose of what would normally be genera-
1519 * ted.
1520 * AFORM = 'N', a random submatrix is generated.
1521 *
1522 * A (local output) DOUBLE PRECISION array
1523 * On entry, A is an array of dimension (LLD_A, *). On exit,
1524 * this array contains the local entries of the randomly genera-
1525 * ted submatrix sub( A ).
1526 *
1527 * LDA (local input) INTEGER
1528 * On entry, LDA specifies the local leading dimension of the
1529 * array A. LDA must be at least one.
1530 *
1531 * LCMT00 (global input) INTEGER
1532 * On entry, LCMT00 is the LCM value specifying the off-diagonal
1533 * of the underlying matrix of interest. LCMT00=0 specifies the
1534 * main diagonal, LCMT00 > 0 specifies a subdiagonal, LCMT00 < 0
1535 * specifies superdiagonals.
1536 *
1537 * IRAN (local input) INTEGER array
1538 * On entry, IRAN is an array of dimension 2 containing respec-
1539 * tively the 16-lower and 16-higher bits of the encoding of the
1540 * entry of the random sequence corresponding locally to the
1541 * first local array entry to generate. Usually, this array is
1542 * computed by PB_SETLOCRAN.
1543 *
1544 * MBLKS (local input) INTEGER
1545 * On entry, MBLKS specifies the local number of blocks of rows.
1546 * MBLKS is at least zero.
1547 *
1548 * IMBLOC (local input) INTEGER
1549 * On entry, IMBLOC specifies the number of rows (size) of the
1550 * local uppest blocks. IMBLOC is at least zero.
1551 *
1552 * MB (global input) INTEGER
1553 * On entry, MB specifies the blocking factor used to partition
1554 * the rows of the matrix. MB must be at least one.
1555 *
1556 * LMBLOC (local input) INTEGER
1557 * On entry, LMBLOC specifies the number of rows (size) of the
1558 * local lowest blocks. LMBLOC is at least zero.
1559 *
1560 * NBLKS (local input) INTEGER
1561 * On entry, NBLKS specifies the local number of blocks of co-
1562 * lumns. NBLKS is at least zero.
1563 *
1564 * INBLOC (local input) INTEGER
1565 * On entry, INBLOC specifies the number of columns (size) of
1566 * the local leftmost blocks. INBLOC is at least zero.
1567 *
1568 * NB (global input) INTEGER
1569 * On entry, NB specifies the blocking factor used to partition
1570 * the the columns of the matrix. NB must be at least one.
1571 *
1572 * LNBLOC (local input) INTEGER
1573 * On entry, LNBLOC specifies the number of columns (size) of
1574 * the local rightmost blocks. LNBLOC is at least zero.
1575 *
1576 * JMP (local input) INTEGER array
1577 * On entry, JMP is an array of dimension JMP_LEN containing the
1578 * different jump values used by the random matrix generator.
1579 *
1580 * IMULADD (local input) INTEGER array
1581 * On entry, IMULADD is an array of dimension (4, JMP_LEN). The
1582 * jth column of this array contains the encoded initial cons-
1583 * tants a_j and c_j to jump from X( n ) to X( n + JMP( j ) )
1584 * (= a_j * X( n ) + c_j) in the random sequence. IMULADD(1:2,j)
1585 * contains respectively the 16-lower and 16-higher bits of the
1586 * constant a_j, and IMULADD(3:4,j) contains the 16-lower and
1587 * 16-higher bits of the constant c_j.
1588 *
1589 * -- Written on April 1, 1998 by
1590 * Antoine Petitet, University of Tennessee, Knoxville 37996, USA.
1591 *
1592 * =====================================================================
1593 *
1594 * .. Parameters ..
1595  INTEGER JMP_1, JMP_COL, JMP_IMBV, JMP_INBV, JMP_LEN,
1596  $ JMP_MB, JMP_NB, JMP_NPIMBLOC, JMP_NPMB,
1597  $ JMP_NQINBLOC, JMP_NQNB, JMP_ROW
1598  parameter( jmp_1 = 1, jmp_row = 2, jmp_col = 3,
1599  $ jmp_mb = 4, jmp_imbv = 5, jmp_npmb = 6,
1600  $ jmp_npimbloc = 7, jmp_nb = 8, jmp_inbv = 9,
1601  $ jmp_nqnb = 10, jmp_nqinbloc = 11,
1602  $ jmp_len = 11 )
1603 * ..
1604 * .. Local Scalars ..
1605  INTEGER I, IB, IBLK, II, IK, ITMP, JB, JBLK, JJ, JK,
1606  $ JTMP, LCMTC, LCMTR, LOW, MNB, UPP
1607  DOUBLE PRECISION DUMMY
1608 * ..
1609 * .. Local Arrays ..
1610  INTEGER IB0( 2 ), IB1( 2 ), IB2( 2 ), IB3( 2 )
1611 * ..
1612 * .. External Subroutines ..
1613  EXTERNAL pb_jumpit
1614 * ..
1615 * .. External Functions ..
1616  LOGICAL LSAME
1617  DOUBLE PRECISION PB_DRAND
1618  EXTERNAL LSAME, PB_DRAND
1619 * ..
1620 * .. Intrinsic Functions ..
1621  INTRINSIC max, min
1622 * ..
1623 * .. Executable Statements ..
1624 *
1625  DO 10 i = 1, 2
1626  ib1( i ) = iran( i )
1627  ib2( i ) = iran( i )
1628  ib3( i ) = iran( i )
1629  10 CONTINUE
1630 *
1631  IF( lsame( aform, 'N' ) ) THEN
1632 *
1633 * Generate random matrix
1634 *
1635  jj = 1
1636 *
1637  DO 50 jblk = 1, nblks
1638 *
1639  IF( jblk.EQ.1 ) THEN
1640  jb = inbloc
1641  ELSE IF( jblk.EQ.nblks ) THEN
1642  jb = lnbloc
1643  ELSE
1644  jb = nb
1645  END IF
1646 *
1647  DO 40 jk = jj, jj + jb - 1
1648 *
1649  ii = 1
1650 *
1651  DO 30 iblk = 1, mblks
1652 *
1653  IF( iblk.EQ.1 ) THEN
1654  ib = imbloc
1655  ELSE IF( iblk.EQ.mblks ) THEN
1656  ib = lmbloc
1657  ELSE
1658  ib = mb
1659  END IF
1660 *
1661 * Blocks are IB by JB
1662 *
1663  DO 20 ik = ii, ii + ib - 1
1664  a( ik, jk ) = pb_drand( 0 )
1665  20 CONTINUE
1666 *
1667  ii = ii + ib
1668 *
1669  IF( iblk.EQ.1 ) THEN
1670 *
1671 * Jump IMBLOC + ( NPROW - 1 ) * MB rows
1672 *
1673  CALL pb_jumpit( imuladd( 1, jmp_npimbloc ), ib1,
1674  $ ib0 )
1675 *
1676  ELSE
1677 *
1678 * Jump NPROW * MB rows
1679 *
1680  CALL pb_jumpit( imuladd( 1, jmp_npmb ), ib1, ib0 )
1681 *
1682  END IF
1683 *
1684  ib1( 1 ) = ib0( 1 )
1685  ib1( 2 ) = ib0( 2 )
1686 *
1687  30 CONTINUE
1688 *
1689 * Jump one column
1690 *
1691  CALL pb_jumpit( imuladd( 1, jmp_col ), ib2, ib0 )
1692 *
1693  ib1( 1 ) = ib0( 1 )
1694  ib1( 2 ) = ib0( 2 )
1695  ib2( 1 ) = ib0( 1 )
1696  ib2( 2 ) = ib0( 2 )
1697 *
1698  40 CONTINUE
1699 *
1700  jj = jj + jb
1701 *
1702  IF( jblk.EQ.1 ) THEN
1703 *
1704 * Jump INBLOC + ( NPCOL - 1 ) * NB columns
1705 *
1706  CALL pb_jumpit( imuladd( 1, jmp_nqinbloc ), ib3, ib0 )
1707 *
1708  ELSE
1709 *
1710 * Jump NPCOL * NB columns
1711 *
1712  CALL pb_jumpit( imuladd( 1, jmp_nqnb ), ib3, ib0 )
1713 *
1714  END IF
1715 *
1716  ib1( 1 ) = ib0( 1 )
1717  ib1( 2 ) = ib0( 2 )
1718  ib2( 1 ) = ib0( 1 )
1719  ib2( 2 ) = ib0( 2 )
1720  ib3( 1 ) = ib0( 1 )
1721  ib3( 2 ) = ib0( 2 )
1722 *
1723  50 CONTINUE
1724 *
1725  ELSE IF( lsame( aform, 'T' ) .OR. lsame( aform, 'C' ) ) THEN
1726 *
1727 * Generate the transpose of the matrix that would be normally
1728 * generated.
1729 *
1730  ii = 1
1731 *
1732  DO 90 iblk = 1, mblks
1733 *
1734  IF( iblk.EQ.1 ) THEN
1735  ib = imbloc
1736  ELSE IF( iblk.EQ.mblks ) THEN
1737  ib = lmbloc
1738  ELSE
1739  ib = mb
1740  END IF
1741 *
1742  DO 80 ik = ii, ii + ib - 1
1743 *
1744  jj = 1
1745 *
1746  DO 70 jblk = 1, nblks
1747 *
1748  IF( jblk.EQ.1 ) THEN
1749  jb = inbloc
1750  ELSE IF( jblk.EQ.nblks ) THEN
1751  jb = lnbloc
1752  ELSE
1753  jb = nb
1754  END IF
1755 *
1756 * Blocks are IB by JB
1757 *
1758  DO 60 jk = jj, jj + jb - 1
1759  a( ik, jk ) = pb_drand( 0 )
1760  60 CONTINUE
1761 *
1762  jj = jj + jb
1763 *
1764  IF( jblk.EQ.1 ) THEN
1765 *
1766 * Jump INBLOC + ( NPCOL - 1 ) * NB columns
1767 *
1768  CALL pb_jumpit( imuladd( 1, jmp_nqinbloc ), ib1,
1769  $ ib0 )
1770 *
1771  ELSE
1772 *
1773 * Jump NPCOL * NB columns
1774 *
1775  CALL pb_jumpit( imuladd( 1, jmp_nqnb ), ib1, ib0 )
1776 *
1777  END IF
1778 *
1779  ib1( 1 ) = ib0( 1 )
1780  ib1( 2 ) = ib0( 2 )
1781 *
1782  70 CONTINUE
1783 *
1784 * Jump one row
1785 *
1786  CALL pb_jumpit( imuladd( 1, jmp_row ), ib2, ib0 )
1787 *
1788  ib1( 1 ) = ib0( 1 )
1789  ib1( 2 ) = ib0( 2 )
1790  ib2( 1 ) = ib0( 1 )
1791  ib2( 2 ) = ib0( 2 )
1792 *
1793  80 CONTINUE
1794 *
1795  ii = ii + ib
1796 *
1797  IF( iblk.EQ.1 ) THEN
1798 *
1799 * Jump IMBLOC + ( NPROW - 1 ) * MB rows
1800 *
1801  CALL pb_jumpit( imuladd( 1, jmp_npimbloc ), ib3, ib0 )
1802 *
1803  ELSE
1804 *
1805 * Jump NPROW * MB rows
1806 *
1807  CALL pb_jumpit( imuladd( 1, jmp_npmb ), ib3, ib0 )
1808 *
1809  END IF
1810 *
1811  ib1( 1 ) = ib0( 1 )
1812  ib1( 2 ) = ib0( 2 )
1813  ib2( 1 ) = ib0( 1 )
1814  ib2( 2 ) = ib0( 2 )
1815  ib3( 1 ) = ib0( 1 )
1816  ib3( 2 ) = ib0( 2 )
1817 *
1818  90 CONTINUE
1819 *
1820  ELSE IF( ( lsame( aform, 'S' ) ).OR.( lsame( aform, 'H' ) ) ) THEN
1821 *
1822 * Generate a symmetric matrix
1823 *
1824  IF( lsame( uplo, 'L' ) ) THEN
1825 *
1826 * generate lower trapezoidal part
1827 *
1828  jj = 1
1829  lcmtc = lcmt00
1830 *
1831  DO 170 jblk = 1, nblks
1832 *
1833  IF( jblk.EQ.1 ) THEN
1834  jb = inbloc
1835  low = 1 - inbloc
1836  ELSE IF( jblk.EQ.nblks ) THEN
1837  jb = lnbloc
1838  low = 1 - nb
1839  ELSE
1840  jb = nb
1841  low = 1 - nb
1842  END IF
1843 *
1844  DO 160 jk = jj, jj + jb - 1
1845 *
1846  ii = 1
1847  lcmtr = lcmtc
1848 *
1849  DO 150 iblk = 1, mblks
1850 *
1851  IF( iblk.EQ.1 ) THEN
1852  ib = imbloc
1853  upp = imbloc - 1
1854  ELSE IF( iblk.EQ.mblks ) THEN
1855  ib = lmbloc
1856  upp = mb - 1
1857  ELSE
1858  ib = mb
1859  upp = mb - 1
1860  END IF
1861 *
1862 * Blocks are IB by JB
1863 *
1864  IF( lcmtr.GT.upp ) THEN
1865 *
1866  DO 100 ik = ii, ii + ib - 1
1867  dummy = pb_drand( 0 )
1868  100 CONTINUE
1869 *
1870  ELSE IF( lcmtr.GE.low ) THEN
1871 *
1872  jtmp = jk - jj + 1
1873  mnb = max( 0, -lcmtr )
1874 *
1875  IF( jtmp.LE.min( mnb, jb ) ) THEN
1876 *
1877  DO 110 ik = ii, ii + ib - 1
1878  a( ik, jk ) = pb_drand( 0 )
1879  110 CONTINUE
1880 *
1881  ELSE IF( ( jtmp.GE.( mnb + 1 ) ) .AND.
1882  $ ( jtmp.LE.min( ib-lcmtr, jb ) ) ) THEN
1883 *
1884  itmp = ii + jtmp + lcmtr - 1
1885 *
1886  DO 120 ik = ii, itmp - 1
1887  dummy = pb_drand( 0 )
1888  120 CONTINUE
1889 *
1890  DO 130 ik = itmp, ii + ib - 1
1891  a( ik, jk ) = pb_drand( 0 )
1892  130 CONTINUE
1893 *
1894  END IF
1895 *
1896  ELSE
1897 *
1898  DO 140 ik = ii, ii + ib - 1
1899  a( ik, jk ) = pb_drand( 0 )
1900  140 CONTINUE
1901 *
1902  END IF
1903 *
1904  ii = ii + ib
1905 *
1906  IF( iblk.EQ.1 ) THEN
1907 *
1908 * Jump IMBLOC + ( NPROW - 1 ) * MB rows
1909 *
1910  lcmtr = lcmtr - jmp( jmp_npimbloc )
1911  CALL pb_jumpit( imuladd( 1, jmp_npimbloc ), ib1,
1912  $ ib0 )
1913 *
1914  ELSE
1915 *
1916 * Jump NPROW * MB rows
1917 *
1918  lcmtr = lcmtr - jmp( jmp_npmb )
1919  CALL pb_jumpit( imuladd( 1, jmp_npmb ), ib1,
1920  $ ib0 )
1921 *
1922  END IF
1923 *
1924  ib1( 1 ) = ib0( 1 )
1925  ib1( 2 ) = ib0( 2 )
1926 *
1927  150 CONTINUE
1928 *
1929 * Jump one column
1930 *
1931  CALL pb_jumpit( imuladd( 1, jmp_col ), ib2, ib0 )
1932 *
1933  ib1( 1 ) = ib0( 1 )
1934  ib1( 2 ) = ib0( 2 )
1935  ib2( 1 ) = ib0( 1 )
1936  ib2( 2 ) = ib0( 2 )
1937 *
1938  160 CONTINUE
1939 *
1940  jj = jj + jb
1941 *
1942  IF( jblk.EQ.1 ) THEN
1943 *
1944 * Jump INBLOC + ( NPCOL - 1 ) * NB columns
1945 *
1946  lcmtc = lcmtc + jmp( jmp_nqinbloc )
1947  CALL pb_jumpit( imuladd( 1, jmp_nqinbloc ), ib3, ib0 )
1948 *
1949  ELSE
1950 *
1951 * Jump NPCOL * NB columns
1952 *
1953  lcmtc = lcmtc + jmp( jmp_nqnb )
1954  CALL pb_jumpit( imuladd( 1, jmp_nqnb ), ib3, ib0 )
1955 *
1956  END IF
1957 *
1958  ib1( 1 ) = ib0( 1 )
1959  ib1( 2 ) = ib0( 2 )
1960  ib2( 1 ) = ib0( 1 )
1961  ib2( 2 ) = ib0( 2 )
1962  ib3( 1 ) = ib0( 1 )
1963  ib3( 2 ) = ib0( 2 )
1964 *
1965  170 CONTINUE
1966 *
1967  ELSE
1968 *
1969 * generate upper trapezoidal part
1970 *
1971  ii = 1
1972  lcmtr = lcmt00
1973 *
1974  DO 250 iblk = 1, mblks
1975 *
1976  IF( iblk.EQ.1 ) THEN
1977  ib = imbloc
1978  upp = imbloc - 1
1979  ELSE IF( iblk.EQ.mblks ) THEN
1980  ib = lmbloc
1981  upp = mb - 1
1982  ELSE
1983  ib = mb
1984  upp = mb - 1
1985  END IF
1986 *
1987  DO 240 ik = ii, ii + ib - 1
1988 *
1989  jj = 1
1990  lcmtc = lcmtr
1991 *
1992  DO 230 jblk = 1, nblks
1993 *
1994  IF( jblk.EQ.1 ) THEN
1995  jb = inbloc
1996  low = 1 - inbloc
1997  ELSE IF( jblk.EQ.nblks ) THEN
1998  jb = lnbloc
1999  low = 1 - nb
2000  ELSE
2001  jb = nb
2002  low = 1 - nb
2003  END IF
2004 *
2005 * Blocks are IB by JB
2006 *
2007  IF( lcmtc.LT.low ) THEN
2008 *
2009  DO 180 jk = jj, jj + jb - 1
2010  dummy = pb_drand( 0 )
2011  180 CONTINUE
2012 *
2013  ELSE IF( lcmtc.LE.upp ) THEN
2014 *
2015  itmp = ik - ii + 1
2016  mnb = max( 0, lcmtc )
2017 *
2018  IF( itmp.LE.min( mnb, ib ) ) THEN
2019 *
2020  DO 190 jk = jj, jj + jb - 1
2021  a( ik, jk ) = pb_drand( 0 )
2022  190 CONTINUE
2023 *
2024  ELSE IF( ( itmp.GE.( mnb + 1 ) ) .AND.
2025  $ ( itmp.LE.min( jb+lcmtc, ib ) ) ) THEN
2026 *
2027  jtmp = jj + itmp - lcmtc - 1
2028 *
2029  DO 200 jk = jj, jtmp - 1
2030  dummy = pb_drand( 0 )
2031  200 CONTINUE
2032 *
2033  DO 210 jk = jtmp, jj + jb - 1
2034  a( ik, jk ) = pb_drand( 0 )
2035  210 CONTINUE
2036 *
2037  END IF
2038 *
2039  ELSE
2040 *
2041  DO 220 jk = jj, jj + jb - 1
2042  a( ik, jk ) = pb_drand( 0 )
2043  220 CONTINUE
2044 *
2045  END IF
2046 *
2047  jj = jj + jb
2048 *
2049  IF( jblk.EQ.1 ) THEN
2050 *
2051 * Jump INBLOC + ( NPCOL - 1 ) * NB columns
2052 *
2053  lcmtc = lcmtc + jmp( jmp_nqinbloc )
2054  CALL pb_jumpit( imuladd( 1, jmp_nqinbloc ), ib1,
2055  $ ib0 )
2056 *
2057  ELSE
2058 *
2059 * Jump NPCOL * NB columns
2060 *
2061  lcmtc = lcmtc + jmp( jmp_nqnb )
2062  CALL pb_jumpit( imuladd( 1, jmp_nqnb ), ib1,
2063  $ ib0 )
2064 *
2065  END IF
2066 *
2067  ib1( 1 ) = ib0( 1 )
2068  ib1( 2 ) = ib0( 2 )
2069 *
2070  230 CONTINUE
2071 *
2072 * Jump one row
2073 *
2074  CALL pb_jumpit( imuladd( 1, jmp_row ), ib2, ib0 )
2075 *
2076  ib1( 1 ) = ib0( 1 )
2077  ib1( 2 ) = ib0( 2 )
2078  ib2( 1 ) = ib0( 1 )
2079  ib2( 2 ) = ib0( 2 )
2080 *
2081  240 CONTINUE
2082 *
2083  ii = ii + ib
2084 *
2085  IF( iblk.EQ.1 ) THEN
2086 *
2087 * Jump IMBLOC + ( NPROW - 1 ) * MB rows
2088 *
2089  lcmtr = lcmtr - jmp( jmp_npimbloc )
2090  CALL pb_jumpit( imuladd( 1, jmp_npimbloc ), ib3, ib0 )
2091 *
2092  ELSE
2093 *
2094 * Jump NPROW * MB rows
2095 *
2096  lcmtr = lcmtr - jmp( jmp_npmb )
2097  CALL pb_jumpit( imuladd( 1, jmp_npmb ), ib3, ib0 )
2098 *
2099  END IF
2100 *
2101  ib1( 1 ) = ib0( 1 )
2102  ib1( 2 ) = ib0( 2 )
2103  ib2( 1 ) = ib0( 1 )
2104  ib2( 2 ) = ib0( 2 )
2105  ib3( 1 ) = ib0( 1 )
2106  ib3( 2 ) = ib0( 2 )
2107 *
2108  250 CONTINUE
2109 *
2110  END IF
2111 *
2112  END IF
2113 *
2114  RETURN
2115 *
2116 * End of PB_DLAGEN
2117 *
2118  END
2119  DOUBLE PRECISION FUNCTION pb_drand( IDUMM )
2120 *
2121 * -- PBLAS test routine (version 2.0) --
2122 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
2123 * and University of California, Berkeley.
2124 * April 1, 1998
2125 *
2126 * .. Scalar Arguments ..
2127  INTEGER idumm
2128 * ..
2129 *
2130 * Purpose
2131 * =======
2132 *
2133 * PB_DRAND generates the next number in the random sequence. This func-
2134 * tion ensures that this number will be in the interval ( -1.0, 1.0 ).
2135 *
2136 * Arguments
2137 * =========
2138 *
2139 * IDUMM (local input) INTEGER
2140 * This argument is ignored, but necessary to a FORTRAN 77 func-
2141 * tion.
2142 *
2143 * Further Details
2144 * ===============
2145 *
2146 * On entry, the array IRAND stored in the common block RANCOM contains
2147 * the information (2 integers) required to generate the next number in
2148 * the sequence X( n ). This number is computed as
2149 *
2150 * X( n ) = ( 2^16 * IRAND( 2 ) + IRAND( 1 ) ) / d,
2151 *
2152 * where the constant d is the largest 32 bit positive integer. The
2153 * array IRAND is then updated for the generation of the next number
2154 * X( n+1 ) in the random sequence as follows X( n+1 ) = a * X( n ) + c.
2155 * The constants a and c should have been preliminarily stored in the
2156 * array IACS as 2 pairs of integers. The initial set up of IRAND and
2157 * IACS is performed by the routine PB_SETRAN.
2158 *
2159 * -- Written on April 1, 1998 by
2160 * Antoine Petitet, University of Tennessee, Knoxville 37996, USA.
2161 *
2162 * =====================================================================
2163 *
2164 * .. Parameters ..
2165  DOUBLE PRECISION one, two
2166  PARAMETER ( one = 1.0d+0, two = 2.0d+0 )
2167 * ..
2168 * .. External Functions ..
2169  DOUBLE PRECISION pb_dran
2170  EXTERNAL pb_dran
2171 * ..
2172 * .. Executable Statements ..
2173 *
2174  pb_drand = one - two * pb_dran( idumm )
2175 *
2176  RETURN
2177 *
2178 * End of PB_DRAND
2179 *
2180  END
2181  DOUBLE PRECISION FUNCTION pb_dran( IDUMM )
2182 *
2183 * -- PBLAS test routine (version 2.0) --
2184 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
2185 * and University of California, Berkeley.
2186 * April 1, 1998
2187 *
2188 * .. Scalar Arguments ..
2189  INTEGER idumm
2190 * ..
2191 *
2192 * Purpose
2193 * =======
2194 *
2195 * PB_DRAN generates the next number in the random sequence.
2196 *
2197 * Arguments
2198 * =========
2199 *
2200 * IDUMM (local input) INTEGER
2201 * This argument is ignored, but necessary to a FORTRAN 77 func-
2202 * tion.
2203 *
2204 * Further Details
2205 * ===============
2206 *
2207 * On entry, the array IRAND stored in the common block RANCOM contains
2208 * the information (2 integers) required to generate the next number in
2209 * the sequence X( n ). This number is computed as
2210 *
2211 * X( n ) = ( 2^16 * IRAND( 2 ) + IRAND( 1 ) ) / d,
2212 *
2213 * where the constant d is the largest 32 bit positive integer. The
2214 * array IRAND is then updated for the generation of the next number
2215 * X( n+1 ) in the random sequence as follows X( n+1 ) = a * X( n ) + c.
2216 * The constants a and c should have been preliminarily stored in the
2217 * array IACS as 2 pairs of integers. The initial set up of IRAND and
2218 * IACS is performed by the routine PB_SETRAN.
2219 *
2220 * -- Written on April 1, 1998 by
2221 * Antoine Petitet, University of Tennessee, Knoxville 37996, USA.
2222 *
2223 * =====================================================================
2224 *
2225 * .. Parameters ..
2226  DOUBLE PRECISION divfac, pow16
2227  PARAMETER ( divfac = 2.147483648d+9,
2228  $ pow16 = 6.5536d+4 )
2229 * ..
2230 * .. Local Arrays ..
2231  INTEGER j( 2 )
2232 * ..
2233 * .. External Subroutines ..
2234  EXTERNAL pb_ladd, pb_lmul
2235 * ..
2236 * .. Intrinsic Functions ..
2237  INTRINSIC dble
2238 * ..
2239 * .. Common Blocks ..
2240  INTEGER iacs( 4 ), irand( 2 )
2241  common /rancom/ irand, iacs
2242 * ..
2243 * .. Save Statements ..
2244  SAVE /rancom/
2245 * ..
2246 * .. Executable Statements ..
2247 *
2248  pb_dran = ( dble( irand( 1 ) ) + pow16 * dble( irand( 2 ) ) ) /
2249  $ divfac
2250 *
2251  CALL pb_lmul( irand, iacs, j )
2252  CALL pb_ladd( j, iacs( 3 ), irand )
2253 *
2254  RETURN
2255 *
2256 * End of PB_DRAN
2257 *
2258  END
pdlascal
subroutine pdlascal(TYPE, M, N, ALPHA, A, IA, JA, DESCA)
Definition: pdblastst.f:7337
pb_ladd
subroutine pb_ladd(J, K, I)
Definition: pblastst.f:4480
max
#define max(A, B)
Definition: pcgemr.c:180
pb_setlocran
subroutine pb_setlocran(SEED, ILOCBLK, JLOCBLK, ILOCOFF, JLOCOFF, MYRDIST, MYCDIST, NPROW, NPCOL, JMP, IMULADD, IRAN)
Definition: pblastst.f:4302
pb_setran
subroutine pb_setran(IRAN, IAC)
Definition: pblastst.f:4759
pb_lmul
subroutine pb_lmul(K, J, I)
Definition: pblastst.f:4559
pdladom
subroutine pdladom(INPLACE, N, ALPHA, A, IA, JA, DESCA)
Definition: pdblastst.f:8242
pb_desctrans
subroutine pb_desctrans(DESCIN, DESCOUT)
Definition: pblastst.f:2964
pb_dlagen
subroutine pb_dlagen(UPLO, AFORM, A, LDA, LCMT00, IRAN, MBLKS, IMBLOC, MB, LMBLOC, NBLKS, INBLOC, NB, LNBLOC, JMP, IMULADD)
Definition: pdblastst.f:9737
pb_dran
double precision function pb_dran(IDUMM)
Definition: pdblastst.f:10440
pb_jumpit
subroutine pb_jumpit(MULADD, IRANN, IRANM)
Definition: pblastst.f:4822
pb_infog2l
subroutine pb_infog2l(I, J, DESC, NPROW, NPCOL, MYROW, MYCOL, II, JJ, PROW, PCOL)
Definition: pblastst.f:1673
pb_drand
double precision function pb_drand(IDUMM)
Definition: pdblastst.f:10378
pb_initmuladd
subroutine pb_initmuladd(MULADD0, JMP, IMULADD)
Definition: pblastst.f:4196
pb_ainfog2l
subroutine pb_ainfog2l(M, N, I, J, DESC, NPROW, NPCOL, MYROW, MYCOL, IMB1, INB1, MP, NQ, II, JJ, PROW, PCOL, RPROW, RPCOL)
Definition: pblastst.f:2023
pb_initjmp
subroutine pb_initjmp(COLMAJ, NVIR, IMBVIR, INBVIR, IMBLOC, INBLOC, MB, NB, RSRC, CSRC, NPROW, NPCOL, STRIDE, JMP)
Definition: pblastst.f:4045
pdlagen
subroutine pdlagen(INPLACE, AFORM, DIAG, OFFA, M, N, IA, JA, DESCA, IASEED, A, LDA)
Definition: pdblastst.f:7845
pb_jump
subroutine pb_jump(K, MULADD, IRANN, IRANM, IMA)
Definition: pblastst.f:4648
pxerbla
subroutine pxerbla(ICTXT, SRNAME, INFO)
Definition: pxerbla.f:2
pb_binfo
subroutine pb_binfo(OFFD, M, N, IMB1, INB1, MB, NB, MRROW, MRCOL, LCMT00, MBLKS, NBLKS, IMBLOC, INBLOC, LMBLOC, LNBLOC, ILOW, LOW, IUPP, UPP)
Definition: pblastst.f:3577
pb_chkmat
subroutine pb_chkmat(ICTXT, M, MPOS0, N, NPOS0, IA, JA, DESCA, DPOS0, INFO)
Definition: pblastst.f:2742
pb_dlascal
subroutine pb_dlascal(UPLO, M, N, IOFFD, ALPHA, A, LDA)
Definition: pdblastst.f:9556
pb_locinfo
subroutine pb_locinfo(I, INB, NB, MYROC, SRCPROC, NPROCS, ILOCBLK, ILOCOFF, MYDIST)
Definition: pblastst.f:3910
min
#define min(A, B)
Definition: pcgemr.c:181