SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
pdlatrd.f
Go to the documentation of this file.
1 SUBROUTINE pdlatrd( UPLO, N, NB, A, IA, JA, DESCA, D, E, TAU, W,
2 $ IW, JW, DESCW, WORK )
3*
4* -- ScaLAPACK auxiliary 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 CHARACTER UPLO
11 INTEGER IA, IW, JA, JW, N, NB
12* ..
13* .. Array Arguments ..
14 INTEGER DESCA( * ), DESCW( * )
15 DOUBLE PRECISION A( * ), D( * ), E( * ), TAU( * ), W( * ),
16 $ work( * )
17* ..
18*
19* Purpose
20* =======
21*
22* PDLATRD reduces NB rows and columns of a real symmetric distributed
23* matrix sub( A ) = A(IA:IA+N-1,JA:JA+N-1) to symmetric tridiagonal
24* form by an orthogonal similarity transformation Q' * sub( A ) * Q,
25* and returns the matrices V and W which are needed to apply the
26* transformation to the unreduced part of sub( A ).
27*
28* If UPLO = 'U', PDLATRD reduces the last NB rows and columns of a
29* matrix, of which the upper triangle is supplied;
30* if UPLO = 'L', PDLATRD reduces the first NB rows and columns of a
31* matrix, of which the lower triangle is supplied.
32*
33* This is an auxiliary routine called by PDSYTRD.
34*
35* Notes
36* =====
37*
38* Each global data object is described by an associated description
39* vector. This vector stores the information required to establish
40* the mapping between an object element and its corresponding process
41* and memory location.
42*
43* Let A be a generic term for any 2D block cyclicly distributed array.
44* Such a global array has an associated description vector DESCA.
45* In the following comments, the character _ should be read as
46* "of the global array".
47*
48* NOTATION STORED IN EXPLANATION
49* --------------- -------------- --------------------------------------
50* DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
51* DTYPE_A = 1.
52* CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
53* the BLACS process grid A is distribu-
54* ted over. The context itself is glo-
55* bal, but the handle (the integer
56* value) may vary.
57* M_A (global) DESCA( M_ ) The number of rows in the global
58* array A.
59* N_A (global) DESCA( N_ ) The number of columns in the global
60* array A.
61* MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
62* the rows of the array.
63* NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
64* the columns of the array.
65* RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
66* row of the array A is distributed.
67* CSRC_A (global) DESCA( CSRC_ ) The process column over which the
68* first column of the array A is
69* distributed.
70* LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
71* array. LLD_A >= MAX(1,LOCr(M_A)).
72*
73* Let K be the number of rows or columns of a distributed matrix,
74* and assume that its process grid has dimension p x q.
75* LOCr( K ) denotes the number of elements of K that a process
76* would receive if K were distributed over the p processes of its
77* process column.
78* Similarly, LOCc( K ) denotes the number of elements of K that a
79* process would receive if K were distributed over the q processes of
80* its process row.
81* The values of LOCr() and LOCc() may be determined via a call to the
82* ScaLAPACK tool function, NUMROC:
83* LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
84* LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
85* An upper bound for these quantities may be computed by:
86* LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
87* LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
88*
89* Arguments
90* =========
91*
92* UPLO (global input) CHARACTER
93* Specifies whether the upper or lower triangular part of the
94* symmetric matrix sub( A ) is stored:
95* = 'U': Upper triangular
96* = 'L': Lower triangular
97*
98* N (global input) INTEGER
99* The number of rows and columns to be operated on, i.e. the
100* order of the distributed submatrix sub( A ). N >= 0.
101*
102* NB (global input) INTEGER
103* The number of rows and columns to be reduced.
104*
105* A (local input/local output) DOUBLE PRECISION pointer into the
106* local memory to an array of dimension (LLD_A,LOCc(JA+N-1)).
107* On entry, this array contains the local pieces of the
108* symmetric distributed matrix sub( A ). If UPLO = 'U', the
109* leading N-by-N upper triangular part of sub( A ) contains
110* the upper triangular part of the matrix, and its strictly
111* lower triangular part is not referenced. If UPLO = 'L', the
112* leading N-by-N lower triangular part of sub( A ) contains the
113* lower triangular part of the matrix, and its strictly upper
114* triangular part is not referenced.
115* On exit, if UPLO = 'U', the last NB columns have been reduced
116* to tridiagonal form, with the diagonal elements overwriting
117* the diagonal elements of sub( A ); the elements above the
118* diagonal with the array TAU, represent the orthogonal matrix
119* Q as a product of elementary reflectors. If UPLO = 'L', the
120* first NB columns have been reduced to tridiagonal form, with
121* the diagonal elements overwriting the diagonal elements of
122* sub( A ); the elements below the diagonal with the array TAU,
123* represent the orthogonal matrix Q as a product of elementary
124* reflectors; See Further Details.
125*
126* IA (global input) INTEGER
127* The row index in the global array A indicating the first
128* row of sub( A ).
129*
130* JA (global input) INTEGER
131* The column index in the global array A indicating the
132* first column of sub( A ).
133*
134* DESCA (global and local input) INTEGER array of dimension DLEN_.
135* The array descriptor for the distributed matrix A.
136*
137* D (local output) DOUBLE PRECISION array, dimension LOCc(JA+N-1)
138* The diagonal elements of the tridiagonal matrix T:
139* D(i) = A(i,i). D is tied to the distributed matrix A.
140*
141* E (local output) DOUBLE PRECISION array, dimension LOCc(JA+N-1)
142* if UPLO = 'U', LOCc(JA+N-2) otherwise. The off-diagonal
143* elements of the tridiagonal matrix T: E(i) = A(i,i+1) if
144* UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'. E is tied to the
145* distributed matrix A.
146*
147* TAU (local output) DOUBLE PRECISION array, dimension
148* LOCc(JA+N-1). This array contains the scalar factors TAU of
149* the elementary reflectors. TAU is tied to the distributed
150* matrix A.
151*
152* W (local output) DOUBLE PRECISION pointer into the local memory
153* to an array of dimension (LLD_W,NB_W), This array contains
154* the local pieces of the N-by-NB_W matrix W required to
155* update the unreduced part of sub( A ).
156*
157* IW (global input) INTEGER
158* The row index in the global array W indicating the first
159* row of sub( W ).
160*
161* JW (global input) INTEGER
162* The column index in the global array W indicating the
163* first column of sub( W ).
164*
165* DESCW (global and local input) INTEGER array of dimension DLEN_.
166* The array descriptor for the distributed matrix W.
167*
168* WORK (local workspace) DOUBLE PRECISION array, dimension (NB_A)
169*
170* Further Details
171* ===============
172*
173* If UPLO = 'U', the matrix Q is represented as a product of elementary
174* reflectors
175*
176* Q = H(n) H(n-1) . . . H(n-nb+1).
177*
178* Each H(i) has the form
179*
180* H(i) = I - tau * v * v'
181*
182* where tau is a real scalar, and v is a real vector with
183* v(i:n) = 0 and v(i-1) = 1; v(1:i-1) is stored on exit in
184* A(ia:ia+i-2,ja+i), and tau in TAU(ja+i-1).
185*
186* If UPLO = 'L', the matrix Q is represented as a product of elementary
187* reflectors
188*
189* Q = H(1) H(2) . . . H(nb).
190*
191* Each H(i) has the form
192*
193* H(i) = I - tau * v * v'
194*
195* where tau is a real scalar, and v is a real vector with
196* v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in
197* A(ia+i+1:ia+n-1,ja+i-1), and tau in TAU(ja+i-1).
198*
199* The elements of the vectors v together form the N-by-NB matrix V
200* which is needed, with W, to apply the transformation to the unreduced
201* part of the matrix, using a symmetric rank-2k update of the form:
202* sub( A ) := sub( A ) - V*W' - W*V'.
203*
204* The contents of A on exit are illustrated by the following examples
205* with n = 5 and nb = 2:
206*
207* if UPLO = 'U': if UPLO = 'L':
208*
209* ( a a a v4 v5 ) ( d )
210* ( a a v4 v5 ) ( 1 d )
211* ( a 1 v5 ) ( v1 1 a )
212* ( d 1 ) ( v1 v2 a a )
213* ( d ) ( v1 v2 a a a )
214*
215* where d denotes a diagonal element of the reduced matrix, a denotes
216* an element of the original matrix that is unchanged, and vi denotes
217* an element of the vector defining H(i).
218*
219* =====================================================================
220*
221* .. Parameters ..
222 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
223 $ lld_, mb_, m_, nb_, n_, rsrc_
224 parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
225 $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
226 $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
227 DOUBLE PRECISION HALF, ONE, ZERO
228 parameter( half = 0.5d+0, one = 1.0d+0, zero = 0.0d+0 )
229* ..
230* .. Local Scalars ..
231 INTEGER I, IACOL, IAROW, ICTXT, II, J, JJ, JP, JWK, K,
232 $ kw, mycol, myrow, npcol, nprow, nq
233 DOUBLE PRECISION ALPHA
234* ..
235* .. Local Arrays ..
236 INTEGER DESCD( DLEN_ ), DESCE( DLEN_ ), DESCWK( DLEN_ )
237* ..
238* .. External Subroutines ..
239 EXTERNAL blacs_gridinfo, descset, dgebr2d, dgebs2d,
240 $ infog2l, pdaxpy, pddot, pdelget,
241 $ pdelset, pdgemv, pdlarfg, pdscal,
242 $ pdsymv
243* ..
244* .. External Functions ..
245 LOGICAL LSAME
246 INTEGER NUMROC
247 EXTERNAL lsame, numroc
248* ..
249* .. Intrinsic Functions ..
250 INTRINSIC min
251* ..
252* .. Executable Statements ..
253*
254* Quick return if possible
255*
256 IF( n.LE.0 )
257 $ RETURN
258*
259 ictxt = desca( ctxt_ )
260 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
261 nq = max( 1, numroc( ja+n-1, desca( nb_ ), mycol, desca( csrc_ ),
262 $ npcol ) )
263 CALL descset( descd, 1, ja+n-1, 1, desca( nb_ ), myrow,
264 $ desca( csrc_ ), desca( ctxt_ ), 1 )
265*
266 IF( lsame( uplo, 'U' ) ) THEN
267*
268 CALL infog2l( n+ia-nb, n+ja-nb, desca, nprow, npcol, myrow,
269 $ mycol, ii, jj, iarow, iacol )
270 CALL descset( descwk, 1, descw( nb_ ), 1, descw( nb_ ), iarow,
271 $ iacol, ictxt, 1 )
272 CALL descset( desce, 1, ja+n-1, 1, desca( nb_ ), myrow,
273 $ desca( csrc_ ), desca( ctxt_ ), 1 )
274*
275* Reduce last NB columns of upper triangle
276*
277 DO 10 j = ja+n-1, ja+n-nb, -1
278 i = ia + j - ja
279 k = j - ja + 1
280 kw = mod( k-1, desca( mb_ ) ) + 1
281*
282* Update A(IA:I,I)
283*
284 CALL pdgemv( 'No transpose', k, n-k, -one, a, ia, j+1,
285 $ desca, w, iw+k-1, jw+kw, descw, descw( m_ ),
286 $ one, a, ia, j, desca, 1 )
287 CALL pdgemv( 'No transpose', k, n-k, -one, w, iw, jw+kw,
288 $ descw, a, i, j+1, desca, desca( m_ ), one, a,
289 $ ia, j, desca, 1 )
290 IF( n-k.GT.0 )
291 $ CALL pdelset( a, i, j+1, desca, e( jp ) )
292*
293* Generate elementary reflector H(i) to annihilate
294* A(IA:I-2,I)
295*
296 jp = min( jj+kw-1, nq )
297 CALL pdlarfg( k-1, e( jp ), i-1, j, a, ia, j, desca, 1,
298 $ tau )
299 CALL pdelset( a, i-1, j, desca, one )
300*
301* Compute W(IW:IW+K-2,JW+KW-1)
302*
303 CALL pdsymv( 'Upper', k-1, one, a, ia, ja, desca, a, ia, j,
304 $ desca, 1, zero, w, iw, jw+kw-1, descw, 1 )
305*
306 jwk = mod( k-1, descwk( nb_ ) ) + 2
307 CALL pdgemv( 'Transpose', k-1, n-k, one, w, iw, jw+kw,
308 $ descw, a, ia, j, desca, 1, zero, work, 1, jwk,
309 $ descwk, descwk( m_ ) )
310 CALL pdgemv( 'No transpose', k-1, n-k, -one, a, ia, j+1,
311 $ desca, work, 1, jwk, descwk, descwk( m_ ), one,
312 $ w, iw, jw+kw-1, descw, 1 )
313 CALL pdgemv( 'Transpose', k-1, n-k, one, a, ia, j+1, desca,
314 $ a, ia, j, desca, 1, zero, work, 1, jwk, descwk,
315 $ descwk( m_ ) )
316 CALL pdgemv( 'No transpose', k-1, n-k, -one, w, iw, jw+kw,
317 $ descw, work, 1, jwk, descwk, descwk( m_ ), one,
318 $ w, iw, jw+kw-1, descw, 1 )
319 CALL pdscal( k-1, tau( jp ), w, iw, jw+kw-1, descw, 1 )
320*
321 CALL pddot( k-1, alpha, w, iw, jw+kw-1, descw, 1, a, ia, j,
322 $ desca, 1 )
323 IF( mycol.EQ.iacol )
324 $ alpha = -half*tau( jp )*alpha
325 CALL pdaxpy( k-1, alpha, a, ia, j, desca, 1, w, iw, jw+kw-1,
326 $ descw, 1 )
327 IF( mycol.EQ.iacol ) THEN
328 CALL pdelget( 'E', ' ', d( jp ), a, i, j, desca )
329 END IF
330*
331 10 CONTINUE
332*
333 ELSE
334*
335 CALL infog2l( ia, ja, desca, nprow, npcol, myrow, mycol, ii,
336 $ jj, iarow, iacol )
337 CALL descset( descwk, 1, descw( nb_ ), 1, descw( nb_ ), iarow,
338 $ iacol, ictxt, 1 )
339 CALL descset( desce, 1, ja+n-2, 1, desca( nb_ ), myrow,
340 $ desca( csrc_ ), desca( ctxt_ ), 1 )
341*
342* Reduce first NB columns of lower triangle
343*
344 DO 20 j = ja, ja+nb-1
345 i = ia + j - ja
346 k = j - ja + 1
347*
348* Update A(J:JA+N-1,J)
349*
350 CALL pdgemv( 'No transpose', n-k+1, k-1, -one, a, i, ja,
351 $ desca, w, iw+k-1, jw, descw, descw( m_ ), one,
352 $ a, i, j, desca, 1 )
353 CALL pdgemv( 'No transpose', n-k+1, k-1, -one, w, iw+k-1,
354 $ jw, descw, a, i, ja, desca, desca( m_ ), one,
355 $ a, i, j, desca, 1 )
356 IF( k.GT.1 )
357 $ CALL pdelset( a, i, j-1, desca, e( jp ) )
358*
359*
360* Generate elementary reflector H(i) to annihilate
361* A(I+2:IA+N-1,I)
362*
363 jp = min( jj+k-1, nq )
364 CALL pdlarfg( n-k, e( jp ), i+1, j, a, i+2, j, desca, 1,
365 $ tau )
366 CALL pdelset( a, i+1, j, desca, one )
367*
368* Compute W(IW+K:IW+N-1,JW+K-1)
369*
370 CALL pdsymv( 'Lower', n-k, one, a, i+1, j+1, desca, a, i+1,
371 $ j, desca, 1, zero, w, iw+k, jw+k-1, descw, 1 )
372*
373 CALL pdgemv( 'Transpose', n-k, k-1, one, w, iw+k, jw, descw,
374 $ a, i+1, j, desca, 1, zero, work, 1, 1, descwk,
375 $ descwk( m_ ) )
376 CALL pdgemv( 'No transpose', n-k, k-1, -one, a, i+1, ja,
377 $ desca, work, 1, 1, descwk, descwk( m_ ), one, w,
378 $ iw+k, jw+k-1, descw, 1 )
379 CALL pdgemv( 'Transpose', n-k, k-1, one, a, i+1, ja, desca,
380 $ a, i+1, j, desca, 1, zero, work, 1, 1, descwk,
381 $ descwk( m_ ) )
382 CALL pdgemv( 'No transpose', n-k, k-1, -one, w, iw+k, jw,
383 $ descw, work, 1, 1, descwk, descwk( m_ ), one, w,
384 $ iw+k, jw+k-1, descw, 1 )
385 CALL pdscal( n-k, tau( jp ), w, iw+k, jw+k-1, descw, 1 )
386 CALL pddot( n-k, alpha, w, iw+k, jw+k-1, descw, 1, a, i+1,
387 $ j, desca, 1 )
388 IF( mycol.EQ.iacol )
389 $ alpha = -half*tau( jp )*alpha
390 CALL pdaxpy( n-k, alpha, a, i+1, j, desca, 1, w, iw+k,
391 $ jw+k-1, descw, 1 )
392 IF( mycol.EQ.iacol ) THEN
393 CALL pdelget( 'E', ' ', d( jp ), a, i, j, desca )
394 END IF
395*
396 20 CONTINUE
397*
398 END IF
399*
400* Broadcast columnwise the diagonal elements into D.
401*
402 IF( mycol.EQ.iacol ) THEN
403 IF( myrow.EQ.iarow ) THEN
404 CALL dgebs2d( ictxt, 'Columnwise', ' ', 1, nb, d( jj ), 1 )
405 ELSE
406 CALL dgebr2d( ictxt, 'Columnwise', ' ', 1, nb, d( jj ), 1,
407 $ iarow, mycol )
408 END IF
409 END IF
410*
411 RETURN
412*
413* End of PDLATRD
414*
415 END
subroutine descset(desc, m, n, mb, nb, irsrc, icsrc, ictxt, lld)
Definition descset.f:3
subroutine infog2l(grindx, gcindx, desc, nprow, npcol, myrow, mycol, lrindx, lcindx, rsrc, csrc)
Definition infog2l.f:3
#define max(A, B)
Definition pcgemr.c:180
#define min(A, B)
Definition pcgemr.c:181
subroutine pdelget(scope, top, alpha, a, ia, ja, desca)
Definition pdelget.f:2
subroutine pdelset(a, ia, ja, desca, alpha)
Definition pdelset.f:2
subroutine pdlarfg(n, alpha, iax, jax, x, ix, jx, descx, incx, tau)
Definition pdlarfg.f:3
subroutine pdlatrd(uplo, n, nb, a, ia, ja, desca, d, e, tau, w, iw, jw, descw, work)
Definition pdlatrd.f:3