ScaLAPACK 2.1  2.1 ScaLAPACK: Scalable Linear Algebra PACKage
pdsyevd.f
Go to the documentation of this file.
1  SUBROUTINE pdsyevd( 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  DOUBLE PRECISION A( * ), W( * ), WORK( * ), Z( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * PDSYEVD 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, PDSYEVD 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 DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension (N)
75 * If INFO=0, the eigenvalues in ascending order.
76 *
77 * Z (local output) DOUBLE PRECISION 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) DOUBLE PRECISION 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.
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  DOUBLE PRECISION ZERO, ONE
161  parameter( zero = 0.0d+0, one = 1.0d+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  DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
171  \$ smlnum
172 * ..
173 * .. Local Arrays ..
174 * ..
175  INTEGER IDUM1( 2 ), IDUM2( 2 )
176 * ..
177 * .. External Functions ..
178  LOGICAL LSAME
179  INTEGER INDXG2P, NUMROC
180  DOUBLE PRECISION PDLAMCH, PDLANSY
181  EXTERNAL lsame, indxg2p, numroc, pdlamch, pdlansy
182 * ..
183 * .. External Subroutines ..
184  EXTERNAL blacs_gridinfo, chk1mat, dscal, pchk1mat,
186  \$ pdsytrd, pxerbla
187 * ..
188 * .. Intrinsic Functions ..
189  INTRINSIC dble, ichar, max, min, mod, sqrt
190 * ..
191 * .. Executable Statements ..
192 * This is just to keep ftnchek and toolpack/1 happy
193  IF( block_cyclic_2d*csrc_*ctxt_*dlen_*dtype_*lld_*mb_*m_*nb_*n_*
194  \$ rsrc_.LT.0 )RETURN
195 *
196 * Quick return
197 *
198  IF( n.EQ.0 )
199  \$ RETURN
200 *
201 * Test the input arguments.
202 *
203  ictxt = descz( ctxt_ )
204  CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
205 *
206  info = 0
207  IF( nprow.EQ.-1 ) THEN
208  info = -( 600+ctxt_ )
209  ELSE
210  CALL chk1mat( n, 3, n, 3, ia, ja, desca, 7, info )
211  CALL chk1mat( n, 3, n, 3, iz, jz, descz, 12, info )
212  IF( info.EQ.0 ) THEN
213  upper = lsame( uplo, 'U' )
214  nb = desca( nb_ )
215  iroffa = mod( ia-1, desca( mb_ ) )
216  icoffa = mod( ja-1, desca( nb_ ) )
217  iroffz = mod( iz-1, descz( mb_ ) )
218  icoffz = mod( jz-1, descz( nb_ ) )
219  iarow = indxg2p( ia, nb, myrow, desca( rsrc_ ), nprow )
220  iacol = indxg2p( ja, nb, mycol, desca( csrc_ ), npcol )
221  np = numroc( n, nb, myrow, iarow, nprow )
222  nq = numroc( n, nb, mycol, iacol, npcol )
223 *
224  lquery = ( lwork.EQ.-1 )
225  trilwmin = 3*n + max( nb*( np+1 ), 3*nb )
226  lwmin = max( 1+6*n+2*np*nq, trilwmin ) + 2*n
227  liwmin = 7*n + 8*npcol + 2
228  work( 1 ) = dble( lwmin )
229  iwork( 1 ) = liwmin
230  IF( .NOT.lsame( jobz, 'V' ) ) THEN
231  info = -1
232  ELSE IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
233  info = -2
234  ELSE IF( iroffa.NE.icoffa .OR. icoffa.NE.0 ) THEN
235  info = -6
236  ELSE IF( iroffa.NE.iroffz .OR. icoffa.NE.icoffz ) THEN
237  info = -10
238  ELSE IF( desca( m_ ).NE.descz( m_ ) ) THEN
239  info = -( 1200+m_ )
240  ELSE IF( desca( mb_ ).NE.desca( nb_ ) ) THEN
241  info = -( 700+nb_ )
242  ELSE IF( descz( mb_ ).NE.descz( nb_ ) ) THEN
243  info = -( 1200+nb_ )
244  ELSE IF( desca( mb_ ).NE.descz( mb_ ) ) THEN
245  info = -( 1200+mb_ )
246  ELSE IF( desca( ctxt_ ).NE.descz( ctxt_ ) ) THEN
247  info = -( 1200+ctxt_ )
248  ELSE IF( desca( rsrc_ ).NE.descz( rsrc_ ) ) THEN
249  info = -( 1200+rsrc_ )
250  ELSE IF( desca( csrc_ ).NE.descz( csrc_ ) ) THEN
251  info = -( 1200+csrc_ )
252  ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
253  info = -14
254  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
255  info = -16
256  END IF
257  END IF
258  IF( upper ) THEN
259  idum1( 1 ) = ichar( 'U' )
260  ELSE
261  idum1( 1 ) = ichar( 'L' )
262  END IF
263  idum2( 1 ) = 2
264  IF( lwork.EQ.-1 ) THEN
265  idum1( 2 ) = -1
266  ELSE
267  idum1( 2 ) = 1
268  END IF
269  idum2( 2 ) = 14
270  CALL pchk1mat( n, 3, n, 3, ia, ja, desca, 7, 2, idum1, idum2,
271  \$ info )
272  END IF
273  IF( info.NE.0 ) THEN
274  CALL pxerbla( ictxt, 'PDSYEVD', -info )
275  RETURN
276  ELSE IF( lquery ) THEN
277  RETURN
278  END IF
279 *
280 * Set up pointers into the WORK array
281 *
282  indtau = 1
283  inde = indtau + n
284  indd = inde + n
285  inde2 = indd + n
286  indwork = inde2 + n
287  llwork = lwork - indwork + 1
288  indwork2 = indd
289  llwork2 = lwork - indwork2 + 1
290 *
291 * Scale matrix to allowable range, if necessary.
292 *
293  iscale = 0
294  safmin = pdlamch( desca( ctxt_ ), 'Safe minimum' )
295  eps = pdlamch( desca( ctxt_ ), 'Precision' )
296  smlnum = safmin / eps
297  bignum = one / smlnum
298  rmin = sqrt( smlnum )
299  rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
300  anrm = pdlansy( 'M', uplo, n, a, ia, ja, desca, work( indwork ) )
301 *
302  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
303  iscale = 1
304  sigma = rmin / anrm
305  ELSE IF( anrm.GT.rmax ) THEN
306  iscale = 1
307  sigma = rmax / anrm
308  END IF
309 *
310  IF( iscale.EQ.1 ) THEN
311  CALL pdlascl( uplo, one, sigma, n, n, a, ia, ja, desca, iinfo )
312  END IF
313 *
314 * Reduce symmetric matrix to tridiagonal form.
315 *
316 *
317  CALL pdsytrd( uplo, n, a, ia, ja, desca, work( indd ),
318  \$ work( inde2 ), work( indtau ), work( indwork ),
319  \$ llwork, iinfo )
320 *
321 * Copy the values of D, E to all processes.
322 *
323  CALL pdlared1d( n, ia, ja, desca, work( indd ), w,
324  \$ work( indwork ), llwork )
325 *
326  CALL pdlared1d( n, ia, ja, desca, work( inde2 ), work( inde ),
327  \$ work( indwork ), llwork )
328 *
329  CALL pdlaset( 'Full', n, n, zero, one, z, 1, 1, descz )
330 *
331  IF( upper ) THEN
332  offset = 1
333  ELSE
334  offset = 0
335  END IF
336  CALL pdstedc( 'I', n, w, work( inde+offset ), z, iz, jz, descz,
337  \$ work( indwork2 ), llwork2, iwork, liwork, info )
338 *
339  CALL pdormtr( 'L', uplo, 'N', n, n, a, ia, ja, desca,
340  \$ work( indtau ), z, iz, jz, descz, work( indwork2 ),
341  \$ llwork2, iinfo )
342 *
343 * If matrix was scaled, then rescale eigenvalues appropriately.
344 *
345  IF( iscale.EQ.1 ) THEN
346  CALL dscal( n, one / sigma, w, 1 )
347  END IF
348 *
349  RETURN
350 *
351 * End of PDSYEVD
352 *
353  END
max
#define max(A, B)
Definition: pcgemr.c:180
pdsyevd
subroutine pdsyevd(JOBZ, UPLO, N, A, IA, JA, DESCA, W, Z, IZ, JZ, DESCZ, WORK, LWORK, IWORK, LIWORK, INFO)
Definition: pdsyevd.f:3
pdstedc
subroutine pdstedc(COMPZ, N, D, E, Q, IQ, JQ, DESCQ, WORK, LWORK, IWORK, LIWORK, INFO)
Definition: pdstedc.f:3
pchk1mat
subroutine pchk1mat(MA, MAPOS0, NA, NAPOS0, IA, JA, DESCA, DESCAPOS0, NEXTRA, EX, EXPOS, INFO)
Definition: pchkxmat.f:3
pdsytrd
subroutine pdsytrd(UPLO, N, A, IA, JA, DESCA, D, E, TAU, WORK, LWORK, INFO)
Definition: pdsytrd.f:3
pdlascl
subroutine pdlascl(TYPE, CFROM, CTO, M, N, A, IA, JA, DESCA, INFO)
Definition: pdlascl.f:3
pdlaset
subroutine pdlaset(UPLO, M, N, ALPHA, BETA, A, IA, JA, DESCA)
Definition: pdblastst.f:6862
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
pdlared1d
subroutine pdlared1d(N, IA, JA, DESC, BYCOL, BYALL, WORK, LWORK)
Definition: pdlared1d.f:2
pdormtr
subroutine pdormtr(SIDE, UPLO, TRANS, M, N, A, IA, JA, DESCA, TAU, C, IC, JC, DESCC, WORK, LWORK, INFO)
Definition: pdormtr.f:3
min
#define min(A, B)
Definition: pcgemr.c:181