ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pssyevd.f
Go to the documentation of this file.
1  SUBROUTINE pssyevd( JOBZ, UPLO, N, A, IA, JA, DESCA, W, Z, IZ, JZ,
2  $ DESCZ, WORK, LWORK, 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 14, 2000
8 *
9 * .. Scalar Arguments ..
10  CHARACTER JOBZ, UPLO
11  INTEGER IA, INFO, IZ, JA, JZ, LIWORK, LWORK, N
12 * ..
13 * .. Array Arguments ..
14  INTEGER DESCA( * ), DESCZ( * ), IWORK( * )
15  REAL A( * ), W( * ), WORK( * ), Z( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * PSSYEVD computes all the eigenvalues and eigenvectors
22 * of a real symmetric matrix A by calling the recommended sequence
23 * of ScaLAPACK routines.
24 *
25 * In its present form, PSSYEVD assumes a homogeneous system and makes
26 * no checks for consistency of the eigenvalues or eigenvectors across
27 * the different processes. Because of this, it is possible that a
28 * heterogeneous system may return incorrect results without any error
29 * messages.
30 *
31 * Arguments
32 * =========
33 *
34 * NP = the number of rows local to a given process.
35 * NQ = the number of columns local to a given process.
36 *
37 * JOBZ (input) CHARACTER*1
38 * = 'N': Compute eigenvalues only; (NOT IMPLEMENTED YET)
39 * = 'V': Compute eigenvalues and eigenvectors.
40 *
41 * UPLO (global input) CHARACTER*1
42 * Specifies whether the upper or lower triangular part of the
43 * symmetric matrix A is stored:
44 * = 'U': Upper triangular
45 * = 'L': Lower triangular
46 *
47 * N (global input) INTEGER
48 * The number of rows and columns to be operated on, i.e. the
49 * order of the distributed submatrix sub( A ). N >= 0.
50 *
51 * A (local input/workspace) block cyclic REAL array,
52 * global dimension (N, N), local dimension ( LLD_A,
53 * LOCc(JA+N-1) )
54 * On entry, the symmetric matrix A. If UPLO = 'U', only the
55 * upper triangular part of A is used to define the elements of
56 * the symmetric matrix. If UPLO = 'L', only the lower
57 * triangular part of A is used to define the elements of the
58 * symmetric matrix.
59 * On exit, the lower triangle (if UPLO='L') or the upper
60 * triangle (if UPLO='U') of A, including the diagonal, is
61 * destroyed.
62 *
63 * IA (global input) INTEGER
64 * A's global row index, which points to the beginning of the
65 * submatrix which is to be operated on.
66 *
67 * JA (global input) INTEGER
68 * A's global column index, which points to the beginning of
69 * the submatrix which is to be operated on.
70 *
71 * DESCA (global and local input) INTEGER array of dimension DLEN_.
72 * The array descriptor for the distributed matrix A.
73 *
74 * W (global output) REAL array, dimension (N)
75 * If INFO=0, the eigenvalues in ascending order.
76 *
77 * Z (local output) REAL array,
78 * global dimension (N, N),
79 * local dimension ( LLD_Z, LOCc(JZ+N-1) )
80 * Z contains the orthonormal eigenvectors
81 * of the symmetric matrix A.
82 *
83 * IZ (global input) INTEGER
84 * Z's global row index, which points to the beginning of the
85 * submatrix which is to be operated on.
86 *
87 * JZ (global input) INTEGER
88 * Z's global column index, which points to the beginning of
89 * the submatrix which is to be operated on.
90 *
91 * DESCZ (global and local input) INTEGER array of dimension DLEN_.
92 * The array descriptor for the distributed matrix Z.
93 * DESCZ( CTXT_ ) must equal DESCA( CTXT_ )
94 *
95 * WORK (local workspace/output) REAL array,
96 * dimension (LWORK)
97 * On output, WORK(1) returns the workspace required.
98 *
99 * LWORK (local input) INTEGER
100 * LWORK >= MAX( 1+6*N+2*NP*NQ, TRILWMIN ) + 2*N
101 * TRILWMIN = 3*N + MAX( NB*( NP+1 ), 3*NB )
102 * NP = NUMROC( N, NB, MYROW, IAROW, NPROW )
103 * NQ = NUMROC( N, NB, MYCOL, IACOL, NPCOL )
104 *
105 * If LWORK = -1, the LWORK is global input and a workspace
106 * query is assumed; the routine only calculates the minimum
107 * size for the WORK array. The required workspace is returned
108 * as the first element of WORK and no error message is issued
109 * by PXERBLA.
110 *
111 * IWORK (local workspace/output) INTEGER array, dimension (LIWORK)
112 * On exit, if LIWORK > 0, IWORK(1) returns the optimal LIWORK.
113 *
114 * LIWORK (input) INTEGER
115 * The dimension of the array IWORK.
116 * LIWORK = 7*N + 8*NPCOL + 2
117 *
118 * INFO (global output) INTEGER
119 * = 0: successful exit
120 * < 0: If the i-th argument is an array and the j-entry had
121 * an illegal value, then INFO = -(i*100+j), if the i-th
122 * argument is a scalar and had an illegal value, then
123 * INFO = -i.
124 * > 0: The algorithm failed to compute the INFO/(N+1) th
125 * eigenvalue while working on the submatrix lying in
126 * global rows and columns mod(INFO,N+1).
127 *
128 * Alignment requirements
129 * ======================
130 *
131 * The distributed submatrices sub( A ), sub( Z ) must verify
132 * some alignment properties, namely the following expression
133 * should be true:
134 * ( MB_A.EQ.NB_A.EQ.MB_Z.EQ.NB_Z .AND. IROFFA.EQ.ICOFFA .AND.
135 * IROFFA.EQ.0 .AND.IROFFA.EQ.IROFFZ. AND. IAROW.EQ.IZROW)
136 * with IROFFA = MOD( IA-1, MB_A )
137 * and ICOFFA = MOD( JA-1, NB_A ).
138 *
139 * Further Details
140 * ======= =======
141 *
142 * Contributed by Francoise Tisseur, University of Manchester.
143 *
144 * Reference: F. Tisseur and J. Dongarra, "A Parallel Divide and
145 * Conquer Algorithm for the Symmetric Eigenvalue Problem
146 * on Distributed Memory Architectures",
147 * SIAM J. Sci. Comput., 6:20 (1999), pp. 2223--2236.
148 * (see also LAPACK Working Note 132)
149 * http://www.netlib.org/lapack/lawns/lawn132.ps
150 *
151 * =====================================================================
152 *
153 * .. Parameters ..
154 *
155  INTEGER BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_,
156  $ mb_, nb_, rsrc_, csrc_, lld_
157  parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
158  $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
159  $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
160  REAL ZERO, ONE
161  parameter( zero = 0.0e+0, one = 1.0e+0 )
162 * ..
163 * .. Local Scalars ..
164  LOGICAL LQUERY, UPPER
165  INTEGER IACOL, IAROW, ICOFFA, ICOFFZ, ICTXT, IINFO,
166  $ indd, inde, inde2, indtau, indwork, indwork2,
167  $ iroffa, iroffz, iscale, liwmin, llwork,
168  $ llwork2, lwmin, mycol, myrow, nb, np, npcol,
169  $ nprow, nq, offset, trilwmin
170  REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
171  $ smlnum
172 * ..
173 * .. Local Arrays ..
174  INTEGER IDUM1( 2 ), IDUM2( 2 )
175 * ..
176 * .. External Functions ..
177  LOGICAL LSAME
178  INTEGER INDXG2P, NUMROC
179  REAL PSLAMCH, PSLANSY
180  EXTERNAL lsame, indxg2p, numroc, pslamch, pslansy
181 * ..
182 * .. External Subroutines ..
183  EXTERNAL blacs_gridinfo, chk1mat, pchk1mat, pslared1d,
184  $ pslascl, pslaset, psormtr, psstedc, pssytrd,
185  $ pxerbla, sscal
186 * ..
187 * .. Intrinsic Functions ..
188  INTRINSIC ichar, max, min, mod, real, sqrt
189 * ..
190 * .. Executable Statements ..
191 * This is just to keep ftnchek and toolpack/1 happy
192  IF( block_cyclic_2d*csrc_*ctxt_*dlen_*dtype_*lld_*mb_*m_*nb_*n_*
193  $ rsrc_.LT.0 )RETURN
194 *
195 * Quick return
196 *
197  IF( n.EQ.0 )
198  $ RETURN
199 *
200 * Test the input arguments.
201 *
202  ictxt = descz( ctxt_ )
203  CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
204 *
205  info = 0
206  IF( nprow.EQ.-1 ) THEN
207  info = -( 600+ctxt_ )
208  ELSE
209  CALL chk1mat( n, 3, n, 3, ia, ja, desca, 7, info )
210  CALL chk1mat( n, 3, n, 3, iz, jz, descz, 12, info )
211  IF( info.EQ.0 ) THEN
212  upper = lsame( uplo, 'U' )
213  nb = desca( nb_ )
214  iroffa = mod( ia-1, desca( mb_ ) )
215  icoffa = mod( ja-1, desca( nb_ ) )
216  iroffz = mod( iz-1, descz( mb_ ) )
217  icoffz = mod( jz-1, descz( nb_ ) )
218  iarow = indxg2p( ia, nb, myrow, desca( rsrc_ ), nprow )
219  iacol = indxg2p( ja, nb, mycol, desca( csrc_ ), npcol )
220  np = numroc( n, nb, myrow, iarow, nprow )
221  nq = numroc( n, nb, mycol, iacol, npcol )
222 *
223  lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
224  trilwmin = 3*n + max( nb*( np+1 ), 3*nb )
225  lwmin = max( 1+6*n+2*np*nq, trilwmin ) + 2*n
226  liwmin = 7*n + 8*npcol + 2
227  work( 1 ) = real( lwmin )
228  iwork( 1 ) = liwmin
229  IF( .NOT.lsame( jobz, 'V' ) ) THEN
230  info = -1
231  ELSE IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
232  info = -2
233  ELSE IF( iroffa.NE.icoffa .OR. icoffa.NE.0 ) THEN
234  info = -6
235  ELSE IF( iroffa.NE.iroffz .OR. icoffa.NE.icoffz ) THEN
236  info = -10
237  ELSE IF( desca( m_ ).NE.descz( m_ ) ) THEN
238  info = -( 1200+m_ )
239  ELSE IF( desca( mb_ ).NE.desca( nb_ ) ) THEN
240  info = -( 700+nb_ )
241  ELSE IF( descz( mb_ ).NE.descz( nb_ ) ) THEN
242  info = -( 1200+nb_ )
243  ELSE IF( desca( mb_ ).NE.descz( mb_ ) ) THEN
244  info = -( 1200+mb_ )
245  ELSE IF( desca( ctxt_ ).NE.descz( ctxt_ ) ) THEN
246  info = -( 1200+ctxt_ )
247  ELSE IF( desca( rsrc_ ).NE.descz( rsrc_ ) ) THEN
248  info = -( 1200+rsrc_ )
249  ELSE IF( desca( csrc_ ).NE.descz( csrc_ ) ) THEN
250  info = -( 1200+csrc_ )
251  ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
252  info = -14
253  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
254  info = -16
255  END IF
256  END IF
257  IF( upper ) THEN
258  idum1( 1 ) = ichar( 'U' )
259  ELSE
260  idum1( 1 ) = ichar( 'L' )
261  END IF
262  idum2( 1 ) = 2
263  IF( lwork.EQ.-1 ) THEN
264  idum1( 2 ) = -1
265  ELSE
266  idum1( 2 ) = 1
267  END IF
268  idum2( 2 ) = 14
269  CALL pchk1mat( n, 3, n, 3, ia, ja, desca, 7, 2, idum1, idum2,
270  $ info )
271  END IF
272  IF( info.NE.0 ) THEN
273  CALL pxerbla( ictxt, 'PSSYEVD', -info )
274  RETURN
275  ELSE IF( lquery ) THEN
276  RETURN
277  END IF
278 *
279 * Set up pointers into the WORK array
280 *
281  indtau = 1
282  inde = indtau + n
283  indd = inde + n
284  inde2 = indd + n
285  indwork = inde2 + n
286  llwork = lwork - indwork + 1
287  indwork2 = indd
288  llwork2 = lwork - indwork2 + 1
289 *
290 * Scale matrix to allowable range, if necessary.
291 *
292  iscale = 0
293  safmin = pslamch( desca( ctxt_ ), 'Safe minimum' )
294  eps = pslamch( desca( ctxt_ ), 'Precision' )
295  smlnum = safmin / eps
296  bignum = one / smlnum
297  rmin = sqrt( smlnum )
298  rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
299  anrm = pslansy( 'M', uplo, n, a, ia, ja, desca, work( indwork ) )
300 *
301  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
302  iscale = 1
303  sigma = rmin / anrm
304  ELSE IF( anrm.GT.rmax ) THEN
305  iscale = 1
306  sigma = rmax / anrm
307  END IF
308 *
309  IF( iscale.EQ.1 ) THEN
310  CALL pslascl( uplo, one, sigma, n, n, a, ia, ja, desca, iinfo )
311  END IF
312 *
313 * Reduce symmetric matrix to tridiagonal form.
314 *
315 *
316  CALL pssytrd( uplo, n, a, ia, ja, desca, work( indd ),
317  $ work( inde2 ), work( indtau ), work( indwork ),
318  $ llwork, iinfo )
319 *
320 * Copy the values of D, E to all processes.
321 *
322  CALL pslared1d( n, ia, ja, desca, work( indd ), w,
323  $ work( indwork ), llwork )
324 *
325  CALL pslared1d( n, ia, ja, desca, work( inde2 ), work( inde ),
326  $ work( indwork ), llwork )
327 *
328  CALL pslaset( 'Full', n, n, zero, one, z, 1, 1, descz )
329 *
330  IF( upper ) THEN
331  offset = 1
332  ELSE
333  offset = 0
334  END IF
335  CALL psstedc( 'I', n, w, work( inde+offset ), z, iz, jz, descz,
336  $ work( indwork2 ), llwork2, iwork, liwork, info )
337 *
338  CALL psormtr( 'L', uplo, 'N', n, n, a, ia, ja, desca,
339  $ work( indtau ), z, iz, jz, descz, work( indwork2 ),
340  $ llwork2, iinfo )
341 *
342 * If matrix was scaled, then rescale eigenvalues appropriately.
343 *
344  IF( iscale.EQ.1 ) THEN
345  CALL sscal( n, one / sigma, w, 1 )
346  END IF
347 *
348  RETURN
349 *
350 * End of PSSYEVD
351 *
352  END
max
#define max(A, B)
Definition: pcgemr.c:180
psstedc
subroutine psstedc(COMPZ, N, D, E, Q, IQ, JQ, DESCQ, WORK, LWORK, IWORK, LIWORK, INFO)
Definition: psstedc.f:3
pslascl
subroutine pslascl(TYPE, CFROM, CTO, M, N, A, IA, JA, DESCA, INFO)
Definition: pslascl.f:3
pssytrd
subroutine pssytrd(UPLO, N, A, IA, JA, DESCA, D, E, TAU, WORK, LWORK, INFO)
Definition: pssytrd.f:3
pssyevd
subroutine pssyevd(JOBZ, UPLO, N, A, IA, JA, DESCA, W, Z, IZ, JZ, DESCZ, WORK, LWORK, IWORK, LIWORK, INFO)
Definition: pssyevd.f:3
pchk1mat
subroutine pchk1mat(MA, MAPOS0, NA, NAPOS0, IA, JA, DESCA, DESCAPOS0, NEXTRA, EX, EXPOS, INFO)
Definition: pchkxmat.f:3
pslared1d
subroutine pslared1d(N, IA, JA, DESC, BYCOL, BYALL, WORK, LWORK)
Definition: pslared1d.f:2
psormtr
subroutine psormtr(SIDE, UPLO, TRANS, M, N, A, IA, JA, DESCA, TAU, C, IC, JC, DESCC, WORK, LWORK, INFO)
Definition: psormtr.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
pslaset
subroutine pslaset(UPLO, M, N, ALPHA, BETA, A, IA, JA, DESCA)
Definition: psblastst.f:6863
min
#define min(A, B)
Definition: pcgemr.c:181