ScaLAPACK 2.1  2.1 ScaLAPACK: Scalable Linear Algebra PACKage
psstedc.f
Go to the documentation of this file.
1  SUBROUTINE psstedc( 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  REAL D( * ), E( * ), Q( * ), WORK( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 * PSSTEDC 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 SLAED3 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) REAL 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) REAL 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) REAL 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) REAL 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.
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  REAL ZERO, ONE
126  parameter( zero = 0.0e+0, one = 1.0e+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  REAL ORGNRM
134 * ..
135 * .. External Functions ..
136  LOGICAL LSAME
137  INTEGER INDXG2P, NUMROC
138  REAL SLANST
139  EXTERNAL indxg2p, lsame, numroc, slanst
140 * ..
141 * .. External Subroutines ..
142  EXTERNAL blacs_gridinfo, chk1mat, infog2l, pslaed0,
143  \$ pslasrt, pxerbla, slascl, sstedc
144 * ..
145 * .. Intrinsic Functions ..
146  INTRINSIC mod, real
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  work( 1 ) = real( lwmin )
175  iwork( 1 ) = liwmin
176  lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
177  IF( .NOT.lsame( compz, 'I' ) ) THEN
178  info = -1
179  ELSE IF( n.LT.0 ) THEN
180  info = -2
181  ELSE IF( iroffq.NE.icoffq .OR. icoffq.NE.0 ) THEN
182  info = -5
183  ELSE IF( descq( mb_ ).NE.descq( nb_ ) ) THEN
184  info = -( 700+nb_ )
185  ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
186  info = -10
187  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
188  info = -12
189  END IF
190  END IF
191  END IF
192  IF( info.NE.0 ) THEN
193  CALL pxerbla( descq( ctxt_ ), 'PSSTEDC', -info )
194  RETURN
195  ELSE IF( lquery ) THEN
196  RETURN
197  END IF
198 *
199 * Quick return
200 *
201  IF( n.EQ.0 )
202  \$ GO TO 10
203  CALL infog2l( iq, jq, descq, nprow, npcol, myrow, mycol, iiq, jjq,
204  \$ iqrow, iqcol )
205  IF( n.EQ.1 ) THEN
206  IF( myrow.EQ.iqrow .AND. mycol.EQ.iqcol )
207  \$ q( 1 ) = one
208  GO TO 10
209  END IF
210 *
211 * If N is smaller than the minimum divide size NB, then
212 * solve the problem with the serial divide and conquer
213 * code locally.
214 *
215  IF( n.LE.nb ) THEN
216  IF( ( myrow.EQ.iqrow ) .AND. ( mycol.EQ.iqcol ) ) THEN
217  ipq = iiq + ( jjq-1 )*ldq
218  CALL sstedc( 'I', n, d, e, q( ipq ), ldq, work, lwork,
219  \$ iwork, liwork, info )
220  IF( info.NE.0 ) THEN
221  info = ( n+1 ) + n
222  GO TO 10
223  END IF
224  END IF
225  GO TO 10
226  END IF
227 *
228 * If P=NPROW*NPCOL=1, solve the problem with SSTEDC.
229 *
230  IF( npcol*nprow.EQ.1 ) THEN
231  ipq = iiq + ( jjq-1 )*ldq
232  CALL sstedc( 'I', n, d, e, q( ipq ), ldq, work, lwork, iwork,
233  \$ liwork, info )
234  GO TO 10
235  END IF
236 *
237 * Scale matrix to allowable range, if necessary.
238 *
239  orgnrm = slanst( 'M', n, d, e )
240  IF( orgnrm.NE.zero ) THEN
241  CALL slascl( 'G', 0, 0, orgnrm, one, n, 1, d, n, info )
242  CALL slascl( 'G', 0, 0, orgnrm, one, n-1, 1, e, n-1, info )
243  END IF
244 *
245  CALL pslaed0( n, d, e, q, iq, jq, descq, work, iwork, info )
246 *
247 * Sort eigenvalues and corresponding eigenvectors
248 *
249  CALL pslasrt( 'I', n, d, q, iq, jq, descq, work, lwork, iwork,
250  \$ liwork, info )
251 *
252 * Scale back.
253 *
254  IF( orgnrm.NE.zero )
255  \$ CALL slascl( 'G', 0, 0, one, orgnrm, n, 1, d, n, info )
256 *
257  10 CONTINUE
258 *
259  IF( lwork.GT.0 )
260  \$ work( 1 ) = real( lwmin )
261  IF( liwork.GT.0 )
262  \$ iwork( 1 ) = liwmin
263  RETURN
264 *
265 * End of PSSTEDC
266 *
267  END
psstedc
subroutine psstedc(COMPZ, N, D, E, Q, IQ, JQ, DESCQ, WORK, LWORK, IWORK, LIWORK, INFO)
Definition: psstedc.f:3
infog2l
subroutine infog2l(GRINDX, GCINDX, DESC, NPROW, NPCOL, MYROW, MYCOL, LRINDX, LCINDX, RSRC, CSRC)
Definition: infog2l.f:3
pslaed0
subroutine pslaed0(N, D, E, Q, IQ, JQ, DESCQ, WORK, IWORK, INFO)
Definition: pslaed0.f:2
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
pslasrt
subroutine pslasrt(ID, N, D, Q, IQ, JQ, DESCQ, WORK, LWORK, IWORK, LIWORK, INFO)
Definition: pslasrt.f:3