LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dsysv.f
Go to the documentation of this file.
1 *> \brief <b> DSYSV computes the solution to system of linear equations A * X = B for SY matrices</b>
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download DSYSV + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsysv.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsysv.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsysv.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DSYSV( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK,
22 * LWORK, INFO )
23 *
24 * .. Scalar Arguments ..
25 * CHARACTER UPLO
26 * INTEGER INFO, LDA, LDB, LWORK, N, NRHS
27 * ..
28 * .. Array Arguments ..
29 * INTEGER IPIV( * )
30 * DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *>
39 *> DSYSV computes the solution to a real system of linear equations
40 *> A * X = B,
41 *> where A is an N-by-N symmetric matrix and X and B are N-by-NRHS
42 *> matrices.
43 *>
44 *> The diagonal pivoting method is used to factor A as
45 *> A = U * D * U**T, if UPLO = 'U', or
46 *> A = L * D * L**T, if UPLO = 'L',
47 *> where U (or L) is a product of permutation and unit upper (lower)
48 *> triangular matrices, and D is symmetric and block diagonal with
49 *> 1-by-1 and 2-by-2 diagonal blocks. The factored form of A is then
50 *> used to solve the system of equations A * X = B.
51 *> \endverbatim
52 *
53 * Arguments:
54 * ==========
55 *
56 *> \param[in] UPLO
57 *> \verbatim
58 *> UPLO is CHARACTER*1
59 *> = 'U': Upper triangle of A is stored;
60 *> = 'L': Lower triangle of A is stored.
61 *> \endverbatim
62 *>
63 *> \param[in] N
64 *> \verbatim
65 *> N is INTEGER
66 *> The number of linear equations, i.e., the order of the
67 *> matrix A. N >= 0.
68 *> \endverbatim
69 *>
70 *> \param[in] NRHS
71 *> \verbatim
72 *> NRHS is INTEGER
73 *> The number of right hand sides, i.e., the number of columns
74 *> of the matrix B. NRHS >= 0.
75 *> \endverbatim
76 *>
77 *> \param[in,out] A
78 *> \verbatim
79 *> A is DOUBLE PRECISION array, dimension (LDA,N)
80 *> On entry, the symmetric matrix A. If UPLO = 'U', the leading
81 *> N-by-N upper triangular part of A contains the upper
82 *> triangular part of the matrix A, and the strictly lower
83 *> triangular part of A is not referenced. If UPLO = 'L', the
84 *> leading N-by-N lower triangular part of A contains the lower
85 *> triangular part of the matrix A, and the strictly upper
86 *> triangular part of A is not referenced.
87 *>
88 *> On exit, if INFO = 0, the block diagonal matrix D and the
89 *> multipliers used to obtain the factor U or L from the
90 *> factorization A = U*D*U**T or A = L*D*L**T as computed by
91 *> DSYTRF.
92 *> \endverbatim
93 *>
94 *> \param[in] LDA
95 *> \verbatim
96 *> LDA is INTEGER
97 *> The leading dimension of the array A. LDA >= max(1,N).
98 *> \endverbatim
99 *>
100 *> \param[out] IPIV
101 *> \verbatim
102 *> IPIV is INTEGER array, dimension (N)
103 *> Details of the interchanges and the block structure of D, as
104 *> determined by DSYTRF. If IPIV(k) > 0, then rows and columns
105 *> k and IPIV(k) were interchanged, and D(k,k) is a 1-by-1
106 *> diagonal block. If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0,
107 *> then rows and columns k-1 and -IPIV(k) were interchanged and
108 *> D(k-1:k,k-1:k) is a 2-by-2 diagonal block. If UPLO = 'L' and
109 *> IPIV(k) = IPIV(k+1) < 0, then rows and columns k+1 and
110 *> -IPIV(k) were interchanged and D(k:k+1,k:k+1) is a 2-by-2
111 *> diagonal block.
112 *> \endverbatim
113 *>
114 *> \param[in,out] B
115 *> \verbatim
116 *> B is DOUBLE PRECISION array, dimension (LDB,NRHS)
117 *> On entry, the N-by-NRHS right hand side matrix B.
118 *> On exit, if INFO = 0, the N-by-NRHS solution matrix X.
119 *> \endverbatim
120 *>
121 *> \param[in] LDB
122 *> \verbatim
123 *> LDB is INTEGER
124 *> The leading dimension of the array B. LDB >= max(1,N).
125 *> \endverbatim
126 *>
127 *> \param[out] WORK
128 *> \verbatim
129 *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
130 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
131 *> \endverbatim
132 *>
133 *> \param[in] LWORK
134 *> \verbatim
135 *> LWORK is INTEGER
136 *> The length of WORK. LWORK >= 1, and for best performance
137 *> LWORK >= max(1,N*NB), where NB is the optimal blocksize for
138 *> DSYTRF.
139 *> for LWORK < N, TRS will be done with Level BLAS 2
140 *> for LWORK >= N, TRS will be done with Level BLAS 3
141 *>
142 *> If LWORK = -1, then a workspace query is assumed; the routine
143 *> only calculates the optimal size of the WORK array, returns
144 *> this value as the first entry of the WORK array, and no error
145 *> message related to LWORK is issued by XERBLA.
146 *> \endverbatim
147 *>
148 *> \param[out] INFO
149 *> \verbatim
150 *> INFO is INTEGER
151 *> = 0: successful exit
152 *> < 0: if INFO = -i, the i-th argument had an illegal value
153 *> > 0: if INFO = i, D(i,i) is exactly zero. The factorization
154 *> has been completed, but the block diagonal matrix D is
155 *> exactly singular, so the solution could not be computed.
156 *> \endverbatim
157 *
158 * Authors:
159 * ========
160 *
161 *> \author Univ. of Tennessee
162 *> \author Univ. of California Berkeley
163 *> \author Univ. of Colorado Denver
164 *> \author NAG Ltd.
165 *
166 *> \date November 2011
167 *
168 *> \ingroup doubleSYsolve
169 *
170 * =====================================================================
171  SUBROUTINE dsysv( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK,
172  $ lwork, info )
173 *
174 * -- LAPACK driver routine (version 3.4.0) --
175 * -- LAPACK is a software package provided by Univ. of Tennessee, --
176 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
177 * November 2011
178 *
179 * .. Scalar Arguments ..
180  CHARACTER uplo
181  INTEGER info, lda, ldb, lwork, n, nrhs
182 * ..
183 * .. Array Arguments ..
184  INTEGER ipiv( * )
185  DOUBLE PRECISION a( lda, * ), b( ldb, * ), work( * )
186 * ..
187 *
188 * =====================================================================
189 *
190 * .. Local Scalars ..
191  LOGICAL lquery
192  INTEGER lwkopt
193 * ..
194 * .. External Functions ..
195  LOGICAL lsame
196  EXTERNAL lsame
197 * ..
198 * .. External Subroutines ..
199  EXTERNAL xerbla, dsytrf, dsytrs, dsytrs2
200 * ..
201 * .. Intrinsic Functions ..
202  INTRINSIC max
203 * ..
204 * .. Executable Statements ..
205 *
206 * Test the input parameters.
207 *
208  info = 0
209  lquery = ( lwork.EQ.-1 )
210  IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lsame( uplo, 'L' ) ) THEN
211  info = -1
212  ELSE IF( n.LT.0 ) THEN
213  info = -2
214  ELSE IF( nrhs.LT.0 ) THEN
215  info = -3
216  ELSE IF( lda.LT.max( 1, n ) ) THEN
217  info = -5
218  ELSE IF( ldb.LT.max( 1, n ) ) THEN
219  info = -8
220  ELSE IF( lwork.LT.1 .AND. .NOT.lquery ) THEN
221  info = -10
222  END IF
223 *
224  IF( info.EQ.0 ) THEN
225  IF( n.EQ.0 ) THEN
226  lwkopt = 1
227  ELSE
228  CALL dsytrf( uplo, n, a, lda, ipiv, work, -1, info )
229  lwkopt = work(1)
230  END IF
231  work( 1 ) = lwkopt
232  END IF
233 *
234  IF( info.NE.0 ) THEN
235  CALL xerbla( 'DSYSV ', -info )
236  return
237  ELSE IF( lquery ) THEN
238  return
239  END IF
240 *
241 * Compute the factorization A = U*D*U**T or A = L*D*L**T.
242 *
243  CALL dsytrf( uplo, n, a, lda, ipiv, work, lwork, info )
244  IF( info.EQ.0 ) THEN
245 *
246 * Solve the system A*X = B, overwriting B with X.
247 *
248  IF ( lwork.LT.n ) THEN
249 *
250 * Solve with TRS ( Use Level BLAS 2)
251 *
252  CALL dsytrs( uplo, n, nrhs, a, lda, ipiv, b, ldb, info )
253 *
254  ELSE
255 *
256 * Solve with TRS2 ( Use Level BLAS 3)
257 *
258  CALL dsytrs2( uplo,n,nrhs,a,lda,ipiv,b,ldb,work,info )
259 *
260  END IF
261 *
262  END IF
263 *
264  work( 1 ) = lwkopt
265 *
266  return
267 *
268 * End of DSYSV
269 *
270  END