OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
complex16 Other Solve Routines
LAPACK » Other Solve Routines

Functions

subroutine zgglse (m, n, p, a, lda, b, ldb, c, d, x, work, lwork, info)
  ZGGLSE solves overdetermined or underdetermined systems for OTHER matrices
subroutine zhpsv (uplo, n, nrhs, ap, ipiv, b, ldb, info)
  ZHPSV computes the solution to system of linear equations A * X = B for OTHER matrices
subroutine zhpsvx (fact, uplo, n, nrhs, ap, afp, ipiv, b, ldb, x, ldx, rcond, ferr, berr, work, rwork, info)
  ZHPSVX computes the solution to system of linear equations A * X = B for OTHER matrices
subroutine zpbsv (uplo, n, kd, nrhs, ab, ldab, b, ldb, info)
  ZPBSV computes the solution to system of linear equations A * X = B for OTHER matrices
subroutine zpbsvx (fact, uplo, n, kd, nrhs, ab, ldab, afb, ldafb, equed, s, b, ldb, x, ldx, rcond, ferr, berr, work, rwork, info)
  ZPBSVX computes the solution to system of linear equations A * X = B for OTHER matrices
subroutine zppsv (uplo, n, nrhs, ap, b, ldb, info)
  ZPPSV computes the solution to system of linear equations A * X = B for OTHER matrices
subroutine zppsvx (fact, uplo, n, nrhs, ap, afp, equed, s, b, ldb, x, ldx, rcond, ferr, berr, work, rwork, info)
  ZPPSVX computes the solution to system of linear equations A * X = B for OTHER matrices
subroutine zspsv (uplo, n, nrhs, ap, ipiv, b, ldb, info)
  ZSPSV computes the solution to system of linear equations A * X = B for OTHER matrices
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

Detailed Description

This is the group of complex16 Other Solve routines

Function Documentation

◆ zgglse()

subroutine zgglse ( integer m,
integer n,
integer p,
complex*16, dimension( lda, * ) a,
integer lda,
complex*16, dimension( ldb, * ) b,
integer ldb,
complex*16, dimension( * ) c,
complex*16, dimension( * ) d,
complex*16, dimension( * ) x,
complex*16, dimension( * ) work,
integer lwork,
integer info )

ZGGLSE solves overdetermined or underdetermined systems for OTHER matrices

Download ZGGLSE + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> ZGGLSE solves the linear equality-constrained least squares (LSE)
!> problem:
!>
!>         minimize || c - A*x ||_2   subject to   B*x = d
!>
!> where A is an M-by-N matrix, B is a P-by-N matrix, c is a given
!> M-vector, and d is a given P-vector. It is assumed that
!> P <= N <= M+P, and
!>
!>          rank(B) = P and  rank( (A) ) = N.
!>                               ( (B) )
!>
!> These conditions ensure that the LSE problem has a unique solution,
!> which is obtained using a generalized RQ factorization of the
!> matrices (B, A) given by
!>
!>    B = (0 R)*Q,   A = Z*T*Q.
!>
Parameters
[in]M
!>          M is INTEGER
!>          The number of rows of the matrix A.  M >= 0.
!>
[in]N
!>          N is INTEGER
!>          The number of columns of the matrices A and B. N >= 0.
!>
[in]P
!>          P is INTEGER
!>          The number of rows of the matrix B. 0 <= P <= N <= M+P.
!>
[in,out]A
!>          A is COMPLEX*16 array, dimension (LDA,N)
!>          On entry, the M-by-N matrix A.
!>          On exit, the elements on and above the diagonal of the array
!>          contain the min(M,N)-by-N upper trapezoidal matrix T.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A. LDA >= max(1,M).
!>
[in,out]B
!>          B is COMPLEX*16 array, dimension (LDB,N)
!>          On entry, the P-by-N matrix B.
!>          On exit, the upper triangle of the subarray B(1:P,N-P+1:N)
!>          contains the P-by-P upper triangular matrix R.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B. LDB >= max(1,P).
!>
[in,out]C
!>          C is COMPLEX*16 array, dimension (M)
!>          On entry, C contains the right hand side vector for the
!>          least squares part of the LSE problem.
!>          On exit, the residual sum of squares for the solution
!>          is given by the sum of squares of elements N-P+1 to M of
!>          vector C.
!>
[in,out]D
!>          D is COMPLEX*16 array, dimension (P)
!>          On entry, D contains the right hand side vector for the
!>          constrained equation.
!>          On exit, D is destroyed.
!>
[out]X
!>          X is COMPLEX*16 array, dimension (N)
!>          On exit, X is the solution of the LSE problem.
!>
[out]WORK
!>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK. LWORK >= max(1,M+N+P).
!>          For optimum performance LWORK >= P+min(M,N)+max(M,N)*NB,
!>          where NB is an upper bound for the optimal blocksizes for
!>          ZGEQRF, CGERQF, ZUNMQR and CUNMRQ.
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal size of the WORK array, returns
!>          this value as the first entry of the WORK array, and no error
!>          message related to LWORK is issued by XERBLA.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit.
!>          < 0:  if INFO = -i, the i-th argument had an illegal value.
!>          = 1:  the upper triangular factor R associated with B in the
!>                generalized RQ factorization of the pair (B, A) is
!>                singular, so that rank(B) < P; the least squares
!>                solution could not be computed.
!>          = 2:  the (N-P) by (N-P) part of the upper trapezoidal factor
!>                T associated with A in the generalized RQ factorization
!>                of the pair (B, A) is singular, so that
!>                rank( (A) ) < N; the least squares solution could not
!>                    ( (B) )
!>                be computed.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 178 of file zgglse.f.

