ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pdlagsy.f
Go to the documentation of this file.
1 *
2 *
3  SUBROUTINE pdlagsy( N, K, D, A, IA, JA, DESCA, ISEED, ORDER, WORK,
4  $ LWORK, INFO )
5 *
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  INTEGER IA, INFO, JA, K, LWORK, N, ORDER
14 * ..
15 * .. Array Arguments ..
16  INTEGER DESCA( * ), ISEED( 4 )
17  DOUBLE PRECISION A( * ), D( * ), WORK( * )
18 * ..
19 *
20 * Notes
21 * =====
22 *
23 * Each global data object is described by an associated description
24 * vector. This vector stores the information required to establish
25 * the mapping between an object element and its corresponding process
26 * and memory location.
27 *
28 * Let A be a generic term for any 2D block cyclicly distributed array.
29 * Such a global array has an associated description vector DESCA.
30 * In the following comments, the character _ should be read as
31 * "of the global array".
32 *
33 * NOTATION STORED IN EXPLANATION
34 * --------------- -------------- --------------------------------------
35 * DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
36 * DTYPE_A = 1.
37 * CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
38 * the BLACS process grid A is distribu-
39 * ted over. The context itself is glo-
40 * bal, but the handle (the integer
41 * value) may vary.
42 * M_A (global) DESCA( M_ ) The number of rows in the global
43 * array A.
44 * N_A (global) DESCA( N_ ) The number of columns in the global
45 * array A.
46 * MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
47 * the rows of the array.
48 * NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
49 * the columns of the array.
50 * RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
51 * row of the array A is distributed.
52 * CSRC_A (global) DESCA( CSRC_ ) The process column over which the
53 * first column of the array A is
54 * distributed.
55 * LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
56 * array. LLD_A >= MAX(1,LOCr(M_A)).
57 *
58 * Let K be the number of rows or columns of a distributed matrix,
59 * and assume that its process grid has dimension p x q.
60 * LOCr( K ) denotes the number of elements of K that a process
61 * would receive if K were distributed over the p processes of its
62 * process column.
63 * Similarly, LOCc( K ) denotes the number of elements of K that a
64 * process would receive if K were distributed over the q processes of
65 * its process row.
66 * The values of LOCr() and LOCc() may be determined via a call to the
67 * ScaLAPACK tool function, NUMROC:
68 * LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
69 * LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
70 * An upper bound for these quantities may be computed by:
71 * LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
72 * LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
73 *
74 * Purpose
75 * =======
76 *
77 * PDLAGSY generates a real symmetric matrix A, by pre- and post-
78 * multiplying a real diagonal matrix D with a random orthogonal matrix:
79 * A = U*D*U'.
80 *
81 * This is just a quick implementation which will be replaced in the
82 * future. The random orthogonal matrix is computed by creating a
83 * random matrix and running QR on it. This requires vastly more
84 * computation than necessary, but not significantly more communication
85 * than is used in the rest of this rouinte, and hence is not that much
86 * slower than an efficient solution.
87 *
88 * Arguments
89 * =========
90 *
91 * N (global input) INTEGER
92 * The size of the matrix A. N >= 0.
93 *
94 * K (global input) INTEGER
95 * The number of nonzero subdiagonals within the band of A.
96 * 0 <= K <= N-1.
97 * ### K must be 0 or N-1, 0 < K < N-1 is not supported yet.
98 *
99 * D (global input) DOUBLE PRECISION array, dimension (N)
100 * The diagonal elements of the diagonal matrix D.
101 *
102 * A (local output) DOUBLE PRECISION array
103 * Global dimension (N, N), local dimension (NP, NQ)
104 * The generated n by n symmetric matrix A (the full matrix is
105 * stored).
106 *
107 * IA (global input) INTEGER
108 * A's global row index, which points to the beginning of the
109 * submatrix which is to be operated on.
110 *
111 * JA (global input) INTEGER
112 * A's global column index, which points to the beginning of
113 * the submatrix which is to be operated on.
114 *
115 * DESCA (global and local input) INTEGER array of dimension DLEN_.
116 * The array descriptor for the distributed matrix A.
117 *
118 * ISEED (global input/output) INTEGER array, dimension (4)
119 * On entry, the seed of the random number generator; the array
120 * elements must be between 0 and 4095, and ISEED(4) must be
121 * odd.
122 * On exit, the seed is updated and will remain identical on
123 * all processes in the context.
124 *
125 * ORDER (global input) INTEGER
126 * Number of reflectors in the matrix Q
127 * At present, ORDER .NE. N is not supported
128 *
129 * WORK (local workspace) DOUBLE PRECISION array, dimension (LWORK)
130 *
131 * LWORK (local input) INTEGER dimension of WORK
132 * LWORK >= SIZETMS as returned by PDLASIZESEP
133 *
134 *
135 * INFO (local output) INTEGER
136 * = 0: successful exit
137 * < 0: If the i-th argument is an array and the j-entry had
138 * an illegal value, then INFO = -(i*100+j), if the i-th
139 * argument is a scalar and had an illegal value, then
140 * INFO = -i.
141 *
142 *
143 * .. Parameters ..
144  INTEGER BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_,
145  $ mb_, nb_, rsrc_, csrc_, lld_
146  parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
147  $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
148  $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
149  DOUBLE PRECISION ZERO
150  parameter( zero = 0.0d+0 )
151 * ..
152 * .. Local Scalars ..
153  INTEGER CSRC_A, I, IACOL, IAROW, ICOFFA, II, IIROW,
154  $ indaa, indtau, indwork, ipostpad, iprepad,
155  $ iroffa, isizesubtst, isizesyevx, isizetst,
156  $ jjcol, ldaa, lii, liii, ljj, ljjj, lwmin, maxi,
157  $ mb_a, mycol, myrow, nb_a, np, npcol, nprow, nq,
158  $ rsrc_a, sizechk, sizemqrleft, sizemqrright,
159  $ sizeqrf, sizeqtq, sizesubtst, sizesyevx,
160  $ sizetms, sizetst
161 * ..
162 * .. External Subroutines ..
163  EXTERNAL blacs_gridinfo, chk1mat, dlaset, pdgeqrf,
165 * ..
166 * .. External Functions ..
167  INTEGER INDXG2P, NUMROC
168  EXTERNAL indxg2p, numroc
169 * ..
170 * .. Intrinsic Functions ..
171 *
172  INTRINSIC max, min, mod
173 * ..
174 * .. Executable Statements ..
175 * This is just to keep ftnchek happy
176  IF( block_cyclic_2d*csrc_*ctxt_*dlen_*dtype_*lld_*mb_*m_*nb_*n_*
177  $ rsrc_.LT.0 )RETURN
178 *
179 * Initialize grid information
180 *
181  CALL blacs_gridinfo( desca( ctxt_ ), nprow, npcol, myrow, mycol )
182 *
183 * Check LWORK
184 *
185  info = 0
186  IF( nprow.EQ.-1 ) THEN
187  info = -( 700+ctxt_ )
188  ELSE
189  CALL chk1mat( n, 1, n, 1, ia, ja, desca, 7, info )
190  END IF
191 *
192  ldaa = desca( lld_ )
193  mb_a = desca( mb_ )
194  nb_a = desca( nb_ )
195  rsrc_a = desca( rsrc_ )
196  csrc_a = desca( csrc_ )
197  iarow = indxg2p( ia, mb_a, myrow, rsrc_a, nprow )
198  iacol = indxg2p( ja, nb_a, mycol, csrc_a, npcol )
199  iroffa = mod( ia-1, mb_a )
200  icoffa = mod( ja-1, nb_a )
201  np = numroc( n+iroffa, mb_a, myrow, iarow, nprow )
202  nq = numroc( n+icoffa, nb_a, mycol, iacol, npcol )
203  iprepad = 0
204  ipostpad = 0
205  CALL pdlasizesep( desca, iprepad, ipostpad, sizemqrleft,
206  $ sizemqrright, sizeqrf, sizetms, sizeqtq,
207  $ sizechk, sizesyevx, isizesyevx, sizesubtst,
208  $ isizesubtst, sizetst, isizetst )
209  lwmin = sizetms
210 *
211 * Test the input arguments
212 *
213  IF( info.EQ.0 ) THEN
214  IF( k.LT.0 .OR. k.GT.n-1 ) THEN
215  info = -2
216  ELSE IF( n.NE.order ) THEN
217  info = -9
218  ELSE IF( lwork.LT.lwmin ) THEN
219  info = -11
220  END IF
221  END IF
222  IF( info.LT.0 ) THEN
223  CALL pxerbla( desca( ctxt_ ), 'PDLAGSY', -info )
224  RETURN
225  END IF
226 *
227  indaa = 1
228  indtau = indaa + ldaa*max( 1, nq )
229  indwork = indtau + max( 1, nq )
230 *
231  IF( k.NE.0 ) THEN
232  CALL dlaset( 'A', ldaa, nq, zero, zero, work( indaa ), ldaa )
233 *
234 *
235 * Build a random matrix
236 *
237 *
238  CALL pdmatgen( desca( ctxt_ ), 'N', 'N', n, order,
239  $ desca( mb_ ), desca( nb_ ), work( indaa ),
240  $ desca( lld_ ), desca( rsrc_ ), desca( csrc_ ),
241  $ iseed( 1 ), 0, np, 0, nq, myrow, mycol, nprow,
242  $ npcol )
243  CALL pdgeqrf( n, order, work( indaa ), ia, ja, desca,
244  $ work( indtau ), work( indwork ), sizeqrf, info )
245 *
246  END IF
247 *
248 * Build a diagonal matrix A with the eigenvalues specified in D
249 *
250  CALL dlaset( 'A', np, nq, zero, zero, a, desca( lld_ ) )
251 *
252  iirow = 0
253  jjcol = 0
254  lii = 1
255  ljj = 1
256 *
257  DO 20 ii = 1, n, desca( mb_ )
258  maxi = min( n, ii+desca( mb_ )-1 )
259  IF( ( myrow.EQ.iirow ) .AND. ( mycol.EQ.jjcol ) ) THEN
260  liii = lii
261  ljjj = ljj
262  DO 10 i = ii, maxi
263  a( liii+( ljjj-1 )*desca( lld_ ) ) = d( i )
264  liii = liii + 1
265  ljjj = ljjj + 1
266  10 CONTINUE
267  END IF
268  IF( myrow.EQ.iirow )
269  $ lii = lii + desca( mb_ )
270  IF( mycol.EQ.jjcol )
271  $ ljj = ljj + desca( mb_ )
272  iirow = mod( iirow+1, nprow )
273  jjcol = mod( jjcol+1, npcol )
274  20 CONTINUE
275 *
276 * A = Q * A
277 *
278  IF( k.NE.0 ) THEN
279 *
280  CALL pdormqr( 'L', 'Transpose', n, n, order, work( indaa ), ia,
281  $ ja, desca, work( indtau ), a, ia, ja, desca,
282  $ work( indwork ), sizemqrleft, info )
283 *
284 *
285 * A = A * Q'
286 *
287 *
288  CALL pdormqr( 'R', 'N', n, n, order, work( indaa ), ia, ja,
289  $ desca, work( indtau ), a, ia, ja, desca,
290  $ work( indwork ), sizemqrright, info )
291 *
292  END IF
293 *
294 * End of PDLAGSY
295 *
296  END
max
#define max(A, B)
Definition: pcgemr.c:180
pdlasizesep
subroutine pdlasizesep(DESCA, IPREPAD, IPOSTPAD, SIZEMQRLEFT, SIZEMQRRIGHT, SIZEQRF, SIZETMS, SIZEQTQ, SIZECHK, SIZESYEVX, ISIZESYEVX, SIZESUBTST, ISIZESUBTST, SIZETST, ISIZETST)
Definition: pdlasizesep.f:8
pdlagsy
subroutine pdlagsy(N, K, D, A, IA, JA, DESCA, ISEED, ORDER, WORK, LWORK, INFO)
Definition: pdlagsy.f:5
pdormqr
subroutine pdormqr(SIDE, TRANS, M, N, K, A, IA, JA, DESCA, TAU, C, IC, JC, DESCC, WORK, LWORK, INFO)
Definition: pdormqr.f:3
pdmatgen
subroutine pdmatgen(ICTXT, AFORM, DIAG, M, N, MB, NB, A, LDA, IAROW, IACOL, ISEED, IROFF, IRNUM, ICOFF, ICNUM, MYROW, MYCOL, NPROW, NPCOL)
Definition: pdmatgen.f:4
pdgeqrf
subroutine pdgeqrf(M, N, A, IA, JA, DESCA, TAU, WORK, LWORK, INFO)
Definition: pdgeqrf.f:3
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