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