180*
181* -- LAPACK driver routine --
182* -- LAPACK is a software package provided by Univ. of Tennessee, --
183* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
184*
185* .. Scalar Arguments ..
186 INTEGER INFO, LDA, LDB, LWORK, M, N, P
187* ..
188* .. Array Arguments ..
189 COMPLEX*16 A( LDA, * ), B( LDB, * ), C( * ), D( * ),
190 $ WORK( * ), X( * )
191* ..
192*
193* =====================================================================
194*
195* .. Parameters ..
196 COMPLEX*16 CONE
197 parameter( cone = ( 1.0d+0, 0.0d+0 ) )
198* ..
199* .. Local Scalars ..
200 LOGICAL LQUERY
201 INTEGER LOPT, LWKMIN, LWKOPT, MN, NB, NB1, NB2, NB3,
202 $ NB4, NR
203* ..
204* .. External Subroutines ..
205 EXTERNAL xerbla, zaxpy, zcopy, zgemv, zggrqf, ztrmv,
206 $ ztrtrs, zunmqr, zunmrq
207* ..
208* .. External Functions ..
209 INTEGER ILAENV
210 EXTERNAL ilaenv
211* ..
212* .. Intrinsic Functions ..
213 INTRINSIC int, max, min
214* ..
215* .. Executable Statements ..
216*
217* Test the input parameters
218*
219 info = 0
220 mn = min( m, n )
221 lquery = ( lwork.EQ.-1 )
222 IF( m.LT.0 ) THEN
223 info = -1
224 ELSE IF( n.LT.0 ) THEN
225 info = -2
226 ELSE IF( p.LT.0 .OR. p.GT.n .OR. p.LT.n-m ) THEN
227 info = -3
228 ELSE IF( lda.LT.max( 1, m ) ) THEN
229 info = -5
230 ELSE IF( ldb.LT.max( 1, p ) ) THEN
231 info = -7
232 END IF
233*
234* Calculate workspace
235*
236 IF( info.EQ.0) THEN
237 IF( n.EQ.0 ) THEN
238 lwkmin = 1
239 lwkopt = 1
240 ELSE
241 nb1 = ilaenv( 1, 'ZGEQRF', ' ', m, n, -1, -1 )
242 nb2 = ilaenv( 1, 'ZGERQF', ' ', m, n, -1, -1 )
243 nb3 = ilaenv( 1, 'ZUNMQR', ' ', m, n, p, -1 )
244 nb4 = ilaenv( 1, 'ZUNMRQ', ' ', m, n, p, -1 )
245 nb = max( nb1, nb2, nb3, nb4 )
246 lwkmin = m + n + p
247 lwkopt = p + mn + max( m, n )*nb
248 END IF
249 work( 1 ) = lwkopt
250*
251 IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
252 info = -12
253 END IF
254 END IF
255*
256 IF( info.NE.0 ) THEN
257 CALL xerbla( 'ZGGLSE', -info )
258 RETURN
259 ELSE IF( lquery ) THEN
260 RETURN
261 END IF
262*
263* Quick return if possible
264*
265 IF( n.EQ.0 )
266 $ RETURN
267*
268* Compute the GRQ factorization of matrices B and A:
269*
270* B*Q**H = ( 0 T12 ) P Z**H*A*Q**H = ( R11 R12 ) N-P
271* N-P P ( 0 R22 ) M+P-N
272* N-P P
273*
274* where T12 and R11 are upper triangular, and Q and Z are
275* unitary.
276*
277 CALL zggrqf( p, m, n, b, ldb, work, a, lda, work( p+1 ),
278 $ work( p+mn+1 ), lwork-p-mn, info )
279 lopt = dble( work( p+mn+1 ) )
280*
281* Update c = Z**H *c = ( c1 ) N-P
282* ( c2 ) M+P-N
283*
284 CALL zunmqr( 'Left', 'Conjugate Transpose', m, 1, mn, a, lda,
285 $ work( p+1 ), c, max( 1, m ), work( p+mn+1 ),
286 $ lwork-p-mn, info )
287 lopt = max( lopt, int( work( p+mn+1 ) ) )
288*
289* Solve T12*x2 = d for x2
290*
291 IF( p.GT.0 ) THEN
292 CALL ztrtrs( 'Upper', 'No transpose', 'Non-unit', p, 1,
293 $ b( 1, n-p+1 ), ldb, d, p, info )
294*
295 IF( info.GT.0 ) THEN
296 info = 1
297 RETURN
298 END IF
299*
300* Put the solution in X
301*
302 CALL zcopy( p, d, 1, x( n-p+1 ), 1 )
303*
304* Update c1
305*
306 CALL zgemv( 'No transpose', n-p, p, -cone, a( 1, n-p+1 ), lda,
307 $ d, 1, cone, c, 1 )
308 END IF
309*
310* Solve R11*x1 = c1 for x1
311*
312 IF( n.GT.p ) THEN
313 CALL ztrtrs( 'Upper', 'No transpose', 'Non-unit', n-p, 1,
314 $ a, lda, c, n-p, info )
315*
316 IF( info.GT.0 ) THEN
317 info = 2
318 RETURN
319 END IF
320*
321* Put the solutions in X
322*
323 CALL zcopy( n-p, c, 1, x, 1 )
324 END IF
325*
326* Compute the residual vector:
327*
328 IF( m.LT.n ) THEN
329 nr = m + p - n
330 IF( nr.GT.0 )
331 $ CALL zgemv( 'No transpose', nr, n-m, -cone, a( n-p+1, m+1 ),
332 $ lda, d( nr+1 ), 1, cone, c( n-p+1 ), 1 )
333 ELSE
334 nr = p
335 END IF
336 IF( nr.GT.0 ) THEN
337 CALL ztrmv( 'Upper', 'No transpose', 'Non unit', nr,
338 $ a( n-p+1, n-p+1 ), lda, d, 1 )
339 CALL zaxpy( nr, -cone, d, 1, c( n-p+1 ), 1 )
340 END IF
341*
342* Backward transformation x = Q**H*x
343*
344 CALL zunmrq( 'Left', 'Conjugate Transpose', n, 1, p, b, ldb,
345 $ work( 1 ), x, n, work( p+mn+1 ), lwork-p-mn, info )
346 work( 1 ) = p + mn + max( lopt, int( work( p+mn+1 ) ) )
347*
348 RETURN
349*
350* End of ZGGLSE
351*
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
zunmqr
subroutine zunmqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
ZUNMQR
Definition zunmqr.f:167
ztrtrs
subroutine ztrtrs(uplo, trans, diag, n, nrhs, a, lda, b, ldb, info)
ZTRTRS
Definition ztrtrs.f:140
zunmrq
subroutine zunmrq(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
ZUNMRQ
Definition zunmrq.f:167
zggrqf
subroutine zggrqf(m, p, n, a, lda, taua, b, ldb, taub, work, lwork, info)
ZGGRQF
Definition zggrqf.f:214
zcopy
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
Definition zcopy.f:81
zaxpy
subroutine zaxpy(n, za, zx, incx, zy, incy)
ZAXPY
Definition zaxpy.f:88
ztrmv
subroutine ztrmv(uplo, trans, diag, n, a, lda, x, incx)
ZTRMV
Definition ztrmv.f:147
zgemv
subroutine zgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
ZGEMV
Definition zgemv.f:158
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21

◆ zhpsv()

subroutine zhpsv ( character uplo,
integer n,
integer nrhs,
complex*16, dimension( * ) ap,
integer, dimension( * ) ipiv,
complex*16, dimension( ldb, * ) b,
integer ldb,
integer info )

ZHPSV computes the solution to system of linear equations A * X = B for OTHER matrices

Download ZHPSV + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> ZHPSV computes the solution to a complex system of linear equations
!>    A * X = B,
!> where A is an N-by-N Hermitian matrix stored in packed format and X
!> and B are N-by-NRHS matrices.
!>
!> The diagonal pivoting method is used to factor A as
!>    A = U * D * U**H,  if UPLO = 'U', or
!>    A = L * D * L**H,  if UPLO = 'L',
!> where U (or L) is a product of permutation and unit upper (lower)
!> triangular matrices, D is Hermitian and block diagonal with 1-by-1
!> and 2-by-2 diagonal blocks.  The factored form of A is then used to
!> solve the system of equations A * X = B.
!>
Parameters
[in]UPLO
!>          UPLO is CHARACTER*1
!>          = 'U':  Upper triangle of A is stored;
!>          = 'L':  Lower triangle of A is stored.
!>
[in]N
!>          N is INTEGER
!>          The number of linear equations, i.e., the order of the
!>          matrix A.  N >= 0.
!>
[in]NRHS
!>          NRHS is INTEGER
!>          The number of right hand sides, i.e., the number of columns
!>          of the matrix B.  NRHS >= 0.
!>
[in,out]AP
!>          AP is COMPLEX*16 array, dimension (N*(N+1)/2)
!>          On entry, the upper or lower triangle of the Hermitian matrix
!>          A, packed columnwise in a linear array.  The j-th column of A
!>          is stored in the array AP as follows:
!>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
!>          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
!>          See below for further details.
!>
!>          On exit, the block diagonal matrix D and the multipliers used
!>          to obtain the factor U or L from the factorization
!>          A = U*D*U**H or A = L*D*L**H as computed by ZHPTRF, stored as
!>          a packed triangular matrix in the same storage format as A.
!>
[out]IPIV
!>          IPIV is INTEGER array, dimension (N)
!>          Details of the interchanges and the block structure of D, as
!>          determined by ZHPTRF.  If IPIV(k) > 0, then rows and columns
!>          k and IPIV(k) were interchanged, and D(k,k) is a 1-by-1
!>          diagonal block.  If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0,
!>          then rows and columns k-1 and -IPIV(k) were interchanged and
!>          D(k-1:k,k-1:k) is a 2-by-2 diagonal block.  If UPLO = 'L' and
!>          IPIV(k) = IPIV(k+1) < 0, then rows and columns k+1 and
!>          -IPIV(k) were interchanged and D(k:k+1,k:k+1) is a 2-by-2
!>          diagonal block.
!>
[in,out]B
!>          B is COMPLEX*16 array, dimension (LDB,NRHS)
!>          On entry, the N-by-NRHS right hand side matrix B.
!>          On exit, if INFO = 0, the N-by-NRHS solution matrix X.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= max(1,N).
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO = i, D(i,i) is exactly zero.  The factorization
!>                has been completed, but the block diagonal matrix D is
!>                exactly singular, so the solution could not be
!>                computed.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  The packed storage scheme is illustrated by the following example
!>  when N = 4, UPLO = 'U':
!>
!>  Two-dimensional storage of the Hermitian matrix A:
!>
!>     a11 a12 a13 a14
!>         a22 a23 a24
!>             a33 a34     (aij = conjg(aji))
!>                 a44
!>
!>  Packed storage of the upper triangle of A:
!>
!>  AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]
!>

Definition at line 161 of file zhpsv.f.

162*
163* -- LAPACK driver routine --
164* -- LAPACK is a software package provided by Univ. of Tennessee, --
165* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
166*
167* .. Scalar Arguments ..
168 CHARACTER UPLO
169 INTEGER INFO, LDB, N, NRHS
170* ..
171* .. Array Arguments ..
172 INTEGER IPIV( * )
173 COMPLEX*16 AP( * ), B( LDB, * )
174* ..
175*
176* =====================================================================
177*
178* .. External Functions ..
179 LOGICAL LSAME
180 EXTERNAL lsame
181* ..
182* .. External Subroutines ..
183 EXTERNAL xerbla, zhptrf, zhptrs
184* ..
185* .. Intrinsic Functions ..
186 INTRINSIC max
187* ..
188* .. Executable Statements ..
189*
190* Test the input parameters.
191*
192 info = 0
193 IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lsame( uplo, 'L' ) ) THEN
194 info = -1
195 ELSE IF( n.LT.0 ) THEN
196 info = -2
197 ELSE IF( nrhs.LT.0 ) THEN
198 info = -3
199 ELSE IF( ldb.LT.max( 1, n ) ) THEN
200 info = -7
201 END IF
202 IF( info.NE.0 ) THEN
203 CALL xerbla( 'ZHPSV ', -info )
204 RETURN
205 END IF
206*
207* Compute the factorization A = U*D*U**H or A = L*D*L**H.
208*
209 CALL zhptrf( uplo, n, ap, ipiv, info )
210 IF( info.EQ.0 ) THEN
211*
212* Solve the system A*X = B, overwriting B with X.
213*
214 CALL zhptrs( uplo, n, nrhs, ap, ipiv, b, ldb, info )
215*
216 END IF
217 RETURN
218*
219* End of ZHPSV
220*
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:53
zhptrs
subroutine zhptrs(uplo, n, nrhs, ap, ipiv, b, ldb, info)
ZHPTRS
Definition zhptrs.f:115
zhptrf
subroutine zhptrf(uplo, n, ap, ipiv, info)
ZHPTRF
Definition zhptrf.f:159

◆ zhpsvx()

subroutine zhpsvx ( character fact,
character uplo,
integer n,
integer nrhs,
complex*16, dimension( * ) ap,
complex*16, dimension( * ) afp,
integer, dimension( * ) ipiv,
complex*16, dimension( ldb, * ) b,
integer ldb,
complex*16, dimension( ldx, * ) x,
integer ldx,
double precision rcond,
double precision, dimension( * ) ferr,
double precision, dimension( * ) berr,
complex*16, dimension( * ) work,
double precision, dimension( * ) rwork,
integer info )

ZHPSVX computes the solution to system of linear equations A * X = B for OTHER matrices

