ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
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
pzelget
subroutine pzelget(SCOPE, TOP, ALPHA, A, IA, JA, DESCA)
Definition: pzelget.f:2
max
#define max(A, B)
Definition: pcgemr.c:180
infog2l
subroutine infog2l(GRINDX, GCINDX, DESC, NPROW, NPCOL, MYROW, MYCOL, LRINDX, LCINDX, RSRC, CSRC)
Definition: infog2l.f:3
pzlatrd
subroutine pzlatrd(UPLO, N, NB, A, IA, JA, DESCA, D, E, TAU, W, IW, JW, DESCW, WORK)
Definition: pzlatrd.f:3
pzlarfg
subroutine pzlarfg(N, ALPHA, IAX, JAX, X, IX, JX, DESCX, INCX, TAU)
Definition: pzlarfg.f:3
descset
subroutine descset(DESC, M, N, MB, NB, IRSRC, ICSRC, ICTXT, LLD)
Definition: descset.f:3
pzlacgv
subroutine pzlacgv(N, X, IX, JX, DESCX, INCX)
Definition: pzlacgv.f:2
pzelset
subroutine pzelset(A, IA, JA, DESCA, ALPHA)
Definition: pzelset.f:2
pdelset
subroutine pdelset(A, IA, JA, DESCA, ALPHA)
Definition: pdelset.f:2
min
#define min(A, B)
Definition: pcgemr.c:181