LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
zspsvx.f
Go to the documentation of this file.
1*> \brief <b> ZSPSVX 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 ZSPSVX + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zspsvx.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zspsvx.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zspsvx.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE ZSPSVX( 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* DOUBLE PRECISION RCOND
28* ..
29* .. Array Arguments ..
30* INTEGER IPIV( * )
31* DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * )
32* COMPLEX*16 AFP( * ), AP( * ), B( LDB, * ), WORK( * ),
33* $ X( LDX, * )
34* ..
35*
36*
37*> \par Purpose:
38* =============
39*>
40*> \verbatim
41*>
42*> ZSPSVX 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*16 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*16 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 ZSPTRF, 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 ZSPTRF, 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 ZSPTRF.
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 ZSPTRF.
158*> \endverbatim
159*>
160*> \param[in] B
161*> \verbatim
162*> B is COMPLEX*16 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*16 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 DOUBLE PRECISION
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 DOUBLE PRECISION 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 DOUBLE PRECISION 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*16 array, dimension (2*N)
218*> \endverbatim
219*>
220*> \param[out] RWORK
221*> \verbatim
222*> RWORK is DOUBLE PRECISION 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 hpsvx
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 zspsvx( 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 DOUBLE PRECISION RCOND
286* ..
287* .. Array Arguments ..
288 INTEGER IPIV( * )
289 DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * )
290 COMPLEX*16 AFP( * ), AP( * ), B( LDB, * ), WORK( * ),
291 $ x( ldx, * )
292* ..
293*
294* =====================================================================
295*
296* .. Parameters ..
297 DOUBLE PRECISION ZERO
298 parameter( zero = 0.0d+0 )
299* ..
300* .. Local Scalars ..
301 LOGICAL NOFACT
302 DOUBLE PRECISION ANORM
303* ..
304* .. External Functions ..
305 LOGICAL LSAME
306 DOUBLE PRECISION DLAMCH, ZLANSP
307 EXTERNAL lsame, dlamch, zlansp
308* ..
309* .. External Subroutines ..
310 EXTERNAL xerbla, zcopy, zlacpy, zspcon, zsprfs, zsptrf,
311 $ zsptrs
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( 'ZSPSVX', -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 zcopy( n*( n+1 ) / 2, ap, 1, afp, 1 )
346 CALL zsptrf( 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 = zlansp( 'I', uplo, n, ap, rwork )
359*
360* Compute the reciprocal of the condition number of A.
361*
362 CALL zspcon( uplo, n, afp, ipiv, anorm, rcond, work, info )
363*
364* Compute the solution vectors X.
365*
366 CALL zlacpy( 'Full', n, nrhs, b, ldb, x, ldx )
367 CALL zsptrs( 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 zsprfs( 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.dlamch( 'Epsilon' ) )
378 $ info = n + 1
379*
380 RETURN
381*
382* End of ZSPSVX
383*
384 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
Definition zcopy.f:81
subroutine zspcon(uplo, n, ap, ipiv, anorm, rcond, work, info)
ZSPCON
Definition zspcon.f:118
subroutine zsprfs(uplo, n, nrhs, ap, afp, ipiv, b, ldb, x, ldx, ferr, berr, work, rwork, info)
ZSPRFS
Definition zsprfs.f:180
subroutine zspsvx(fact, uplo, n, nrhs, ap, afp, ipiv, b, ldb, x, ldx, rcond, ferr, berr, work, rwork, info)
ZSPSVX computes the solution to system of linear equations A * X = B for OTHER matrices
Definition zspsvx.f:277
subroutine zsptrf(uplo, n, ap, ipiv, info)
ZSPTRF
Definition zsptrf.f:158
subroutine zsptrs(uplo, n, nrhs, ap, ipiv, b, ldb, info)
ZSPTRS
Definition zsptrs.f:115
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
Definition zlacpy.f:103