Download ZHPSVX + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> ZHPSVX uses the diagonal pivoting factorization A = U*D*U**H or
!> A = L*D*L**H to compute the solution to a complex system of linear
!> equations A * X = B, where A is an N-by-N Hermitian matrix stored
!> in packed format and X and B are N-by-NRHS matrices.
!>
!> Error bounds on the solution and a condition estimate are also
!> provided.
!>
Description:
!>
!> The following steps are performed:
!>
!> 1. If FACT = 'N', the diagonal pivoting method is used to factor A as
!>       A = U * D * U**H,  if UPLO = 'U', or
!>       A = L * D * L**H,  if UPLO = 'L',
!>    where U (or L) is a product of permutation and unit upper (lower)
!>    triangular matrices and D is Hermitian and block diagonal with
!>    1-by-1 and 2-by-2 diagonal blocks.
!>
!> 2. If some D(i,i)=0, so that D is exactly singular, then the routine
!>    returns with INFO = i. Otherwise, the factored form of A is used
!>    to estimate the condition number of the matrix A.  If the
!>    reciprocal of the condition number is less than machine precision,
!>    INFO = N+1 is returned as a warning, but the routine still goes on
!>    to solve for X and compute error bounds as described below.
!>
!> 3. The system of equations is solved for X using the factored form
!>    of A.
!>
!> 4. Iterative refinement is applied to improve the computed solution
!>    matrix and calculate error bounds and backward error estimates
!>    for it.
!>
Parameters
[in]FACT
!>          FACT is CHARACTER*1
!>          Specifies whether or not the factored form of A has been
!>          supplied on entry.
!>          = 'F':  On entry, AFP and IPIV contain the factored form of
!>                  A.  AFP and IPIV will not be modified.
!>          = 'N':  The matrix A will be copied to AFP and factored.
!>
[in]UPLO
!>          UPLO is CHARACTER*1
!>          = 'U':  Upper triangle of A is stored;
!>          = 'L':  Lower triangle of A is stored.
!>
[in]N
!>          N is INTEGER
!>          The number of linear equations, i.e., the order of the
!>          matrix A.  N >= 0.
!>
[in]NRHS
!>          NRHS is INTEGER
!>          The number of right hand sides, i.e., the number of columns
!>          of the matrices B and X.  NRHS >= 0.
!>
[in]AP
!>          AP is COMPLEX*16 array, dimension (N*(N+1)/2)
!>          The upper or lower triangle of the Hermitian matrix A, packed
!>          columnwise in a linear array.  The j-th column of A is stored
!>          in the array AP as follows:
!>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
!>          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
!>          See below for further details.
!>
[in,out]AFP
!>          AFP is COMPLEX*16 array, dimension (N*(N+1)/2)
!>          If FACT = 'F', then AFP is an input argument and on entry
!>          contains the block diagonal matrix D and the multipliers used
!>          to obtain the factor U or L from the factorization
!>          A = U*D*U**H or A = L*D*L**H as computed by ZHPTRF, stored as
!>          a packed triangular matrix in the same storage format as A.
!>
!>          If FACT = 'N', then AFP is an output argument and on exit
!>          contains the block diagonal matrix D and the multipliers used
!>          to obtain the factor U or L from the factorization
!>          A = U*D*U**H or A = L*D*L**H as computed by ZHPTRF, stored as
!>          a packed triangular matrix in the same storage format as A.
!>
[in,out]IPIV
!>          IPIV is INTEGER array, dimension (N)
!>          If FACT = 'F', then IPIV is an input argument and on entry
!>          contains details of the interchanges and the block structure
!>          of D, as determined by ZHPTRF.
!>          If IPIV(k) > 0, then rows and columns k and IPIV(k) were
!>          interchanged and D(k,k) is a 1-by-1 diagonal block.
!>          If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
!>          columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
!>          is a 2-by-2 diagonal block.  If UPLO = 'L' and IPIV(k) =
!>          IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
!>          interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
!>
!>          If FACT = 'N', then IPIV is an output argument and on exit
!>          contains details of the interchanges and the block structure
!>          of D, as determined by ZHPTRF.
!>
[in]B
!>          B is COMPLEX*16 array, dimension (LDB,NRHS)
!>          The N-by-NRHS right hand side matrix B.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= max(1,N).
!>
[out]X
!>          X is COMPLEX*16 array, dimension (LDX,NRHS)
!>          If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.
!>
[in]LDX
!>          LDX is INTEGER
!>          The leading dimension of the array X.  LDX >= max(1,N).
!>
[out]RCOND
!>          RCOND is DOUBLE PRECISION
!>          The estimate of the reciprocal condition number of the matrix
!>          A.  If RCOND is less than the machine precision (in
!>          particular, if RCOND = 0), the matrix is singular to working
!>          precision.  This condition is indicated by a return code of
!>          INFO > 0.
!>
[out]FERR
!>          FERR is DOUBLE PRECISION array, dimension (NRHS)
!>          The estimated forward error bound for each solution vector
!>          X(j) (the j-th column of the solution matrix X).
!>          If XTRUE is the true solution corresponding to X(j), FERR(j)
!>          is an estimated upper bound for the magnitude of the largest
!>          element in (X(j) - XTRUE) divided by the magnitude of the
!>          largest element in X(j).  The estimate is as reliable as
!>          the estimate for RCOND, and is almost always a slight
!>          overestimate of the true error.
!>
[out]BERR
!>          BERR is DOUBLE PRECISION array, dimension (NRHS)
!>          The componentwise relative backward error of each solution
!>          vector X(j) (i.e., the smallest relative change in
!>          any element of A or B that makes X(j) an exact solution).
!>
[out]WORK
!>          WORK is COMPLEX*16 array, dimension (2*N)
!>
[out]RWORK
!>          RWORK is DOUBLE PRECISION array, dimension (N)
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0: successful exit
!>          < 0: if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO = i, and i is
!>                <= N:  D(i,i) is exactly zero.  The factorization
!>                       has been completed but the factor D is exactly
!>                       singular, so the solution and error bounds could
!>                       not be computed. RCOND = 0 is returned.
!>                = N+1: D is nonsingular, but RCOND is less than machine
!>                       precision, meaning that the matrix is singular
!>                       to working precision.  Nevertheless, the
!>                       solution and error bounds are computed because
!>                       there are a number of situations where the
!>                       computed solution can be more accurate than the
!>                       value of RCOND would suggest.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  The packed storage scheme is illustrated by the following example
!>  when N = 4, UPLO = 'U':
!>
!>  Two-dimensional storage of the Hermitian matrix A:
!>
!>     a11 a12 a13 a14
!>         a22 a23 a24
!>             a33 a34     (aij = conjg(aji))
!>                 a44
!>
!>  Packed storage of the upper triangle of A:
!>
!>  AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]
!>

