LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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