1*> \brief <b> CSPSVX computes the solution to system of linear equations A * X = B for OTHER 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 CSPSVX + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cspsvx.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cspsvx.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cspsvx.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE CSPSVX( FACT, UPLO, N, NRHS, AP, AFP, IPIV, B, LDB, X, 22* LDX, RCOND, FERR, BERR, WORK, RWORK, INFO ) 23* 24* .. Scalar Arguments .. 25* CHARACTER FACT, UPLO 26* INTEGER INFO, LDB, LDX, N, NRHS 27* REAL RCOND 28* .. 29* .. Array Arguments .. 30* INTEGER IPIV( * ) 31* REAL BERR( * ), FERR( * ), RWORK( * ) 32* COMPLEX AFP( * ), AP( * ), B( LDB, * ), WORK( * ), 33* $ X( LDX, * ) 34* .. 35* 36* 37*> \par Purpose: 38* ============= 39*> 40*> \verbatim 41*> 42*> CSPSVX uses the diagonal pivoting factorization A = U*D*U**T or 43*> A = L*D*L**T to compute the solution to a complex system of linear 44*> equations A * X = B, where A is an N-by-N symmetric matrix stored 45*> in packed format and X and B are N-by-NRHS matrices. 46*> 47*> Error bounds on the solution and a condition estimate are also 48*> provided. 49*> \endverbatim 50* 51*> \par Description: 52* ================= 53*> 54*> \verbatim 55*> 56*> The following steps are performed: 57*> 58*> 1. If FACT = 'N', the diagonal pivoting method is used to factor A as 59*> A = U * D * U**T, if UPLO = 'U', or 60*> A = L * D * L**T, if UPLO = 'L', 61*> where U (or L) is a product of permutation and unit upper (lower) 62*> triangular matrices and D is symmetric and block diagonal with 63*> 1-by-1 and 2-by-2 diagonal blocks. 64*> 65*> 2. If some D(i,i)=0, so that D is exactly singular, then the routine 66*> returns with INFO = i. Otherwise, the factored form of A is used 67*> to estimate the condition number of the matrix A. If the 68*> reciprocal of the condition number is less than machine precision, 69*> INFO = N+1 is returned as a warning, but the routine still goes on 70*> to solve for X and compute error bounds as described below. 71*> 72*> 3. The system of equations is solved for X using the factored form 73*> of A. 74*> 75*> 4. Iterative refinement is applied to improve the computed solution 76*> matrix and calculate error bounds and backward error estimates 77*> for it. 78*> \endverbatim 79* 80* Arguments: 81* ========== 82* 83*> \param[in] FACT 84*> \verbatim 85*> FACT is CHARACTER*1 86*> Specifies whether or not the factored form of A has been 87*> supplied on entry. 88*> = 'F': On entry, AFP and IPIV contain the factored form 89*> of A. AP, AFP and IPIV will not be modified. 90*> = 'N': The matrix A will be copied to AFP and factored. 91*> \endverbatim 92*> 93*> \param[in] UPLO 94*> \verbatim 95*> UPLO is CHARACTER*1 96*> = 'U': Upper triangle of A is stored; 97*> = 'L': Lower triangle of A is stored. 98*> \endverbatim 99*> 100*> \param[in] N 101*> \verbatim 102*> N is INTEGER 103*> The number of linear equations, i.e., the order of the 104*> matrix A. N >= 0. 105*> \endverbatim 106*> 107*> \param[in] NRHS 108*> \verbatim 109*> NRHS is INTEGER 110*> The number of right hand sides, i.e., the number of columns 111*> of the matrices B and X. NRHS >= 0. 112*> \endverbatim 113*> 114*> \param[in] AP 115*> \verbatim 116*> AP is COMPLEX array, dimension (N*(N+1)/2) 117*> The upper or lower triangle of the symmetric matrix A, packed 118*> columnwise in a linear array. The j-th column of A is stored 119*> in the array AP as follows: 120*> if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; 121*> if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n. 122*> See below for further details. 123*> \endverbatim 124*> 125*> \param[in,out] AFP 126*> \verbatim 127*> AFP is COMPLEX array, dimension (N*(N+1)/2) 128*> If FACT = 'F', then AFP is an input argument and on entry 129*> contains the block diagonal matrix D and the multipliers used 130*> to obtain the factor U or L from the factorization 131*> A = U*D*U**T or A = L*D*L**T as computed by CSPTRF, stored as 132*> a packed triangular matrix in the same storage format as A. 133*> 134*> If FACT = 'N', then AFP is an output argument and on exit 135*> contains the block diagonal matrix D and the multipliers used 136*> to obtain the factor U or L from the factorization 137*> A = U*D*U**T or A = L*D*L**T as computed by CSPTRF, stored as 138*> a packed triangular matrix in the same storage format as A. 139*> \endverbatim 140*> 141*> \param[in,out] IPIV 142*> \verbatim 143*> IPIV is INTEGER array, dimension (N) 144*> If FACT = 'F', then IPIV is an input argument and on entry 145*> contains details of the interchanges and the block structure 146*> of D, as determined by CSPTRF. 147*> If IPIV(k) > 0, then rows and columns k and IPIV(k) were 148*> interchanged and D(k,k) is a 1-by-1 diagonal block. 149*> If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and 150*> columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k) 151*> is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) = 152*> IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were 153*> interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block. 154*> 155*> If FACT = 'N', then IPIV is an output argument and on exit 156*> contains details of the interchanges and the block structure 157*> of D, as determined by CSPTRF. 158*> \endverbatim 159*> 160*> \param[in] B 161*> \verbatim 162*> B is COMPLEX array, dimension (LDB,NRHS) 163*> The N-by-NRHS right hand side matrix B. 164*> \endverbatim 165*> 166*> \param[in] LDB 167*> \verbatim 168*> LDB is INTEGER 169*> The leading dimension of the array B. LDB >= max(1,N). 170*> \endverbatim 171*> 172*> \param[out] X 173*> \verbatim 174*> X is COMPLEX array, dimension (LDX,NRHS) 175*> If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X. 176*> \endverbatim 177*> 178*> \param[in] LDX 179*> \verbatim 180*> LDX is INTEGER 181*> The leading dimension of the array X. LDX >= max(1,N). 182*> \endverbatim 183*> 184*> \param[out] RCOND 185*> \verbatim 186*> RCOND is REAL 187*> The estimate of the reciprocal condition number of the matrix 188*> A. If RCOND is less than the machine precision (in 189*> particular, if RCOND = 0), the matrix is singular to working 190*> precision. This condition is indicated by a return code of 191*> INFO > 0. 192*> \endverbatim 193*> 194*> \param[out] FERR 195*> \verbatim 196*> FERR is REAL array, dimension (NRHS) 197*> The estimated forward error bound for each solution vector 198*> X(j) (the j-th column of the solution matrix X). 199*> If XTRUE is the true solution corresponding to X(j), FERR(j) 200*> is an estimated upper bound for the magnitude of the largest 201*> element in (X(j) - XTRUE) divided by the magnitude of the 202*> largest element in X(j). The estimate is as reliable as 203*> the estimate for RCOND, and is almost always a slight 204*> overestimate of the true error. 205*> \endverbatim 206*> 207*> \param[out] BERR 208*> \verbatim 209*> BERR is REAL array, dimension (NRHS) 210*> The componentwise relative backward error of each solution 211*> vector X(j) (i.e., the smallest relative change in 212*> any element of A or B that makes X(j) an exact solution). 213*> \endverbatim 214*> 215*> \param[out] WORK 216*> \verbatim 217*> WORK is COMPLEX array, dimension (2*N) 218*> \endverbatim 219*> 220*> \param[out] RWORK 221*> \verbatim 222*> RWORK is REAL array, dimension (N) 223*> \endverbatim 224*> 225*> \param[out] INFO 226*> \verbatim 227*> INFO is INTEGER 228*> = 0: successful exit 229*> < 0: if INFO = -i, the i-th argument had an illegal value 230*> > 0: if INFO = i, and i is 231*> <= N: D(i,i) is exactly zero. The factorization 232*> has been completed but the factor D is exactly 233*> singular, so the solution and error bounds could 234*> not be computed. RCOND = 0 is returned. 235*> = N+1: D is nonsingular, but RCOND is less than machine 236*> precision, meaning that the matrix is singular 237*> to working precision. Nevertheless, the 238*> solution and error bounds are computed because 239*> there are a number of situations where the 240*> computed solution can be more accurate than the 241*> value of RCOND would suggest. 242*> \endverbatim 243* 244* Authors: 245* ======== 246* 247*> \author Univ. of Tennessee 248*> \author Univ. of California Berkeley 249*> \author Univ. of Colorado Denver 250*> \author NAG Ltd. 251* 252*> \ingroup complexOTHERsolve 253* 254*> \par Further Details: 255* ===================== 256*> 257*> \verbatim 258*> 259*> The packed storage scheme is illustrated by the following example 260*> when N = 4, UPLO = 'U': 261*> 262*> Two-dimensional storage of the symmetric matrix A: 263*> 264*> a11 a12 a13 a14 265*> a22 a23 a24 266*> a33 a34 (aij = aji) 267*> a44 268*> 269*> Packed storage of the upper triangle of A: 270*> 271*> AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ] 272*> \endverbatim 273*> 274* ===================================================================== 275 SUBROUTINE CSPSVX( FACT, UPLO, N, NRHS, AP, AFP, IPIV, B, LDB, X, 276 $ LDX, RCOND, FERR, BERR, WORK, RWORK, INFO ) 277* 278* -- LAPACK driver routine -- 279* -- LAPACK is a software package provided by Univ. of Tennessee, -- 280* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 281* 282* .. Scalar Arguments .. 283 CHARACTER FACT, UPLO 284 INTEGER INFO, LDB, LDX, N, NRHS 285 REAL RCOND 286* .. 287* .. Array Arguments .. 288 INTEGER IPIV( * ) 289 REAL BERR( * ), FERR( * ), RWORK( * ) 290 COMPLEX AFP( * ), AP( * ), B( LDB, * ), WORK( * ), 291 $ X( LDX, * ) 292* .. 293* 294* ===================================================================== 295* 296* .. Parameters .. 297 REAL ZERO 298 PARAMETER ( ZERO = 0.0E+0 ) 299* .. 300* .. Local Scalars .. 301 LOGICAL NOFACT 302 REAL ANORM 303* .. 304* .. External Functions .. 305 LOGICAL LSAME 306 REAL CLANSP, SLAMCH 307 EXTERNAL LSAME, CLANSP, SLAMCH 308* .. 309* .. External Subroutines .. 310 EXTERNAL CCOPY, CLACPY, CSPCON, CSPRFS, CSPTRF, CSPTRS, 311 $ XERBLA 312* .. 313* .. Intrinsic Functions .. 314 INTRINSIC MAX 315* .. 316* .. Executable Statements .. 317* 318* Test the input parameters. 319* 320 INFO = 0 321 NOFACT = LSAME( FACT, 'N' ) 322 IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN 323 INFO = -1 324 ELSE IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) 325 $ THEN 326 INFO = -2 327 ELSE IF( N.LT.0 ) THEN 328 INFO = -3 329 ELSE IF( NRHS.LT.0 ) THEN 330 INFO = -4 331 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 332 INFO = -9 333 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN 334 INFO = -11 335 END IF 336 IF( INFO.NE.0 ) THEN 337 CALL XERBLA( 'CSPSVX', -INFO ) 338 RETURN 339 END IF 340* 341 IF( NOFACT ) THEN 342* 343* Compute the factorization A = U*D*U**T or A = L*D*L**T. 344* 345 CALL CCOPY( N*( N+1 ) / 2, AP, 1, AFP, 1 ) 346 CALL CSPTRF( UPLO, N, AFP, IPIV, INFO ) 347* 348* Return if INFO is non-zero. 349* 350 IF( INFO.GT.0 )THEN 351 RCOND = ZERO 352 RETURN 353 END IF 354 END IF 355* 356* Compute the norm of the matrix A. 357* 358 ANORM = CLANSP( 'I', UPLO, N, AP, RWORK ) 359* 360* Compute the reciprocal of the condition number of A. 361* 362 CALL CSPCON( UPLO, N, AFP, IPIV, ANORM, RCOND, WORK, INFO ) 363* 364* Compute the solution vectors X. 365* 366 CALL CLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 367 CALL CSPTRS( UPLO, N, NRHS, AFP, IPIV, X, LDX, INFO ) 368* 369* Use iterative refinement to improve the computed solutions and 370* compute error bounds and backward error estimates for them. 371* 372 CALL CSPRFS( UPLO, N, NRHS, AP, AFP, IPIV, B, LDB, X, LDX, FERR, 373 $ BERR, WORK, RWORK, INFO ) 374* 375* Set INFO = N+1 if the matrix is singular to working precision. 376* 377 IF( RCOND.LT.SLAMCH( 'Epsilon' ) ) 378 $ INFO = N + 1 379* 380 RETURN 381* 382* End of CSPSVX 383* 384 END 385