Definition at line 275 of file zhpsvx.f.

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, ZLANHP
307 EXTERNAL lsame, dlamch, zlanhp
308* ..
309* .. External Subroutines ..
310 EXTERNAL xerbla, zcopy, zhpcon, zhprfs, zhptrf, zhptrs,
311 $ zlacpy
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( 'ZHPSVX', -info )
338 RETURN
339 END IF
340*
341 IF( nofact ) THEN
342*
343* Compute the factorization A = U*D*U**H or A = L*D*L**H.
344*
345 CALL zcopy( n*( n+1 ) / 2, ap, 1, afp, 1 )
346 CALL zhptrf( 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 = zlanhp( 'I', uplo, n, ap, rwork )
359*
360* Compute the reciprocal of the condition number of A.
361*
362 CALL zhpcon( 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 zhptrs( 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 zhprfs( 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 ZHPSVX
383*
zlacpy
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
zlanhp
double precision function zlanhp(norm, uplo, n, ap, work)
ZLANHP returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition zlanhp.f:117
zhpcon
subroutine zhpcon(uplo, n, ap, ipiv, anorm, rcond, work, info)
ZHPCON
Definition zhpcon.f:118
zhprfs
subroutine zhprfs(uplo, n, nrhs, ap, afp, ipiv, b, ldb, x, ldx, ferr, berr, work, rwork, info)
ZHPRFS
Definition zhprfs.f:180
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69

◆ zpbsv()

subroutine zpbsv ( character uplo,
integer n,
integer kd,
integer nrhs,
complex*16, dimension( ldab, * ) ab,
integer ldab,
complex*16, dimension( ldb, * ) b,
integer ldb,
integer info )

ZPBSV computes the solution to system of linear equations A * X = B for OTHER matrices

Download ZPBSV + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> ZPBSV computes the solution to a complex system of linear equations
!>    A * X = B,
!> where A is an N-by-N Hermitian positive definite band matrix and X
!> and B are N-by-NRHS matrices.
!>
!> The Cholesky decomposition is used to factor A as
!>    A = U**H * U,  if UPLO = 'U', or
!>    A = L * L**H,  if UPLO = 'L',
!> where U is an upper triangular band matrix, and L is a lower
!> triangular band matrix, with the same number of superdiagonals or
!> subdiagonals as A.  The factored form of A is then used to solve the
!> system of equations A * X = B.
!>
Parameters
[in]UPLO
!>          UPLO is CHARACTER*1
!>          = 'U':  Upper triangle of A is stored;
!>          = 'L':  Lower triangle of A is stored.
!>
[in]N
!>          N is INTEGER
!>          The number of linear equations, i.e., the order of the
!>          matrix A.  N >= 0.
!>
[in]KD
!>          KD is INTEGER
!>          The number of superdiagonals of the matrix A if UPLO = 'U',
!>          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
!>
[in]NRHS
!>          NRHS is INTEGER
!>          The number of right hand sides, i.e., the number of columns
!>          of the matrix B.  NRHS >= 0.
!>
[in,out]AB
!>          AB is COMPLEX*16 array, dimension (LDAB,N)
!>          On entry, the upper or lower triangle of the Hermitian band
!>          matrix A, stored in the first KD+1 rows of the array.  The
!>          j-th column of A is stored in the j-th column of the array AB
!>          as follows:
!>          if UPLO = 'U', AB(KD+1+i-j,j) = A(i,j) for max(1,j-KD)<=i<=j;
!>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(N,j+KD).
!>          See below for further details.
!>
!>          On exit, if INFO = 0, the triangular factor U or L from the
!>          Cholesky factorization A = U**H *U or A = L*L**H of the band
!>          matrix A, in the same storage format as A.
!>
[in]LDAB
!>          LDAB is INTEGER
!>          The leading dimension of the array AB.  LDAB >= KD+1.
!>
[in,out]B
!>          B is COMPLEX*16 array, dimension (LDB,NRHS)
!>          On entry, the N-by-NRHS right hand side matrix B.
!>          On exit, if INFO = 0, the N-by-NRHS solution matrix X.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= max(1,N).
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO = i, the leading minor of order i of A is not
!>                positive definite, so the factorization could not be
!>                completed, and the solution has not been computed.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  The band storage scheme is illustrated by the following example, when
!>  N = 6, KD = 2, and UPLO = 'U':
!>
!>  On entry:                       On exit:
!>
!>      *    *   a13  a24  a35  a46      *    *   u13  u24  u35  u46
!>      *   a12  a23  a34  a45  a56      *   u12  u23  u34  u45  u56
!>     a11  a22  a33  a44  a55  a66     u11  u22  u33  u44  u55  u66
!>
!>  Similarly, if UPLO = 'L' the format of A is as follows:
!>
!>  On entry:                       On exit:
!>
!>     a11  a22  a33  a44  a55  a66     l11  l22  l33  l44  l55  l66
!>     a21  a32  a43  a54  a65   *      l21  l32  l43  l54  l65   *
!>     a31  a42  a53  a64   *    *      l31  l42  l53  l64   *    *
!>
!>  Array elements marked * are not used by the routine.
!>

Definition at line 163 of file zpbsv.f.

164*
165* -- LAPACK driver routine --
166* -- LAPACK is a software package provided by Univ. of Tennessee, --
167* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
168*
169* .. Scalar Arguments ..
170 CHARACTER UPLO
171 INTEGER INFO, KD, LDAB, LDB, N, NRHS
172* ..
173* .. Array Arguments ..
174 COMPLEX*16 AB( LDAB, * ), B( LDB, * )
175* ..
176*
177* =====================================================================
178*
179* .. External Functions ..
180 LOGICAL LSAME
181 EXTERNAL lsame
182* ..
183* .. External Subroutines ..
184 EXTERNAL xerbla, zpbtrf, zpbtrs
185* ..
186* .. Intrinsic Functions ..
187 INTRINSIC max
188* ..
189* .. Executable Statements ..
190*
191* Test the input parameters.
192*
193 info = 0
194 IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lsame( uplo, 'L' ) ) THEN
195 info = -1
196 ELSE IF( n.LT.0 ) THEN
197 info = -2
198 ELSE IF( kd.LT.0 ) THEN
199 info = -3
200 ELSE IF( nrhs.LT.0 ) THEN
201 info = -4
202 ELSE IF( ldab.LT.kd+1 ) THEN
203 info = -6
204 ELSE IF( ldb.LT.max( 1, n ) ) THEN
205 info = -8
206 END IF
207 IF( info.NE.0 ) THEN
208 CALL xerbla( 'ZPBSV ', -info )
209 RETURN
210 END IF
211*
212* Compute the Cholesky factorization A = U**H *U or A = L*L**H.
213*
214 CALL zpbtrf( uplo, n, kd, ab, ldab, info )
215 IF( info.EQ.0 ) THEN
216*
217* Solve the system A*X = B, overwriting B with X.
218*
219 CALL zpbtrs( uplo, n, kd, nrhs, ab, ldab, b, ldb, info )
220*
221 END IF
222 RETURN
223*
224* End of ZPBSV
225*
zpbtrs
subroutine zpbtrs(uplo, n, kd, nrhs, ab, ldab, b, ldb, info)
ZPBTRS
Definition zpbtrs.f:121
zpbtrf
subroutine zpbtrf(uplo, n, kd, ab, ldab, info)
ZPBTRF
Definition zpbtrf.f:142

◆ zpbsvx()

subroutine zpbsvx ( character fact,
character uplo,
integer n,
integer kd,
integer nrhs,
complex*16, dimension( ldab, * ) ab,
integer ldab,
complex*16, dimension( ldafb, * ) afb,
integer ldafb,
character equed,
double precision, dimension( * ) s,
complex*16, dimension( ldb, * ) b,
integer ldb,
complex*16, dimension( ldx, * ) x,
integer ldx,
double precision rcond,
double precision, dimension( * ) ferr,
double precision, dimension( * ) berr,
complex*16, dimension( * ) work,
double precision, dimension( * ) rwork,
integer info )

ZPBSVX computes the solution to system of linear equations A * X = B for OTHER matrices

Download ZPBSVX + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> ZPBSVX uses the Cholesky factorization A = U**H*U or A = L*L**H to
!> compute the solution to a complex system of linear equations
!>    A * X = B,
!> where A is an N-by-N Hermitian positive definite band matrix and X
!> and B are N-by-NRHS matrices.
!>
!> Error bounds on the solution and a condition estimate are also
!> provided.
!>
Description:
!>
!> The following steps are performed:
!>
!> 1. If FACT = 'E', real scaling factors are computed to equilibrate
!>    the system:
!>       diag(S) * A * diag(S) * inv(diag(S)) * X = diag(S) * B
!>    Whether or not the system will be equilibrated depends on the
!>    scaling of the matrix A, but if equilibration is used, A is
!>    overwritten by diag(S)*A*diag(S) and B by diag(S)*B.
!>
!> 2. If FACT = 'N' or 'E', the Cholesky decomposition is used to
!>    factor the matrix A (after equilibration if FACT = 'E') as
!>       A = U**H * U,  if UPLO = 'U', or
!>       A = L * L**H,  if UPLO = 'L',
!>    where U is an upper triangular band matrix, and L is a lower
!>    triangular band matrix.
!>
!> 3. If the leading i-by-i principal minor is not positive definite,
!>    then the routine returns with INFO = i. Otherwise, the factored
!>    form of A is used to estimate the condition number of the matrix
!>    A.  If the reciprocal of the condition number is less than machine
!>    precision, INFO = N+1 is returned as a warning, but the routine
!>    still goes on to solve for X and compute error bounds as
!>    described below.
!>
!> 4. The system of equations is solved for X using the factored form
!>    of A.
!>
!> 5. Iterative refinement is applied to improve the computed solution
!>    matrix and calculate error bounds and backward error estimates
!>    for it.
!>
!> 6. If equilibration was used, the matrix X is premultiplied by
!>    diag(S) so that it solves the original system before
!>    equilibration.
!>
Parameters
[in]FACT
!>          FACT is CHARACTER*1
!>          Specifies whether or not the factored form of the matrix A is
!>          supplied on entry, and if not, whether the matrix A should be
!>          equilibrated before it is factored.
!>          = 'F':  On entry, AFB contains the factored form of A.
!>                  If EQUED = 'Y', the matrix A has been equilibrated
!>                  with scaling factors given by S.  AB and AFB will not
!>                  be modified.
!>          = 'N':  The matrix A will be copied to AFB and factored.
!>          = 'E':  The matrix A will be equilibrated if necessary, then
!>                  copied to AFB and factored.
!>
[in]UPLO
!>          UPLO is CHARACTER*1
!>          = 'U':  Upper triangle of A is stored;
!>          = 'L':  Lower triangle of A is stored.
!>
[in]N
!>          N is INTEGER
!>          The number of linear equations, i.e., the order of the
!>          matrix A.  N >= 0.
!>
[in]KD
!>          KD is INTEGER
!>          The number of superdiagonals of the matrix A if UPLO = 'U',
!>          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
!>
[in]NRHS
!>          NRHS is INTEGER
!>          The number of right-hand sides, i.e., the number of columns
!>          of the matrices B and X.  NRHS >= 0.
!>
[in,out]AB
!>          AB is COMPLEX*16 array, dimension (LDAB,N)
!>          On entry, the upper or lower triangle of the Hermitian band
!>          matrix A, stored in the first KD+1 rows of the array, except
!>          if FACT = 'F' and EQUED = 'Y', then A must contain the
!>          equilibrated matrix diag(S)*A*diag(S).  The j-th column of A
!>          is stored in the j-th column of the array AB as follows:
!>          if UPLO = 'U', AB(KD+1+i-j,j) = A(i,j) for max(1,j-KD)<=i<=j;
!>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(N,j+KD).
!>          See below for further details.
!>
!>          On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by
!>          diag(S)*A*diag(S).
!>
[in]LDAB
!>          LDAB is INTEGER
!>          The leading dimension of the array A.  LDAB >= KD+1.
!>
[in,out]AFB
!>          AFB is COMPLEX*16 array, dimension (LDAFB,N)
!>          If FACT = 'F', then AFB is an input argument and on entry
!>          contains the triangular factor U or L from the Cholesky
!>          factorization A = U**H *U or A = L*L**H of the band matrix
!>          A, in the same storage format as A (see AB).  If EQUED = 'Y',
!>          then AFB is the factored form of the equilibrated matrix A.
!>
!>          If FACT = 'N', then AFB is an output argument and on exit
!>          returns the triangular factor U or L from the Cholesky
!>          factorization A = U**H *U or A = L*L**H.
!>
!>          If FACT = 'E', then AFB is an output argument and on exit
!>          returns the triangular factor U or L from the Cholesky
!>          factorization A = U**H *U or A = L*L**H of the equilibrated
!>          matrix A (see the description of A for the form of the
!>          equilibrated matrix).
!>
[in]LDAFB
!>          LDAFB is INTEGER
!>          The leading dimension of the array AFB.  LDAFB >= KD+1.
!>
[in,out]EQUED
!>          EQUED is CHARACTER*1
!>          Specifies the form of equilibration that was done.
!>          = 'N':  No equilibration (always true if FACT = 'N').
!>          = 'Y':  Equilibration was done, i.e., A has been replaced by
!>                  diag(S) * A * diag(S).
!>          EQUED is an input argument if FACT = 'F'; otherwise, it is an
!>          output argument.
!>
[in,out]S
!>          S is DOUBLE PRECISION array, dimension (N)
!>          The scale factors for A; not accessed if EQUED = 'N'.  S is
!>          an input argument if FACT = 'F'; otherwise, S is an output
!>          argument.  If FACT = 'F' and EQUED = 'Y', each element of S
!>          must be positive.
!>
[in,out]B
!>          B is COMPLEX*16 array, dimension (LDB,NRHS)
!>          On entry, the N-by-NRHS right hand side matrix B.
!>          On exit, if EQUED = 'N', B is not modified; if EQUED = 'Y',
!>          B is overwritten by diag(S) * B.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= max(1,N).
!>
[out]X
!>          X is COMPLEX*16 array, dimension (LDX,NRHS)
!>          If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to
!>          the original system of equations.  Note that if EQUED = 'Y',
!>          A and B are modified on exit, and the solution to the
!>          equilibrated system is inv(diag(S))*X.
!>
[in]LDX
!>          LDX is INTEGER
!>          The leading dimension of the array X.  LDX >= max(1,N).
!>
[out]RCOND
!>          RCOND is DOUBLE PRECISION
!>          The estimate of the reciprocal condition number of the matrix
!>          A after equilibration (if done).  If RCOND is less than the
!>          machine precision (in particular, if RCOND = 0), the matrix
!>          is singular to working precision.  This condition is
!>          indicated by a return code of INFO > 0.
!>
[out]FERR
!>          FERR is DOUBLE PRECISION array, dimension (NRHS)
!>          The estimated forward error bound for each solution vector
!>          X(j) (the j-th column of the solution matrix X).
!>          If XTRUE is the true solution corresponding to X(j), FERR(j)
!>          is an estimated upper bound for the magnitude of the largest
!>          element in (X(j) - XTRUE) divided by the magnitude of the
!>          largest element in X(j).  The estimate is as reliable as
!>          the estimate for RCOND, and is almost always a slight
!>          overestimate of the true error.
!>
[out]BERR
!>          BERR is DOUBLE PRECISION array, dimension (NRHS)
!>          The componentwise relative backward error of each solution
!>          vector X(j) (i.e., the smallest relative change in
!>          any element of A or B that makes X(j) an exact solution).
!>
[out]WORK
!>          WORK is COMPLEX*16 array, dimension (2*N)
!>
[out]RWORK
!>          RWORK is DOUBLE PRECISION array, dimension (N)
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0: successful exit
!>          < 0: if INFO = -i, the i-th argument had an illegal value
!>          > 0: if INFO = i, and i is
!>                <= N:  the leading minor of order i of A is
!>                       not positive definite, so the factorization
!>                       could not be completed, and the solution has not
!>                       been computed. RCOND = 0 is returned.
!>                = N+1: U is nonsingular, but RCOND is less than machine
!>                       precision, meaning that the matrix is singular
!>                       to working precision.  Nevertheless, the
!>                       solution and error bounds are computed because
!>                       there are a number of situations where the
!>                       computed solution can be more accurate than the
!>                       value of RCOND would suggest.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  The band storage scheme is illustrated by the following example, when
!>  N = 6, KD = 2, and UPLO = 'U':
!>
!>  Two-dimensional storage of the Hermitian matrix A:
!>
!>     a11  a12  a13
!>          a22  a23  a24
!>               a33  a34  a35
!>                    a44  a45  a46
!>                         a55  a56
!>     (aij=conjg(aji))         a66
!>
!>  Band storage of the upper triangle of A:
!>
!>      *    *   a13  a24  a35  a46
!>      *   a12  a23  a34  a45  a56
!>     a11  a22  a33  a44  a55  a66
!>
!>  Similarly, if UPLO = 'L' the format of A is as follows:
!>
!>     a11  a22  a33  a44  a55  a66
!>     a21  a32  a43  a54  a65   *
!>     a31  a42  a53  a64   *    *
!>
!>  Array elements marked * are not used by the routine.
!>

Definition at line 339 of file zpbsvx.f.

342*
343* -- LAPACK driver routine --
344* -- LAPACK is a software package provided by Univ. of Tennessee, --
345* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
346*
347* .. Scalar Arguments ..
348 CHARACTER EQUED, FACT, UPLO
349 INTEGER INFO, KD, LDAB, LDAFB, LDB, LDX, N, NRHS
350 DOUBLE PRECISION RCOND
351* ..
352* .. Array Arguments ..
353 DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * ), S( * )
354 COMPLEX*16 AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ),
355 $ WORK( * ), X( LDX, * )
356* ..
357*
358* =====================================================================
359*
360* .. Parameters ..
361 DOUBLE PRECISION ZERO, ONE
362 parameter( zero = 0.0d+0, one = 1.0d+0 )
363* ..
364* .. Local Scalars ..
365 LOGICAL EQUIL, NOFACT, RCEQU, UPPER
366 INTEGER I, INFEQU, J, J1, J2
367 DOUBLE PRECISION AMAX, ANORM, BIGNUM, SCOND, SMAX, SMIN, SMLNUM
368* ..
369* .. External Functions ..
370 LOGICAL LSAME
371 DOUBLE PRECISION DLAMCH, ZLANHB
372 EXTERNAL lsame, dlamch, zlanhb
373* ..
374* .. External Subroutines ..
375 EXTERNAL xerbla, zcopy, zlacpy, zlaqhb, zpbcon, zpbequ,
376 $ zpbrfs, zpbtrf, zpbtrs
377* ..
378* .. Intrinsic Functions ..
379 INTRINSIC max, min
380* ..
381* .. Executable Statements ..
382*
383 info = 0
384 nofact = lsame( fact, 'N' )
385 equil = lsame( fact, 'E' )
386 upper = lsame( uplo, 'U' )
387 IF( nofact .OR. equil ) THEN
388 equed = 'N'
389 rcequ = .false.
390 ELSE
391 rcequ = lsame( equed, 'Y' )
392 smlnum = dlamch( 'Safe minimum' )
393 bignum = one / smlnum
394 END IF
395*
396* Test the input parameters.
397*
398 IF( .NOT.nofact .AND. .NOT.equil .AND. .NOT.lsame( fact, 'F' ) )
399 $ THEN
400 info = -1
401 ELSE IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
402 info = -2
403 ELSE IF( n.LT.0 ) THEN
404 info = -3
405 ELSE IF( kd.LT.0 ) THEN
406 info = -4
407 ELSE IF( nrhs.LT.0 ) THEN
408 info = -5
409 ELSE IF( ldab.LT.kd+1 ) THEN
410 info = -7
411 ELSE IF( ldafb.LT.kd+1 ) THEN
412 info = -9
413 ELSE IF( lsame( fact, 'F' ) .AND. .NOT.
414 $ ( rcequ .OR. lsame( equed, 'N' ) ) ) THEN
415 info = -10
416 ELSE
417 IF( rcequ ) THEN
418 smin = bignum
419 smax = zero
420 DO 10 j = 1, n
421 smin = min( smin, s( j ) )
422 smax = max( smax, s( j ) )
423 10 CONTINUE
424 IF( smin.LE.zero ) THEN
425 info = -11
426 ELSE IF( n.GT.0 ) THEN
427 scond = max( smin, smlnum ) / min( smax, bignum )
428 ELSE
429 scond = one
430 END IF
431 END IF
432 IF( info.EQ.0 ) THEN
433 IF( ldb.LT.max( 1, n ) ) THEN
434 info = -13
435 ELSE IF( ldx.LT.max( 1, n ) ) THEN
436 info = -15
437 END IF
438 END IF
439 END IF
440*
441 IF( info.NE.0 ) THEN
442 CALL xerbla( 'ZPBSVX', -info )
443 RETURN
444 END IF
445*
446 IF( equil ) THEN
447*
448* Compute row and column scalings to equilibrate the matrix A.
449*
450 CALL zpbequ( uplo, n, kd, ab, ldab, s, scond, amax, infequ )
451 IF( infequ.EQ.0 ) THEN
452*
453* Equilibrate the matrix.
454*
455 CALL zlaqhb( uplo, n, kd, ab, ldab, s, scond, amax, equed )
456 rcequ = lsame( equed, 'Y' )
457 END IF
458 END IF
459*
460* Scale the right-hand side.
461*
462 IF( rcequ ) THEN
463 DO 30 j = 1, nrhs
464 DO 20 i = 1, n
465 b( i, j ) = s( i )*b( i, j )
466 20 CONTINUE
467 30 CONTINUE
468 END IF
469*
470 IF( nofact .OR. equil ) THEN
471*
472* Compute the Cholesky factorization A = U**H *U or A = L*L**H.
473*
474 IF( upper ) THEN
475 DO 40 j = 1, n
476 j1 = max( j-kd, 1 )
477 CALL zcopy( j-j1+1, ab( kd+1-j+j1, j ), 1,
478 $ afb( kd+1-j+j1, j ), 1 )
479 40 CONTINUE
480 ELSE
481 DO 50 j = 1, n
482 j2 = min( j+kd, n )
483 CALL zcopy( j2-j+1, ab( 1, j ), 1, afb( 1, j ), 1 )
484 50 CONTINUE
485 END IF
486*
487 CALL zpbtrf( uplo, n, kd, afb, ldafb, info )
488*
489* Return if INFO is non-zero.
490*
491 IF( info.GT.0 )THEN
492 rcond = zero
493 RETURN
494 END IF
495 END IF
496*
497* Compute the norm of the matrix A.
498*
499 anorm = zlanhb( '1', uplo, n, kd, ab, ldab, rwork )
500*
501* Compute the reciprocal of the condition number of A.
502*
503 CALL zpbcon( uplo, n, kd, afb, ldafb, anorm, rcond, work, rwork,
504 $ info )
505*
506* Compute the solution matrix X.
507*
508 CALL zlacpy( 'Full', n, nrhs, b, ldb, x, ldx )
509 CALL zpbtrs( uplo, n, kd, nrhs, afb, ldafb, x, ldx, info )
510*
511* Use iterative refinement to improve the computed solution and
512* compute error bounds and backward error estimates for it.
513*
514 CALL zpbrfs( uplo, n, kd, nrhs, ab, ldab, afb, ldafb, b, ldb, x,
515 $ ldx, ferr, berr, work, rwork, info )
516*
517* Transform the solution matrix X to a solution of the original
518* system.
519*
520 IF( rcequ ) THEN
521 DO 70 j = 1, nrhs
522 DO 60 i = 1, n
523 x( i, j ) = s( i )*x( i, j )
524 60 CONTINUE
525 70 CONTINUE
526 DO 80 j = 1, nrhs
527 ferr( j ) = ferr( j ) / scond
528 80 CONTINUE
529 END IF
530*
531* Set INFO = N+1 if the matrix is singular to working precision.
532*
533 IF( rcond.LT.dlamch( 'Epsilon' ) )
534 $ info = n + 1
535*
536 RETURN
537*
538* End of ZPBSVX
539*
zlanhb
double precision function zlanhb(norm, uplo, n, k, ab, ldab, work)
ZLANHB returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition zlanhb.f:132
zlaqhb
subroutine zlaqhb(uplo, n, kd, ab, ldab, s, scond, amax, equed)
ZLAQHB scales a Hermitian band matrix, using scaling factors computed by cpbequ.
Definition zlaqhb.f:141
zpbcon
subroutine zpbcon(uplo, n, kd, ab, ldab, anorm, rcond, work, rwork, info)
ZPBCON
Definition zpbcon.f:133
zpbequ
subroutine zpbequ(uplo, n, kd, ab, ldab, s, scond, amax, info)
ZPBEQU
Definition zpbequ.f:130
zpbrfs
subroutine zpbrfs(uplo, n, kd, nrhs, ab, ldab, afb, ldafb, b, ldb, x, ldx, ferr, berr, work, rwork, info)
ZPBRFS
Definition zpbrfs.f:189

◆ zppsv()

subroutine zppsv ( character uplo,
integer n,
integer nrhs,
complex*16, dimension( * ) ap,
complex*16, dimension( ldb, * ) b,
integer ldb,
integer info )

ZPPSV computes the solution to system of linear equations A * X = B for OTHER matrices

Download ZPPSV + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> ZPPSV computes the solution to a complex system of linear equations
!>    A * X = B,
!> where A is an N-by-N Hermitian positive definite matrix stored in
!> packed format and X and B are N-by-NRHS matrices.
!>
!> The Cholesky decomposition is used to factor A as
!>    A = U**H * U,  if UPLO = 'U', or
!>    A = L * L**H,  if UPLO = 'L',
!> where U is an upper triangular matrix and L is a lower triangular
!> matrix.  The factored form of A is then used to solve the system of
!> equations A * X = B.
!>
Parameters
[in]UPLO
!>          UPLO is CHARACTER*1
!>          = 'U':  Upper triangle of A is stored;
!>          = 'L':  Lower triangle of A is stored.
!>
[in]N
!>          N is INTEGER
!>          The number of linear equations, i.e., the order of the
!>          matrix A.  N >= 0.
!>
[in]NRHS
!>          NRHS is INTEGER
!>          The number of right hand sides, i.e., the number of columns
!>          of the matrix B.  NRHS >= 0.
!>
[in,out]AP
!>          AP is COMPLEX*16 array, dimension (N*(N+1)/2)
!>          On entry, the upper or lower triangle of the Hermitian matrix
!>          A, packed columnwise in a linear array.  The j-th column of A
!>          is stored in the array AP as follows:
!>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
!>          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
!>          See below for further details.
!>
!>          On exit, if INFO = 0, the factor U or L from the Cholesky
!>          factorization A = U**H*U or A = L*L**H, in the same storage
!>          format as A.
!>
[in,out]B
!>          B is COMPLEX*16 array, dimension (LDB,NRHS)
!>          On entry, the N-by-NRHS right hand side matrix B.
!>          On exit, if INFO = 0, the N-by-NRHS solution matrix X.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= max(1,N).
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO = i, the leading minor of order i of A is not
!>                positive definite, so the factorization could not be
!>                completed, and the solution has not been computed.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  The packed storage scheme is illustrated by the following example
!>  when N = 4, UPLO = 'U':
!>
!>  Two-dimensional storage of the Hermitian matrix A:
!>
!>     a11 a12 a13 a14
!>         a22 a23 a24
!>             a33 a34     (aij = conjg(aji))
!>                 a44
!>
!>  Packed storage of the upper triangle of A:
!>
!>  AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]
!>

