LAPACK  3.6.0
LAPACK: Linear Algebra PACKage
dposvx.f
Go to the documentation of this file.
1 *> \brief <b> DPOSVX computes the solution to system of linear equations A * X = B for PO 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 DPOSVX + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dposvx.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dposvx.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dposvx.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DPOSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, EQUED,
22 * S, B, LDB, X, LDX, RCOND, FERR, BERR, WORK,
23 * IWORK, INFO )
24 *
25 * .. Scalar Arguments ..
26 * CHARACTER EQUED, FACT, UPLO
27 * INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS
28 * DOUBLE PRECISION RCOND
29 * ..
30 * .. Array Arguments ..
31 * INTEGER IWORK( * )
32 * DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
33 * $ BERR( * ), FERR( * ), S( * ), WORK( * ),
34 * $ X( LDX, * )
35 * ..
36 *
37 *
38 *> \par Purpose:
39 * =============
40 *>
41 *> \verbatim
42 *>
43 *> DPOSVX uses the Cholesky factorization A = U**T*U or A = L*L**T to
44 *> compute the solution to a real system of linear equations
45 *> A * X = B,
46 *> where A is an N-by-N symmetric positive definite matrix and X and B
47 *> are N-by-NRHS matrices.
48 *>
49 *> Error bounds on the solution and a condition estimate are also
50 *> provided.
51 *> \endverbatim
52 *
53 *> \par Description:
54 * =================
55 *>
56 *> \verbatim
57 *>
58 *> The following steps are performed:
59 *>
60 *> 1. If FACT = 'E', real scaling factors are computed to equilibrate
61 *> the system:
62 *> diag(S) * A * diag(S) * inv(diag(S)) * X = diag(S) * B
63 *> Whether or not the system will be equilibrated depends on the
64 *> scaling of the matrix A, but if equilibration is used, A is
65 *> overwritten by diag(S)*A*diag(S) and B by diag(S)*B.
66 *>
67 *> 2. If FACT = 'N' or 'E', the Cholesky decomposition is used to
68 *> factor the matrix A (after equilibration if FACT = 'E') as
69 *> A = U**T* U, if UPLO = 'U', or
70 *> A = L * L**T, if UPLO = 'L',
71 *> where U is an upper triangular matrix and L is a lower triangular
72 *> matrix.
73 *>
74 *> 3. If the leading i-by-i principal minor is not positive definite,
75 *> then the routine returns with INFO = i. Otherwise, the factored
76 *> form of A is used to estimate the condition number of the matrix
77 *> A. If the reciprocal of the condition number is less than machine
78 *> precision, INFO = N+1 is returned as a warning, but the routine
79 *> still goes on to solve for X and compute error bounds as
80 *> described below.
81 *>
82 *> 4. The system of equations is solved for X using the factored form
83 *> of A.
84 *>
85 *> 5. Iterative refinement is applied to improve the computed solution
86 *> matrix and calculate error bounds and backward error estimates
87 *> for it.
88 *>
89 *> 6. If equilibration was used, the matrix X is premultiplied by
90 *> diag(S) so that it solves the original system before
91 *> equilibration.
92 *> \endverbatim
93 *
94 * Arguments:
95 * ==========
96 *
97 *> \param[in] FACT
98 *> \verbatim
99 *> FACT is CHARACTER*1
100 *> Specifies whether or not the factored form of the matrix A is
101 *> supplied on entry, and if not, whether the matrix A should be
102 *> equilibrated before it is factored.
103 *> = 'F': On entry, AF contains the factored form of A.
104 *> If EQUED = 'Y', the matrix A has been equilibrated
105 *> with scaling factors given by S. A and AF will not
106 *> be modified.
107 *> = 'N': The matrix A will be copied to AF and factored.
108 *> = 'E': The matrix A will be equilibrated if necessary, then
109 *> copied to AF and factored.
110 *> \endverbatim
111 *>
112 *> \param[in] UPLO
113 *> \verbatim
114 *> UPLO is CHARACTER*1
115 *> = 'U': Upper triangle of A is stored;
116 *> = 'L': Lower triangle of A is stored.
117 *> \endverbatim
118 *>
119 *> \param[in] N
120 *> \verbatim
121 *> N is INTEGER
122 *> The number of linear equations, i.e., the order of the
123 *> matrix A. N >= 0.
124 *> \endverbatim
125 *>
126 *> \param[in] NRHS
127 *> \verbatim
128 *> NRHS is INTEGER
129 *> The number of right hand sides, i.e., the number of columns
130 *> of the matrices B and X. NRHS >= 0.
131 *> \endverbatim
132 *>
133 *> \param[in,out] A
134 *> \verbatim
135 *> A is DOUBLE PRECISION array, dimension (LDA,N)
136 *> On entry, the symmetric matrix A, except if FACT = 'F' and
137 *> EQUED = 'Y', then A must contain the equilibrated matrix
138 *> diag(S)*A*diag(S). If UPLO = 'U', the leading
139 *> N-by-N upper triangular part of A contains the upper
140 *> triangular part of the matrix A, and the strictly lower
141 *> triangular part of A is not referenced. If UPLO = 'L', the
142 *> leading N-by-N lower triangular part of A contains the lower
143 *> triangular part of the matrix A, and the strictly upper
144 *> triangular part of A is not referenced. A is not modified if
145 *> FACT = 'F' or 'N', or if FACT = 'E' and EQUED = 'N' on exit.
146 *>
147 *> On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by
148 *> diag(S)*A*diag(S).
149 *> \endverbatim
150 *>
151 *> \param[in] LDA
152 *> \verbatim
153 *> LDA is INTEGER
154 *> The leading dimension of the array A. LDA >= max(1,N).
155 *> \endverbatim
156 *>
157 *> \param[in,out] AF
158 *> \verbatim
159 *> AF is DOUBLE PRECISION array, dimension (LDAF,N)
160 *> If FACT = 'F', then AF is an input argument and on entry
161 *> contains the triangular factor U or L from the Cholesky
162 *> factorization A = U**T*U or A = L*L**T, in the same storage
163 *> format as A. If EQUED .ne. 'N', then AF is the factored form
164 *> of the equilibrated matrix diag(S)*A*diag(S).
165 *>
166 *> If FACT = 'N', then AF is an output argument and on exit
167 *> returns the triangular factor U or L from the Cholesky
168 *> factorization A = U**T*U or A = L*L**T of the original
169 *> matrix A.
170 *>
171 *> If FACT = 'E', then AF is an output argument and on exit
172 *> returns the triangular factor U or L from the Cholesky
173 *> factorization A = U**T*U or A = L*L**T of the equilibrated
174 *> matrix A (see the description of A for the form of the
175 *> equilibrated matrix).
176 *> \endverbatim
177 *>
178 *> \param[in] LDAF
179 *> \verbatim
180 *> LDAF is INTEGER
181 *> The leading dimension of the array AF. LDAF >= max(1,N).
182 *> \endverbatim
183 *>
184 *> \param[in,out] EQUED
185 *> \verbatim
186 *> EQUED is CHARACTER*1
187 *> Specifies the form of equilibration that was done.
188 *> = 'N': No equilibration (always true if FACT = 'N').
189 *> = 'Y': Equilibration was done, i.e., A has been replaced by
190 *> diag(S) * A * diag(S).
191 *> EQUED is an input argument if FACT = 'F'; otherwise, it is an
192 *> output argument.
193 *> \endverbatim
194 *>
195 *> \param[in,out] S
196 *> \verbatim
197 *> S is DOUBLE PRECISION array, dimension (N)
198 *> The scale factors for A; not accessed if EQUED = 'N'. S is
199 *> an input argument if FACT = 'F'; otherwise, S is an output
200 *> argument. If FACT = 'F' and EQUED = 'Y', each element of S
201 *> must be positive.
202 *> \endverbatim
203 *>
204 *> \param[in,out] B
205 *> \verbatim
206 *> B is DOUBLE PRECISION array, dimension (LDB,NRHS)
207 *> On entry, the N-by-NRHS right hand side matrix B.
208 *> On exit, if EQUED = 'N', B is not modified; if EQUED = 'Y',
209 *> B is overwritten by diag(S) * B.
210 *> \endverbatim
211 *>
212 *> \param[in] LDB
213 *> \verbatim
214 *> LDB is INTEGER
215 *> The leading dimension of the array B. LDB >= max(1,N).
216 *> \endverbatim
217 *>
218 *> \param[out] X
219 *> \verbatim
220 *> X is DOUBLE PRECISION array, dimension (LDX,NRHS)
221 *> If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to
222 *> the original system of equations. Note that if EQUED = 'Y',
223 *> A and B are modified on exit, and the solution to the
224 *> equilibrated system is inv(diag(S))*X.
225 *> \endverbatim
226 *>
227 *> \param[in] LDX
228 *> \verbatim
229 *> LDX is INTEGER
230 *> The leading dimension of the array X. LDX >= max(1,N).
231 *> \endverbatim
232 *>
233 *> \param[out] RCOND
234 *> \verbatim
235 *> RCOND is DOUBLE PRECISION
236 *> The estimate of the reciprocal condition number of the matrix
237 *> A after equilibration (if done). If RCOND is less than the
238 *> machine precision (in particular, if RCOND = 0), the matrix
239 *> is singular to working precision. This condition is
240 *> indicated by a return code of INFO > 0.
241 *> \endverbatim
242 *>
243 *> \param[out] FERR
244 *> \verbatim
245 *> FERR is DOUBLE PRECISION array, dimension (NRHS)
246 *> The estimated forward error bound for each solution vector
247 *> X(j) (the j-th column of the solution matrix X).
248 *> If XTRUE is the true solution corresponding to X(j), FERR(j)
249 *> is an estimated upper bound for the magnitude of the largest
250 *> element in (X(j) - XTRUE) divided by the magnitude of the
251 *> largest element in X(j). The estimate is as reliable as
252 *> the estimate for RCOND, and is almost always a slight
253 *> overestimate of the true error.
254 *> \endverbatim
255 *>
256 *> \param[out] BERR
257 *> \verbatim
258 *> BERR is DOUBLE PRECISION array, dimension (NRHS)
259 *> The componentwise relative backward error of each solution
260 *> vector X(j) (i.e., the smallest relative change in
261 *> any element of A or B that makes X(j) an exact solution).
262 *> \endverbatim
263 *>
264 *> \param[out] WORK
265 *> \verbatim
266 *> WORK is DOUBLE PRECISION array, dimension (3*N)
267 *> \endverbatim
268 *>
269 *> \param[out] IWORK
270 *> \verbatim
271 *> IWORK is INTEGER array, dimension (N)
272 *> \endverbatim
273 *>
274 *> \param[out] INFO
275 *> \verbatim
276 *> INFO is INTEGER
277 *> = 0: successful exit
278 *> < 0: if INFO = -i, the i-th argument had an illegal value
279 *> > 0: if INFO = i, and i is
280 *> <= N: the leading minor of order i of A is
281 *> not positive definite, so the factorization
282 *> could not be completed, and the solution has not
283 *> been computed. RCOND = 0 is returned.
284 *> = N+1: U is nonsingular, but RCOND is less than machine
285 *> precision, meaning that the matrix is singular
286 *> to working precision. Nevertheless, the
287 *> solution and error bounds are computed because
288 *> there are a number of situations where the
289 *> computed solution can be more accurate than the
290 *> value of RCOND would suggest.
291 *> \endverbatim
292 *
293 * Authors:
294 * ========
295 *
296 *> \author Univ. of Tennessee
297 *> \author Univ. of California Berkeley
298 *> \author Univ. of Colorado Denver
299 *> \author NAG Ltd.
300 *
301 *> \date April 2012
302 *
303 *> \ingroup doublePOsolve
304 *
305 * =====================================================================
306  SUBROUTINE dposvx( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, EQUED,
307  $ s, b, ldb, x, ldx, rcond, ferr, berr, work,
308  $ iwork, info )
309 *
310 * -- LAPACK driver routine (version 3.4.1) --
311 * -- LAPACK is a software package provided by Univ. of Tennessee, --
312 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
313 * April 2012
314 *
315 * .. Scalar Arguments ..
316  CHARACTER EQUED, FACT, UPLO
317  INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS
318  DOUBLE PRECISION RCOND
319 * ..
320 * .. Array Arguments ..
321  INTEGER IWORK( * )
322  DOUBLE PRECISION A( lda, * ), AF( ldaf, * ), B( ldb, * ),
323  $ berr( * ), ferr( * ), s( * ), work( * ),
324  $ x( ldx, * )
325 * ..
326 *
327 * =====================================================================
328 *
329 * .. Parameters ..
330  DOUBLE PRECISION ZERO, ONE
331  parameter( zero = 0.0d+0, one = 1.0d+0 )
332 * ..
333 * .. Local Scalars ..
334  LOGICAL EQUIL, NOFACT, RCEQU
335  INTEGER I, INFEQU, J
336  DOUBLE PRECISION AMAX, ANORM, BIGNUM, SCOND, SMAX, SMIN, SMLNUM
337 * ..
338 * .. External Functions ..
339  LOGICAL LSAME
340  DOUBLE PRECISION DLAMCH, DLANSY
341  EXTERNAL lsame, dlamch, dlansy
342 * ..
343 * .. External Subroutines ..
344  EXTERNAL dlacpy, dlaqsy, dpocon, dpoequ, dporfs, dpotrf,
345  $ dpotrs, xerbla
346 * ..
347 * .. Intrinsic Functions ..
348  INTRINSIC max, min
349 * ..
350 * .. Executable Statements ..
351 *
352  info = 0
353  nofact = lsame( fact, 'N' )
354  equil = lsame( fact, 'E' )
355  IF( nofact .OR. equil ) THEN
356  equed = 'N'
357  rcequ = .false.
358  ELSE
359  rcequ = lsame( equed, 'Y' )
360  smlnum = dlamch( 'Safe minimum' )
361  bignum = one / smlnum
362  END IF
363 *
364 * Test the input parameters.
365 *
366  IF( .NOT.nofact .AND. .NOT.equil .AND. .NOT.lsame( fact, 'F' ) )
367  $ THEN
368  info = -1
369  ELSE IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lsame( uplo, 'L' ) )
370  $ THEN
371  info = -2
372  ELSE IF( n.LT.0 ) THEN
373  info = -3
374  ELSE IF( nrhs.LT.0 ) THEN
375  info = -4
376  ELSE IF( lda.LT.max( 1, n ) ) THEN
377  info = -6
378  ELSE IF( ldaf.LT.max( 1, n ) ) THEN
379  info = -8
380  ELSE IF( lsame( fact, 'F' ) .AND. .NOT.
381  $ ( rcequ .OR. lsame( equed, 'N' ) ) ) THEN
382  info = -9
383  ELSE
384  IF( rcequ ) THEN
385  smin = bignum
386  smax = zero
387  DO 10 j = 1, n
388  smin = min( smin, s( j ) )
389  smax = max( smax, s( j ) )
390  10 CONTINUE
391  IF( smin.LE.zero ) THEN
392  info = -10
393  ELSE IF( n.GT.0 ) THEN
394  scond = max( smin, smlnum ) / min( smax, bignum )
395  ELSE
396  scond = one
397  END IF
398  END IF
399  IF( info.EQ.0 ) THEN
400  IF( ldb.LT.max( 1, n ) ) THEN
401  info = -12
402  ELSE IF( ldx.LT.max( 1, n ) ) THEN
403  info = -14
404  END IF
405  END IF
406  END IF
407 *
408  IF( info.NE.0 ) THEN
409  CALL xerbla( 'DPOSVX', -info )
410  RETURN
411  END IF
412 *
413  IF( equil ) THEN
414 *
415 * Compute row and column scalings to equilibrate the matrix A.
416 *
417  CALL dpoequ( n, a, lda, s, scond, amax, infequ )
418  IF( infequ.EQ.0 ) THEN
419 *
420 * Equilibrate the matrix.
421 *
422  CALL dlaqsy( uplo, n, a, lda, s, scond, amax, equed )
423  rcequ = lsame( equed, 'Y' )
424  END IF
425  END IF
426 *
427 * Scale the right hand side.
428 *
429  IF( rcequ ) THEN
430  DO 30 j = 1, nrhs
431  DO 20 i = 1, n
432  b( i, j ) = s( i )*b( i, j )
433  20 CONTINUE
434  30 CONTINUE
435  END IF
436 *
437  IF( nofact .OR. equil ) THEN
438 *
439 * Compute the Cholesky factorization A = U**T *U or A = L*L**T.
440 *
441  CALL dlacpy( uplo, n, n, a, lda, af, ldaf )
442  CALL dpotrf( uplo, n, af, ldaf, info )
443 *
444 * Return if INFO is non-zero.
445 *
446  IF( info.GT.0 )THEN
447  rcond = zero
448  RETURN
449  END IF
450  END IF
451 *
452 * Compute the norm of the matrix A.
453 *
454  anorm = dlansy( '1', uplo, n, a, lda, work )
455 *
456 * Compute the reciprocal of the condition number of A.
457 *
458  CALL dpocon( uplo, n, af, ldaf, anorm, rcond, work, iwork, info )
459 *
460 * Compute the solution matrix X.
461 *
462  CALL dlacpy( 'Full', n, nrhs, b, ldb, x, ldx )
463  CALL dpotrs( uplo, n, nrhs, af, ldaf, x, ldx, info )
464 *
465 * Use iterative refinement to improve the computed solution and
466 * compute error bounds and backward error estimates for it.
467 *
468  CALL dporfs( uplo, n, nrhs, a, lda, af, ldaf, b, ldb, x, ldx,
469  $ ferr, berr, work, iwork, info )
470 *
471 * Transform the solution matrix X to a solution of the original
472 * system.
473 *
474  IF( rcequ ) THEN
475  DO 50 j = 1, nrhs
476  DO 40 i = 1, n
477  x( i, j ) = s( i )*x( i, j )
478  40 CONTINUE
479  50 CONTINUE
480  DO 60 j = 1, nrhs
481  ferr( j ) = ferr( j ) / scond
482  60 CONTINUE
483  END IF
484 *
485 * Set INFO = N+1 if the matrix is singular to working precision.
486 *
487  IF( rcond.LT.dlamch( 'Epsilon' ) )
488  $ info = n + 1
489 *
490  RETURN
491 *
492 * End of DPOSVX
493 *
494  END
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:105
subroutine dlaqsy(UPLO, N, A, LDA, S, SCOND, AMAX, EQUED)
DLAQSY scales a symmetric/Hermitian matrix, using scaling factors computed by spoequ.
Definition: dlaqsy.f:135
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dpoequ(N, A, LDA, S, SCOND, AMAX, INFO)
DPOEQU
Definition: dpoequ.f:114
subroutine dpotrf(UPLO, N, A, LDA, INFO)
DPOTRF
Definition: dpotrf.f:109
subroutine dposvx(FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, EQUED, S, B, LDB, X, LDX, RCOND, FERR, BERR, WORK, IWORK, INFO)
DPOSVX computes the solution to system of linear equations A * X = B for PO matrices ...
Definition: dposvx.f:309
subroutine dpotrs(UPLO, N, NRHS, A, LDA, B, LDB, INFO)
DPOTRS
Definition: dpotrs.f:112
subroutine dpocon(UPLO, N, A, LDA, ANORM, RCOND, WORK, IWORK, INFO)
DPOCON
Definition: dpocon.f:123
subroutine dporfs(UPLO, N, NRHS, A, LDA, AF, LDAF, B, LDB, X, LDX, FERR, BERR, WORK, IWORK, INFO)
DPORFS
Definition: dporfs.f:185