DPBSVX(l)	     LAPACK driver routine (version 1.1)	    DPBSVX(l)

NAME
  DPBSVX - use the Cholesky factorization A = U**T*U or	A = L*L**T to compute
  the solution to a real system	of linear equations  A * X = B,

SYNOPSIS

  SUBROUTINE DPBSVX( FACT, UPLO, N, KD,	NRHS, AB, LDAB,	AFB, LDAFB, EQUED, S,
		     B,	LDB, X,	LDX, RCOND, FERR, BERR,	WORK, IWORK, INFO )

      CHARACTER	     EQUED, FACT, UPLO

      INTEGER	     INFO, KD, LDAB, LDAFB, LDB, LDX, N, NRHS

      DOUBLE	     PRECISION RCOND

      INTEGER	     IWORK( * )

      DOUBLE	     PRECISION AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ),
		     BERR( * ),	FERR( *	), S( *	), WORK( * ), X( LDX, *	)

PURPOSE
  DPBSVX uses the Cholesky factorization A = U**T*U or A = L*L**T to compute
  the solution to a real system	of linear equations
     A * X = B,	where A	is an N-by-N symmetric 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**T * U,  if UPLO = 'U', or
	A = L *	L**T,  if UPLO = 'L',
     where U is	an upper triangular band matrix, and L is a lower
     triangular	band matrix.

  3. 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, steps	4-6 are	skipped.

  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.

ARGUMENTS

  FACT	  (input) CHARACTER*1
	  Specifies whether or not the factored	form of	the matrix A is	sup-
	  plied	on entry, and if not, whether the matrix A should be equili-
	  brated before	it is factored.	 = 'F':	 On entry, AFB contains	the
	  factored form	of A.  If EQUED	= 'Y', the matrix A has	been equili-
	  brated 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.

  UPLO	  (input) CHARACTER*1
	  = 'U':  Upper	triangle of A is stored;
	  = 'L':  Lower	triangle of A is stored.

  N	  (input) INTEGER
	  The number of	linear equations, i.e.,	the order of the matrix	A.  N
	  >= 0.

  KD	  (input) INTEGER
	  The number of	superdiagonals of the matrix A if UPLO = 'U', or the
	  number of subdiagonals if UPLO = 'L'.	 KD >= 0.

  NRHS	  (input) INTEGER
	  The number of	right-hand sides, i.e.,	the number of columns of the
	  matrices B and X.  NRHS >= 0.

  AB	  (input/output) DOUBLE	PRECISION array, dimension (LDAB,N)
	  On entry, the	upper or lower triangle	of the symmetric 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).

  LDAB	  (input) INTEGER
	  The leading dimension	of the array A.	 LDAB >= KD+1.

  AFB	  (input or output) DOUBLE PRECISION 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**T*U or A =	L*L**T of the band matrix A, in	the same storage for-
	  mat 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**T*U or A =	L*L**T.

	  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**T*U or A =	L*L**T of the equilibrated matrix A (see the
	  description of A for the form	of the equilibrated matrix).

  LDAFB	  (input) INTEGER
	  The leading dimension	of the array AFB.  LDAFB >= KD+1.

  EQUED	  (input/output) 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 variable if	FACT = 'F';
	  otherwise, it	is an output variable.

  S	  (input/output) DOUBLE	PRECISION array, dimension (N)
	  The scale factors for	A; not accessed	if EQUED = 'N'.	 S is an
	  input	variable if FACT = 'F';	otherwise, S is	an output variable.
	  If FACT = 'F'	and EQUED = 'Y', each element of S must	be positive.

  B	  (input/output) DOUBLE	PRECISION 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.

  LDB	  (input) INTEGER
	  The leading dimension	of the array B.	 LDB >=	max(1,N).

  X	  (output) DOUBLE PRECISION array, dimension (LDX,NRHS)
	  If INFO = 0, 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.

  LDX	  (input) INTEGER
	  The leading dimension	of the array X.	 LDX >=	max(1,N).

  RCOND	  (output) 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, and	the solution and error bounds are not computed.

  FERR	  (output) DOUBLE PRECISION array, dimension (NRHS)
	  The estimated	forward	error bounds for each solution vector X(j)
	  (the j-th column of the solution matrix X).  If XTRUE	is the true
	  solution, FERR(j) bounds the magnitude of the	largest	entry in
	  (X(j)	- XTRUE) divided by the	magnitude of the largest entry in
	  X(j).	 The quality of	the error bound	depends	on the quality of the
	  estimate of norm(inv(A)) computed in the code; if the	estimate of
	  norm(inv(A)) is accurate, the	error bound is guaranteed.

  BERR	  (output) DOUBLE PRECISION array, dimension (NRHS)
	  The componentwise relative backward error of each solution vector
	  X(j) (i.e., the smallest relative change in any entry	of A or	B
	  that makes X(j) an exact solution).

  WORK	  (workspace) DOUBLE PRECISION array, dimension	(3*N)

  IWORK	  (workspace) INTEGER array, dimension (N)

  INFO	  (output) INTEGER
	  = 0:	successful exit
	  < 0:	if INFO	= -i, the i-th argument	had an illegal value
	  > 0:	if INFO	= i, and i is
	  <= N:	if INFO	= i, the leading minor of order	i of A is not posi-
	  tive definite, so the	factorization could not	be completed, and the
	  solution has not been	computed.  = N+1: RCOND	is less	than machine
	  precision.  The factorization	has been completed, but	the matrix is
	  singular to working precision, and the solution and error bounds
	  have not been	computed.

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 symmetric 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.


Back to the listing of expert driver routines