Definition at line 143 of file zppsv.f.

144*
145* -- LAPACK driver routine --
146* -- LAPACK is a software package provided by Univ. of Tennessee, --
147* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
148*
149* .. Scalar Arguments ..
150 CHARACTER UPLO
151 INTEGER INFO, LDB, N, NRHS
152* ..
153* .. Array Arguments ..
154 COMPLEX*16 AP( * ), B( LDB, * )
155* ..
156*
157* =====================================================================
158*
159* .. External Functions ..
160 LOGICAL LSAME
161 EXTERNAL lsame
162* ..
163* .. External Subroutines ..
164 EXTERNAL xerbla, zpptrf, zpptrs
165* ..
166* .. Intrinsic Functions ..
167 INTRINSIC max
168* ..
169* .. Executable Statements ..
170*
171* Test the input parameters.
172*
173 info = 0
174 IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lsame( uplo, 'L' ) ) THEN
175 info = -1
176 ELSE IF( n.LT.0 ) THEN
177 info = -2
178 ELSE IF( nrhs.LT.0 ) THEN
179 info = -3
180 ELSE IF( ldb.LT.max( 1, n ) ) THEN
181 info = -6
182 END IF
183 IF( info.NE.0 ) THEN
184 CALL xerbla( 'ZPPSV ', -info )
185 RETURN
186 END IF
187*
188* Compute the Cholesky factorization A = U**H *U or A = L*L**H.
189*
190 CALL zpptrf( uplo, n, ap, info )
191 IF( info.EQ.0 ) THEN
192*
193* Solve the system A*X = B, overwriting B with X.
194*
195 CALL zpptrs( uplo, n, nrhs, ap, b, ldb, info )
196*
197 END IF
198 RETURN
199*
200* End of ZPPSV
201*
zpptrs
subroutine zpptrs(uplo, n, nrhs, ap, b, ldb, info)
ZPPTRS
Definition zpptrs.f:108
zpptrf
subroutine zpptrf(uplo, n, ap, info)
ZPPTRF
Definition zpptrf.f:119

