ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pdstedc.f
Go to the documentation of this file.
1  SUBROUTINE pdstedc( COMPZ, N, D, E, Q, IQ, JQ, DESCQ, WORK, LWORK,
2  $ IWORK, LIWORK, INFO )
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 *
268  END
infog2l
subroutine infog2l(GRINDX, GCINDX, DESC, NPROW, NPCOL, MYROW, MYCOL, LRINDX, LCINDX, RSRC, CSRC)
Definition: infog2l.f:3
pdstedc
subroutine pdstedc(COMPZ, N, D, E, Q, IQ, JQ, DESCQ, WORK, LWORK, IWORK, LIWORK, INFO)
Definition: pdstedc.f:3
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
chk1mat
subroutine chk1mat(MA, MAPOS0, NA, NAPOS0, IA, JA, DESCA, DESCAPOS0, INFO)
Definition: chk1mat.f:3
pxerbla
subroutine pxerbla(ICTXT, SRNAME, INFO)
Definition: pxerbla.f:2