LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
ssysvx.f
Go to the documentation of this file.
1 *> \brief <b> SSYSVX 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 SSYSVX + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssysvx.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssysvx.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssysvx.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SSYSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B,
22 * LDB, X, LDX, RCOND, FERR, BERR, WORK, LWORK,
23 * IWORK, INFO )
24 *
25 * .. Scalar Arguments ..
26 * CHARACTER FACT, UPLO
27 * INTEGER INFO, LDA, LDAF, LDB, LDX, LWORK, N, NRHS
28 * REAL RCOND
29 * ..
30 * .. Array Arguments ..
31 * INTEGER IPIV( * ), IWORK( * )
32 * REAL A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
33 * $ BERR( * ), FERR( * ), WORK( * ), X( LDX, * )
34 * ..
35 *
36 *
37 *> \par Purpose:
38 * =============
39 *>
40 *> \verbatim
41 *>
42 *> SSYSVX uses the diagonal pivoting factorization to compute the
43 *> solution to a real system of linear equations A * X = B,
44 *> where A is an N-by-N symmetric matrix and X and B are N-by-NRHS
45 *> matrices.
46 *>
47 *> Error bounds on the solution and a condition estimate are also
48 *> provided.
49 *> \endverbatim
50 *
51 *> \par Description:
52 * =================
53 *>
54 *> \verbatim
55 *>
56 *> The following steps are performed:
57 *>
58 *> 1. If FACT = 'N', the diagonal pivoting method is used to factor A.
59 *> The form of the factorization is
60 *> A = U * D * U**T, if UPLO = 'U', or
61 *> A = L * D * L**T, if UPLO = 'L',
62 *> where U (or L) is a product of permutation and unit upper (lower)
63 *> triangular matrices, and D is symmetric and block diagonal with
64 *> 1-by-1 and 2-by-2 diagonal blocks.
65 *>
66 *> 2. If some D(i,i)=0, so that D is exactly singular, then the routine
67 *> returns with INFO = i. Otherwise, the factored form of A is used
68 *> to estimate the condition number of the matrix A. If the
69 *> reciprocal of the condition number is less than machine precision,
70 *> INFO = N+1 is returned as a warning, but the routine still goes on
71 *> to solve for X and compute error bounds as described below.
72 *>
73 *> 3. The system of equations is solved for X using the factored form
74 *> of A.
75 *>
76 *> 4. Iterative refinement is applied to improve the computed solution
77 *> matrix and calculate error bounds and backward error estimates
78 *> for it.
79 *> \endverbatim
80 *
81 * Arguments:
82 * ==========
83 *
84 *> \param[in] FACT
85 *> \verbatim
86 *> FACT is CHARACTER*1
87 *> Specifies whether or not the factored form of A has been
88 *> supplied on entry.
89 *> = 'F': On entry, AF and IPIV contain the factored form of
90 *> A. AF and IPIV will not be modified.
91 *> = 'N': The matrix A will be copied to AF and factored.
92 *> \endverbatim
93 *>
94 *> \param[in] UPLO
95 *> \verbatim
96 *> UPLO is CHARACTER*1
97 *> = 'U': Upper triangle of A is stored;
98 *> = 'L': Lower triangle of A is stored.
99 *> \endverbatim
100 *>
101 *> \param[in] N
102 *> \verbatim
103 *> N is INTEGER
104 *> The number of linear equations, i.e., the order of the
105 *> matrix A. N >= 0.
106 *> \endverbatim
107 *>
108 *> \param[in] NRHS
109 *> \verbatim
110 *> NRHS is INTEGER
111 *> The number of right hand sides, i.e., the number of columns
112 *> of the matrices B and X. NRHS >= 0.
113 *> \endverbatim
114 *>
115 *> \param[in] A
116 *> \verbatim
117 *> A is REAL array, dimension (LDA,N)
118 *> The symmetric matrix A. If UPLO = 'U', the leading N-by-N
119 *> upper triangular part of A contains the upper triangular part
120 *> of the matrix A, and the strictly lower triangular part of A
121 *> is not referenced. If UPLO = 'L', the leading N-by-N lower
122 *> triangular part of A contains the lower triangular part of
123 *> the matrix A, and the strictly upper triangular part of A is
124 *> not referenced.
125 *> \endverbatim
126 *>
127 *> \param[in] LDA
128 *> \verbatim
129 *> LDA is INTEGER
130 *> The leading dimension of the array A. LDA >= max(1,N).
131 *> \endverbatim
132 *>
133 *> \param[in,out] AF
134 *> \verbatim
135 *> AF is REAL array, dimension (LDAF,N)
136 *> If FACT = 'F', then AF is an input argument and on entry
137 *> contains the block diagonal matrix D and the multipliers used
138 *> to obtain the factor U or L from the factorization
139 *> A = U*D*U**T or A = L*D*L**T as computed by SSYTRF.
140 *>
141 *> If FACT = 'N', then AF is an output argument and on exit
142 *> returns the block diagonal matrix D and the multipliers used
143 *> to obtain the factor U or L from the factorization
144 *> A = U*D*U**T or A = L*D*L**T.
145 *> \endverbatim
146 *>
147 *> \param[in] LDAF
148 *> \verbatim
149 *> LDAF is INTEGER
150 *> The leading dimension of the array AF. LDAF >= max(1,N).
151 *> \endverbatim
152 *>
153 *> \param[in,out] IPIV
154 *> \verbatim
155 *> IPIV is INTEGER array, dimension (N)
156 *> If FACT = 'F', then IPIV is an input argument and on entry
157 *> contains details of the interchanges and the block structure
158 *> of D, as determined by SSYTRF.
159 *> If IPIV(k) > 0, then rows and columns k and IPIV(k) were
160 *> interchanged and D(k,k) is a 1-by-1 diagonal block.
161 *> If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
162 *> columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
163 *> is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) =
164 *> IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
165 *> interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
166 *>
167 *> If FACT = 'N', then IPIV is an output argument and on exit
168 *> contains details of the interchanges and the block structure
169 *> of D, as determined by SSYTRF.
170 *> \endverbatim
171 *>
172 *> \param[in] B
173 *> \verbatim
174 *> B is REAL array, dimension (LDB,NRHS)
175 *> The N-by-NRHS right hand side matrix B.
176 *> \endverbatim
177 *>
178 *> \param[in] LDB
179 *> \verbatim
180 *> LDB is INTEGER
181 *> The leading dimension of the array B. LDB >= max(1,N).
182 *> \endverbatim
183 *>
184 *> \param[out] X
185 *> \verbatim
186 *> X is REAL array, dimension (LDX,NRHS)
187 *> If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.
188 *> \endverbatim
189 *>
190 *> \param[in] LDX
191 *> \verbatim
192 *> LDX is INTEGER
193 *> The leading dimension of the array X. LDX >= max(1,N).
194 *> \endverbatim
195 *>
196 *> \param[out] RCOND
197 *> \verbatim
198 *> RCOND is REAL
199 *> The estimate of the reciprocal condition number of the matrix
200 *> A. If RCOND is less than the machine precision (in
201 *> particular, if RCOND = 0), the matrix is singular to working
202 *> precision. This condition is indicated by a return code of
203 *> INFO > 0.
204 *> \endverbatim
205 *>
206 *> \param[out] FERR
207 *> \verbatim
208 *> FERR is REAL array, dimension (NRHS)
209 *> The estimated forward error bound for each solution vector
210 *> X(j) (the j-th column of the solution matrix X).
211 *> If XTRUE is the true solution corresponding to X(j), FERR(j)
212 *> is an estimated upper bound for the magnitude of the largest
213 *> element in (X(j) - XTRUE) divided by the magnitude of the
214 *> largest element in X(j). The estimate is as reliable as
215 *> the estimate for RCOND, and is almost always a slight
216 *> overestimate of the true error.
217 *> \endverbatim
218 *>
219 *> \param[out] BERR
220 *> \verbatim
221 *> BERR is REAL array, dimension (NRHS)
222 *> The componentwise relative backward error of each solution
223 *> vector X(j) (i.e., the smallest relative change in
224 *> any element of A or B that makes X(j) an exact solution).
225 *> \endverbatim
226 *>
227 *> \param[out] WORK
228 *> \verbatim
229 *> WORK is REAL array, dimension (MAX(1,LWORK))
230 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
231 *> \endverbatim
232 *>
233 *> \param[in] LWORK
234 *> \verbatim
235 *> LWORK is INTEGER
236 *> The length of WORK. LWORK >= max(1,3*N), and for best
237 *> performance, when FACT = 'N', LWORK >= max(1,3*N,N*NB), where
238 *> NB is the optimal blocksize for SSYTRF.
239 *>
240 *> If LWORK = -1, then a workspace query is assumed; the routine
241 *> only calculates the optimal size of the WORK array, returns
242 *> this value as the first entry of the WORK array, and no error
243 *> message related to LWORK is issued by XERBLA.
244 *> \endverbatim
245 *>
246 *> \param[out] IWORK
247 *> \verbatim
248 *> IWORK is INTEGER array, dimension (N)
249 *> \endverbatim
250 *>
251 *> \param[out] INFO
252 *> \verbatim
253 *> INFO is INTEGER
254 *> = 0: successful exit
255 *> < 0: if INFO = -i, the i-th argument had an illegal value
256 *> > 0: if INFO = i, and i is
257 *> <= N: D(i,i) is exactly zero. The factorization
258 *> has been completed but the factor D is exactly
259 *> singular, so the solution and error bounds could
260 *> not be computed. RCOND = 0 is returned.
261 *> = N+1: D is nonsingular, but RCOND is less than machine
262 *> precision, meaning that the matrix is singular
263 *> to working precision. Nevertheless, the
264 *> solution and error bounds are computed because
265 *> there are a number of situations where the
266 *> computed solution can be more accurate than the
267 *> value of RCOND would suggest.
268 *> \endverbatim
269 *
270 * Authors:
271 * ========
272 *
273 *> \author Univ. of Tennessee
274 *> \author Univ. of California Berkeley
275 *> \author Univ. of Colorado Denver
276 *> \author NAG Ltd.
277 *
278 *> \date April 2012
279 *
280 *> \ingroup realSYsolve
281 *
282 * =====================================================================
283  SUBROUTINE ssysvx( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B,
284  $ ldb, x, ldx, rcond, ferr, berr, work, lwork,
285  $ iwork, info )
286 *
287 * -- LAPACK driver routine (version 3.4.1) --
288 * -- LAPACK is a software package provided by Univ. of Tennessee, --
289 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
290 * April 2012
291 *
292 * .. Scalar Arguments ..
293  CHARACTER fact, uplo
294  INTEGER info, lda, ldaf, ldb, ldx, lwork, n, nrhs
295  REAL rcond
296 * ..
297 * .. Array Arguments ..
298  INTEGER ipiv( * ), iwork( * )
299  REAL a( lda, * ), af( ldaf, * ), b( ldb, * ),
300  $ berr( * ), ferr( * ), work( * ), x( ldx, * )
301 * ..
302 *
303 * =====================================================================
304 *
305 * .. Parameters ..
306  REAL zero
307  parameter( zero = 0.0e+0 )
308 * ..
309 * .. Local Scalars ..
310  LOGICAL lquery, nofact
311  INTEGER lwkopt, nb
312  REAL anorm
313 * ..
314 * .. External Functions ..
315  LOGICAL lsame
316  INTEGER ilaenv
317  REAL slamch, slansy
318  EXTERNAL ilaenv, lsame, slamch, slansy
319 * ..
320 * .. External Subroutines ..
321  EXTERNAL slacpy, ssycon, ssyrfs, ssytrf, ssytrs, xerbla
322 * ..
323 * .. Intrinsic Functions ..
324  INTRINSIC max
325 * ..
326 * .. Executable Statements ..
327 *
328 * Test the input parameters.
329 *
330  info = 0
331  nofact = lsame( fact, 'N' )
332  lquery = ( lwork.EQ.-1 )
333  IF( .NOT.nofact .AND. .NOT.lsame( fact, 'F' ) ) THEN
334  info = -1
335  ELSE IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lsame( uplo, 'L' ) )
336  $ THEN
337  info = -2
338  ELSE IF( n.LT.0 ) THEN
339  info = -3
340  ELSE IF( nrhs.LT.0 ) THEN
341  info = -4
342  ELSE IF( lda.LT.max( 1, n ) ) THEN
343  info = -6
344  ELSE IF( ldaf.LT.max( 1, n ) ) THEN
345  info = -8
346  ELSE IF( ldb.LT.max( 1, n ) ) THEN
347  info = -11
348  ELSE IF( ldx.LT.max( 1, n ) ) THEN
349  info = -13
350  ELSE IF( lwork.LT.max( 1, 3*n ) .AND. .NOT.lquery ) THEN
351  info = -18
352  END IF
353 *
354  IF( info.EQ.0 ) THEN
355  lwkopt = max( 1, 3*n )
356  IF( nofact ) THEN
357  nb = ilaenv( 1, 'SSYTRF', uplo, n, -1, -1, -1 )
358  lwkopt = max( lwkopt, n*nb )
359  END IF
360  work( 1 ) = lwkopt
361  END IF
362 *
363  IF( info.NE.0 ) THEN
364  CALL xerbla( 'SSYSVX', -info )
365  return
366  ELSE IF( lquery ) THEN
367  return
368  END IF
369 *
370  IF( nofact ) THEN
371 *
372 * Compute the factorization A = U*D*U**T or A = L*D*L**T.
373 *
374  CALL slacpy( uplo, n, n, a, lda, af, ldaf )
375  CALL ssytrf( uplo, n, af, ldaf, ipiv, work, lwork, info )
376 *
377 * Return if INFO is non-zero.
378 *
379  IF( info.GT.0 )THEN
380  rcond = zero
381  return
382  END IF
383  END IF
384 *
385 * Compute the norm of the matrix A.
386 *
387  anorm = slansy( 'I', uplo, n, a, lda, work )
388 *
389 * Compute the reciprocal of the condition number of A.
390 *
391  CALL ssycon( uplo, n, af, ldaf, ipiv, anorm, rcond, work, iwork,
392  $ info )
393 *
394 * Compute the solution vectors X.
395 *
396  CALL slacpy( 'Full', n, nrhs, b, ldb, x, ldx )
397  CALL ssytrs( uplo, n, nrhs, af, ldaf, ipiv, x, ldx, info )
398 *
399 * Use iterative refinement to improve the computed solutions and
400 * compute error bounds and backward error estimates for them.
401 *
402  CALL ssyrfs( uplo, n, nrhs, a, lda, af, ldaf, ipiv, b, ldb, x,
403  $ ldx, ferr, berr, work, iwork, info )
404 *
405 * Set INFO = N+1 if the matrix is singular to working precision.
406 *
407  IF( rcond.LT.slamch( 'Epsilon' ) )
408  $ info = n + 1
409 *
410  work( 1 ) = lwkopt
411 *
412  return
413 *
414 * End of SSYSVX
415 *
416  END