◆ zppsvx()

subroutine zppsvx ( character fact,
character uplo,
integer n,
integer nrhs,
complex*16, dimension( * ) ap,
complex*16, dimension( * ) afp,
character equed,
double precision, dimension( * ) s,
complex*16, dimension( ldb, * ) b,
integer ldb,
complex*16, dimension( ldx, * ) x,
integer ldx,
double precision rcond,
double precision, dimension( * ) ferr,
double precision, dimension( * ) berr,
complex*16, dimension( * ) work,
double precision, dimension( * ) rwork,
integer info )

ZPPSVX computes the solution to system of linear equations A * X = B for OTHER matrices

Download ZPPSVX + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> ZPPSVX uses the Cholesky factorization A = U**H * U or A = L * L**H to
!> compute the solution to a complex system of linear equations
!>    A * X = B,
!> where A is an N-by-N Hermitian positive definite matrix stored in
!> packed format and X and B are N-by-NRHS matrices.
!>
!> Error bounds on the solution and a condition estimate are also
!> provided.
!>
Description:
!>
!> The following steps are performed:
!>
!> 1. If FACT = 'E', real scaling factors are computed to equilibrate
!>    the system:
!>       diag(S) * A * diag(S) * inv(diag(S)) * X = diag(S) * B
!>    Whether or not the system will be equilibrated depends on the
!>    scaling of the matrix A, but if equilibration is used, A is
!>    overwritten by diag(S)*A*diag(S) and B by diag(S)*B.
!>
!> 2. If FACT = 'N' or 'E', the Cholesky decomposition is used to
!>    factor the matrix A (after equilibration if FACT = 'E') as
!>       A = U**H * U ,  if UPLO = 'U', or
!>       A = L * L**H,  if UPLO = 'L',
!>    where U is an upper triangular matrix, L is a lower triangular
!>    matrix, and **H indicates conjugate transpose.
!>
!> 3. If the leading i-by-i principal minor is not positive definite,
!>    then the routine returns with INFO = i. Otherwise, the factored
!>    form of A is used to estimate the condition number of the matrix
!>    A.  If the reciprocal of the condition number is less than machine
!>    precision, INFO = N+1 is returned as a warning, but the routine
!>    still goes on to solve for X and compute error bounds as
!>    described below.
!>
!> 4. The system of equations is solved for X using the factored form
!>    of A.
!>
!> 5. Iterative refinement is applied to improve the computed solution
!>    matrix and calculate error bounds and backward error estimates
!>    for it.
!>
!> 6. If equilibration was used, the matrix X is premultiplied by
!>    diag(S) so that it solves the original system before
!>    equilibration.
!>
Parameters
[in]FACT
!>          FACT is CHARACTER*1
!>          Specifies whether or not the factored form of the matrix A is
!>          supplied on entry, and if not, whether the matrix A should be
!>          equilibrated before it is factored.
!>          = 'F':  On entry, AFP contains the factored form of A.
!>                  If EQUED = 'Y', the matrix A has been equilibrated
!>                  with scaling factors given by S.  AP and AFP will not
!>                  be modified.
!>          = 'N':  The matrix A will be copied to AFP and factored.
!>          = 'E':  The matrix A will be equilibrated if necessary, then
!>                  copied to AFP and factored.
!>
[in]UPLO
!>          UPLO is CHARACTER*1
!>          = 'U':  Upper triangle of A is stored;
!>          = 'L':  Lower triangle of A is stored.
!>
[in]N
!>          N is INTEGER
!>          The number of linear equations, i.e., the order of the
!>          matrix A.  N >= 0.
!>
[in]NRHS
!>          NRHS is INTEGER
!>          The number of right hand sides, i.e., the number of columns
!>          of the matrices B and X.  NRHS >= 0.
!>
[in,out]AP
!>          AP is COMPLEX*16 array, dimension (N*(N+1)/2)
!>          On entry, the upper or lower triangle of the Hermitian matrix
!>          A, packed columnwise in a linear array, except if FACT = 'F'
!>          and EQUED = 'Y', then A must contain the equilibrated matrix
!>          diag(S)*A*diag(S).  The j-th column of A is stored in the
!>          array AP as follows:
!>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
!>          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
!>          See below for further details.  A is not modified if
!>          FACT = 'F' or 'N', or if FACT = 'E' and EQUED = 'N' on exit.
!>
!>          On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by
!>          diag(S)*A*diag(S).
!>
[in,out]AFP
!>          AFP is COMPLEX*16 array, dimension (N*(N+1)/2)
!>          If FACT = 'F', then AFP is an input argument and on entry
!>          contains the triangular factor U or L from the Cholesky
!>          factorization A = U**H*U or A = L*L**H, in the same storage
!>          format as A.  If EQUED .ne. 'N', then AFP is the factored
!>          form of the equilibrated matrix A.
!>
!>          If FACT = 'N', then AFP is an output argument and on exit
!>          returns the triangular factor U or L from the Cholesky
!>          factorization A = U**H * U or A = L * L**H of the original
!>          matrix A.
!>
!>          If FACT = 'E', then AFP is an output argument and on exit
!>          returns the triangular factor U or L from the Cholesky
!>          factorization A = U**H * U or A = L * L**H of the equilibrated
!>          matrix A (see the description of AP for the form of the
!>          equilibrated matrix).
!>
[in,out]EQUED
!>          EQUED is CHARACTER*1
!>          Specifies the form of equilibration that was done.
!>          = 'N':  No equilibration (always true if FACT = 'N').
!>          = 'Y':  Equilibration was done, i.e., A has been replaced by
!>                  diag(S) * A * diag(S).
!>          EQUED is an input argument if FACT = 'F'; otherwise, it is an
!>          output argument.
!>
[in,out]S
!>          S is DOUBLE PRECISION array, dimension (N)
!>          The scale factors for A; not accessed if EQUED = 'N'.  S is
!>          an input argument if FACT = 'F'; otherwise, S is an output
!>          argument.  If FACT = 'F' and EQUED = 'Y', each element of S
!>          must be positive.
!>
[in,out]B
!>          B is COMPLEX*16 array, dimension (LDB,NRHS)
!>          On entry, the N-by-NRHS right hand side matrix B.
!>          On exit, if EQUED = 'N', B is not modified; if EQUED = 'Y',
!>          B is overwritten by diag(S) * B.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= max(1,N).
!>
[out]X
!>          X is COMPLEX*16 array, dimension (LDX,NRHS)
!>          If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to
!>          the original system of equations.  Note that if EQUED = 'Y',
!>          A and B are modified on exit, and the solution to the
!>          equilibrated system is inv(diag(S))*X.
!>
[in]LDX
!>          LDX is INTEGER
!>          The leading dimension of the array X.  LDX >= max(1,N).
!>
[out]RCOND
!>          RCOND is DOUBLE PRECISION
!>          The estimate of the reciprocal condition number of the matrix
!>          A after equilibration (if done).  If RCOND is less than the
!>          machine precision (in particular, if RCOND = 0), the matrix
!>          is singular to working precision.  This condition is
!>          indicated by a return code of INFO > 0.
!>
[out]FERR
!>          FERR is DOUBLE PRECISION array, dimension (NRHS)
!>          The estimated forward error bound for each solution vector
!>          X(j) (the j-th column of the solution matrix X).
!>          If XTRUE is the true solution corresponding to X(j), FERR(j)
!>          is an estimated upper bound for the magnitude of the largest
!>          element in (X(j) - XTRUE) divided by the magnitude of the
!>          largest element in X(j).  The estimate is as reliable as
!>          the estimate for RCOND, and is almost always a slight
!>          overestimate of the true error.
!>
[out]BERR
!>          BERR is DOUBLE PRECISION array, dimension (NRHS)
!>          The componentwise relative backward error of each solution
!>          vector X(j) (i.e., the smallest relative change in
!>          any element of A or B that makes X(j) an exact solution).
!>
[out]WORK
!>          WORK is COMPLEX*16 array, dimension (2*N)
!>
[out]RWORK
!>          RWORK is DOUBLE PRECISION array, dimension (N)
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO = i, and i is
!>                <= N:  the leading minor of order i of A is
!>                       not positive definite, so the factorization
!>                       could not be completed, and the solution has not
!>                       been computed. RCOND = 0 is returned.
!>                = N+1: U is nonsingular, but RCOND is less than machine
!>                       precision, meaning that the matrix is singular
!>                       to working precision.  Nevertheless, the
!>                       solution and error bounds are computed because
!>                       there are a number of situations where the
!>                       computed solution can be more accurate than the
!>                       value of RCOND would suggest.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  The packed storage scheme is illustrated by the following example
!>  when N = 4, UPLO = 'U':
!>
!>  Two-dimensional storage of the Hermitian matrix A:
!>
!>     a11 a12 a13 a14
!>         a22 a23 a24
!>             a33 a34     (aij = conjg(aji))
!>                 a44
!>
!>  Packed storage of the upper triangle of A:
!>
!>  AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]
!>

