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