1: SUBROUTINE ZSPSVX( FACT, UPLO, N, NRHS, AP, AFP, IPIV, B, LDB, X,
2: $ LDX, RCOND, FERR, BERR, WORK, RWORK, INFO )
3: *
4: * -- LAPACK driver routine (version 3.2) --
5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7: * November 2006
8: *
9: * .. Scalar Arguments ..
10: CHARACTER FACT, UPLO
11: INTEGER INFO, LDB, LDX, N, NRHS
12: DOUBLE PRECISION RCOND
13: * ..
14: * .. Array Arguments ..
15: INTEGER IPIV( * )
16: DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * )
17: COMPLEX*16 AFP( * ), AP( * ), B( LDB, * ), WORK( * ),
18: $ X( LDX, * )
19: * ..
20: *
21: * Purpose
22: * =======
23: *
24: * ZSPSVX uses the diagonal pivoting factorization A = U*D*U**T or
25: * A = L*D*L**T to compute the solution to a complex system of linear
26: * equations A * X = B, where A is an N-by-N symmetric matrix stored
27: * in packed format and X and B are N-by-NRHS matrices.
28: *
29: * Error bounds on the solution and a condition estimate are also
30: * provided.
31: *
32: * Description
33: * ===========
34: *
35: * The following steps are performed:
36: *
37: * 1. If FACT = 'N', the diagonal pivoting method is used to factor A as
38: * A = U * D * U**T, if UPLO = 'U', or
39: * A = L * D * L**T, if UPLO = 'L',
40: * where U (or L) is a product of permutation and unit upper (lower)
41: * triangular matrices and D is symmetric and block diagonal with
42: * 1-by-1 and 2-by-2 diagonal blocks.
43: *
44: * 2. If some D(i,i)=0, so that D is exactly singular, then the routine
45: * returns with INFO = i. Otherwise, the factored form of A is used
46: * to estimate the condition number of the matrix A. If the
47: * reciprocal of the condition number is less than machine precision,
48: * INFO = N+1 is returned as a warning, but the routine still goes on
49: * to solve for X and compute error bounds as described below.
50: *
51: * 3. The system of equations is solved for X using the factored form
52: * of A.
53: *
54: * 4. Iterative refinement is applied to improve the computed solution
55: * matrix and calculate error bounds and backward error estimates
56: * for it.
57: *
58: * Arguments
59: * =========
60: *
61: * FACT (input) CHARACTER*1
62: * Specifies whether or not the factored form of A has been
63: * supplied on entry.
64: * = 'F': On entry, AFP and IPIV contain the factored form
65: * of A. AP, AFP and IPIV will not be modified.
66: * = 'N': The matrix A will be copied to AFP and factored.
67: *
68: * UPLO (input) CHARACTER*1
69: * = 'U': Upper triangle of A is stored;
70: * = 'L': Lower triangle of A is stored.
71: *
72: * N (input) INTEGER
73: * The number of linear equations, i.e., the order of the
74: * matrix A. N >= 0.
75: *
76: * NRHS (input) INTEGER
77: * The number of right hand sides, i.e., the number of columns
78: * of the matrices B and X. NRHS >= 0.
79: *
80: * AP (input) COMPLEX*16 array, dimension (N*(N+1)/2)
81: * The upper or lower triangle of the symmetric matrix A, packed
82: * columnwise in a linear array. The j-th column of A is stored
83: * in the array AP as follows:
84: * if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
85: * if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
86: * See below for further details.
87: *
88: * AFP (input or output) COMPLEX*16 array, dimension (N*(N+1)/2)
89: * If FACT = 'F', then AFP is an input argument and on entry
90: * contains the block diagonal matrix D and the multipliers used
91: * to obtain the factor U or L from the factorization
92: * A = U*D*U**T or A = L*D*L**T as computed by ZSPTRF, stored as
93: * a packed triangular matrix in the same storage format as A.
94: *
95: * If FACT = 'N', then AFP is an output argument and on exit
96: * contains the block diagonal matrix D and the multipliers used
97: * to obtain the factor U or L from the factorization
98: * A = U*D*U**T or A = L*D*L**T as computed by ZSPTRF, stored as
99: * a packed triangular matrix in the same storage format as A.
100: *
101: * IPIV (input or output) INTEGER array, dimension (N)
102: * If FACT = 'F', then IPIV is an input argument and on entry
103: * contains details of the interchanges and the block structure
104: * of D, as determined by ZSPTRF.
105: * If IPIV(k) > 0, then rows and columns k and IPIV(k) were
106: * interchanged and D(k,k) is a 1-by-1 diagonal block.
107: * If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
108: * columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
109: * is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) =
110: * IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
111: * interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
112: *
113: * If FACT = 'N', then IPIV is an output argument and on exit
114: * contains details of the interchanges and the block structure
115: * of D, as determined by ZSPTRF.
116: *
117: * B (input) COMPLEX*16 array, dimension (LDB,NRHS)
118: * The N-by-NRHS right hand side matrix B.
119: *
120: * LDB (input) INTEGER
121: * The leading dimension of the array B. LDB >= max(1,N).
122: *
123: * X (output) COMPLEX*16 array, dimension (LDX,NRHS)
124: * If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.
125: *
126: * LDX (input) INTEGER
127: * The leading dimension of the array X. LDX >= max(1,N).
128: *
129: * RCOND (output) DOUBLE PRECISION
130: * The estimate of the reciprocal condition number of the matrix
131: * A. If RCOND is less than the machine precision (in
132: * particular, if RCOND = 0), the matrix is singular to working
133: * precision. This condition is indicated by a return code of
134: * INFO > 0.
135: *
136: * FERR (output) DOUBLE PRECISION array, dimension (NRHS)
137: * The estimated forward error bound for each solution vector
138: * X(j) (the j-th column of the solution matrix X).
139: * If XTRUE is the true solution corresponding to X(j), FERR(j)
140: * is an estimated upper bound for the magnitude of the largest
141: * element in (X(j) - XTRUE) divided by the magnitude of the
142: * largest element in X(j). The estimate is as reliable as
143: * the estimate for RCOND, and is almost always a slight
144: * overestimate of the true error.
145: *
146: * BERR (output) DOUBLE PRECISION array, dimension (NRHS)
147: * The componentwise relative backward error of each solution
148: * vector X(j) (i.e., the smallest relative change in
149: * any element of A or B that makes X(j) an exact solution).
150: *
151: * WORK (workspace) COMPLEX*16 array, dimension (2*N)
152: *
153: * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
154: *
155: * INFO (output) INTEGER
156: * = 0: successful exit
157: * < 0: if INFO = -i, the i-th argument had an illegal value
158: * > 0: if INFO = i, and i is
159: * <= N: D(i,i) is exactly zero. The factorization
160: * has been completed but the factor D is exactly
161: * singular, so the solution and error bounds could
162: * not be computed. RCOND = 0 is returned.
163: * = N+1: D is nonsingular, but RCOND is less than machine
164: * precision, meaning that the matrix is singular
165: * to working precision. Nevertheless, the
166: * solution and error bounds are computed because
167: * there are a number of situations where the
168: * computed solution can be more accurate than the
169: * value of RCOND would suggest.
170: *
171: * Further Details
172: * ===============
173: *
174: * The packed storage scheme is illustrated by the following example
175: * when N = 4, UPLO = 'U':
176: *
177: * Two-dimensional storage of the symmetric matrix A:
178: *
179: * a11 a12 a13 a14
180: * a22 a23 a24
181: * a33 a34 (aij = aji)
182: * a44
183: *
184: * Packed storage of the upper triangle of A:
185: *
186: * AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]
187: *
188: * =====================================================================
189: *
190: * .. Parameters ..
191: DOUBLE PRECISION ZERO
192: PARAMETER ( ZERO = 0.0D+0 )
193: * ..
194: * .. Local Scalars ..
195: LOGICAL NOFACT
196: DOUBLE PRECISION ANORM
197: * ..
198: * .. External Functions ..
199: LOGICAL LSAME
200: DOUBLE PRECISION DLAMCH, ZLANSP
201: EXTERNAL LSAME, DLAMCH, ZLANSP
202: * ..
203: * .. External Subroutines ..
204: EXTERNAL XERBLA, ZCOPY, ZLACPY, ZSPCON, ZSPRFS, ZSPTRF,
205: $ ZSPTRS
206: * ..
207: * .. Intrinsic Functions ..
208: INTRINSIC MAX
209: * ..
210: * .. Executable Statements ..
211: *
212: * Test the input parameters.
213: *
214: INFO = 0
215: NOFACT = LSAME( FACT, 'N' )
216: IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN
217: INFO = -1
218: ELSE IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) )
219: $ THEN
220: INFO = -2
221: ELSE IF( N.LT.0 ) THEN
222: INFO = -3
223: ELSE IF( NRHS.LT.0 ) THEN
224: INFO = -4
225: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
226: INFO = -9
227: ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
228: INFO = -11
229: END IF
230: IF( INFO.NE.0 ) THEN
231: CALL XERBLA( 'ZSPSVX', -INFO )
232: RETURN
233: END IF
234: *
235: IF( NOFACT ) THEN
236: *
237: * Compute the factorization A = U*D*U' or A = L*D*L'.
238: *
239: CALL ZCOPY( N*( N+1 ) / 2, AP, 1, AFP, 1 )
240: CALL ZSPTRF( UPLO, N, AFP, IPIV, INFO )
241: *
242: * Return if INFO is non-zero.
243: *
244: IF( INFO.GT.0 )THEN
245: RCOND = ZERO
246: RETURN
247: END IF
248: END IF
249: *
250: * Compute the norm of the matrix A.
251: *
252: ANORM = ZLANSP( 'I', UPLO, N, AP, RWORK )
253: *
254: * Compute the reciprocal of the condition number of A.
255: *
256: CALL ZSPCON( UPLO, N, AFP, IPIV, ANORM, RCOND, WORK, INFO )
257: *
258: * Compute the solution vectors X.
259: *
260: CALL ZLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
261: CALL ZSPTRS( UPLO, N, NRHS, AFP, IPIV, X, LDX, INFO )
262: *
263: * Use iterative refinement to improve the computed solutions and
264: * compute error bounds and backward error estimates for them.
265: *
266: CALL ZSPRFS( UPLO, N, NRHS, AP, AFP, IPIV, B, LDB, X, LDX, FERR,
267: $ BERR, WORK, RWORK, INFO )
268: *
269: * Set INFO = N+1 if the matrix is singular to working precision.
270: *
271: IF( RCOND.LT.DLAMCH( 'Epsilon' ) )
272: $ INFO = N + 1
273: *
274: RETURN
275: *
276: * End of ZSPSVX
277: *
278: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zspsvx.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack