LAPACK  3.4.2 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
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