![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief <b> ZHESVX computes the solution to system of linear equations A * X = B for HE matrices</b> 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download ZHESVX + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhesvx.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhesvx.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhesvx.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZHESVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B, 00022 * LDB, X, LDX, RCOND, FERR, BERR, WORK, LWORK, 00023 * RWORK, INFO ) 00024 * 00025 * .. Scalar Arguments .. 00026 * CHARACTER FACT, UPLO 00027 * INTEGER INFO, LDA, LDAF, LDB, LDX, LWORK, N, NRHS 00028 * DOUBLE PRECISION RCOND 00029 * .. 00030 * .. Array Arguments .. 00031 * INTEGER IPIV( * ) 00032 * DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * ) 00033 * COMPLEX*16 A( LDA, * ), AF( LDAF, * ), B( LDB, * ), 00034 * $ WORK( * ), X( LDX, * ) 00035 * .. 00036 * 00037 * 00038 *> \par Purpose: 00039 * ============= 00040 *> 00041 *> \verbatim 00042 *> 00043 *> ZHESVX uses the diagonal pivoting factorization to compute the 00044 *> solution to a complex system of linear equations A * X = B, 00045 *> where A is an N-by-N Hermitian matrix and X and B are N-by-NRHS 00046 *> matrices. 00047 *> 00048 *> Error bounds on the solution and a condition estimate are also 00049 *> provided. 00050 *> \endverbatim 00051 * 00052 *> \par Description: 00053 * ================= 00054 *> 00055 *> \verbatim 00056 *> 00057 *> The following steps are performed: 00058 *> 00059 *> 1. If FACT = 'N', the diagonal pivoting method is used to factor A. 00060 *> The form of the factorization is 00061 *> A = U * D * U**H, if UPLO = 'U', or 00062 *> A = L * D * L**H, if UPLO = 'L', 00063 *> where U (or L) is a product of permutation and unit upper (lower) 00064 *> triangular matrices, and D is Hermitian and block diagonal with 00065 *> 1-by-1 and 2-by-2 diagonal blocks. 00066 *> 00067 *> 2. If some D(i,i)=0, so that D is exactly singular, then the routine 00068 *> returns with INFO = i. Otherwise, the factored form of A is used 00069 *> to estimate the condition number of the matrix A. If the 00070 *> reciprocal of the condition number is less than machine precision, 00071 *> INFO = N+1 is returned as a warning, but the routine still goes on 00072 *> to solve for X and compute error bounds as described below. 00073 *> 00074 *> 3. The system of equations is solved for X using the factored form 00075 *> of A. 00076 *> 00077 *> 4. Iterative refinement is applied to improve the computed solution 00078 *> matrix and calculate error bounds and backward error estimates 00079 *> for it. 00080 *> \endverbatim 00081 * 00082 * Arguments: 00083 * ========== 00084 * 00085 *> \param[in] FACT 00086 *> \verbatim 00087 *> FACT is CHARACTER*1 00088 *> Specifies whether or not the factored form of A has been 00089 *> supplied on entry. 00090 *> = 'F': On entry, AF and IPIV contain the factored form 00091 *> of A. A, AF and IPIV will not be modified. 00092 *> = 'N': The matrix A will be copied to AF and factored. 00093 *> \endverbatim 00094 *> 00095 *> \param[in] UPLO 00096 *> \verbatim 00097 *> UPLO is CHARACTER*1 00098 *> = 'U': Upper triangle of A is stored; 00099 *> = 'L': Lower triangle of A is stored. 00100 *> \endverbatim 00101 *> 00102 *> \param[in] N 00103 *> \verbatim 00104 *> N is INTEGER 00105 *> The number of linear equations, i.e., the order of the 00106 *> matrix A. N >= 0. 00107 *> \endverbatim 00108 *> 00109 *> \param[in] NRHS 00110 *> \verbatim 00111 *> NRHS is INTEGER 00112 *> The number of right hand sides, i.e., the number of columns 00113 *> of the matrices B and X. NRHS >= 0. 00114 *> \endverbatim 00115 *> 00116 *> \param[in] A 00117 *> \verbatim 00118 *> A is COMPLEX*16 array, dimension (LDA,N) 00119 *> The Hermitian matrix A. If UPLO = 'U', the leading N-by-N 00120 *> upper triangular part of A contains the upper triangular part 00121 *> of the matrix A, and the strictly lower triangular part of A 00122 *> is not referenced. If UPLO = 'L', the leading N-by-N lower 00123 *> triangular part of A contains the lower triangular part of 00124 *> the matrix A, and the strictly upper triangular part of A is 00125 *> not referenced. 00126 *> \endverbatim 00127 *> 00128 *> \param[in] LDA 00129 *> \verbatim 00130 *> LDA is INTEGER 00131 *> The leading dimension of the array A. LDA >= max(1,N). 00132 *> \endverbatim 00133 *> 00134 *> \param[in,out] AF 00135 *> \verbatim 00136 *> AF is COMPLEX*16 array, dimension (LDAF,N) 00137 *> If FACT = 'F', then AF is an input argument and on entry 00138 *> contains the block diagonal matrix D and the multipliers used 00139 *> to obtain the factor U or L from the factorization 00140 *> A = U*D*U**H or A = L*D*L**H as computed by ZHETRF. 00141 *> 00142 *> If FACT = 'N', then AF is an output argument and on exit 00143 *> returns the block diagonal matrix D and the multipliers used 00144 *> to obtain the factor U or L from the factorization 00145 *> A = U*D*U**H or A = L*D*L**H. 00146 *> \endverbatim 00147 *> 00148 *> \param[in] LDAF 00149 *> \verbatim 00150 *> LDAF is INTEGER 00151 *> The leading dimension of the array AF. LDAF >= max(1,N). 00152 *> \endverbatim 00153 *> 00154 *> \param[in,out] IPIV 00155 *> \verbatim 00156 *> IPIV is INTEGER array, dimension (N) 00157 *> If FACT = 'F', then IPIV is an input argument and on entry 00158 *> contains details of the interchanges and the block structure 00159 *> of D, as determined by ZHETRF. 00160 *> If IPIV(k) > 0, then rows and columns k and IPIV(k) were 00161 *> interchanged and D(k,k) is a 1-by-1 diagonal block. 00162 *> If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and 00163 *> columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k) 00164 *> is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) = 00165 *> IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were 00166 *> interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block. 00167 *> 00168 *> If FACT = 'N', then IPIV is an output argument and on exit 00169 *> contains details of the interchanges and the block structure 00170 *> of D, as determined by ZHETRF. 00171 *> \endverbatim 00172 *> 00173 *> \param[in] B 00174 *> \verbatim 00175 *> B is COMPLEX*16 array, dimension (LDB,NRHS) 00176 *> The N-by-NRHS right hand side matrix B. 00177 *> \endverbatim 00178 *> 00179 *> \param[in] LDB 00180 *> \verbatim 00181 *> LDB is INTEGER 00182 *> The leading dimension of the array B. LDB >= max(1,N). 00183 *> \endverbatim 00184 *> 00185 *> \param[out] X 00186 *> \verbatim 00187 *> X is COMPLEX*16 array, dimension (LDX,NRHS) 00188 *> If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X. 00189 *> \endverbatim 00190 *> 00191 *> \param[in] LDX 00192 *> \verbatim 00193 *> LDX is INTEGER 00194 *> The leading dimension of the array X. LDX >= max(1,N). 00195 *> \endverbatim 00196 *> 00197 *> \param[out] RCOND 00198 *> \verbatim 00199 *> RCOND is DOUBLE PRECISION 00200 *> The estimate of the reciprocal condition number of the matrix 00201 *> A. If RCOND is less than the machine precision (in 00202 *> particular, if RCOND = 0), the matrix is singular to working 00203 *> precision. This condition is indicated by a return code of 00204 *> INFO > 0. 00205 *> \endverbatim 00206 *> 00207 *> \param[out] FERR 00208 *> \verbatim 00209 *> FERR is DOUBLE PRECISION array, dimension (NRHS) 00210 *> The estimated forward error bound for each solution vector 00211 *> X(j) (the j-th column of the solution matrix X). 00212 *> If XTRUE is the true solution corresponding to X(j), FERR(j) 00213 *> is an estimated upper bound for the magnitude of the largest 00214 *> element in (X(j) - XTRUE) divided by the magnitude of the 00215 *> largest element in X(j). The estimate is as reliable as 00216 *> the estimate for RCOND, and is almost always a slight 00217 *> overestimate of the true error. 00218 *> \endverbatim 00219 *> 00220 *> \param[out] BERR 00221 *> \verbatim 00222 *> BERR is DOUBLE PRECISION array, dimension (NRHS) 00223 *> The componentwise relative backward error of each solution 00224 *> vector X(j) (i.e., the smallest relative change in 00225 *> any element of A or B that makes X(j) an exact solution). 00226 *> \endverbatim 00227 *> 00228 *> \param[out] WORK 00229 *> \verbatim 00230 *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK)) 00231 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00232 *> \endverbatim 00233 *> 00234 *> \param[in] LWORK 00235 *> \verbatim 00236 *> LWORK is INTEGER 00237 *> The length of WORK. LWORK >= max(1,2*N), and for best 00238 *> performance, when FACT = 'N', LWORK >= max(1,2*N,N*NB), where 00239 *> NB is the optimal blocksize for ZHETRF. 00240 *> 00241 *> If LWORK = -1, then a workspace query is assumed; the routine 00242 *> only calculates the optimal size of the WORK array, returns 00243 *> this value as the first entry of the WORK array, and no error 00244 *> message related to LWORK is issued by XERBLA. 00245 *> \endverbatim 00246 *> 00247 *> \param[out] RWORK 00248 *> \verbatim 00249 *> RWORK is DOUBLE PRECISION array, dimension (N) 00250 *> \endverbatim 00251 *> 00252 *> \param[out] INFO 00253 *> \verbatim 00254 *> INFO is INTEGER 00255 *> = 0: successful exit 00256 *> < 0: if INFO = -i, the i-th argument had an illegal value 00257 *> > 0: if INFO = i, and i is 00258 *> <= N: D(i,i) is exactly zero. The factorization 00259 *> has been completed but the factor D is exactly 00260 *> singular, so the solution and error bounds could 00261 *> not be computed. RCOND = 0 is returned. 00262 *> = N+1: D is nonsingular, but RCOND is less than machine 00263 *> precision, meaning that the matrix is singular 00264 *> to working precision. Nevertheless, the 00265 *> solution and error bounds are computed because 00266 *> there are a number of situations where the 00267 *> computed solution can be more accurate than the 00268 *> value of RCOND would suggest. 00269 *> \endverbatim 00270 * 00271 * Authors: 00272 * ======== 00273 * 00274 *> \author Univ. of Tennessee 00275 *> \author Univ. of California Berkeley 00276 *> \author Univ. of Colorado Denver 00277 *> \author NAG Ltd. 00278 * 00279 *> \date April 2012 00280 * 00281 *> \ingroup complex16HEsolve 00282 * 00283 * ===================================================================== 00284 SUBROUTINE ZHESVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B, 00285 $ LDB, X, LDX, RCOND, FERR, BERR, WORK, LWORK, 00286 $ RWORK, INFO ) 00287 * 00288 * -- LAPACK driver routine (version 3.4.1) -- 00289 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00290 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00291 * April 2012 00292 * 00293 * .. Scalar Arguments .. 00294 CHARACTER FACT, UPLO 00295 INTEGER INFO, LDA, LDAF, LDB, LDX, LWORK, N, NRHS 00296 DOUBLE PRECISION RCOND 00297 * .. 00298 * .. Array Arguments .. 00299 INTEGER IPIV( * ) 00300 DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * ) 00301 COMPLEX*16 A( LDA, * ), AF( LDAF, * ), B( LDB, * ), 00302 $ WORK( * ), X( LDX, * ) 00303 * .. 00304 * 00305 * ===================================================================== 00306 * 00307 * .. Parameters .. 00308 DOUBLE PRECISION ZERO 00309 PARAMETER ( ZERO = 0.0D+0 ) 00310 * .. 00311 * .. Local Scalars .. 00312 LOGICAL LQUERY, NOFACT 00313 INTEGER LWKOPT, NB 00314 DOUBLE PRECISION ANORM 00315 * .. 00316 * .. External Functions .. 00317 LOGICAL LSAME 00318 INTEGER ILAENV 00319 DOUBLE PRECISION DLAMCH, ZLANHE 00320 EXTERNAL LSAME, ILAENV, DLAMCH, ZLANHE 00321 * .. 00322 * .. External Subroutines .. 00323 EXTERNAL XERBLA, ZHECON, ZHERFS, ZHETRF, ZHETRS, ZLACPY 00324 * .. 00325 * .. Intrinsic Functions .. 00326 INTRINSIC MAX 00327 * .. 00328 * .. Executable Statements .. 00329 * 00330 * Test the input parameters. 00331 * 00332 INFO = 0 00333 NOFACT = LSAME( FACT, 'N' ) 00334 LQUERY = ( LWORK.EQ.-1 ) 00335 IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN 00336 INFO = -1 00337 ELSE IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) 00338 $ THEN 00339 INFO = -2 00340 ELSE IF( N.LT.0 ) THEN 00341 INFO = -3 00342 ELSE IF( NRHS.LT.0 ) THEN 00343 INFO = -4 00344 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00345 INFO = -6 00346 ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN 00347 INFO = -8 00348 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00349 INFO = -11 00350 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN 00351 INFO = -13 00352 ELSE IF( LWORK.LT.MAX( 1, 2*N ) .AND. .NOT.LQUERY ) THEN 00353 INFO = -18 00354 END IF 00355 * 00356 IF( INFO.EQ.0 ) THEN 00357 LWKOPT = MAX( 1, 2*N ) 00358 IF( NOFACT ) THEN 00359 NB = ILAENV( 1, 'ZHETRF', UPLO, N, -1, -1, -1 ) 00360 LWKOPT = MAX( LWKOPT, N*NB ) 00361 END IF 00362 WORK( 1 ) = LWKOPT 00363 END IF 00364 * 00365 IF( INFO.NE.0 ) THEN 00366 CALL XERBLA( 'ZHESVX', -INFO ) 00367 RETURN 00368 ELSE IF( LQUERY ) THEN 00369 RETURN 00370 END IF 00371 * 00372 IF( NOFACT ) THEN 00373 * 00374 * Compute the factorization A = U*D*U**H or A = L*D*L**H. 00375 * 00376 CALL ZLACPY( UPLO, N, N, A, LDA, AF, LDAF ) 00377 CALL ZHETRF( UPLO, N, AF, LDAF, IPIV, WORK, LWORK, INFO ) 00378 * 00379 * Return if INFO is non-zero. 00380 * 00381 IF( INFO.GT.0 )THEN 00382 RCOND = ZERO 00383 RETURN 00384 END IF 00385 END IF 00386 * 00387 * Compute the norm of the matrix A. 00388 * 00389 ANORM = ZLANHE( 'I', UPLO, N, A, LDA, RWORK ) 00390 * 00391 * Compute the reciprocal of the condition number of A. 00392 * 00393 CALL ZHECON( UPLO, N, AF, LDAF, IPIV, ANORM, RCOND, WORK, INFO ) 00394 * 00395 * Compute the solution vectors X. 00396 * 00397 CALL ZLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 00398 CALL ZHETRS( UPLO, N, NRHS, AF, LDAF, IPIV, X, LDX, INFO ) 00399 * 00400 * Use iterative refinement to improve the computed solutions and 00401 * compute error bounds and backward error estimates for them. 00402 * 00403 CALL ZHERFS( UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB, X, 00404 $ LDX, FERR, BERR, WORK, RWORK, INFO ) 00405 * 00406 * Set INFO = N+1 if the matrix is singular to working precision. 00407 * 00408 IF( RCOND.LT.DLAMCH( 'Epsilon' ) ) 00409 $ INFO = N + 1 00410 * 00411 WORK( 1 ) = LWKOPT 00412 * 00413 RETURN 00414 * 00415 * End of ZHESVX 00416 * 00417 END