Definition at line 309 of file zppsvx.f.

311*
312* -- LAPACK driver routine --
313* -- LAPACK is a software package provided by Univ. of Tennessee, --
314* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
315*
316* .. Scalar Arguments ..
317 CHARACTER EQUED, FACT, UPLO
318 INTEGER INFO, LDB, LDX, N, NRHS
319 DOUBLE PRECISION RCOND
320* ..
321* .. Array Arguments ..
322 DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * ), S( * )
323 COMPLEX*16 AFP( * ), AP( * ), B( LDB, * ), WORK( * ),
324 $ X( LDX, * )
325* ..
326*
327* =====================================================================
328*
329* .. Parameters ..
330 DOUBLE PRECISION ZERO, ONE
331 parameter( zero = 0.0d+0, one = 1.0d+0 )
332* ..
333* .. Local Scalars ..
334 LOGICAL EQUIL, NOFACT, RCEQU
335 INTEGER I, INFEQU, J
336 DOUBLE PRECISION AMAX, ANORM, BIGNUM, SCOND, SMAX, SMIN, SMLNUM
337* ..
338* .. External Functions ..
339 LOGICAL LSAME
340 DOUBLE PRECISION DLAMCH, ZLANHP
341 EXTERNAL lsame, dlamch, zlanhp
342* ..
343* .. External Subroutines ..
344 EXTERNAL xerbla, zcopy, zlacpy, zlaqhp, zppcon, zppequ,
345 $ zpprfs, zpptrf, zpptrs
346* ..
347* .. Intrinsic Functions ..
348 INTRINSIC max, min
349* ..
350* .. Executable Statements ..
351*
352 info = 0
353 nofact = lsame( fact, 'N' )
354 equil = lsame( fact, 'E' )
355 IF( nofact .OR. equil ) THEN
356 equed = 'N'
357 rcequ = .false.
358 ELSE
359 rcequ = lsame( equed, 'Y' )
360 smlnum = dlamch( 'Safe minimum' )
361 bignum = one / smlnum
362 END IF
363*
364* Test the input parameters.
365*
366 IF( .NOT.nofact .AND. .NOT.equil .AND. .NOT.lsame( fact, 'F' ) )
367 $ THEN
368 info = -1
369 ELSE IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lsame( uplo, 'L' ) )
370 $ THEN
371 info = -2
372 ELSE IF( n.LT.0 ) THEN
373 info = -3
374 ELSE IF( nrhs.LT.0 ) THEN
375 info = -4
376 ELSE IF( lsame( fact, 'F' ) .AND. .NOT.
377 $ ( rcequ .OR. lsame( equed, 'N' ) ) ) THEN
378 info = -7
379 ELSE
380 IF( rcequ ) THEN
381 smin = bignum
382 smax = zero
383 DO 10 j = 1, n
384 smin = min( smin, s( j ) )
385 smax = max( smax, s( j ) )
386 10 CONTINUE
387 IF( smin.LE.zero ) THEN
388 info = -8
389 ELSE IF( n.GT.0 ) THEN
390 scond = max( smin, smlnum ) / min( smax, bignum )
391 ELSE
392 scond = one
393 END IF
394 END IF
395 IF( info.EQ.0 ) THEN
396 IF( ldb.LT.max( 1, n ) ) THEN
397 info = -10
398 ELSE IF( ldx.LT.max( 1, n ) ) THEN
399 info = -12
400 END IF
401 END IF
402 END IF
403*
404 IF( info.NE.0 ) THEN
405 CALL xerbla( 'ZPPSVX', -info )
406 RETURN
407 END IF
408*
409 IF( equil ) THEN
410*
411* Compute row and column scalings to equilibrate the matrix A.
412*
413 CALL zppequ( uplo, n, ap, s, scond, amax, infequ )
414 IF( infequ.EQ.0 ) THEN
415*
416* Equilibrate the matrix.
417*
418 CALL zlaqhp( uplo, n, ap, s, scond, amax, equed )
419 rcequ = lsame( equed, 'Y' )
420 END IF
421 END IF
422*
423* Scale the right-hand side.
424*
425 IF( rcequ ) THEN
426 DO 30 j = 1, nrhs
427 DO 20 i = 1, n
428 b( i, j ) = s( i )*b( i, j )
429 20 CONTINUE
430 30 CONTINUE
431 END IF
432*
433 IF( nofact .OR. equil ) THEN
434*
435* Compute the Cholesky factorization A = U**H * U or A = L * L**H.
436*
437 CALL zcopy( n*( n+1 ) / 2, ap, 1, afp, 1 )
438 CALL zpptrf( uplo, n, afp, info )
439*
440* Return if INFO is non-zero.
441*
442 IF( info.GT.0 )THEN
443 rcond = zero
444 RETURN
445 END IF
446 END IF
447*
448* Compute the norm of the matrix A.
449*
450 anorm = zlanhp( 'I', uplo, n, ap, rwork )
451*
452* Compute the reciprocal of the condition number of A.
453*
454 CALL zppcon( uplo, n, afp, anorm, rcond, work, rwork, info )
455*
456* Compute the solution matrix X.
457*
458 CALL zlacpy( 'Full', n, nrhs, b, ldb, x, ldx )
459 CALL zpptrs( uplo, n, nrhs, afp, x, ldx, info )
460*
461* Use iterative refinement to improve the computed solution and
462* compute error bounds and backward error estimates for it.
463*
464 CALL zpprfs( uplo, n, nrhs, ap, afp, b, ldb, x, ldx, ferr, berr,
465 $ work, rwork, info )
466*
467* Transform the solution matrix X to a solution of the original
468* system.
469*
470 IF( rcequ ) THEN
471 DO 50 j = 1, nrhs
472 DO 40 i = 1, n
473 x( i, j ) = s( i )*x( i, j )
474 40 CONTINUE
475 50 CONTINUE
476 DO 60 j = 1, nrhs
477 ferr( j ) = ferr( j ) / scond
478 60 CONTINUE
479 END IF
480*
481* Set INFO = N+1 if the matrix is singular to working precision.
482*
483 IF( rcond.LT.dlamch( 'Epsilon' ) )
484 $ info = n + 1
485*
486 RETURN
487*
488* End of ZPPSVX
489*
zlaqhp
subroutine zlaqhp(uplo, n, ap, s, scond, amax, equed)
ZLAQHP scales a Hermitian matrix stored in packed form.
Definition zlaqhp.f:126
zppequ
subroutine zppequ(uplo, n, ap, s, scond, amax, info)
ZPPEQU
Definition zppequ.f:117
zpprfs
subroutine zpprfs(uplo, n, nrhs, ap, afp, b, ldb, x, ldx, ferr, berr, work, rwork, info)
ZPPRFS
Definition zpprfs.f:171
zppcon
subroutine zppcon(uplo, n, ap, anorm, rcond, work, rwork, info)
ZPPCON
Definition zppcon.f:118

◆ zspsv()

subroutine zspsv ( character uplo,
integer n,
integer nrhs,
complex*16, dimension( * ) ap,
integer, dimension( * ) ipiv,
complex*16, dimension( ldb, * ) b,
integer ldb,
integer info )

ZSPSV computes the solution to system of linear equations A * X = B for OTHER matrices

Download ZSPSV + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> ZSPSV computes the solution to a complex system of linear equations
!>    A * X = B,
!> where A is an N-by-N symmetric matrix stored in packed format and X
!> and B are N-by-NRHS matrices.
!>
!> The diagonal pivoting method is used to factor A as
!>    A = U * D * U**T,  if UPLO = 'U', or
!>    A = L * D * L**T,  if UPLO = 'L',
!> where U (or L) is a product of permutation and unit upper (lower)
!> triangular matrices, D is symmetric and block diagonal with 1-by-1
!> and 2-by-2 diagonal blocks.  The factored form of A is then used to
!> solve the system of equations A * X = B.
!>
Parameters
[in]UPLO
!>          UPLO is CHARACTER*1
!>          = 'U':  Upper triangle of A is stored;
!>          = 'L':  Lower triangle of A is stored.
!>
[in]N
!>          N is INTEGER
!>          The number of linear equations, i.e., the order of the
!>          matrix A.  N >= 0.
!>
[in]NRHS
!>          NRHS is INTEGER
!>          The number of right hand sides, i.e., the number of columns
!>          of the matrix B.  NRHS >= 0.
!>
[in,out]AP
!>          AP is COMPLEX*16 array, dimension (N*(N+1)/2)
!>          On entry, the upper or lower triangle of the symmetric matrix
!>          A, packed columnwise in a linear array.  The j-th column of A
!>          is stored in the array AP as follows:
!>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
!>          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
!>          See below for further details.
!>
!>          On exit, the block diagonal matrix D and the multipliers used
!>          to obtain the factor U or L from the factorization
!>          A = U*D*U**T or A = L*D*L**T as computed by ZSPTRF, stored as
!>          a packed triangular matrix in the same storage format as A.
!>
[out]IPIV
!>          IPIV is INTEGER array, dimension (N)
!>          Details of the interchanges and the block structure of D, as
!>          determined by ZSPTRF.  If IPIV(k) > 0, then rows and columns
!>          k and IPIV(k) were interchanged, and D(k,k) is a 1-by-1
!>          diagonal block.  If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0,
!>          then rows and columns k-1 and -IPIV(k) were interchanged and
!>          D(k-1:k,k-1:k) is a 2-by-2 diagonal block.  If UPLO = 'L' and
!>          IPIV(k) = IPIV(k+1) < 0, then rows and columns k+1 and
!>          -IPIV(k) were interchanged and D(k:k+1,k:k+1) is a 2-by-2
!>          diagonal block.
!>
[in,out]B
!>          B is COMPLEX*16 array, dimension (LDB,NRHS)
!>          On entry, the N-by-NRHS right hand side matrix B.
!>          On exit, if INFO = 0, the N-by-NRHS solution matrix X.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= max(1,N).
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO = i, D(i,i) is exactly zero.  The factorization
!>                has been completed, but the block diagonal matrix D is
!>                exactly singular, so the solution could not be
!>                computed.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  The packed storage scheme is illustrated by the following example
!>  when N = 4, UPLO = 'U':
!>
!>  Two-dimensional storage of the symmetric matrix A:
!>
!>     a11 a12 a13 a14
!>         a22 a23 a24
!>             a33 a34     (aij = aji)
!>                 a44
!>
!>  Packed storage of the upper triangle of A:
!>
!>  AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]
!>

