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

◆ cptsvx()

subroutine cptsvx ( character fact,
integer n,
integer nrhs,
real, dimension( * ) d,
complex, dimension( * ) e,
real, dimension( * ) df,
complex, dimension( * ) ef,
complex, dimension( ldb, * ) b,
integer ldb,
complex, dimension( ldx, * ) x,
integer ldx,
real rcond,
real, dimension( * ) ferr,
real, dimension( * ) berr,
complex, dimension( * ) work,
real, dimension( * ) rwork,
integer info )

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

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

Purpose:
!>
!> CPTSVX 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 REAL array, dimension (N)
!>          The n diagonal elements of the tridiagonal matrix A.
!>
[in]E
!>          E is COMPLEX array, dimension (N-1)
!>          The (n-1) subdiagonal elements of the tridiagonal matrix A.
!>
[in,out]DF
!>          DF is REAL 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 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 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 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 REAL
!>          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 REAL 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 REAL 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 array, dimension (N)
!>
[out]RWORK
!>          RWORK is REAL 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 cptsvx.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 REAL RCOND
241* ..
242* .. Array Arguments ..
243 REAL BERR( * ), D( * ), DF( * ), FERR( * ),
244 $ RWORK( * )
245 COMPLEX B( LDB, * ), E( * ), EF( * ), WORK( * ),
246 $ X( LDX, * )
247* ..
248*
249* =====================================================================
250*
251* .. Parameters ..
252 REAL ZERO
253 parameter( zero = 0.0e+0 )
254* ..
255* .. Local Scalars ..
256 LOGICAL NOFACT
257 REAL ANORM
258* ..
259* .. External Functions ..
260 LOGICAL LSAME
261 REAL CLANHT, SLAMCH
262 EXTERNAL lsame, clanht, slamch
263* ..
264* .. External Subroutines ..
265 EXTERNAL ccopy, clacpy, cptcon, cptrfs, cpttrf,
266 $ cpttrs,
267 $ scopy, xerbla
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( 'CPTSVX', -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 scopy( n, d, 1, df, 1 )
299 IF( n.GT.1 )
300 $ CALL ccopy( n-1, e, 1, ef, 1 )
301 CALL cpttrf( 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 = clanht( '1', n, d, e )
314*
315* Compute the reciprocal of the condition number of A.
316*
317 CALL cptcon( n, df, ef, anorm, rcond, rwork, info )
318*
319* Compute the solution vectors X.
320*
321 CALL clacpy( 'Full', n, nrhs, b, ldb, x, ldx )
322 CALL cpttrs( '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 cptrfs( '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.slamch( 'Epsilon' ) )
334 $ info = n + 1
335*
336 RETURN
337*
338* End of CPTSVX
339*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
scopy
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
ccopy
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
clacpy
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
Definition clacpy.f:101
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
clanht
real function clanht(norm, n, d, e)
CLANHT returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition clanht.f:99
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
cptcon
subroutine cptcon(n, d, e, anorm, rcond, rwork, info)
CPTCON
Definition cptcon.f:117
cptrfs
subroutine cptrfs(uplo, n, nrhs, d, e, df, ef, b, ldb, x, ldx, ferr, berr, work, rwork, info)
CPTRFS
Definition cptrfs.f:181
cpttrf
subroutine cpttrf(n, d, e, info)
CPTTRF
Definition cpttrf.f:90
cpttrs
subroutine cpttrs(uplo, n, nrhs, d, e, b, ldb, info)
CPTTRS
Definition cpttrs.f:119
Here is the call graph for this function:
Here is the caller graph for this function: