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