LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ zptsvx()

subroutine zptsvx ( character fact,
integer n,
integer nrhs,
double precision, dimension( * ) d,
complex*16, dimension( * ) e,
double precision, dimension( * ) df,
complex*16, dimension( * ) ef,
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 )

ZPTSVX computes the solution to system of linear equations A * X = B for PT matrices

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

Purpose:
!>
!> ZPTSVX uses the factorization 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 positive definite tridiagonal 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 = 'N', the matrix A is factored as A = L*D*L**H, where L
!>    is a unit lower bidiagonal matrix and D is diagonal.  The
!>    factorization can also be regarded as having the form
!>    A = U**H*D*U.
!>
!> 2. If the leading principal minor of order i is not positive,
!>    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 the matrix
!>          A is supplied on entry.
!>          = 'F':  On entry, DF and EF contain the factored form of A.
!>                  D, E, DF, and EF will not be modified.
!>          = 'N':  The matrix A will be copied to DF and EF and
!>                  factored.
!> 
[in]N
!>          N is INTEGER
!>          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]D
!>          D is DOUBLE PRECISION array, dimension (N)
!>          The n diagonal elements of the tridiagonal matrix A.
!> 
[in]E
!>          E is COMPLEX*16 array, dimension (N-1)
!>          The (n-1) subdiagonal elements of the tridiagonal matrix A.
!> 
[in,out]DF
!>          DF is DOUBLE PRECISION array, dimension (N)
!>          If FACT = 'F', then DF is an input argument and on entry
!>          contains the n diagonal elements of the diagonal matrix D
!>          from the L*D*L**H factorization of A.
!>          If FACT = 'N', then DF is an output argument and on exit
!>          contains the n diagonal elements of the diagonal matrix D
!>          from the L*D*L**H factorization of A.
!> 
[in,out]EF
!>          EF is COMPLEX*16 array, dimension (N-1)
!>          If FACT = 'F', then EF is an input argument and on entry
!>          contains the (n-1) subdiagonal elements of the unit
!>          bidiagonal factor L from the L*D*L**H factorization of A.
!>          If FACT = 'N', then EF is an output argument and on exit
!>          contains the (n-1) subdiagonal elements of the unit
!>          bidiagonal factor L from the L*D*L**H factorization of A.
!> 
[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 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 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).
!> 
[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 (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 principal minor of order i of A
!>                       is not positive, 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.

Definition at line 230 of file zptsvx.f.

232*
233* -- LAPACK driver routine --
234* -- LAPACK is a software package provided by Univ. of Tennessee, --
235* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
236*
237* .. Scalar Arguments ..
238 CHARACTER FACT
239 INTEGER INFO, LDB, LDX, N, NRHS
240 DOUBLE PRECISION RCOND
241* ..
242* .. Array Arguments ..
243 DOUBLE PRECISION BERR( * ), D( * ), DF( * ), FERR( * ),
244 $ RWORK( * )
245 COMPLEX*16 B( LDB, * ), E( * ), EF( * ), WORK( * ),
246 $ X( LDX, * )
247* ..
248*
249* =====================================================================
250*
251* .. Parameters ..
252 DOUBLE PRECISION ZERO
253 parameter( zero = 0.0d+0 )
254* ..
255* .. Local Scalars ..
256 LOGICAL NOFACT
257 DOUBLE PRECISION ANORM
258* ..
259* .. External Functions ..
260 LOGICAL LSAME
261 DOUBLE PRECISION DLAMCH, ZLANHT
262 EXTERNAL lsame, dlamch, zlanht
263* ..
264* .. External Subroutines ..
265 EXTERNAL dcopy, xerbla, zcopy, zlacpy, zptcon,
266 $ zptrfs,
267 $ zpttrf, zpttrs
268* ..
269* .. Intrinsic Functions ..
270 INTRINSIC max
271* ..
272* .. Executable Statements ..
273*
274* Test the input parameters.
275*
276 info = 0
277 nofact = lsame( fact, 'N' )
278 IF( .NOT.nofact .AND. .NOT.lsame( fact, 'F' ) ) THEN
279 info = -1
280 ELSE IF( n.LT.0 ) THEN
281 info = -2
282 ELSE IF( nrhs.LT.0 ) THEN
283 info = -3
284 ELSE IF( ldb.LT.max( 1, n ) ) THEN
285 info = -9
286 ELSE IF( ldx.LT.max( 1, n ) ) THEN
287 info = -11
288 END IF
289 IF( info.NE.0 ) THEN
290 CALL xerbla( 'ZPTSVX', -info )
291 RETURN
292 END IF
293*
294 IF( nofact ) THEN
295*
296* Compute the L*D*L**H (or U**H*D*U) factorization of A.
297*
298 CALL dcopy( n, d, 1, df, 1 )
299 IF( n.GT.1 )
300 $ CALL zcopy( n-1, e, 1, ef, 1 )
301 CALL zpttrf( n, df, ef, info )
302*
303* Return if INFO is non-zero.
304*
305 IF( info.GT.0 )THEN
306 rcond = zero
307 RETURN
308 END IF
309 END IF
310*
311* Compute the norm of the matrix A.
312*
313 anorm = zlanht( '1', n, d, e )
314*
315* Compute the reciprocal of the condition number of A.
316*
317 CALL zptcon( n, df, ef, anorm, rcond, rwork, info )
318*
319* Compute the solution vectors X.
320*
321 CALL zlacpy( 'Full', n, nrhs, b, ldb, x, ldx )
322 CALL zpttrs( 'Lower', n, nrhs, df, ef, x, ldx, info )
323*
324* Use iterative refinement to improve the computed solutions and
325* compute error bounds and backward error estimates for them.
326*
327 CALL zptrfs( 'Lower', n, nrhs, d, e, df, ef, b, ldb, x, ldx,
328 $ ferr,
329 $ berr, work, rwork, info )
330*
331* Set INFO = N+1 if the matrix is singular to working precision.
332*
333 IF( rcond.LT.dlamch( 'Epsilon' ) )
334 $ info = n + 1
335*
336 RETURN
337*
338* End of ZPTSVX
339*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
Definition zcopy.f:81
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
Definition zlacpy.f:101
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
double precision function zlanht(norm, n, d, e)
ZLANHT returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition zlanht.f:99
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine zptcon(n, d, e, anorm, rcond, rwork, info)
ZPTCON
Definition zptcon.f:117
subroutine zptrfs(uplo, n, nrhs, d, e, df, ef, b, ldb, x, ldx, ferr, berr, work, rwork, info)
ZPTRFS
Definition zptrfs.f:181
subroutine zpttrf(n, d, e, info)
ZPTTRF
Definition zpttrf.f:90
subroutine zpttrs(uplo, n, nrhs, d, e, b, ldb, info)
ZPTTRS
Definition zpttrs.f:119
Here is the call graph for this function:
Here is the caller graph for this function: