KTH framework for Nek5000 toolboxes; testing version  0.0.1
dsygs2.f
Go to the documentation of this file.
1  SUBROUTINE dsygs2( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
2 *
3 * -- LAPACK routine (version 3.0) --
4 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
5 * Courant Institute, Argonne National Lab, and Rice University
6 * February 29, 1992
7 *
8 * .. Scalar Arguments ..
9  CHARACTER UPLO
10  INTEGER INFO, ITYPE, LDA, LDB, N
11 * ..
12 * .. Array Arguments ..
13  DOUBLE PRECISION A( LDA, * ), B( LDB, * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DSYGS2 reduces a real symmetric-definite generalized eigenproblem
20 * to standard form.
21 *
22 * If ITYPE = 1, the problem is A*x = lambda*B*x,
23 * and A is overwritten by inv(U')*A*inv(U) or inv(L)*A*inv(L')
24 *
25 * If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or
26 * B*A*x = lambda*x, and A is overwritten by U*A*U` or L'*A*L.
27 *
28 * B must have been previously factorized as U'*U or L*L' by DPOTRF.
29 *
30 * Arguments
31 * =========
32 *
33 * ITYPE (input) INTEGER
34 * = 1: compute inv(U')*A*inv(U) or inv(L)*A*inv(L');
35 * = 2 or 3: compute U*A*U' or L'*A*L.
36 *
37 * UPLO (input) CHARACTER
38 * Specifies whether the upper or lower triangular part of the
39 * symmetric matrix A is stored, and how B has been factorized.
40 * = 'U': Upper triangular
41 * = 'L': Lower triangular
42 *
43 * N (input) INTEGER
44 * The order of the matrices A and B. N >= 0.
45 *
46 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
47 * On entry, the symmetric matrix A. If UPLO = 'U', the leading
48 * n by n upper triangular part of A contains the upper
49 * triangular part of the matrix A, and the strictly lower
50 * triangular part of A is not referenced. If UPLO = 'L', the
51 * leading n by n lower triangular part of A contains the lower
52 * triangular part of the matrix A, and the strictly upper
53 * triangular part of A is not referenced.
54 *
55 * On exit, if INFO = 0, the transformed matrix, stored in the
56 * same format as A.
57 *
58 * LDA (input) INTEGER
59 * The leading dimension of the array A. LDA >= max(1,N).
60 *
61 * B (input) DOUBLE PRECISION array, dimension (LDB,N)
62 * The triangular factor from the Cholesky factorization of B,
63 * as returned by DPOTRF.
64 *
65 * LDB (input) INTEGER
66 * The leading dimension of the array B. LDB >= max(1,N).
67 *
68 * INFO (output) INTEGER
69 * = 0: successful exit.
70 * < 0: if INFO = -i, the i-th argument had an illegal value.
71 *
72 * =====================================================================
73 *
74 * .. Parameters ..
75  DOUBLE PRECISION ONE, HALF
76  parameter( one = 1.0d0, half = 0.5d0 )
77 * ..
78 * .. Local Scalars ..
79  LOGICAL UPPER
80  INTEGER K
81  DOUBLE PRECISION AKK, BKK, CT
82 * ..
83 * .. External Subroutines ..
84  EXTERNAL daxpy, dscal, dsyr2, dtrmv, dtrsv, xerbla
85 * ..
86 * .. Intrinsic Functions ..
87  INTRINSIC max
88 * ..
89 * .. External Functions ..
90  LOGICAL LSAME
91  EXTERNAL lsame
92 * ..
93 * .. Executable Statements ..
94 *
95 * Test the input parameters.
96 *
97  info = 0
98  upper = lsame( uplo, 'U' )
99  IF( itype.LT.1 .OR. itype.GT.3 ) THEN
100  info = -1
101  ELSE IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
102  info = -2
103  ELSE IF( n.LT.0 ) THEN
104  info = -3
105  ELSE IF( lda.LT.max( 1, n ) ) THEN
106  info = -5
107  ELSE IF( ldb.LT.max( 1, n ) ) THEN
108  info = -7
109  END IF
110  IF( info.NE.0 ) THEN
111  CALL xerbla( 'DSYGS2', -info )
112  RETURN
113  END IF
114 *
115  IF( itype.EQ.1 ) THEN
116  IF( upper ) THEN
117 *
118 * Compute inv(U')*A*inv(U)
119 *
120  DO 10 k = 1, n
121 *
122 * Update the upper triangle of A(k:n,k:n)
123 *
124  akk = a( k, k )
125  bkk = b( k, k )
126  akk = akk / bkk**2
127  a( k, k ) = akk
128  IF( k.LT.n ) THEN
129  CALL dscal( n-k, one / bkk, a( k, k+1 ), lda )
130  ct = -half*akk
131  CALL daxpy( n-k, ct, b( k, k+1 ), ldb, a( k, k+1 ),
132  $ lda )
133  CALL dsyr2( uplo, n-k, -one, a( k, k+1 ), lda,
134  $ b( k, k+1 ), ldb, a( k+1, k+1 ), lda )
135  CALL daxpy( n-k, ct, b( k, k+1 ), ldb, a( k, k+1 ),
136  $ lda )
137  CALL dtrsv( uplo, 'Transpose', 'Non-unit', n-k,
138  $ b( k+1, k+1 ), ldb, a( k, k+1 ), lda )
139  END IF
140  10 CONTINUE
141  ELSE
142 *
143 * Compute inv(L)*A*inv(L')
144 *
145  DO 20 k = 1, n
146 *
147 * Update the lower triangle of A(k:n,k:n)
148 *
149  akk = a( k, k )
150  bkk = b( k, k )
151  akk = akk / bkk**2
152  a( k, k ) = akk
153  IF( k.LT.n ) THEN
154  CALL dscal( n-k, one / bkk, a( k+1, k ), 1 )
155  ct = -half*akk
156  CALL daxpy( n-k, ct, b( k+1, k ), 1, a( k+1, k ), 1 )
157  CALL dsyr2( uplo, n-k, -one, a( k+1, k ), 1,
158  $ b( k+1, k ), 1, a( k+1, k+1 ), lda )
159  CALL daxpy( n-k, ct, b( k+1, k ), 1, a( k+1, k ), 1 )
160  CALL dtrsv( uplo, 'No transpose', 'Non-unit', n-k,
161  $ b( k+1, k+1 ), ldb, a( k+1, k ), 1 )
162  END IF
163  20 CONTINUE
164  END IF
165  ELSE
166  IF( upper ) THEN
167 *
168 * Compute U*A*U'
169 *
170  DO 30 k = 1, n
171 *
172 * Update the upper triangle of A(1:k,1:k)
173 *
174  akk = a( k, k )
175  bkk = b( k, k )
176  CALL dtrmv( uplo, 'No transpose', 'Non-unit', k-1, b,
177  $ ldb, a( 1, k ), 1 )
178  ct = half*akk
179  CALL daxpy( k-1, ct, b( 1, k ), 1, a( 1, k ), 1 )
180  CALL dsyr2( uplo, k-1, one, a( 1, k ), 1, b( 1, k ), 1,
181  $ a, lda )
182  CALL daxpy( k-1, ct, b( 1, k ), 1, a( 1, k ), 1 )
183  CALL dscal( k-1, bkk, a( 1, k ), 1 )
184  a( k, k ) = akk*bkk**2
185  30 CONTINUE
186  ELSE
187 *
188 * Compute L'*A*L
189 *
190  DO 40 k = 1, n
191 *
192 * Update the lower triangle of A(1:k,1:k)
193 *
194  akk = a( k, k )
195  bkk = b( k, k )
196  CALL dtrmv( uplo, 'Transpose', 'Non-unit', k-1, b, ldb,
197  $ a( k, 1 ), lda )
198  ct = half*akk
199  CALL daxpy( k-1, ct, b( k, 1 ), ldb, a( k, 1 ), lda )
200  CALL dsyr2( uplo, k-1, one, a( k, 1 ), lda, b( k, 1 ),
201  $ ldb, a, lda )
202  CALL daxpy( k-1, ct, b( k, 1 ), ldb, a( k, 1 ), lda )
203  CALL dscal( k-1, bkk, a( k, 1 ), lda )
204  a( k, k ) = akk*bkk**2
205  40 CONTINUE
206  END IF
207  END IF
208  RETURN
209 *
210 * End of DSYGS2
211 *
212  END
subroutine daxpy(n, da, dx, incx, dy, incy)
Definition: daxpy.f:2
subroutine dscal(n, da, dx, incx)
Definition: dscal.f:2
subroutine dsygs2(ITYPE, UPLO, N, A, LDA, B, LDB, INFO)
Definition: dsygs2.f:2
subroutine dsyr2(UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA)
Definition: dsyr2.f:2
subroutine dtrmv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
Definition: dtrmv.f:2
subroutine dtrsv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
Definition: dtrsv.f:2
subroutine xerbla(SRNAME, INFO)
Definition: xerbla.f:2