ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
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
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
pdelget
subroutine pdelget(SCOPE, TOP, ALPHA, A, IA, JA, DESCA)
Definition: pdelget.f:2
descset
subroutine descset(DESC, M, N, MB, NB, IRSRC, ICSRC, ICTXT, LLD)
Definition: descset.f:3
pdlatrd
subroutine pdlatrd(UPLO, N, NB, A, IA, JA, DESCA, D, E, TAU, W, IW, JW, DESCW, WORK)
Definition: pdlatrd.f:3
pdlarfg
subroutine pdlarfg(N, ALPHA, IAX, JAX, X, IX, JX, DESCX, INCX, TAU)
Definition: pdlarfg.f:3
pdelset
subroutine pdelset(A, IA, JA, DESCA, ALPHA)
Definition: pdelset.f:2
min
#define min(A, B)
Definition: pcgemr.c:181