OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
pdsytdrv.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine pdsytdrv (uplo, n, a, ia, ja, desca, d, e, tau, work, info)

Function/Subroutine Documentation

◆ pdsytdrv()

subroutine pdsytdrv ( character uplo,
integer n,
double precision, dimension( * ) a,
integer ia,
integer ja,
integer, dimension( * ) desca,
double precision, dimension( * ) d,
double precision, dimension( * ) e,
double precision, dimension( * ) tau,
double precision, dimension( * ) work,
integer info )

Definition at line 1 of file pdsytdrv.f.

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 CHARACTER UPLO
11 INTEGER IA, INFO, JA, N
12* ..
13* .. Array Arguments ..
14 INTEGER DESCA( * )
15 DOUBLE PRECISION A( * ), D( * ), E( * ), TAU( * ), WORK( * )
16* ..
17*
18* Purpose
19* =======
20*
21* PDSYTDRV computes sub( A ) = A(IA:IA+N-1,JA:JA+N-1) from Q, the
22* symmetric tridiagonal matrix T (or D and E), and TAU, which were
23* computed by PDSYTRD: sub( A ) := Q * T * Q'.
24*
25* Notes
26* =====
27*
28* Each global data object is described by an associated description
29* vector. This vector stores the information required to establish
30* the mapping between an object element and its corresponding process
31* and memory location.
32*
33* Let A be a generic term for any 2D block cyclicly distributed array.
34* Such a global array has an associated description vector DESCA.
35* In the following comments, the character _ should be read as
36* "of the global array".
37*
38* NOTATION STORED IN EXPLANATION
39* --------------- -------------- --------------------------------------
40* DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
41* DTYPE_A = 1.
42* CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
43* the BLACS process grid A is distribu-
44* ted over. The context itself is glo-
45* bal, but the handle (the integer
46* value) may vary.
47* M_A (global) DESCA( M_ ) The number of rows in the global
48* array A.
49* N_A (global) DESCA( N_ ) The number of columns in the global
50* array A.
51* MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
52* the rows of the array.
53* NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
54* the columns of the array.
55* RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
56* row of the array A is distributed.
57* CSRC_A (global) DESCA( CSRC_ ) The process column over which the
58* first column of the array A is
59* distributed.
60* LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
61* array. LLD_A >= MAX(1,LOCr(M_A)).
62*
63* Let K be the number of rows or columns of a distributed matrix,
64* and assume that its process grid has dimension p x q.
65* LOCr( K ) denotes the number of elements of K that a process
66* would receive if K were distributed over the p processes of its
67* process column.
68* Similarly, LOCc( K ) denotes the number of elements of K that a
69* process would receive if K were distributed over the q processes of
70* its process row.
71* The values of LOCr() and LOCc() may be determined via a call to the
72* ScaLAPACK tool function, NUMROC:
73* LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
74* LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
75* An upper bound for these quantities may be computed by:
76* LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
77* LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
78*
79* Arguments
80* =========
81*
82* UPLO (global input) CHARACTER
83* Specifies whether the upper or lower triangular part of the
84* symmetric matrix sub( A ) is stored:
85* = 'U': Upper triangular
86* = 'L': Lower triangular
87*
88* N (global input) INTEGER
89* The number of rows and columns to be operated on, i.e. the
90* order of the distributed submatrix sub( A ). N >= 0.
91*
92* A (local input/local output) DOUBLE PRECISION pointer into the
93* local memory to an array of dimension (LLD_A,LOCc(JA+N-1)).
94* This array contains the local pieces of sub( A ). On entry,
95* if UPLO='U', the diagonal and first superdiagonal of sub( A )
96* have the corresponding elements of the tridiagonal matrix T,
97* and the elements above the first superdiagonal, with the
98* array TAU, represent the orthogonal matrix Q as a product of
99* elementary reflectors, and the strictly lower triangular part
100* of sub( A ) is not referenced. If UPLO='L', the diagonal and
101* first subdiagonal of sub( A ) have the corresponding elements
102* of the tridiagonal matrix T, and the elements below the first
103* subdiagonal, with the array TAU, represent the orthogonal
104* matrix Q as a product of elementary reflectors, and the
105* strictly upper triangular part of sub( A ) is not referenced.
106* On exit, if UPLO = 'U', the upper triangular part of the
107* distributed symmetric matrix sub( A ) is recovered.
108* If UPLO='L', the lower triangular part of the distributed
109* symmetric matrix sub( A ) is recovered.
110*
111* IA (global input) INTEGER
112* The row index in the global array A indicating the first
113* row of sub( A ).
114*
115* JA (global input) INTEGER
116* The column index in the global array A indicating the
117* first column of sub( A ).
118*
119* DESCA (global and local input) INTEGER array of dimension DLEN_.
120* The array descriptor for the distributed matrix A.
121*
122* D (local input) DOUBLE PRECISION array, dimension LOCc(JA+N-1)
123* The diagonal elements of the tridiagonal matrix T:
124* D(i) = A(i,i). D is tied to the distributed matrix A.
125*
126* E (local input) DOUBLE PRECISION array, dimension LOCc(JA+N-1)
127* if UPLO = 'U', LOCc(JA+N-2) otherwise. The off-diagonal
128* elements of the tridiagonal matrix T: E(i) = A(i,i+1) if
129* UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'. E is tied to the
130* distributed matrix A.
131*
132* TAU (local input) DOUBLE PRECISION, array, dimension
133* LOCc(JA+N-1). This array contains the scalar factors TAU of
134* the elementary reflectors. TAU is tied to the distributed
135* matrix A.
136*
137* WORK (local workspace) DOUBLE PRECISION array, dimension (LWORK)
138* LWORK >= 2 * NB *( NB + NP )
139*
140* where NB = MB_A = NB_A,
141* NP = NUMROC( N, NB, MYROW, IAROW, NPROW ),
142* IAROW = INDXG2P( IA, NB, MYROW, RSRC_A, NPROW ).
143*
144* INDXG2P and NUMROC are ScaLAPACK tool functions;
145* MYROW, MYCOL, NPROW and NPCOL can be determined by calling
146* the subroutine BLACS_GRIDINFO.
147*
148* INFO (global output) INTEGER
149* On exit, if INFO <> 0, a discrepancy has been found between
150* the diagonal and off-diagonal elements of A and the copies
151* contained in the arrays D and E.
152*
153* =====================================================================
154*
155* .. Parameters ..
156 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
157 $ LLD_, MB_, M_, NB_, N_, RSRC_
158 parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
159 $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
160 $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
161 DOUBLE PRECISION EIGHT, HALF, ONE, ZERO
162 parameter( eight = 8.0d+0, half = 0.5d+0, one = 1.0d+0,
163 $ zero = 0.0d+0 )
164* ..
165* .. Local Scalars ..
166 LOGICAL UPPER
167 INTEGER I, IACOL, IAROW, ICTXT, II, IPT, IPV, IPX,
168 $ IPY, J, JB, JJ, JL, K, MYCOL, MYROW, NB, NP,
169 $ NPCOL, NPROW
170 DOUBLE PRECISION ADDBND, D1, D2, E1, E2
171* ..
172* .. Local Arrays ..
173 INTEGER DESCD( DLEN_ ), DESCE( DLEN_ ), DESCV( DLEN_ ),
174 $ DESCT( DLEN_ )
175* ..
176* .. External Functions ..
177 LOGICAL LSAME
178 INTEGER INDXG2P, NUMROC
179 DOUBLE PRECISION PDLAMCH
180 EXTERNAL indxg2p, lsame, numroc, pdlamch
181* ..
182* .. External Subroutines ..
183 EXTERNAL blacs_gridinfo, descset, infog2l, igsum2d,
184 $ pdelget, pdgemm, pdlacpy,
185 $ pdlarft, pdlaset, pdsymm,
186 $ pdsyr2k, pdtrmm
187* ..
188* .. Intrinsic Functions ..
189 INTRINSIC abs, max, min, mod
190* ..
191* .. Executable statements ..
192*
193 ictxt = desca( ctxt_ )
194 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
195*
196 info = 0
197 nb = desca( mb_ )
198 upper = lsame( uplo, 'U' )
199 CALL infog2l( ia, ja, desca, nprow, npcol, myrow, mycol, ii, jj,
200 $ iarow, iacol )
201 np = numroc( n, nb, myrow, iarow, nprow )
202*
203 ipt = 1
204 ipv = nb * nb + ipt
205 ipx = nb * np + ipv
206 ipy = nb * np + ipx
207*
208 CALL descset( descd, 1, ja+n-1, 1, desca( nb_ ), myrow,
209 $ desca( csrc_ ), desca( ctxt_ ), 1 )
210*
211 addbnd = eight * pdlamch( ictxt, 'eps' )
212*
213 IF( upper ) THEN
214*
215 CALL descset( desce, 1, ja+n-1, 1, desca( nb_ ), myrow,
216 $ desca( csrc_ ), desca( ctxt_ ), 1 )
217*
218 DO 10 j = 0, n-1
219 d1 = zero
220 e1 = zero
221 d2 = zero
222 e2 = zero
223 CALL pdelget( ' ', ' ', d2, d, 1, ja+j, descd )
224 CALL pdelget( 'Columnwise', ' ', d1, a, ia+j, ja+j, desca )
225 IF( j.LT.(n-1) ) THEN
226 CALL pdelget( ' ', ' ', e2, e, 1, ja+j+1, desce )
227 CALL pdelget( 'Columnwise', ' ', e1, a, ia+j, ja+j+1,
228 $ desca )
229 END IF
230*
231 IF( ( abs( d1 - d2 ).GT.( abs( d2 ) * addbnd ) ) .OR.
232 $ ( abs( e1 - e2 ).GT.( abs( e2 ) * addbnd ) ) )
233 $ info = info + 1
234 10 CONTINUE
235*
236* Compute the upper triangle of sub( A ).
237*
238 CALL descset( descv, n, nb, nb, nb, iarow, iacol, ictxt,
239 $ max( 1, np ) )
240 CALL descset( desct, nb, nb, nb, nb, iarow, iacol, ictxt, nb )
241*
242 DO 20 k = 0, n-1, nb
243 jb = min( nb, n-k )
244 i = ia + k
245 j = ja + k
246*
247* Compute the lower triangular matrix T.
248*
249 CALL pdlarft( 'Backward', 'Columnwise', k+jb-1, jb, a, ia,
250 $ j, desca, tau, work( ipt ), work( ipv ) )
251*
252* Copy Householder vectors into WORK( IPV ).
253*
254 CALL pdlacpy( 'All', k+jb-1, jb, a, ia, j, desca,
255 $ work( ipv ), 1, 1, descv )
256*
257 IF( k.GT.0 ) THEN
258 CALL pdlaset( 'Lower', jb+1, jb, zero, one, work( ipv ),
259 $ k, 1, descv )
260 ELSE
261 CALL pdlaset( 'Lower', jb, jb-1, zero, one, work( ipv ),
262 $ 1, 2, descv )
263 CALL pdlaset( 'Ge', jb, 1, zero, zero, work( ipv ), 1,
264 $ 1, descv )
265 END IF
266*
267* Zero out the strict upper triangular part of A.
268*
269 IF( k.GT.0 ) THEN
270 CALL pdlaset( 'Ge', k-1, jb, zero, zero, a, ia, j,
271 $ desca )
272 CALL pdlaset( 'Upper', jb-1, jb-1, zero, zero, a, i-1,
273 $ j+1, desca )
274 ELSE IF( jb.GT.1 ) THEN
275 CALL pdlaset( 'Upper', jb-2, jb-2, zero, zero, a, ia,
276 $ j+2, desca )
277 END IF
278*
279* (1) X := A * V * T'
280*
281 CALL pdsymm( 'Left', 'Upper', k+jb, jb, one, a, ia, ja,
282 $ desca, work( ipv ), 1, 1, descv, zero,
283 $ work( ipx ), 1, 1, descv )
284 CALL pdtrmm( 'Right', 'Lower', 'Transpose', 'Non-Unit',
285 $ k+jb, jb, one, work( ipt ), 1, 1, desct,
286 $ work( ipx ), 1, 1, descv )
287*
288* (2) X := X - 1/2 * V * (T * V' * X)
289*
290 CALL pdgemm( 'transpose', 'no transpose', JB, JB, K+JB, ONE,
291 $ WORK( IPV ), 1, 1, DESCV, WORK( IPX ), 1, 1,
292 $ DESCV, ZERO, WORK( IPY ), 1, 1, DESCT )
293 CALL PDTRMM( 'left', 'lower', 'no transpose', 'non-unit',
294 $ JB, JB, ONE, WORK( IPT ), 1, 1, DESCT,
295 $ WORK( IPY ), 1, 1, DESCT )
296 CALL PDGEMM( 'no tranpose', 'no transpose', K+JB, JB, JB,
297 $ -HALF, WORK( IPV ), 1, 1, DESCV, WORK( IPY ),
298 $ 1, 1, DESCT, ONE, WORK( IPX ), 1, 1, DESCV )
299*
300* (3) A := A - X * V' - V * X'
301*
302 CALL PDSYR2K( 'upper', 'no transpose', K+JB, JB, -ONE,
303 $ WORK( IPV ), 1, 1, DESCV, WORK( IPX ), 1, 1,
304 $ DESCV, ONE, A, IA, JA, DESCA )
305*
306 DESCV( CSRC_ ) = MOD( DESCV( CSRC_ ) + 1, NPCOL )
307 DESCT( CSRC_ ) = MOD( DESCT( CSRC_ ) + 1, NPCOL )
308*
309 20 CONTINUE
310*
311 ELSE
312*
313 CALL DESCSET( DESCE, 1, JA+N-2, 1, DESCA( NB_ ), MYROW,
314 $ DESCA( CSRC_ ), DESCA( CTXT_ ), 1 )
315*
316 DO 30 J = 0, N-1
317 D1 = ZERO
318 E1 = ZERO
319 D2 = ZERO
320 E2 = ZERO
321 CALL PDELGET( ' ', ' ', D2, D, 1, JA+J, DESCD )
322 CALL PDELGET( 'columnwise', ' ', D1, A, IA+J, JA+J, DESCA )
323.LT. IF( J(N-1) ) THEN
324 CALL PDELGET( ' ', ' ', E2, E, 1, JA+J, DESCE )
325 CALL PDELGET( 'columnwise', ' ', E1, A, IA+J+1, JA+J,
326 $ DESCA )
327 END IF
328*
329.GT..OR. IF( ( ABS( D1 - D2 )( ABS( D2 ) * ADDBND ) )
330.GT. $ ( ABS( E1 - E2 )( ABS( E2 ) * ADDBND ) ) )
331 $ INFO = INFO + 1
332 30 CONTINUE
333*
334* Compute the lower triangle of sub( A ).
335*
336 JL = MAX( ( ( JA+N-2 ) / NB ) * NB + 1, JA )
337 IACOL = INDXG2P( JL, NB, MYCOL, DESCA( CSRC_ ), NPCOL )
338 CALL DESCSET( DESCV, N, NB, NB, NB, IAROW, IACOL, ICTXT,
339 $ MAX( 1, NP ) )
340 CALL DESCSET( DESCT, NB, NB, NB, NB, INDXG2P( IA+JL-JA+1, NB,
341 $ MYROW, DESCA( RSRC_ ), NPROW ), IACOL, ICTXT,
342 $ NB )
343*
344 DO 40 J = JL, JA, -NB
345 K = J - JA + 1
346 I = IA + K - 1
347 JB = MIN( N-K+1, NB )
348*
349* Compute upper triangular matrix T from TAU.
350*
351 CALL PDLARFT( 'forward', 'columnwise', N-K, JB, A, I+1, J,
352 $ DESCA, TAU, WORK( IPT ), WORK( IPV ) )
353*
354* Copy Householder vectors into WORK( IPV ).
355*
356 CALL PDLACPY( 'lower', N-K, JB, A, I+1, J, DESCA,
357 $ WORK( IPV ), K+1, 1, DESCV )
358 CALL PDLASET( 'upper', N-K, JB, ZERO, ONE, WORK( IPV ),
359 $ K+1, 1, DESCV )
360 CALL PDLASET( 'ge', 1, JB, ZERO, ZERO, WORK( IPV ), K, 1,
361 $ DESCV )
362*
363* Zero out the strict lower triangular part of A.
364*
365 CALL PDLASET( 'lower', N-K-1, JB, ZERO, ZERO, A, I+2, J,
366 $ DESCA )
367*
368* (1) X := A * V * T'
369*
370 CALL PDSYMM( 'left', 'lower', N-K+1, JB, ONE, A, I, J,
371 $ DESCA, WORK( IPV ), K, 1, DESCV, ZERO,
372 $ WORK( IPX ), K, 1, DESCV )
373 CALL PDTRMM( 'right', 'upper', 'transpose', 'non-unit',
374 $ N-K+1, JB, ONE, WORK( IPT ), 1, 1, DESCT,
375 $ WORK( IPX ), K, 1, DESCV )
376*
377* (2) X := X - 1/2 * V * (T * V' * X)
378*
379 CALL PDGEMM( 'transpose', 'no transpose', JB, JB, N-K+1,
380 $ ONE, WORK( IPV ), K, 1, DESCV, WORK( IPX ),
381 $ K, 1, DESCV, ZERO, WORK( IPY ), 1, 1, DESCT )
382 CALL PDTRMM( 'left', 'upper', 'no transpose', 'non-unit',
383 $ JB, JB, ONE, WORK( IPT ), 1, 1, DESCT,
384 $ WORK( IPY ), 1, 1, DESCT )
385 CALL PDGEMM( 'no transpose', 'no transpose', N-K+1, JB, JB,
386 $ -HALF, WORK( IPV ), K, 1, DESCV, WORK( IPY ),
387 $ 1, 1, DESCT, ONE, WORK( IPX ), K, 1, DESCV )
388*
389* (3) A := A - X * V' - V * X'
390*
391 CALL PDSYR2K( 'lower', 'no tranpose', N-K+1, JB, -ONE,
392 $ WORK( IPV ), K, 1, DESCV, WORK( IPX ), K, 1,
393 $ DESCV, ONE, A, I, J, DESCA )
394*
395 DESCV( CSRC_ ) = MOD( DESCV( CSRC_ ) + NPCOL - 1, NPCOL )
396 DESCT( RSRC_ ) = MOD( DESCT( RSRC_ ) + NPROW - 1, NPROW )
397 DESCT( CSRC_ ) = MOD( DESCT( CSRC_ ) + NPCOL - 1, NPCOL )
398*
399 40 CONTINUE
400*
401 END IF
402*
403 CALL IGSUM2D( ICTXT, 'all', ' ', 1, 1, INFO, 1, -1, 0 )
404*
405 RETURN
406*
407* End of PDSYTDRV
408*
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:53
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
indxg2p
integer function indxg2p(indxglob, nb, iproc, isrcproc, nprocs)
Definition mpi.f:947
descset
subroutine descset(desc, m, n, mb, nb, irsrc, icsrc, ictxt, lld)
Definition mpi.f:1610
infog2l
subroutine infog2l(grindx, gcindx, desc, nprow, npcol, myrow, mycol, lrindx, lcindx, rsrc, csrc)
Definition mpi.f:937
blacs_gridinfo
subroutine blacs_gridinfo(cntxt, nprow, npcol, myrow, mycol)
Definition mpi.f:754
numroc
integer function numroc(n, nb, iproc, isrcproc, nprocs)
Definition mpi.f:786
pdlaset
subroutine pdlaset(uplo, m, n, alpha, beta, a, ia, ja, desca)
Definition pdblastst.f:6862
pdlamch
double precision function pdlamch(ictxt, cmach)
Definition pdblastst.f:6769
pdelget
subroutine pdelget(scope, top, alpha, a, ia, ja, desca)
Definition pdelget.f:2
pdlacpy
subroutine pdlacpy(uplo, m, n, a, ia, ja, desca, b, ib, jb, descb)
Definition pdlacpy.f:3
pdlarft
subroutine pdlarft(direct, storev, n, k, v, iv, jv, descv, tau, t, work)
Definition pdlarft.f:3