SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
pzlatrd.f
Go to the documentation of this file.
1 SUBROUTINE pzlatrd( 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 D( * ), E( * )
16 COMPLEX*16 A( * ), TAU( * ), W( * ), WORK( * )
17* ..
18*
19* Purpose
20* =======
21*
22* PZLATRD reduces NB rows and columns of a complex Hermitian
23* distributed matrix sub( A ) = A(IA:IA+N-1,JA:JA+N-1) to complex
24* tridiagonal form by an unitary similarity transformation
25* Q' * sub( A ) * Q, and returns the matrices V and W which are
26* needed to apply the transformation to the unreduced part of sub( A ).
27*
28* If UPLO = 'U', PZLATRD reduces the last NB rows and columns of a
29* matrix, of which the upper triangle is supplied;
30* if UPLO = 'L', PZLATRD 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 PZHETRD.
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* Hermitian 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) COMPLEX*16 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* Hermitian 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 unitary matrix Q
119* 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 unitary 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) COMPLEX*16, 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) COMPLEX*16 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) COMPLEX*16 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 complex scalar, and v is a complex 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 complex scalar, and v is a complex 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 Hermitian 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 COMPLEX*16 HALF, ONE, ZERO
228 parameter( half = ( 0.5d+0, 0.0d+0 ),
229 $ one = ( 1.0d+0, 0.0d+0 ),
230 $ zero = ( 0.0d+0, 0.0d+0 ) )
231* ..
232* .. Local Scalars ..
233 INTEGER I, IACOL, IAROW, ICTXT, II, J, JJ, JP, JWK, K,
234 $ kw, mycol, myrow, npcol, nprow, nq
235 COMPLEX*16 AII, ALPHA, BETA
236* ..
237* .. Local Arrays ..
238 INTEGER DESCD( DLEN_ ), DESCE( DLEN_ ), DESCWK( DLEN_ )
239* ..
240* .. External Subroutines ..
241 EXTERNAL blacs_gridinfo, descset, dgebr2d, dgebs2d,
242 $ infog2l, pdelset, pzaxpy, pzdotc,
243 $ pzelget, pzelset, pzgemv, pzhemv,
244 $ pzlacgv, pzlarfg, pzscal
245* ..
246* .. External Functions ..
247 LOGICAL LSAME
248 INTEGER NUMROC
249 EXTERNAL lsame, numroc
250* ..
251* .. Intrinsic Functions ..
252 INTRINSIC dble, dcmplx, min
253* ..
254* .. Executable Statements ..
255*
256* Quick return if possible
257*
258 IF( n.LE.0 )
259 $ RETURN
260*
261 ictxt = desca( ctxt_ )
262 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
263 nq = max( 1, numroc( ja+n-1, desca( nb_ ), mycol, desca( csrc_ ),
264 $ npcol ) )
265 CALL descset( descd, 1, ja+n-1, 1, desca( nb_ ), myrow,
266 $ desca( csrc_ ), desca( ctxt_ ), 1 )
267 aii = zero
268 beta = zero
269*
270 IF( lsame( uplo, 'U' ) ) THEN
271*
272 CALL infog2l( n+ia-nb, n+ja-nb, desca, nprow, npcol, myrow,
273 $ mycol, ii, jj, iarow, iacol )
274 CALL descset( descwk, 1, descw( nb_ ), 1, descw( nb_ ), iarow,
275 $ iacol, ictxt, 1 )
276 CALL descset( desce, 1, ja+n-1, 1, desca( nb_ ), myrow,
277 $ desca( csrc_ ), desca( ctxt_ ), 1 )
278*
279* Reduce last NB columns of upper triangle
280*
281 DO 10 j = ja+n-1, ja+n-nb, -1
282 i = ia + j - ja
283 k = j - ja + 1
284 kw = mod( k-1, desca( mb_ ) ) + 1
285*
286* Update A(IA:I,I)
287*
288 CALL pzelget( 'E', ' ', aii, a, i, j, desca )
289 CALL pzelset( a, i, j, desca, dcmplx( dble( aii ) ) )
290 CALL pzlacgv( n-k, w, iw+k-1, jw+kw, descw, descw( m_ ) )
291 CALL pzgemv( 'No transpose', k, n-k, -one, a, ia, j+1,
292 $ desca, w, iw+k-1, jw+kw, descw, descw( m_ ),
293 $ one, a, ia, j, desca, 1 )
294 CALL pzlacgv( n-k, w, iw+k-1, jw+kw, descw, descw( m_ ) )
295 CALL pzlacgv( n-k, a, i, j+1, desca, desca( m_ ) )
296 CALL pzgemv( 'No transpose', k, n-k, -one, w, iw, jw+kw,
297 $ descw, a, i, j+1, desca, desca( m_ ), one, a,
298 $ ia, j, desca, 1 )
299 CALL pzlacgv( n-k, a, i, j+1, desca, desca( m_ ) )
300 CALL pzelget( 'E', ' ', aii, a, i, j, desca )
301 CALL pzelset( a, i, j, desca, dcmplx( dble( aii ) ) )
302 IF( n-k.GT.0 )
303 $ CALL pzelset( a, i, j+1, desca, dcmplx( e( jp ) ) )
304*
305* Generate elementary reflector H(i) to annihilate
306* A(IA:I-2,I)
307*
308 jp = min( jj+kw-1, nq )
309 CALL pzlarfg( k-1, beta, i-1, j, a, ia, j, desca, 1,
310 $ tau )
311 CALL pdelset( e, 1, j, desce, dble( beta ) )
312 CALL pzelset( a, i-1, j, desca, one )
313*
314* Compute W(IW:IW+K-2,JW+KW-1)
315*
316 CALL pzhemv( 'Upper', k-1, one, a, ia, ja, desca, a, ia, j,
317 $ desca, 1, zero, w, iw, jw+kw-1, descw, 1 )
318*
319 jwk = mod( k-1, descwk( nb_ ) ) + 2
320 CALL pzgemv( 'Conjugate transpose', k-1, n-k, one, w, iw,
321 $ jw+kw, descw, a, ia, j, desca, 1, zero, work,
322 $ 1, jwk, descwk, descwk( m_ ) )
323 CALL pzgemv( 'No transpose', k-1, n-k, -one, a, ia, j+1,
324 $ desca, work, 1, jwk, descwk, descwk( m_ ), one,
325 $ w, iw, jw+kw-1, descw, 1 )
326 CALL pzgemv( 'Conjugate transpose', k-1, n-k, one, a, ia,
327 $ j+1, desca, a, ia, j, desca, 1, zero, work, 1,
328 $ jwk, descwk, descwk( m_ ) )
329 CALL pzgemv( 'No transpose', k-1, n-k, -one, w, iw, jw+kw,
330 $ descw, work, 1, jwk, descwk, descwk( m_ ), one,
331 $ w, iw, jw+kw-1, descw, 1 )
332 CALL pzscal( k-1, tau( jp ), w, iw, jw+kw-1, descw, 1 )
333*
334 CALL pzdotc( k-1, alpha, w, iw, jw+kw-1, descw, 1, a, ia, j,
335 $ desca, 1 )
336 IF( mycol.EQ.iacol )
337 $ alpha = -half*tau( jp )*alpha
338 CALL pzaxpy( k-1, alpha, a, ia, j, desca, 1, w, iw, jw+kw-1,
339 $ descw, 1 )
340 CALL pzelget( 'E', ' ', beta, a, i, j, desca )
341 CALL pdelset( d, 1, j, descd, dble( beta ) )
342*
343 10 CONTINUE
344*
345 ELSE
346*
347 CALL infog2l( ia, ja, desca, nprow, npcol, myrow, mycol, ii,
348 $ jj, iarow, iacol )
349 CALL descset( descwk, 1, descw( nb_ ), 1, descw( nb_ ), iarow,
350 $ iacol, ictxt, 1 )
351 CALL descset( desce, 1, ja+n-2, 1, desca( nb_ ), myrow,
352 $ desca( csrc_ ), desca( ctxt_ ), 1 )
353*
354* Reduce first NB columns of lower triangle
355*
356 DO 20 j = ja, ja+nb-1
357 i = ia + j - ja
358 k = j - ja + 1
359*
360* Update A(J:JA+N-1,J)
361*
362 CALL pzelget( 'E', ' ', aii, a, i, j, desca )
363 CALL pzelset( a, i, j, desca, dcmplx( dble( aii ) ) )
364 CALL pzlacgv( k-1, w, iw+k-1, jw, descw, descw( m_ ) )
365 CALL pzgemv( 'No transpose', n-k+1, k-1, -one, a, i, ja,
366 $ desca, w, iw+k-1, jw, descw, descw( m_ ), one,
367 $ a, i, j, desca, 1 )
368 CALL pzlacgv( k-1, w, iw+k-1, jw, descw, descw( m_ ) )
369 CALL pzlacgv( k-1, a, i, ja, desca, desca( m_ ) )
370 CALL pzgemv( 'No transpose', n-k+1, k-1, -one, w, iw+k-1,
371 $ jw, descw, a, i, ja, desca, desca( m_ ), one,
372 $ a, i, j, desca, 1 )
373 CALL pzlacgv( k-1, a, i, ja, desca, desca( m_ ) )
374 CALL pzelget( 'E', ' ', aii, a, i, j, desca )
375 CALL pzelset( a, i, j, desca, dcmplx( dble( aii ) ) )
376 IF( k.GT.1 )
377 $ CALL pzelset( a, i, j-1, desca, dcmplx( e( jp ) ) )
378*
379*
380* Generate elementary reflector H(i) to annihilate
381* A(I+2:IA+N-1,I)
382*
383 jp = min( jj+k-1, nq )
384 CALL pzlarfg( n-k, beta, i+1, j, a, i+2, j, desca, 1,
385 $ tau )
386 CALL pdelset( e, 1, j, desce, dble( beta ) )
387 CALL pzelset( a, i+1, j, desca, one )
388*
389* Compute W(IW+K:IW+N-1,JW+K-1)
390*
391 CALL pzhemv( 'Lower', n-k, one, a, i+1, j+1, desca, a, i+1,
392 $ j, desca, 1, zero, w, iw+k, jw+k-1, descw, 1 )
393*
394 CALL pzgemv( 'Conjugate Transpose', n-k, k-1, one, w, iw+k,
395 $ jw, descw, a, i+1, j, desca, 1, zero, work, 1,
396 $ 1, descwk, descwk( m_ ) )
397 CALL pzgemv( 'No transpose', n-k, k-1, -one, a, i+1, ja,
398 $ desca, work, 1, 1, descwk, descwk( m_ ), one, w,
399 $ iw+k, jw+k-1, descw, 1 )
400 CALL pzgemv( 'Conjugate transpose', n-k, k-1, one, a, i+1,
401 $ ja, desca, a, i+1, j, desca, 1, zero, work, 1,
402 $ 1, descwk, descwk( m_ ) )
403 CALL pzgemv( 'No transpose', n-k, k-1, -one, w, iw+k, jw,
404 $ descw, work, 1, 1, descwk, descwk( m_ ), one, w,
405 $ iw+k, jw+k-1, descw, 1 )
406 CALL pzscal( n-k, tau( jp ), w, iw+k, jw+k-1, descw, 1 )
407 CALL pzdotc( n-k, alpha, w, iw+k, jw+k-1, descw, 1, a, i+1,
408 $ j, desca, 1 )
409 IF( mycol.EQ.iacol )
410 $ alpha = -half*tau( jp )*alpha
411 CALL pzaxpy( n-k, alpha, a, i+1, j, desca, 1, w, iw+k,
412 $ jw+k-1, descw, 1 )
413 CALL pzelget( 'E', ' ', beta, a, i, j, desca )
414 CALL pdelset( d, 1, j, descd, dble( beta ) )
415*
416 20 CONTINUE
417*
418 END IF
419*
420* Broadcast columnwise the diagonal elements into D.
421*
422 IF( mycol.EQ.iacol ) THEN
423 IF( myrow.EQ.iarow ) THEN
424 CALL dgebs2d( ictxt, 'Columnwise', ' ', 1, nb, d( jj ), 1 )
425 ELSE
426 CALL dgebr2d( ictxt, 'Columnwise', ' ', 1, nb, d( jj ), 1,
427 $ iarow, mycol )
428 END IF
429 END IF
430*
431 RETURN
432*
433* End of PZLATRD
434*
435 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 pdelset(a, ia, ja, desca, alpha)
Definition pdelset.f:2
subroutine pzelget(scope, top, alpha, a, ia, ja, desca)
Definition pzelget.f:2
subroutine pzelset(a, ia, ja, desca, alpha)
Definition pzelset.f:2
subroutine pzlacgv(n, x, ix, jx, descx, incx)
Definition pzlacgv.f:2
subroutine pzlarfg(n, alpha, iax, jax, x, ix, jx, descx, incx, tau)
Definition pzlarfg.f:3
subroutine pzlatrd(uplo, n, nb, a, ia, ja, desca, d, e, tau, w, iw, jw, descw, work)
Definition pzlatrd.f:3