--- rpl/lapack/lapack/zlatrs.f 2010/08/06 15:32:47 1.4
+++ rpl/lapack/lapack/zlatrs.f 2012/07/31 11:06:39 1.10
@@ -1,10 +1,248 @@
+*> \brief \b ZLATRS
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download ZLATRS + dependencies
+*>
+*> [TGZ]
+*>
+*> [ZIP]
+*>
+*> [TXT]
+*> \endhtmlonly
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE ZLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE,
+* CNORM, INFO )
+*
+* .. Scalar Arguments ..
+* CHARACTER DIAG, NORMIN, TRANS, UPLO
+* INTEGER INFO, LDA, N
+* DOUBLE PRECISION SCALE
+* ..
+* .. Array Arguments ..
+* DOUBLE PRECISION CNORM( * )
+* COMPLEX*16 A( LDA, * ), X( * )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> ZLATRS solves one of the triangular systems
+*>
+*> A * x = s*b, A**T * x = s*b, or A**H * x = s*b,
+*>
+*> with scaling to prevent overflow. Here A is an upper or lower
+*> triangular matrix, A**T denotes the transpose of A, A**H denotes the
+*> conjugate transpose of A, x and b are n-element vectors, and s is a
+*> scaling factor, usually less than or equal to 1, chosen so that the
+*> components of x will be less than the overflow threshold. If the
+*> unscaled problem will not cause overflow, the Level 2 BLAS routine
+*> ZTRSV is called. If the matrix A is singular (A(j,j) = 0 for some j),
+*> then s is set to 0 and a non-trivial solution to A*x = 0 is returned.
+*> \endverbatim
+*
+* Arguments:
+* ==========
+*
+*> \param[in] UPLO
+*> \verbatim
+*> UPLO is CHARACTER*1
+*> Specifies whether the matrix A is upper or lower triangular.
+*> = 'U': Upper triangular
+*> = 'L': Lower triangular
+*> \endverbatim
+*>
+*> \param[in] TRANS
+*> \verbatim
+*> TRANS is CHARACTER*1
+*> Specifies the operation applied to A.
+*> = 'N': Solve A * x = s*b (No transpose)
+*> = 'T': Solve A**T * x = s*b (Transpose)
+*> = 'C': Solve A**H * x = s*b (Conjugate transpose)
+*> \endverbatim
+*>
+*> \param[in] DIAG
+*> \verbatim
+*> DIAG is CHARACTER*1
+*> Specifies whether or not the matrix A is unit triangular.
+*> = 'N': Non-unit triangular
+*> = 'U': Unit triangular
+*> \endverbatim
+*>
+*> \param[in] NORMIN
+*> \verbatim
+*> NORMIN is CHARACTER*1
+*> Specifies whether CNORM has been set or not.
+*> = 'Y': CNORM contains the column norms on entry
+*> = 'N': CNORM is not set on entry. On exit, the norms will
+*> be computed and stored in CNORM.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*> N is INTEGER
+*> The order of the matrix A. N >= 0.
+*> \endverbatim
+*>
+*> \param[in] A
+*> \verbatim
+*> A is COMPLEX*16 array, dimension (LDA,N)
+*> The triangular matrix A. If UPLO = 'U', the leading n by n
+*> upper triangular part of the array A contains the upper
+*> triangular matrix, and the strictly lower triangular part of
+*> A is not referenced. If UPLO = 'L', the leading n by n lower
+*> triangular part of the array A contains the lower triangular
+*> matrix, and the strictly upper triangular part of A is not
+*> referenced. If DIAG = 'U', the diagonal elements of A are
+*> also not referenced and are assumed to be 1.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*> LDA is INTEGER
+*> The leading dimension of the array A. LDA >= max (1,N).
+*> \endverbatim
+*>
+*> \param[in,out] X
+*> \verbatim
+*> X is COMPLEX*16 array, dimension (N)
+*> On entry, the right hand side b of the triangular system.
+*> On exit, X is overwritten by the solution vector x.
+*> \endverbatim
+*>
+*> \param[out] SCALE
+*> \verbatim
+*> SCALE is DOUBLE PRECISION
+*> The scaling factor s for the triangular system
+*> A * x = s*b, A**T * x = s*b, or A**H * x = s*b.
+*> If SCALE = 0, the matrix A is singular or badly scaled, and
+*> the vector x is an exact or approximate solution to A*x = 0.
+*> \endverbatim
+*>
+*> \param[in,out] CNORM
+*> \verbatim
+*> CNORM is DOUBLE PRECISION array, dimension (N)
+*>
+*> If NORMIN = 'Y', CNORM is an input argument and CNORM(j)
+*> contains the norm of the off-diagonal part of the j-th column
+*> of A. If TRANS = 'N', CNORM(j) must be greater than or equal
+*> to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j)
+*> must be greater than or equal to the 1-norm.
+*>
+*> If NORMIN = 'N', CNORM is an output argument and CNORM(j)
+*> returns the 1-norm of the offdiagonal part of the j-th column
+*> of A.
+*> \endverbatim
+*>
+*> \param[out] INFO
+*> \verbatim
+*> INFO is INTEGER
+*> = 0: successful exit
+*> < 0: if INFO = -k, the k-th argument had an illegal value
+*> \endverbatim
+*
+* Authors:
+* ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \date April 2012
+*
+*> \ingroup complex16OTHERauxiliary
+*
+*> \par Further Details:
+* =====================
+*>
+*> \verbatim
+*>
+*> A rough bound on x is computed; if that is less than overflow, ZTRSV
+*> is called, otherwise, specific code is used which checks for possible
+*> overflow or divide-by-zero at every operation.
+*>
+*> A columnwise scheme is used for solving A*x = b. The basic algorithm
+*> if A is lower triangular is
+*>
+*> x[1:n] := b[1:n]
+*> for j = 1, ..., n
+*> x(j) := x(j) / A(j,j)
+*> x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]
+*> end
+*>
+*> Define bounds on the components of x after j iterations of the loop:
+*> M(j) = bound on x[1:j]
+*> G(j) = bound on x[j+1:n]
+*> Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.
+*>
+*> Then for iteration j+1 we have
+*> M(j+1) <= G(j) / | A(j+1,j+1) |
+*> G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |
+*> <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )
+*>
+*> where CNORM(j+1) is greater than or equal to the infinity-norm of
+*> column j+1 of A, not counting the diagonal. Hence
+*>
+*> G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )
+*> 1<=i<=j
+*> and
+*>
+*> |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )
+*> 1<=i< j
+*>
+*> Since |x(j)| <= M(j), we use the Level 2 BLAS routine ZTRSV if the
+*> reciprocal of the largest M(j), j=1,..,n, is larger than
+*> max(underflow, 1/overflow).
+*>
+*> The bound on x(j) is also used to determine when a step in the
+*> columnwise method can be performed without fear of overflow. If
+*> the computed bound is greater than a large constant, x is scaled to
+*> prevent overflow, but if the bound overflows, x is set to 0, x(j) to
+*> 1, and scale to 0, and a non-trivial solution to A*x = 0 is found.
+*>
+*> Similarly, a row-wise scheme is used to solve A**T *x = b or
+*> A**H *x = b. The basic algorithm for A upper triangular is
+*>
+*> for j = 1, ..., n
+*> x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j)
+*> end
+*>
+*> We simultaneously compute two bounds
+*> G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j
+*> M(j) = bound on x(i), 1<=i<=j
+*>
+*> The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we
+*> add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.
+*> Then the bound on x(j) is
+*>
+*> M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |
+*>
+*> <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )
+*> 1<=i<=j
+*>
+*> and we can safely call ZTRSV if 1/M(n) and 1/G(n) are both greater
+*> than max(underflow, 1/overflow).
+*> \endverbatim
+*>
+* =====================================================================
SUBROUTINE ZLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE,
$ CNORM, INFO )
*
-* -- LAPACK auxiliary routine (version 3.2) --
+* -- LAPACK auxiliary routine (version 3.4.1) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-* November 2006
+* April 2012
*
* .. Scalar Arguments ..
CHARACTER DIAG, NORMIN, TRANS, UPLO
@@ -16,158 +254,6 @@
COMPLEX*16 A( LDA, * ), X( * )
* ..
*
-* Purpose
-* =======
-*
-* ZLATRS solves one of the triangular systems
-*
-* A * x = s*b, A**T * x = s*b, or A**H * x = s*b,
-*
-* with scaling to prevent overflow. Here A is an upper or lower
-* triangular matrix, A**T denotes the transpose of A, A**H denotes the
-* conjugate transpose of A, x and b are n-element vectors, and s is a
-* scaling factor, usually less than or equal to 1, chosen so that the
-* components of x will be less than the overflow threshold. If the
-* unscaled problem will not cause overflow, the Level 2 BLAS routine
-* ZTRSV is called. If the matrix A is singular (A(j,j) = 0 for some j),
-* then s is set to 0 and a non-trivial solution to A*x = 0 is returned.
-*
-* Arguments
-* =========
-*
-* UPLO (input) CHARACTER*1
-* Specifies whether the matrix A is upper or lower triangular.
-* = 'U': Upper triangular
-* = 'L': Lower triangular
-*
-* TRANS (input) CHARACTER*1
-* Specifies the operation applied to A.
-* = 'N': Solve A * x = s*b (No transpose)
-* = 'T': Solve A**T * x = s*b (Transpose)
-* = 'C': Solve A**H * x = s*b (Conjugate transpose)
-*
-* DIAG (input) CHARACTER*1
-* Specifies whether or not the matrix A is unit triangular.
-* = 'N': Non-unit triangular
-* = 'U': Unit triangular
-*
-* NORMIN (input) CHARACTER*1
-* Specifies whether CNORM has been set or not.
-* = 'Y': CNORM contains the column norms on entry
-* = 'N': CNORM is not set on entry. On exit, the norms will
-* be computed and stored in CNORM.
-*
-* N (input) INTEGER
-* The order of the matrix A. N >= 0.
-*
-* A (input) COMPLEX*16 array, dimension (LDA,N)
-* The triangular matrix A. If UPLO = 'U', the leading n by n
-* upper triangular part of the array A contains the upper
-* triangular matrix, and the strictly lower triangular part of
-* A is not referenced. If UPLO = 'L', the leading n by n lower
-* triangular part of the array A contains the lower triangular
-* matrix, and the strictly upper triangular part of A is not
-* referenced. If DIAG = 'U', the diagonal elements of A are
-* also not referenced and are assumed to be 1.
-*
-* LDA (input) INTEGER
-* The leading dimension of the array A. LDA >= max (1,N).
-*
-* X (input/output) COMPLEX*16 array, dimension (N)
-* On entry, the right hand side b of the triangular system.
-* On exit, X is overwritten by the solution vector x.
-*
-* SCALE (output) DOUBLE PRECISION
-* The scaling factor s for the triangular system
-* A * x = s*b, A**T * x = s*b, or A**H * x = s*b.
-* If SCALE = 0, the matrix A is singular or badly scaled, and
-* the vector x is an exact or approximate solution to A*x = 0.
-*
-* CNORM (input or output) DOUBLE PRECISION array, dimension (N)
-*
-* If NORMIN = 'Y', CNORM is an input argument and CNORM(j)
-* contains the norm of the off-diagonal part of the j-th column
-* of A. If TRANS = 'N', CNORM(j) must be greater than or equal
-* to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j)
-* must be greater than or equal to the 1-norm.
-*
-* If NORMIN = 'N', CNORM is an output argument and CNORM(j)
-* returns the 1-norm of the offdiagonal part of the j-th column
-* of A.
-*
-* INFO (output) INTEGER
-* = 0: successful exit
-* < 0: if INFO = -k, the k-th argument had an illegal value
-*
-* Further Details
-* ======= =======
-*
-* A rough bound on x is computed; if that is less than overflow, ZTRSV
-* is called, otherwise, specific code is used which checks for possible
-* overflow or divide-by-zero at every operation.
-*
-* A columnwise scheme is used for solving A*x = b. The basic algorithm
-* if A is lower triangular is
-*
-* x[1:n] := b[1:n]
-* for j = 1, ..., n
-* x(j) := x(j) / A(j,j)
-* x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]
-* end
-*
-* Define bounds on the components of x after j iterations of the loop:
-* M(j) = bound on x[1:j]
-* G(j) = bound on x[j+1:n]
-* Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.
-*
-* Then for iteration j+1 we have
-* M(j+1) <= G(j) / | A(j+1,j+1) |
-* G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |
-* <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )
-*
-* where CNORM(j+1) is greater than or equal to the infinity-norm of
-* column j+1 of A, not counting the diagonal. Hence
-*
-* G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )
-* 1<=i<=j
-* and
-*
-* |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )
-* 1<=i< j
-*
-* Since |x(j)| <= M(j), we use the Level 2 BLAS routine ZTRSV if the
-* reciprocal of the largest M(j), j=1,..,n, is larger than
-* max(underflow, 1/overflow).
-*
-* The bound on x(j) is also used to determine when a step in the
-* columnwise method can be performed without fear of overflow. If
-* the computed bound is greater than a large constant, x is scaled to
-* prevent overflow, but if the bound overflows, x is set to 0, x(j) to
-* 1, and scale to 0, and a non-trivial solution to A*x = 0 is found.
-*
-* Similarly, a row-wise scheme is used to solve A**T *x = b or
-* A**H *x = b. The basic algorithm for A upper triangular is
-*
-* for j = 1, ..., n
-* x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j)
-* end
-*
-* We simultaneously compute two bounds
-* G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j
-* M(j) = bound on x(i), 1<=i<=j
-*
-* The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we
-* add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.
-* Then the bound on x(j) is
-*
-* M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |
-*
-* <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )
-* 1<=i<=j
-*
-* and we can safely call ZTRSV if 1/M(n) and 1/G(n) are both greater
-* than max(underflow, 1/overflow).
-*
* =====================================================================
*
* .. Parameters ..