SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ pdstedc()

subroutine pdstedc ( character  compz,
integer  n,
double precision, dimension( * )  d,
double precision, dimension( * )  e,
double precision, dimension( * )  q,
integer  iq,
integer  jq,
integer, dimension( * )  descq,
double precision, dimension( * )  work,
integer  lwork,
integer, dimension( * )  iwork,
integer  liwork,
integer  info 
)

Definition at line 1 of file pdstedc.f.

3*
4* -- ScaLAPACK routine (version 1.7) --
5* University of Tennessee, Knoxville, Oak Ridge National Laboratory,
6* and University of California, Berkeley.
7* March 13, 2000
8*
9* .. Scalar Arguments ..
10 CHARACTER COMPZ
11 INTEGER INFO, IQ, JQ, LIWORK, LWORK, N
12* ..
13* .. Array Arguments ..
14 INTEGER DESCQ( * ), IWORK( * )
15 DOUBLE PRECISION D( * ), E( * ), Q( * ), WORK( * )
16* ..
17*
18* Purpose
19* =======
20* PDSTEDC computes all eigenvalues and eigenvectors of a
21* symmetric tridiagonal matrix in parallel, using the divide and
22* conquer algorithm.
23*
24* This code makes very mild assumptions about floating point
25* arithmetic. It will work on machines with a guard digit in
26* add/subtract, or on those binary machines without guard digits
27* which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
28* It could conceivably fail on hexadecimal or decimal machines
29* without guard digits, but we know of none. See DLAED3 for details.
30*
31* Arguments
32* =========
33*
34* COMPZ (input) CHARACTER*1
35* = 'N': Compute eigenvalues only. (NOT IMPLEMENTED YET)
36* = 'I': Compute eigenvectors of tridiagonal matrix also.
37* = 'V': Compute eigenvectors of original dense symmetric
38* matrix also. On entry, Z contains the orthogonal
39* matrix used to reduce the original matrix to
40* tridiagonal form. (NOT IMPLEMENTED YET)
41*
42* N (global input) INTEGER
43* The order of the tridiagonal matrix T. N >= 0.
44*
45* D (global input/output) DOUBLE PRECISION array, dimension (N)
46* On entry, the diagonal elements of the tridiagonal matrix.
47* On exit, if INFO = 0, the eigenvalues in descending order.
48*
49* E (global input/output) DOUBLE PRECISION array, dimension (N-1)
50* On entry, the subdiagonal elements of the tridiagonal matrix.
51* On exit, E has been destroyed.
52*
53* Q (local output) DOUBLE PRECISION array,
54* local dimension ( LLD_Q, LOCc(JQ+N-1))
55* Q contains the orthonormal eigenvectors of the symmetric
56* tridiagonal matrix.
57* On output, Q is distributed across the P processes in block
58* cyclic format.
59*
60* IQ (global input) INTEGER
61* Q's global row index, which points to the beginning of the
62* submatrix which is to be operated on.
63*
64* JQ (global input) INTEGER
65* Q's global column index, which points to the beginning of
66* the submatrix which is to be operated on.
67*
68* DESCQ (global and local input) INTEGER array of dimension DLEN_.
69* The array descriptor for the distributed matrix Z.
70*
71*
72* WORK (local workspace/output) DOUBLE PRECISION array,
73* dimension (LWORK)
74* On output, WORK(1) returns the workspace needed.
75*
76* LWORK (local input/output) INTEGER,
77* the dimension of the array WORK.
78* LWORK = 6*N + 2*NP*NQ
79* NP = NUMROC( N, NB, MYROW, DESCQ( RSRC_ ), NPROW )
80* NQ = NUMROC( N, NB, MYCOL, DESCQ( CSRC_ ), NPCOL )
81*
82* If LWORK = -1, the LWORK is global input and a workspace
83* query is assumed; the routine only calculates the minimum
84* size for the WORK array. The required workspace is returned
85* as the first element of WORK and no error message is issued
86* by PXERBLA.
87*
88* IWORK (local workspace/output) INTEGER array, dimension (LIWORK)
89* On exit, if LIWORK > 0, IWORK(1) returns the optimal LIWORK.
90*
91* LIWORK (input) INTEGER
92* The dimension of the array IWORK.
93* LIWORK = 2 + 7*N + 8*NPCOL
94*
95* INFO (global output) INTEGER
96* = 0: successful exit
97* < 0: If the i-th argument is an array and the j-entry had
98* an illegal value, then INFO = -(i*100+j), if the i-th
99* argument is a scalar and had an illegal value, then
100* INFO = -i.
101* > 0: The algorithm failed to compute the INFO/(N+1) th
102* eigenvalue while working on the submatrix lying in
103* global rows and columns mod(INFO,N+1).
104*
105* Further Details
106* ======= =======
107*
108* Contributed by Francoise Tisseur, University of Manchester.
109*
110* Reference: F. Tisseur and J. Dongarra, "A Parallel Divide and
111* Conquer Algorithm for the Symmetric Eigenvalue Problem
112* on Distributed Memory Architectures",
113* SIAM J. Sci. Comput., 6:20 (1999), pp. 2223--2236.
114* (see also LAPACK Working Note 132)
115* http://www.netlib.org/lapack/lawns/lawn132.ps
116*
117* =====================================================================
118*
119* .. Parameters ..
120 INTEGER BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_,
121 $ MB_, NB_, RSRC_, CSRC_, LLD_
122 parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
123 $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
124 $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
125 DOUBLE PRECISION ZERO, ONE
126 parameter( zero = 0.0d+0, one = 1.0d+0 )
127* ..
128* .. Local Scalars ..
129 LOGICAL LQUERY
130 INTEGER ICOFFQ, IIQ, IPQ, IQCOL, IQROW, IROFFQ, JJQ,
131 $ LDQ, LIWMIN, LWMIN, MYCOL, MYROW, NB, NP,
132 $ NPCOL, NPROW, NQ
133 DOUBLE PRECISION ORGNRM
134* ..
135* .. External Functions ..
136 LOGICAL LSAME
137 INTEGER INDXG2P, NUMROC
138 DOUBLE PRECISION DLANST
139 EXTERNAL indxg2p, lsame, numroc, dlanst
140* ..
141* .. External Subroutines ..
142 EXTERNAL blacs_gridinfo, chk1mat, dlascl, dstedc,
143 $ infog2l, pdlaed0, pdlasrt, pxerbla
144* ..
145* .. Intrinsic Functions ..
146 INTRINSIC dble, mod
147* ..
148* .. Executable Statements ..
149*
150* This is just to keep ftnchek and toolpack/1 happy
151 IF( block_cyclic_2d*csrc_*ctxt_*dlen_*dtype_*lld_*mb_*m_*nb_*n_*
152 $ rsrc_.LT.0 )RETURN
153*
154* Test the input parameters.
155*
156 CALL blacs_gridinfo( descq( ctxt_ ), nprow, npcol, myrow, mycol )
157 ldq = descq( lld_ )
158 nb = descq( nb_ )
159 np = numroc( n, nb, myrow, descq( rsrc_ ), nprow )
160 nq = numroc( n, nb, mycol, descq( csrc_ ), npcol )
161 info = 0
162 IF( nprow.EQ.-1 ) THEN
163 info = -( 600+ctxt_ )
164 ELSE
165 CALL chk1mat( n, 2, n, 2, iq, jq, descq, 8, info )
166 IF( info.EQ.0 ) THEN
167 nb = descq( nb_ )
168 iroffq = mod( iq-1, descq( mb_ ) )
169 icoffq = mod( jq-1, descq( nb_ ) )
170 iqrow = indxg2p( iq, nb, myrow, descq( rsrc_ ), nprow )
171 iqcol = indxg2p( jq, nb, mycol, descq( csrc_ ), npcol )
172 lwmin = 6*n + 2*np*nq
173 liwmin = 2 + 7*n + 8*npcol
174*
175 work( 1 ) = dble( lwmin )
176 iwork( 1 ) = liwmin
177 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
178 IF( .NOT.lsame( compz, 'I' ) ) THEN
179 info = -1
180 ELSE IF( n.LT.0 ) THEN
181 info = -2
182 ELSE IF( iroffq.NE.icoffq .OR. icoffq.NE.0 ) THEN
183 info = -5
184 ELSE IF( descq( mb_ ).NE.descq( nb_ ) ) THEN
185 info = -( 700+nb_ )
186 ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
187 info = -10
188 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
189 info = -12
190 END IF
191 END IF
192 END IF
193 IF( info.NE.0 ) THEN
194 CALL pxerbla( descq( ctxt_ ), 'PDSTEDC', -info )
195 RETURN
196 ELSE IF( lquery ) THEN
197 RETURN
198 END IF
199*
200* Quick return
201*
202 IF( n.EQ.0 )
203 $ GO TO 10
204 CALL infog2l( iq, jq, descq, nprow, npcol, myrow, mycol, iiq, jjq,
205 $ iqrow, iqcol )
206 IF( n.EQ.1 ) THEN
207 IF( myrow.EQ.iqrow .AND. mycol.EQ.iqcol )
208 $ q( 1 ) = one
209 GO TO 10
210 END IF
211*
212* If N is smaller than the minimum divide size NB, then
213* solve the problem with the serial divide and conquer
214* code locally.
215*
216 IF( n.LE.nb ) THEN
217 IF( ( myrow.EQ.iqrow ) .AND. ( mycol.EQ.iqcol ) ) THEN
218 ipq = iiq + ( jjq-1 )*ldq
219 CALL dstedc( 'I', n, d, e, q( ipq ), ldq, work, lwork,
220 $ iwork, liwork, info )
221 IF( info.NE.0 ) THEN
222 info = ( n+1 ) + n
223 GO TO 10
224 END IF
225 END IF
226 GO TO 10
227 END IF
228*
229* If P=NPROW*NPCOL=1, solve the problem with DSTEDC.
230*
231 IF( npcol*nprow.EQ.1 ) THEN
232 ipq = iiq + ( jjq-1 )*ldq
233 CALL dstedc( 'I', n, d, e, q( ipq ), ldq, work, lwork, iwork,
234 $ liwork, info )
235 GO TO 10
236 END IF
237*
238* Scale matrix to allowable range, if necessary.
239*
240 orgnrm = dlanst( 'M', n, d, e )
241 IF( orgnrm.NE.zero ) THEN
242 CALL dlascl( 'G', 0, 0, orgnrm, one, n, 1, d, n, info )
243 CALL dlascl( 'G', 0, 0, orgnrm, one, n-1, 1, e, n-1, info )
244 END IF
245*
246 CALL pdlaed0( n, d, e, q, iq, jq, descq, work, iwork, info )
247*
248* Sort eigenvalues and corresponding eigenvectors
249*
250 CALL pdlasrt( 'I', n, d, q, iq, jq, descq, work, lwork, iwork,
251 $ liwork, info )
252*
253* Scale back.
254*
255 IF( orgnrm.NE.zero )
256 $ CALL dlascl( 'G', 0, 0, one, orgnrm, n, 1, d, n, info )
257*
258 10 CONTINUE
259*
260 IF( lwork.GT.0 )
261 $ work( 1 ) = dble( lwmin )
262 IF( liwork.GT.0 )
263 $ iwork( 1 ) = liwmin
264 RETURN
265*
266* End of PDSTEDC
267*
chk1mat
subroutine chk1mat(ma, mapos0, na, napos0, ia, ja, desca, descapos0, info)
Definition chk1mat.f:3
indxg2p
integer function indxg2p(indxglob, nb, iproc, isrcproc, nprocs)
Definition indxg2p.f:2
infog2l
subroutine infog2l(grindx, gcindx, desc, nprow, npcol, myrow, mycol, lrindx, lcindx, rsrc, csrc)
Definition infog2l.f:3
numroc
integer function numroc(n, nb, iproc, isrcproc, nprocs)
Definition numroc.f:2
pdlaed0
subroutine pdlaed0(n, d, e, q, iq, jq, descq, work, iwork, info)
Definition pdlaed0.f:2
pdlasrt
subroutine pdlasrt(id, n, d, q, iq, jq, descq, work, lwork, iwork, liwork, info)
Definition pdlasrt.f:3
pxerbla
subroutine pxerbla(ictxt, srname, info)
Definition pxerbla.f:2
lsame
logical function lsame(ca, cb)
Definition tools.f:1724
Here is the call graph for this function:
Here is the caller graph for this function: