001:       SUBROUTINE SSPGVD( ITYPE, JOBZ, UPLO, N, AP, BP, W, Z, LDZ, WORK,
002:      $                   LWORK, IWORK, LIWORK, INFO )
003: *
004: *  -- LAPACK driver routine (version 3.2) --
005: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
006: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       CHARACTER          JOBZ, UPLO
011:       INTEGER            INFO, ITYPE, LDZ, LIWORK, LWORK, N
012: *     ..
013: *     .. Array Arguments ..
014:       INTEGER            IWORK( * )
015:       REAL               AP( * ), BP( * ), W( * ), WORK( * ),
016:      $                   Z( LDZ, * )
017: *     ..
018: *
019: *  Purpose
020: *  =======
021: *
022: *  SSPGVD computes all the eigenvalues, and optionally, the eigenvectors
023: *  of a real generalized symmetric-definite eigenproblem, of the form
024: *  A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.  Here A and
025: *  B are assumed to be symmetric, stored in packed format, and B is also
026: *  positive definite.
027: *  If eigenvectors are desired, it uses a divide and conquer algorithm.
028: *
029: *  The divide and conquer algorithm makes very mild assumptions about
030: *  floating point arithmetic. It will work on machines with a guard
031: *  digit in add/subtract, or on those binary machines without guard
032: *  digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
033: *  Cray-2. It could conceivably fail on hexadecimal or decimal machines
034: *  without guard digits, but we know of none.
035: *
036: *  Arguments
037: *  =========
038: *
039: *  ITYPE   (input) INTEGER
040: *          Specifies the problem type to be solved:
041: *          = 1:  A*x = (lambda)*B*x
042: *          = 2:  A*B*x = (lambda)*x
043: *          = 3:  B*A*x = (lambda)*x
044: *
045: *  JOBZ    (input) CHARACTER*1
046: *          = 'N':  Compute eigenvalues only;
047: *          = 'V':  Compute eigenvalues and eigenvectors.
048: *
049: *  UPLO    (input) CHARACTER*1
050: *          = 'U':  Upper triangles of A and B are stored;
051: *          = 'L':  Lower triangles of A and B are stored.
052: *
053: *  N       (input) INTEGER
054: *          The order of the matrices A and B.  N >= 0.
055: *
056: *  AP      (input/output) REAL array, dimension (N*(N+1)/2)
057: *          On entry, the upper or lower triangle of the symmetric matrix
058: *          A, packed columnwise in a linear array.  The j-th column of A
059: *          is stored in the array AP as follows:
060: *          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
061: *          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
062: *
063: *          On exit, the contents of AP are destroyed.
064: *
065: *  BP      (input/output) REAL array, dimension (N*(N+1)/2)
066: *          On entry, the upper or lower triangle of the symmetric matrix
067: *          B, packed columnwise in a linear array.  The j-th column of B
068: *          is stored in the array BP as follows:
069: *          if UPLO = 'U', BP(i + (j-1)*j/2) = B(i,j) for 1<=i<=j;
070: *          if UPLO = 'L', BP(i + (j-1)*(2*n-j)/2) = B(i,j) for j<=i<=n.
071: *
072: *          On exit, the triangular factor U or L from the Cholesky
073: *          factorization B = U**T*U or B = L*L**T, in the same storage
074: *          format as B.
075: *
076: *  W       (output) REAL array, dimension (N)
077: *          If INFO = 0, the eigenvalues in ascending order.
078: *
079: *  Z       (output) REAL array, dimension (LDZ, N)
080: *          If JOBZ = 'V', then if INFO = 0, Z contains the matrix Z of
081: *          eigenvectors.  The eigenvectors are normalized as follows:
082: *          if ITYPE = 1 or 2, Z**T*B*Z = I;
083: *          if ITYPE = 3, Z**T*inv(B)*Z = I.
084: *          If JOBZ = 'N', then Z is not referenced.
085: *
086: *  LDZ     (input) INTEGER
087: *          The leading dimension of the array Z.  LDZ >= 1, and if
088: *          JOBZ = 'V', LDZ >= max(1,N).
089: *
090: *  WORK    (workspace/output) REAL array, dimension (MAX(1,LWORK))
091: *          On exit, if INFO = 0, WORK(1) returns the required LWORK.
092: *
093: *  LWORK   (input) INTEGER
094: *          The dimension of the array WORK.
095: *          If N <= 1,               LWORK >= 1.
096: *          If JOBZ = 'N' and N > 1, LWORK >= 2*N.
097: *          If JOBZ = 'V' and N > 1, LWORK >= 1 + 6*N + 2*N**2.
098: *
099: *          If LWORK = -1, then a workspace query is assumed; the routine
100: *          only calculates the required sizes of the WORK and IWORK
101: *          arrays, returns these values as the first entries of the WORK
102: *          and IWORK arrays, and no error message related to LWORK or
103: *          LIWORK is issued by XERBLA.
104: *
105: *  IWORK   (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
106: *          On exit, if INFO = 0, IWORK(1) returns the required LIWORK.
107: *
108: *  LIWORK  (input) INTEGER
109: *          The dimension of the array IWORK.
110: *          If JOBZ  = 'N' or N <= 1, LIWORK >= 1.
111: *          If JOBZ  = 'V' and N > 1, LIWORK >= 3 + 5*N.
112: *
113: *          If LIWORK = -1, then a workspace query is assumed; the
114: *          routine only calculates the required sizes of the WORK and
115: *          IWORK arrays, returns these values as the first entries of
116: *          the WORK and IWORK arrays, and no error message related to
117: *          LWORK or LIWORK is issued by XERBLA.
118: *
119: *  INFO    (output) INTEGER
120: *          = 0:  successful exit
121: *          < 0:  if INFO = -i, the i-th argument had an illegal value
122: *          > 0:  SPPTRF or SSPEVD returned an error code:
123: *             <= N:  if INFO = i, SSPEVD failed to converge;
124: *                    i off-diagonal elements of an intermediate
125: *                    tridiagonal form did not converge to zero;
126: *             > N:   if INFO = N + i, for 1 <= i <= N, then the leading
127: *                    minor of order i of B is not positive definite.
128: *                    The factorization of B could not be completed and
129: *                    no eigenvalues or eigenvectors were computed.
130: *
131: *  Further Details
132: *  ===============
133: *
134: *  Based on contributions by
135: *     Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
136: *
137: *  =====================================================================
138: *
139: *     .. Parameters ..
140:       REAL               TWO
141:       PARAMETER          ( TWO = 2.0E+0 )
142: *     ..
143: *     .. Local Scalars ..
144:       LOGICAL            LQUERY, UPPER, WANTZ
145:       CHARACTER          TRANS
146:       INTEGER            J, LIWMIN, LWMIN, NEIG
147: *     ..
148: *     .. External Functions ..
149:       LOGICAL            LSAME
150:       EXTERNAL           LSAME
151: *     ..
152: *     .. External Subroutines ..
153:       EXTERNAL           SPPTRF, SSPEVD, SSPGST, STPMV, STPSV, XERBLA
154: *     ..
155: *     .. Intrinsic Functions ..
156:       INTRINSIC          MAX, REAL
157: *     ..
158: *     .. Executable Statements ..
159: *
160: *     Test the input parameters.
161: *
162:       WANTZ = LSAME( JOBZ, 'V' )
163:       UPPER = LSAME( UPLO, 'U' )
164:       LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
165: *
166:       INFO = 0
167:       IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
168:          INFO = -1
169:       ELSE IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
170:          INFO = -2
171:       ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN
172:          INFO = -3
173:       ELSE IF( N.LT.0 ) THEN
174:          INFO = -4
175:       ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
176:          INFO = -9
177:       END IF
178: *
179:       IF( INFO.EQ.0 ) THEN
180:          IF( N.LE.1 ) THEN
181:             LIWMIN = 1
182:             LWMIN = 1
183:          ELSE
184:             IF( WANTZ ) THEN
185:                LIWMIN = 3 + 5*N
186:                LWMIN = 1 + 6*N + 2*N**2
187:             ELSE
188:                LIWMIN = 1
189:                LWMIN = 2*N
190:             END IF
191:          END IF
192:          WORK( 1 ) = LWMIN
193:          IWORK( 1 ) = LIWMIN
194: *
195:          IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
196:             INFO = -11
197:          ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
198:             INFO = -13
199:          END IF
200:       END IF
201: *
202:       IF( INFO.NE.0 ) THEN
203:          CALL XERBLA( 'SSPGVD', -INFO )
204:          RETURN
205:       ELSE IF( LQUERY ) THEN
206:          RETURN
207:       END IF
208: *
209: *     Quick return if possible
210: *
211:       IF( N.EQ.0 )
212:      $   RETURN
213: *
214: *     Form a Cholesky factorization of BP.
215: *
216:       CALL SPPTRF( UPLO, N, BP, INFO )
217:       IF( INFO.NE.0 ) THEN
218:          INFO = N + INFO
219:          RETURN
220:       END IF
221: *
222: *     Transform problem to standard eigenvalue problem and solve.
223: *
224:       CALL SSPGST( ITYPE, UPLO, N, AP, BP, INFO )
225:       CALL SSPEVD( JOBZ, UPLO, N, AP, W, Z, LDZ, WORK, LWORK, IWORK,
226:      $             LIWORK, INFO )
227:       LWMIN = MAX( REAL( LWMIN ), REAL( WORK( 1 ) ) )
228:       LIWMIN = MAX( REAL( LIWMIN ), REAL( IWORK( 1 ) ) )
229: *
230:       IF( WANTZ ) THEN
231: *
232: *        Backtransform eigenvectors to the original problem.
233: *
234:          NEIG = N
235:          IF( INFO.GT.0 )
236:      $      NEIG = INFO - 1
237:          IF( ITYPE.EQ.1 .OR. ITYPE.EQ.2 ) THEN
238: *
239: *           For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
240: *           backtransform eigenvectors: x = inv(L)'*y or inv(U)*y
241: *
242:             IF( UPPER ) THEN
243:                TRANS = 'N'
244:             ELSE
245:                TRANS = 'T'
246:             END IF
247: *
248:             DO 10 J = 1, NEIG
249:                CALL STPSV( UPLO, TRANS, 'Non-unit', N, BP, Z( 1, J ),
250:      $                     1 )
251:    10       CONTINUE
252: *
253:          ELSE IF( ITYPE.EQ.3 ) THEN
254: *
255: *           For B*A*x=(lambda)*x;
256: *           backtransform eigenvectors: x = L*y or U'*y
257: *
258:             IF( UPPER ) THEN
259:                TRANS = 'T'
260:             ELSE
261:                TRANS = 'N'
262:             END IF
263: *
264:             DO 20 J = 1, NEIG
265:                CALL STPMV( UPLO, TRANS, 'Non-unit', N, BP, Z( 1, J ),
266:      $                     1 )
267:    20       CONTINUE
268:          END IF
269:       END IF
270: *
271:       WORK( 1 ) = LWMIN
272:       IWORK( 1 ) = LIWMIN
273: *
274:       RETURN
275: *
276: *     End of SSPGVD
277: *
278:       END
279: