ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pdlatms.f
Go to the documentation of this file.
1 *
2 *
3  SUBROUTINE pdlatms( M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX,
4  $ KL, KU, PACK, A, IA, JA, DESCA, ORDER, WORK,
5  $ LWORK, INFO )
6 *
7 * -- ScaLAPACK routine (version 1.7) --
8 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
9 * and University of California, Berkeley.
10 * May 1, 1997
11 *
12 * .. Scalar Arguments ..
13  CHARACTER DIST, PACK, SYM
14  INTEGER IA, INFO, JA, KL, KU, LWORK, M, MODE, N, ORDER
15  DOUBLE PRECISION COND, DMAX
16 * ..
17 * .. Array Arguments ..
18  INTEGER DESCA( * ), ISEED( 4 )
19  DOUBLE PRECISION A( * ), D( * ), WORK( * )
20 * ..
21 *
22 * Purpose
23 * =======
24 *
25 * PDLATMS generates random symmetric matrices with specified
26 * eigenvalues for testing SCALAPACK programs.
27 *
28 * PDLATMS operates by applying the following sequence of
29 * operations:
30 *
31 * Set the diagonal to D, where D may be input or
32 * computed according to MODE, COND, DMAX, and SYM
33 * as described below.
34 *
35 * Generate a dense M x N matrix by multiplying D on the left
36 * and the right by random unitary matrices, then:
37 *
38 * Reduce the bandwidth according to KL and KU, using
39 * Householder transformations.
40 * ### bandwidth reduction NOT SUPPORTED ###
41 *
42 * Arguments
43 * =========
44 *
45 * M - (global input) INTEGER
46 * The number of rows of A. Not modified.
47 *
48 * N - (global input) INTEGER
49 * The number of columns of A. Not modified.
50 * ### M .ne. N unsupported
51 *
52 * DIST - (global input) CHARACTER*1
53 * On entry, DIST specifies the type of distribution to be used
54 * to generate the random eigen-/singular values.
55 * 'U' => UNIFORM( 0, 1 ) ( 'U' for uniform )
56 * 'S' => UNIFORM( -1, 1 ) ( 'S' for symmetric )
57 * 'N' => NORMAL( 0, 1 ) ( 'N' for normal )
58 * Not modified.
59 *
60 * ISEED - (global input) INTEGER array, dimension ( 4 )
61 * On entry ISEED specifies the seed of the random number
62 * generator. They should lie between 0 and 4095 inclusive,
63 * and ISEED(4) should be odd. The random number generator
64 * uses a linear congruential sequence limited to small
65 * integers, and so should produce machine independent
66 * random numbers. The values of ISEED are changed on
67 * exit, and can be used in the next call to DLATMS
68 * to continue the same random number sequence.
69 * Changed on exit.
70 *
71 * SYM - (global input) CHARACTER*1
72 * If SYM='S' or 'H', the generated matrix is symmetric, with
73 * eigenvalues specified by D, COND, MODE, and DMAX; they
74 * may be positive, negative, or zero.
75 * If SYM='P', the generated matrix is symmetric, with
76 * eigenvalues (= singular values) specified by D, COND,
77 * MODE, and DMAX; they will not be negative.
78 * If SYM='N', the generated matrix is nonsymmetric, with
79 * singular values specified by D, COND, MODE, and DMAX;
80 * they will not be negative.
81 * ### SYM = 'N' NOT SUPPORTED ###
82 * Not modified.
83 *
84 * D - (local input/output) DOUBLE PRECISION array,
85 * dimension ( MIN( M , N ) )
86 * This array is used to specify the singular values or
87 * eigenvalues of A (see SYM, above.) If MODE=0, then D is
88 * assumed to contain the singular/eigenvalues, otherwise
89 * they will be computed according to MODE, COND, and DMAX,
90 * and placed in D.
91 * Modified if MODE is nonzero.
92 *
93 * MODE - (global input) INTEGER
94 * On entry this describes how the singular/eigenvalues are to
95 * be specified:
96 * MODE = 0 means use D as input
97 * MODE = 1 sets D(1)=1 and D(2:N)=1.0/COND
98 * MODE = 2 sets D(1:N-1)=1 and D(N)=1.0/COND
99 * MODE = 3 sets D(I)=COND**(-(I-1)/(N-1))
100 * MODE = 4 sets D(i)=1 - (i-1)/(N-1)*(1 - 1/COND)
101 * MODE = 5 sets D to random numbers in the range
102 * ( 1/COND , 1 ) such that their logarithms
103 * are uniformly distributed.
104 * MODE = 6 set D to random numbers from same distribution
105 * as the rest of the matrix.
106 * MODE < 0 has the same meaning as ABS(MODE), except that
107 * the order of the elements of D is reversed.
108 * Thus if MODE is positive, D has entries ranging from
109 * 1 to 1/COND, if negative, from 1/COND to 1,
110 * If SYM='S' or 'H', and MODE is neither 0, 6, nor -6, then
111 * the elements of D will also be multiplied by a random
112 * sign (i.e., +1 or -1.)
113 * Not modified.
114 *
115 * COND - (global input) DOUBLE PRECISION
116 * On entry, this is used as described under MODE above.
117 * If used, it must be >= 1. Not modified.
118 *
119 * DMAX - (global input) DOUBLE PRECISION
120 * If MODE is neither -6, 0 nor 6, the contents of D, as
121 * computed according to MODE and COND, will be scaled by
122 * DMAX / max(abs(D(i))); thus, the maximum absolute eigen- or
123 * singular value (which is to say the norm) will be abs(DMAX).
124 * Note that DMAX need not be positive: if DMAX is negative
125 * (or zero), D will be scaled by a negative number (or zero).
126 * Not modified.
127 *
128 * KL - (global input) INTEGER
129 * This specifies the lower bandwidth of the matrix. For
130 * example, KL=0 implies upper triangular, KL=1 implies upper
131 * Hessenberg, and KL being at least M-1 means that the matrix
132 * has full lower bandwidth. KL must equal KU if the matrix
133 * is symmetric.
134 * Not modified.
135 * ### 1 <= KL < N-1 is NOT SUPPORTED ###
136 *
137 * KU - (global input) INTEGER
138 * This specifies the upper bandwidth of the matrix. For
139 * example, KU=0 implies lower triangular, KU=1 implies lower
140 * Hessenberg, and KU being at least N-1 means that the matrix
141 * has full upper bandwidth. KL must equal KU if the matrix
142 * is symmetric.
143 * Not modified.
144 * ### 1 <= KU < N-1 is NOT SUPPORTED ###
145 *
146 * PACK - (global input) CHARACTER*1
147 * This specifies packing of matrix as follows:
148 * 'N' => no packing
149 * ### PACK must be 'N' all other options NOT SUPPORTED ###
150 *
151 * A - (local output) DOUBLE PRECISION array
152 * Global dimension (M, N), local dimension (MP, NQ)
153 * On exit A is the desired test matrix.
154 *
155 * IA (global input) INTEGER
156 * A's global row index, which points to the beginning of the
157 * submatrix which is to be operated on.
158 *
159 * JA (global input) INTEGER
160 * A's global column index, which points to the beginning of
161 * the submatrix which is to be operated on.
162 *
163 * DESCA (global and local input) INTEGER array of dimension DLEN_.
164 * The array descriptor for the distributed matrix A.
165 *
166 * ORDER - (input) INTEGER
167 * The number of reflectors used to define the orthogonal
168 * matrix Q. A = Q * D * Q'
169 * Higher ORDER requires more computation and communication.
170 *
171 * WORK - (local input/output) DOUBLE PRECISION array,
172 * dimension (LWORK)
173 *
174 * LWORK - (local input) INTEGER dimension of WORK
175 * LWORK >= SIZETMS as returned by PDLASIZESEP
176 *
177 * INFO - (global output) INTEGER
178 * Error code. On exit, INFO will be set to one of the
179 * following values:
180 * 0 => normal return
181 * -1 => M negative or unequal to N and SYM='S', 'H', or 'P'
182 * -2 => N negative
183 * -3 => DIST illegal string
184 * -5 => SYM illegal string
185 * -7 => MODE not in range -6 to 6
186 * -8 => COND less than 1.0, and MODE neither -6, 0 nor 6
187 * -10 => KL negative
188 * -11 => KU negative, or SYM='S' or 'H' and KU not equal to KL
189 * -16 => DESCA is inconsistent
190 * -17 => ORDER not in the range 0 to N inclusive
191 * 1 => Error return from DLATM1
192 * 2 => Cannot scale to DMAX (max. sing. value is 0)
193 * 3 => Error return from PDLAGSY
194 *
195 *-----------------------------------------------------------------------
196 *
197 *
198 * .. Parameters ..
199  INTEGER BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_,
200  $ MB_, NB_, RSRC_, CSRC_, LLD_
201  parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
202  $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
203  $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
204  DOUBLE PRECISION ZERO, ONE
205  PARAMETER ( ZERO = 0.0d+0, one = 1.0d+0 )
206 * ..
207 * .. Local Scalars ..
208  INTEGER I, IDIST, IINFO, IPACK, IRSIGN, ISYM, LLB,
209  $ MNMIN, MYCOL, MYROW, NP, NPCOL, NPROW, NQ
210  DOUBLE PRECISION ALPHA, TEMP
211 * ..
212 * .. Local Arrays ..
213  INTEGER IDUM1( 1 ), IDUM2( 1 )
214 * ..
215 * .. External Functions ..
216  LOGICAL LSAME
217  INTEGER NUMROC
218  EXTERNAL lsame, numroc
219 * ..
220 * .. External Subroutines ..
221  EXTERNAL blacs_gridinfo, chk1mat, dlaset, dlatm1, dscal,
222  $ pchk1mat, pdlagsy, pxerbla
223 * ..
224 * .. Intrinsic Functions ..
225  INTRINSIC abs, max, min, mod
226 * ..
227 * .. Executable Statements ..
228 * This is just to keep ftnchek happy
229  IF( block_cyclic_2d*csrc_*ctxt_*dlen_*dtype_*lld_*mb_*m_*nb_*n_*
230  $ rsrc_.LT.0 )RETURN
231 *
232 * 1) Decode and Test the input parameters.
233 * Initialize flags & seed.
234 *
235 *
236  info = 0
237 *
238  CALL blacs_gridinfo( desca( ctxt_ ), nprow, npcol, myrow, mycol )
239  IF( ( myrow.GE.nprow .OR. myrow.LT.0 ) .OR.
240  $ ( mycol.GE.npcol .OR. mycol.LT.0 ) )RETURN
241 *
242  np = numroc( n, desca( mb_ ), myrow, 0, nprow )
243  nq = numroc( n, desca( nb_ ), mycol, 0, npcol )
244 *
245 * Quick return if possible
246 *
247  IF( m.EQ.0 .OR. n.EQ.0 )
248  $ RETURN
249 *
250 * Decode DIST
251 *
252  IF( lsame( dist, 'U' ) ) THEN
253  idist = 1
254  ELSE IF( lsame( dist, 'S' ) ) THEN
255  idist = 2
256  ELSE IF( lsame( dist, 'N' ) ) THEN
257  idist = 3
258  ELSE
259  idist = -1
260  END IF
261 *
262 * Decode SYM
263 *
264  IF( lsame( sym, 'N' ) ) THEN
265  isym = 1
266  irsign = 0
267  ELSE IF( lsame( sym, 'P' ) ) THEN
268  isym = 2
269  irsign = 0
270  ELSE IF( lsame( sym, 'S' ) ) THEN
271  isym = 2
272  irsign = 1
273  ELSE IF( lsame( sym, 'H' ) ) THEN
274  isym = 2
275  irsign = 1
276  ELSE
277  isym = -1
278  END IF
279 *
280 * Decode PACK
281 *
282  IF( lsame( pack, 'N' ) ) THEN
283  ipack = 0
284  ELSE
285  ipack = 1
286  END IF
287 *
288 * Set certain internal parameters
289 *
290  mnmin = min( m, n )
291  llb = min( kl, m-1 )
292 *
293  IF( order.EQ.0 )
294  $ order = n
295 *
296 * Set INFO if an error
297 *
298  IF( nprow.EQ.-1 ) THEN
299  info = -( 1600+ctxt_ )
300  ELSE
301  CALL chk1mat( m, 1, n, 2, ia, ja, desca, 16, info )
302  IF( info.EQ.0 ) THEN
303  IF( m.NE.n .AND. isym.NE.1 ) THEN
304  info = -2
305  ELSE IF( idist.EQ.-1 ) THEN
306  info = -3
307  ELSE IF( isym.EQ.-1 ) THEN
308  info = -5
309  ELSE IF( abs( mode ).GT.6 ) THEN
310  info = -7
311  ELSE IF( ( mode.NE.0 .AND. abs( mode ).NE.6 ) .AND. cond.LT.
312  $ one ) THEN
313  info = -8
314  ELSE IF( kl.LT.0 ) THEN
315  info = -10
316  ELSE IF( ku.LT.0 .OR. ( isym.NE.1 .AND. kl.NE.ku ) ) THEN
317  info = -11
318  ELSE IF( ( order.LT.0 ) .OR. ( order.GT.n ) ) THEN
319  info = -17
320  END IF
321  END IF
322  CALL pchk1mat( m, 1, n, 2, ia, ja, desca, 16, 0, idum1, idum2,
323  $ info )
324  END IF
325 *
326 * Check for unsupported features
327 *
328  IF( isym.NE.2 ) THEN
329  info = -5
330  ELSE IF( ipack.NE.0 ) THEN
331  info = -12
332  ELSE IF( kl.GT.0 .AND. kl.LT.m-1 ) THEN
333  info = -10
334  ELSE IF( ku.GT.0 .AND. ku.LT.n-1 ) THEN
335  info = -11
336  ELSE IF( llb.NE.0 .AND. llb.NE.m-1 ) THEN
337  info = -10
338  END IF
339  IF( info.NE.0 ) THEN
340  CALL pxerbla( desca( ctxt_ ), 'PDLATMS', -info )
341  RETURN
342  END IF
343 *
344 * Initialize random number generator
345 *
346  DO 10 i = 1, 4
347  iseed( i ) = mod( abs( iseed( i ) ), 4096 )
348  10 CONTINUE
349 *
350  IF( mod( iseed( 4 ), 2 ).NE.1 )
351  $ iseed( 4 ) = iseed( 4 ) + 1
352 *
353 * 2) Set up D if indicated.
354 *
355 * Compute D according to COND and MODE
356 *
357  CALL dlatm1( mode, cond, irsign, idist, iseed, d, mnmin, iinfo )
358 *
359  IF( iinfo.NE.0 ) THEN
360  info = 1
361  RETURN
362  END IF
363 *
364 *
365  IF( mode.NE.0 .AND. abs( mode ).NE.6 ) THEN
366 *
367 * Scale by DMAX
368 *
369  temp = abs( d( 1 ) )
370  DO 20 i = 2, mnmin
371  temp = max( temp, abs( d( i ) ) )
372  20 CONTINUE
373 *
374  IF( temp.GT.zero ) THEN
375  alpha = dmax / temp
376  ELSE
377  info = 2
378  RETURN
379  END IF
380 *
381  CALL dscal( mnmin, alpha, d, 1 )
382 *
383  END IF
384 *
385  CALL dlaset( 'A', np, nq, zero, zero, a, desca( lld_ ) )
386 *
387 * symmetric -- A = U D U'
388 *
389  CALL pdlagsy( m, llb, d, a, ia, ja, desca, iseed, order, work,
390  $ lwork, iinfo )
391 *
392  RETURN
393 *
394 * End of PDLATMS
395 *
396  END
max
#define max(A, B)
Definition: pcgemr.c:180
pdlagsy
subroutine pdlagsy(N, K, D, A, IA, JA, DESCA, ISEED, ORDER, WORK, LWORK, INFO)
Definition: pdlagsy.f:5
pchk1mat
subroutine pchk1mat(MA, MAPOS0, NA, NAPOS0, IA, JA, DESCA, DESCAPOS0, NEXTRA, EX, EXPOS, INFO)
Definition: pchkxmat.f:3
pdlatms
subroutine pdlatms(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, IA, JA, DESCA, ORDER, WORK, LWORK, INFO)
Definition: pdlatms.f:6
dlatm1
subroutine dlatm1(MODE, COND, IRSIGN, IDIST, ISEED, D, N, INFO)
Definition: dlatm1.f:2
chk1mat
subroutine chk1mat(MA, MAPOS0, NA, NAPOS0, IA, JA, DESCA, DESCAPOS0, INFO)
Definition: chk1mat.f:3
pxerbla
subroutine pxerbla(ICTXT, SRNAME, INFO)
Definition: pxerbla.f:2
min
#define min(A, B)
Definition: pcgemr.c:181