KTH framework for Nek5000 toolboxes; testing version  0.0.1
dsygst.f
Go to the documentation of this file.
1  SUBROUTINE dsygst( 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 * September 30, 1994
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 * DSYGST 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**T)*A*inv(U) or inv(L)*A*inv(L**T)
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**T or L**T*A*L.
27 *
28 * B must have been previously factorized as U**T*U or L*L**T by DPOTRF.
29 *
30 * Arguments
31 * =========
32 *
33 * ITYPE (input) INTEGER
34 * = 1: compute inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T);
35 * = 2 or 3: compute U*A*U**T or L**T*A*L.
36 *
37 * UPLO (input) CHARACTER
38 * = 'U': Upper triangle of A is stored and B is factored as
39 * U**T*U;
40 * = 'L': Lower triangle of A is stored and B is factored as
41 * L*L**T.
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, KB, NB
81 * ..
82 * .. External Subroutines ..
83  EXTERNAL dsygs2, dsymm, dsyr2k, dtrmm, dtrsm, xerbla
84 * ..
85 * .. Intrinsic Functions ..
86  INTRINSIC max, min
87 * ..
88 * .. External Functions ..
89  LOGICAL LSAME
90  INTEGER ILAENV
91  EXTERNAL lsame, ilaenv
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( 'DSYGST', -info )
112  RETURN
113  END IF
114 *
115 * Quick return if possible
116 *
117  IF( n.EQ.0 )
118  $ RETURN
119 *
120 * Determine the block size for this environment.
121 *
122  nb = ilaenv( 1, 'DSYGST', uplo, n, -1, -1, -1 )
123 *
124  IF( nb.LE.1 .OR. nb.GE.n ) THEN
125 *
126 * Use unblocked code
127 *
128  CALL dsygs2( itype, uplo, n, a, lda, b, ldb, info )
129  ELSE
130 *
131 * Use blocked code
132 *
133  IF( itype.EQ.1 ) THEN
134  IF( upper ) THEN
135 *
136 * Compute inv(U')*A*inv(U)
137 *
138  DO 10 k = 1, n, nb
139  kb = min( n-k+1, nb )
140 *
141 * Update the upper triangle of A(k:n,k:n)
142 *
143  CALL dsygs2( itype, uplo, kb, a( k, k ), lda,
144  $ b( k, k ), ldb, info )
145  IF( k+kb.LE.n ) THEN
146  CALL dtrsm( 'Left', uplo, 'Transpose', 'Non-unit',
147  $ kb, n-k-kb+1, one, b( k, k ), ldb,
148  $ a( k, k+kb ), lda )
149  CALL dsymm( 'Left', uplo, kb, n-k-kb+1, -half,
150  $ a( k, k ), lda, b( k, k+kb ), ldb, one,
151  $ a( k, k+kb ), lda )
152  CALL dsyr2k( uplo, 'Transpose', n-k-kb+1, kb, -one,
153  $ a( k, k+kb ), lda, b( k, k+kb ), ldb,
154  $ one, a( k+kb, k+kb ), lda )
155  CALL dsymm( 'Left', uplo, kb, n-k-kb+1, -half,
156  $ a( k, k ), lda, b( k, k+kb ), ldb, one,
157  $ a( k, k+kb ), lda )
158  CALL dtrsm( 'Right', uplo, 'No transpose',
159  $ 'Non-unit', kb, n-k-kb+1, one,
160  $ b( k+kb, k+kb ), ldb, a( k, k+kb ),
161  $ lda )
162  END IF
163  10 CONTINUE
164  ELSE
165 *
166 * Compute inv(L)*A*inv(L')
167 *
168  DO 20 k = 1, n, nb
169  kb = min( n-k+1, nb )
170 *
171 * Update the lower triangle of A(k:n,k:n)
172 *
173  CALL dsygs2( itype, uplo, kb, a( k, k ), lda,
174  $ b( k, k ), ldb, info )
175  IF( k+kb.LE.n ) THEN
176  CALL dtrsm( 'Right', uplo, 'Transpose', 'Non-unit',
177  $ n-k-kb+1, kb, one, b( k, k ), ldb,
178  $ a( k+kb, k ), lda )
179  CALL dsymm( 'Right', uplo, n-k-kb+1, kb, -half,
180  $ a( k, k ), lda, b( k+kb, k ), ldb, one,
181  $ a( k+kb, k ), lda )
182  CALL dsyr2k( uplo, 'No transpose', n-k-kb+1, kb,
183  $ -one, a( k+kb, k ), lda, b( k+kb, k ),
184  $ ldb, one, a( k+kb, k+kb ), lda )
185  CALL dsymm( 'Right', uplo, n-k-kb+1, kb, -half,
186  $ a( k, k ), lda, b( k+kb, k ), ldb, one,
187  $ a( k+kb, k ), lda )
188  CALL dtrsm( 'Left', uplo, 'No transpose',
189  $ 'Non-unit', n-k-kb+1, kb, one,
190  $ b( k+kb, k+kb ), ldb, a( k+kb, k ),
191  $ lda )
192  END IF
193  20 CONTINUE
194  END IF
195  ELSE
196  IF( upper ) THEN
197 *
198 * Compute U*A*U'
199 *
200  DO 30 k = 1, n, nb
201  kb = min( n-k+1, nb )
202 *
203 * Update the upper triangle of A(1:k+kb-1,1:k+kb-1)
204 *
205  CALL dtrmm( 'Left', uplo, 'No transpose', 'Non-unit',
206  $ k-1, kb, one, b, ldb, a( 1, k ), lda )
207  CALL dsymm( 'Right', uplo, k-1, kb, half, a( k, k ),
208  $ lda, b( 1, k ), ldb, one, a( 1, k ), lda )
209  CALL dsyr2k( uplo, 'No transpose', k-1, kb, one,
210  $ a( 1, k ), lda, b( 1, k ), ldb, one, a,
211  $ lda )
212  CALL dsymm( 'Right', uplo, k-1, kb, half, a( k, k ),
213  $ lda, b( 1, k ), ldb, one, a( 1, k ), lda )
214  CALL dtrmm( 'Right', uplo, 'Transpose', 'Non-unit',
215  $ k-1, kb, one, b( k, k ), ldb, a( 1, k ),
216  $ lda )
217  CALL dsygs2( itype, uplo, kb, a( k, k ), lda,
218  $ b( k, k ), ldb, info )
219  30 CONTINUE
220  ELSE
221 *
222 * Compute L'*A*L
223 *
224  DO 40 k = 1, n, nb
225  kb = min( n-k+1, nb )
226 *
227 * Update the lower triangle of A(1:k+kb-1,1:k+kb-1)
228 *
229  CALL dtrmm( 'Right', uplo, 'No transpose', 'Non-unit',
230  $ kb, k-1, one, b, ldb, a( k, 1 ), lda )
231  CALL dsymm( 'Left', uplo, kb, k-1, half, a( k, k ),
232  $ lda, b( k, 1 ), ldb, one, a( k, 1 ), lda )
233  CALL dsyr2k( uplo, 'Transpose', k-1, kb, one,
234  $ a( k, 1 ), lda, b( k, 1 ), ldb, one, a,
235  $ lda )
236  CALL dsymm( 'Left', uplo, kb, k-1, half, a( k, k ),
237  $ lda, b( k, 1 ), ldb, one, a( k, 1 ), lda )
238  CALL dtrmm( 'Left', uplo, 'Transpose', 'Non-unit', kb,
239  $ k-1, one, b( k, k ), ldb, a( k, 1 ), lda )
240  CALL dsygs2( itype, uplo, kb, a( k, k ), lda,
241  $ b( k, k ), ldb, info )
242  40 CONTINUE
243  END IF
244  END IF
245  END IF
246  RETURN
247 *
248 * End of DSYGST
249 *
250  END
subroutine dsygs2(ITYPE, UPLO, N, A, LDA, B, LDB, INFO)
Definition: dsygs2.f:2
subroutine dsygst(ITYPE, UPLO, N, A, LDA, B, LDB, INFO)
Definition: dsygst.f:2
subroutine dsymm(SIDE, UPLO, M, N, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
Definition: dsymm.f:3
subroutine dsyr2k(UPLO, TRANS, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
Definition: dsyr2k.f:3
subroutine dtrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
Definition: dtrmm.f:3
subroutine dtrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
Definition: dtrsm.f:3
subroutine xerbla(SRNAME, INFO)
Definition: xerbla.f:2