ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pdlagge.f
Go to the documentation of this file.
1  SUBROUTINE pdlagge( M, N, D, A, IA, JA, DESCA, ISEED, ORDER, WORK,
2  $ LWORK, INFO )
3 *
4 * -- ScaLAPACK routine (version 1.7) --
5 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
6 * and University of California, Berkeley.
7 * May 1, 1997
8 *
9 * .. Scalar Arguments ..
10  INTEGER IA, INFO, JA, LWORK, M, N, ORDER
11 * ..
12 * .. Array Arguments ..
13  INTEGER DESCA( * ), ISEED( 4 )
14  DOUBLE PRECISION A( * ), D( * ), WORK( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * PDLAGGE generates a real symmetric matrix A, by pre- and post-
21 * multiplying a real diagonal matrix D with a random orthogonal
22 * matrices: A = U*D*VT.
23 *
24 * This is just a quick implementation which will be replaced in the
25 * future. The random matrix A1(m,n) is generated and random left
26 * orthogonal matrix U(m,m) is obtained by running QR on A1:
27 * A1(m,n) = U(m,m)*R,
28 * where U(m,m) is a product of min(m,n) Householder rotations.
29 * Afterwards the space of A1 is reused for a second random matrix
30 * A2(m,n), which is used to obtain the right orthogonal matrix VT(n,n)
31 * by running LQ on A2:
32 * A2(m,n) = L*VT(n,n).
33 * This requires vastly more computation than necessary, but not
34 * significantly more communication than is used in the rest of this
35 * routine, and hence is not that much slower than an efficient
36 * solution.
37 *
38 * Notes
39 * =====
40 *
41 * Each global data object is described by an associated description
42 * vector. This vector stores the information required to establish
43 * the mapping between an object element and its corresponding process
44 * and memory location.
45 *
46 * Let A be a generic term for any 2D block cyclicly distributed array.
47 * Such a global array has an associated description vector DESCA.
48 * In the following comments, the character _ should be read as
49 * "of the global array".
50 *
51 * NOTATION STORED IN EXPLANATION
52 * --------------- -------------- --------------------------------------
53 * DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
54 * DTYPE_A = 1.
55 * CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
56 * the BLACS process grid A is distribu-
57 * ted over. The context itself is glo-
58 * bal, but the handle (the integer
59 * value) may vary.
60 * M_A (global) DESCA( M_ ) The number of rows in the global
61 * array A.
62 * N_A (global) DESCA( N_ ) The number of columns in the global
63 * array A.
64 * MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
65 * the rows of the array.
66 * NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
67 * the columns of the array.
68 * RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
69 * row of the array A is distributed.
70 * CSRC_A (global) DESCA( CSRC_ ) The process column over which the
71 * first column of the array A is
72 * distributed.
73 * LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
74 * array. LLD_A >= MAX(1,LOCr(M_A)).
75 *
76 * Let K be the number of rows or columns of a distributed matrix,
77 * and assume that its process grid has dimension p x q.
78 * LOCr( K ) denotes the number of elements of K that a process
79 * would receive if K were distributed over the p processes of its
80 * process column.
81 * Similarly, LOCc( K ) denotes the number of elements of K that a
82 * process would receive if K were distributed over the q processes of
83 * its process row.
84 * The values of LOCr() and LOCc() may be determined via a call to the
85 * ScaLAPACK tool function, NUMROC:
86 * LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
87 * LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
88 * An upper bound for these quantities may be computed by:
89 * LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
90 * LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
91 *
92 * Arguments
93 * =========
94 *
95 * M (global input) INTEGER
96 * Number of rows of the matrix A. M >= 0.
97 *
98 * N (global input) INTEGER
99 * Number of columns of matrix A. N >= 0.
100 *
101 * D (local input) DOUBLE PRECISION array, dimension (N)
102 * The diagonal elements of the diagonal matrix D.
103 *
104 * A (local output) DOUBLE PRECISION array
105 * Global dimension (M, N), local dimension (MP, NQ)
106 *
107 * IA (global input) INTEGER
108 * The global row index of the submatrix of the distributed
109 * matrix A to operate on.
110 *
111 * JA (global input) INTEGER
112 * The global column index of the submatrix of the distributed
113 * matrix A to operate 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. On exit, the seed is updated and will remain identical
122 * on all processes in the context.
123 *
124 * ORDER (global input) INTEGER
125 * Number of reflectors in the matrix Q
126 * At present, ORDER .NE. N is not supported
127 *
128 * WORK (local workspace) DOUBLE PRECISION array, dimension (LWORK)
129 *
130 * LWORK (local input) INTEGER dimension of WORK
131 * LWORK >= MAX( QR_WORK, LQ_WORK )
132 * QR_WORK = LDAA*MAX( 1, NQ ) + 200 + MAX( 1, DTAU1 ) +
133 * MAX( SIZEMQRLEFT, SIZEQRF)
134 * LQ_WORK = LDAA*MAX( 1, NQ ) + 200 + MAX( 1, DTAU2) +
135 * MAX( SIZEMLQRIGHT, SIZEQRF )
136 * Where:
137 * LDAA = DESCA( LLD_ )
138 * MB_A = DESCA( MB_ )
139 * NB_A = DESCA( NB_ )
140 * RSRC_A = DESCA( RSRC_ )
141 * CSRC_A = DESCA( CSRC_ )
142 * LCM = ILCM( NPROW, NPCOL )
143 * LCMQ = LCM / NPCOL
144 * IROFFA = MOD( IA-1, MB_A )
145 * ICOFFA = MOD( JA-1, NB_A )
146 * IAROW = INDXG2P( IA, MB_A, MYROW, RSRC_A, NPROW )
147 * IACOL = INDXG2P( JA, NB_A, MYCOL, CSRC_A, NPCOL )
148 * MP = NUMROC( N+IROFFA, MB_A, MYROW, IAROW, NPROW )
149 * NQ = NUMROC( N+ICOFFA, NB_A, MYCOL, IACOL, NPCOL )
150 * DTAU1 = NUMROC( JA + SIZE- 1, NB_A, MYCOL, IACOL, NPROW )
151 * DTAU2 = NUMROC( IA + SIZE- 1, MB_A, MYROW, IAROW, NPROW )
152 * SIZEMQRLEFT = MAX( (MB_A*(MB_A-1))/2, ( MP + NQ ) * MB_A )
153 * + ( MP + NB_A ) * NB_A
154 * SIZEMLQRIGHT = MAX( (MB_A*(MB_A-1))/2, (MP + NQ)*MB_A ) +
155 * MB_A * MB_A
156 * SIZEQRF = NB_A*NP + MB_A*NQ + NB_A*NB_A
157 *
158 * INFO (local output) INTEGER
159 *
160 * = 0: successful exit
161 * < 0: If the i-th argument is an array and the j-entry had
162 * an illegal value, then INFO = -(i*100+j), if the i-th
163 * argument is a scalar and had an illegal value, then
164 * INFO = -i.
165 *
166 * ======================================================================
167 *
168 * .. Parameters ..
169  INTEGER BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_,
170  $ mb_, nb_, rsrc_, csrc_, lld_
171  parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
172  $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
173  $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
174  DOUBLE PRECISION ZERO
175  parameter( zero = 0.0d+0 )
176 * ..
177 * .. Local Scalars ..
178  INTEGER CSRC_A, DTAU1, DTAU2, I, IACOL, IAROW, ICOFFA,
179  $ iroffa, lcm, lcmq, ldaa, lq_work, lwmin, mb_a,
180  $ mp, mycol, myrow, nb_a, npcol, nprow, nq,
181  $ ptr2aa, ptr2tau, ptr2work, qr_work, rsrc_a,
182  $ SIZE, sizelqf, sizemlqright, sizemqrleft,
183  $ sizeqrf
184 * ..
185 * .. External Subroutines ..
186  EXTERNAL blacs_gridinfo, chk1mat, pdelset, pdgelqf,
188  $ pxerbla
189 * ..
190 * .. External Functions ..
191  INTEGER ILCM, INDXG2P, NUMROC
192  EXTERNAL ilcm, indxg2p, numroc
193 * ..
194 * .. Intrinsic Functions ..
195  INTRINSIC max, min, mod
196 * ..
197 * .. Executable Statements ..
198 * This is just to keep ftnchek happy
199  IF( block_cyclic_2d*dlen_*dtype_*m_*n_.LT.0 )RETURN
200 *
201 * Initialize grid information.
202 *
203  CALL blacs_gridinfo( desca( ctxt_ ), nprow, npcol, myrow, mycol )
204 *
205 * Check LWORK.
206 *
207  info = 0
208  SIZE = min( m, n )
209  IF( nprow.EQ.-1 ) THEN
210  info = -607
211  ELSE
212  CALL chk1mat( m, 1, n, 2, ia, ja, desca, 8, info )
213  END IF
214 * Calculation of a minimum workspace.
215  ldaa = desca( lld_ )
216  mb_a = desca( mb_ )
217  nb_a = desca( nb_ )
218  rsrc_a = desca( rsrc_ )
219  csrc_a = desca( csrc_ )
220  lcm = ilcm( nprow, npcol )
221  lcmq = lcm / npcol
222  iroffa = mod( ia-1, mb_a )
223  icoffa = mod( ja-1, nb_a )
224  iarow = indxg2p( ia, mb_a, myrow, rsrc_a, nprow )
225  iacol = indxg2p( ja, nb_a, mycol, csrc_a, npcol )
226  dtau1 = numroc( ja+size-1, nb_a, mycol, iacol, npcol )
227  dtau2 = numroc( ia+size-1, mb_a, myrow, iarow, nprow )
228  mp = numroc( m+iroffa, mb_a, myrow, iarow, nprow )
229  nq = numroc( n+icoffa, nb_a, mycol, iacol, npcol )
230 *
231  sizemqrleft = max( ( mb_a*( mb_a-1 ) ) / 2, ( mp+nq )*mb_a ) +
232  $ ( mp+nb_a )*nb_a
233  sizemlqright = max( ( mb_a*( mb_a-1 ) ) / 2, ( mp+nq )*mb_a ) +
234  $ mb_a*mb_a
235  sizeqrf = nb_a*mp + mb_a*nq + nb_a*nb_a + 100
236  sizelqf = nb_a*( mp+nq+nb_a ) + 100
237 *
238  qr_work = ldaa*max( 1, nq ) + 200 + max( 1, dtau1 ) +
239  $ max( sizemqrleft, sizeqrf )
240  lq_work = ldaa*max( 1, nq ) + 200 + max( 1, dtau2 ) +
241  $ max( sizemlqright, sizelqf )
242  lwmin = max( qr_work, lq_work )
243  work( 1 ) = lwmin
244  IF( lwork.EQ.-1 )
245  $ GO TO 20
246 *
247 * Test the input arguments.
248 *
249  IF( info.EQ.0 ) THEN
250  IF( size.NE.order ) THEN
251  info = -9
252  ELSE IF( lwork.LT.lwmin ) THEN
253  info = -11
254  END IF
255  END IF
256  IF( info.LT.0 ) THEN
257  CALL pxerbla( desca( ctxt_ ), 'PDLAGGE', -info )
258  RETURN
259  END IF
260 *
261 * Build a diagonal matrix A with the eigenvalues specified in D.
262 *
263  CALL pdlaset( 'Full', m, n, zero, zero, a, ia, ja, desca )
264  DO 10 i = 1, SIZE
265  CALL pdelset( a, i, i, desca, d( i ) )
266  10 CONTINUE
267 *
268 * Local dimension of array TAU in tis case is LOCc(JA+MIN(M,N)-1).
269 *
270  ptr2aa = 2
271  ptr2tau = ptr2aa + ldaa*max( 1, nq ) + 100
272  ptr2work = ptr2tau + max( 1, dtau1 ) + 100
273 *
274  CALL pdlaset( 'All', m, n, zero, zero, work( ptr2aa ), ia, ja,
275  $ desca )
276 *
277 * Build a random matrix AA1.
278 *
279  CALL pdmatgen( desca( ctxt_ ), 'N', 'N', m, n, desca( mb_ ),
280  $ desca( nb_ ), work( ptr2aa ), desca( lld_ ),
281  $ desca( rsrc_ ), desca( csrc_ ), iseed( 1 ), 0, mp,
282  $ 0, nq, myrow, mycol, nprow, npcol )
283 *
284 * Produce QR decomposition AA1 -> U*R.
285 *
286  CALL pdgeqrf( m, n, work( ptr2aa ), ia, ja, desca,
287  $ work( ptr2tau ), work( ptr2work ), sizeqrf, info )
288 *
289 * A = U*A.
290 *
291  CALL pdormqr( 'L', 'N', m, n, SIZE, work( ptr2aa ), ia, ja, desca,
292  $ work( ptr2tau ), a, ia, ja, desca, work( ptr2work ),
293  $ sizemqrleft, info )
294 *
295 * Reinitialize pointer to WORK array. Dimension of array TAU in
296 * this case is LOCr(IA+MIN(M,N)-1).
297 *
298  ptr2work = ptr2tau + max( 1, dtau2 ) + 100
299 *
300 * Use the same workspace to generate a random matrix AA2.
301 *
302  CALL pdmatgen( desca( ctxt_ ), 'N', 'N', m, n, desca( mb_ ),
303  $ desca( nb_ ), work( ptr2aa ), desca( lld_ ),
304  $ desca( rsrc_ ), desca( csrc_ ), iseed( 2 ), 0, mp,
305  $ 0, nq, myrow, mycol, nprow, npcol )
306 *
307 * Produce LQ decomposition of random matrix AA2 -> L*VT.
308 *
309  CALL pdgelqf( m, n, work( ptr2aa ), ia, ja, desca,
310  $ work( ptr2tau ), work( ptr2work ), sizelqf, info )
311 *
312 * Calculate A = A*VT.
313 *
314  CALL pdormlq( 'R', 'N', m, n, SIZE, work( ptr2aa ), ia, ja, desca,
315  $ work( ptr2tau ), a, ia, ja, desca, work( ptr2work ),
316  $ sizemlqright, info )
317 *
318 * End of PDLAGGE
319 *
320  20 CONTINUE
321  RETURN
322  END
pdlagge
subroutine pdlagge(M, N, D, A, IA, JA, DESCA, ISEED, ORDER, WORK, LWORK, INFO)
Definition: pdlagge.f:3
max
#define max(A, B)
Definition: pcgemr.c:180
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
pdlaset
subroutine pdlaset(UPLO, M, N, ALPHA, BETA, A, IA, JA, DESCA)
Definition: pdblastst.f:6862
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
pdgelqf
subroutine pdgelqf(M, N, A, IA, JA, DESCA, TAU, WORK, LWORK, INFO)
Definition: pdgelqf.f:3
pdormlq
subroutine pdormlq(SIDE, TRANS, M, N, K, A, IA, JA, DESCA, TAU, C, IC, JC, DESCC, WORK, LWORK, INFO)
Definition: pdormlq.f:3
pdelset
subroutine pdelset(A, IA, JA, DESCA, ALPHA)
Definition: pdelset.f:2
min
#define min(A, B)
Definition: pcgemr.c:181