Definition at line 161 of file zspsv.f.

162*
163* -- LAPACK driver routine --
164* -- LAPACK is a software package provided by Univ. of Tennessee, --
165* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
166*
167* .. Scalar Arguments ..
168 CHARACTER UPLO
169 INTEGER INFO, LDB, N, NRHS
170* ..
171* .. Array Arguments ..
172 INTEGER IPIV( * )
173 COMPLEX*16 AP( * ), B( LDB, * )
174* ..
175*
176* =====================================================================
177*
178* .. External Functions ..
179 LOGICAL LSAME
180 EXTERNAL lsame
181* ..
182* .. External Subroutines ..
183 EXTERNAL xerbla, zsptrf, zsptrs
184* ..
185* .. Intrinsic Functions ..
186 INTRINSIC max
187* ..
188* .. Executable Statements ..
189*
190* Test the input parameters.
191*
192 info = 0
193 IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lsame( uplo, 'L' ) ) THEN
194 info = -1
195 ELSE IF( n.LT.0 ) THEN
196 info = -2
197 ELSE IF( nrhs.LT.0 ) THEN
198 info = -3
199 ELSE IF( ldb.LT.max( 1, n ) ) THEN
200 info = -7
201 END IF
202 IF( info.NE.0 ) THEN
203 CALL xerbla( 'ZSPSV ', -info )
204 RETURN
205 END IF
206*
207* Compute the factorization A = U*D*U**T or A = L*D*L**T.
208*
209 CALL zsptrf( uplo, n, ap, ipiv, info )
210 IF( info.EQ.0 ) THEN
211*
212* Solve the system A*X = B, overwriting B with X.
213*
214 CALL zsptrs( uplo, n, nrhs, ap, ipiv, b, ldb, info )
215*
216 END IF
217 RETURN
218*
219* End of ZSPSV
220*
zsptrs
subroutine zsptrs(uplo, n, nrhs, ap, ipiv, b, ldb, info)
ZSPTRS
Definition zsptrs.f:115
zsptrf
subroutine zsptrf(uplo, n, ap, ipiv, info)
ZSPTRF
Definition zsptrf.f:158

◆ zspsvx()

subroutine zspsvx ( character fact,
character uplo,
integer n,
integer nrhs,
complex*16, dimension( * ) ap,
complex*16, dimension( * ) afp,
integer, dimension( * ) ipiv,
complex*16, dimension( ldb, * ) b,
integer ldb,
complex*16, dimension( ldx, * ) x,
integer ldx,
double precision rcond,
double precision, dimension( * ) ferr,
double precision, dimension( * ) berr,
complex*16, dimension( * ) work,
double precision, dimension( * ) rwork,
integer info )

ZSPSVX computes the solution to system of linear equations A * X = B for OTHER matrices

Download ZSPSVX + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> ZSPSVX uses the diagonal pivoting factorization A = U*D*U**T or
!> A = L*D*L**T to compute the solution to a complex system of linear
!> equations A * X = B, where A is an N-by-N symmetric matrix stored
!> in packed format and X and B are N-by-NRHS matrices.
!>
!> Error bounds on the solution and a condition estimate are also
!> provided.
!>
Description:
!>
!> The following steps are performed:
!>
!> 1. If FACT = 'N', the diagonal pivoting method is used to factor A as
!>       A = U * D * U**T,  if UPLO = 'U', or
!>       A = L * D * L**T,  if UPLO = 'L',
!>    where U (or L) is a product of permutation and unit upper (lower)
!>    triangular matrices and D is symmetric and block diagonal with
!>    1-by-1 and 2-by-2 diagonal blocks.
!>
!> 2. If some D(i,i)=0, so that D is exactly singular, then the routine
!>    returns with INFO = i. Otherwise, the factored form of A is used
!>    to estimate the condition number of the matrix A.  If the
!>    reciprocal of the condition number is less than machine precision,
!>    INFO = N+1 is returned as a warning, but the routine still goes on
!>    to solve for X and compute error bounds as described below.
!>
!> 3. The system of equations is solved for X using the factored form
!>    of A.
!>
!> 4. Iterative refinement is applied to improve the computed solution
!>    matrix and calculate error bounds and backward error estimates
!>    for it.
!>
Parameters
[in]FACT
!>          FACT is CHARACTER*1
!>          Specifies whether or not the factored form of A has been
!>          supplied on entry.
!>          = 'F':  On entry, AFP and IPIV contain the factored form
!>                  of A.  AP, AFP and IPIV will not be modified.
!>          = 'N':  The matrix A will be copied to AFP and factored.
!>
[in]UPLO
!>          UPLO is CHARACTER*1
!>          = 'U':  Upper triangle of A is stored;
!>          = 'L':  Lower triangle of A is stored.
!>
[in]N
!>          N is INTEGER
!>          The number of linear equations, i.e., the order of the
!>          matrix A.  N >= 0.
!>
[in]NRHS
!>          NRHS is INTEGER
!>          The number of right hand sides, i.e., the number of columns
!>          of the matrices B and X.  NRHS >= 0.
!>
[in]AP
!>          AP is COMPLEX*16 array, dimension (N*(N+1)/2)
!>          The upper or lower triangle of the symmetric matrix A, packed
!>          columnwise in a linear array.  The j-th column of A is stored
!>          in the array AP as follows:
!>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
!>          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
!>          See below for further details.
!>
[in,out]AFP
!>          AFP is COMPLEX*16 array, dimension (N*(N+1)/2)
!>          If FACT = 'F', then AFP is an input argument and on entry
!>          contains the block diagonal matrix D and the multipliers used
!>          to obtain the factor U or L from the factorization
!>          A = U*D*U**T or A = L*D*L**T as computed by ZSPTRF, stored as
!>          a packed triangular matrix in the same storage format as A.
!>
!>          If FACT = 'N', then AFP is an output argument and on exit
!>          contains the block diagonal matrix D and the multipliers used
!>          to obtain the factor U or L from the factorization
!>          A = U*D*U**T or A = L*D*L**T as computed by ZSPTRF, stored as
!>          a packed triangular matrix in the same storage format as A.
!>
[in,out]IPIV
!>          IPIV is INTEGER array, dimension (N)
!>          If FACT = 'F', then IPIV is an input argument and on entry
!>          contains details of the interchanges and the block structure
!>          of D, as determined by ZSPTRF.
!>          If IPIV(k) > 0, then rows and columns k and IPIV(k) were
!>          interchanged and D(k,k) is a 1-by-1 diagonal block.
!>          If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
!>          columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
!>          is a 2-by-2 diagonal block.  If UPLO = 'L' and IPIV(k) =
!>          IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
!>          interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
!>
!>          If FACT = 'N', then IPIV is an output argument and on exit
!>          contains details of the interchanges and the block structure
!>          of D, as determined by ZSPTRF.
!>
[in]B
!>          B is COMPLEX*16 array, dimension (LDB,NRHS)
!>          The N-by-NRHS right hand side matrix B.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= max(1,N).
!>
[out]X
!>          X is COMPLEX*16 array, dimension (LDX,NRHS)
!>          If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.
!>
[in]LDX
!>          LDX is INTEGER
!>          The leading dimension of the array X.  LDX >= max(1,N).
!>
[out]RCOND
!>          RCOND is DOUBLE PRECISION
!>          The estimate of the reciprocal condition number of the matrix
!>          A.  If RCOND is less than the machine precision (in
!>          particular, if RCOND = 0), the matrix is singular to working
!>          precision.  This condition is indicated by a return code of
!>          INFO > 0.
!>
[out]FERR
!>          FERR is DOUBLE PRECISION array, dimension (NRHS)
!>          The estimated forward error bound for each solution vector
!>          X(j) (the j-th column of the solution matrix X).
!>          If XTRUE is the true solution corresponding to X(j), FERR(j)
!>          is an estimated upper bound for the magnitude of the largest
!>          element in (X(j) - XTRUE) divided by the magnitude of the
!>          largest element in X(j).  The estimate is as reliable as
!>          the estimate for RCOND, and is almost always a slight
!>          overestimate of the true error.
!>
[out]BERR
!>          BERR is DOUBLE PRECISION array, dimension (NRHS)
!>          The componentwise relative backward error of each solution
!>          vector X(j) (i.e., the smallest relative change in
!>          any element of A or B that makes X(j) an exact solution).
!>
[out]WORK
!>          WORK is COMPLEX*16 array, dimension (2*N)
!>
[out]RWORK
!>          RWORK is DOUBLE PRECISION array, dimension (N)
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0: successful exit
!>          < 0: if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO = i, and i is
!>                <= N:  D(i,i) is exactly zero.  The factorization
!>                       has been completed but the factor D is exactly
!>                       singular, so the solution and error bounds could
!>                       not be computed. RCOND = 0 is returned.
!>                = N+1: D is nonsingular, but RCOND is less than machine
!>                       precision, meaning that the matrix is singular
!>                       to working precision.  Nevertheless, the
!>                       solution and error bounds are computed because
!>                       there are a number of situations where the
!>                       computed solution can be more accurate than the
!>                       value of RCOND would suggest.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  The packed storage scheme is illustrated by the following example
!>  when N = 4, UPLO = 'U':
!>
!>  Two-dimensional storage of the symmetric matrix A:
!>
!>     a11 a12 a13 a14
!>         a22 a23 a24
!>             a33 a34     (aij = aji)
!>                 a44
!>
!>  Packed storage of the upper triangle of A:
!>
!>  AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]
!>

Definition at line 275 of file zspsvx.f.

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*
zlansp
double precision function zlansp(norm, uplo, n, ap, work)
ZLANSP returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition zlansp.f:115
zsprfs
subroutine zsprfs(uplo, n, nrhs, ap, afp, ipiv, b, ldb, x, ldx, ferr, berr, work, rwork, info)
ZSPRFS
Definition zsprfs.f:180
zspcon
subroutine zspcon(uplo, n, ap, ipiv, anorm, rcond, work, info)
ZSPCON
Definition zspcon.f:118