--- rpl/lapack/lapack/dlatps.f 2010/08/07 13:22:21 1.5
+++ rpl/lapack/lapack/dlatps.f 2016/08/27 15:34:32 1.16
@@ -1,10 +1,238 @@
+*> \brief \b DLATPS solves a triangular system of equations with the matrix held in packed storage.
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download DLATPS + dependencies
+*>
+*> [TGZ]
+*>
+*> [ZIP]
+*>
+*> [TXT]
+*> \endhtmlonly
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE DLATPS( UPLO, TRANS, DIAG, NORMIN, N, AP, X, SCALE,
+* CNORM, INFO )
+*
+* .. Scalar Arguments ..
+* CHARACTER DIAG, NORMIN, TRANS, UPLO
+* INTEGER INFO, N
+* DOUBLE PRECISION SCALE
+* ..
+* .. Array Arguments ..
+* DOUBLE PRECISION AP( * ), CNORM( * ), X( * )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> DLATPS solves one of the triangular systems
+*>
+*> A *x = s*b or A**T*x = s*b
+*>
+*> with scaling to prevent overflow, where A is an upper or lower
+*> triangular matrix stored in packed form. Here A**T denotes the
+*> 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
+*> DTPSV 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**T* x = s*b (Conjugate transpose = 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] AP
+*> \verbatim
+*> AP is DOUBLE PRECISION array, dimension (N*(N+1)/2)
+*> The upper or lower triangular 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.
+*> \endverbatim
+*>
+*> \param[in,out] X
+*> \verbatim
+*> X is DOUBLE PRECISION 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 or A**T* 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 September 2012
+*
+*> \ingroup doubleOTHERauxiliary
+*
+*> \par Further Details:
+* =====================
+*>
+*> \verbatim
+*>
+*> A rough bound on x is computed; if that is less than overflow, DTPSV
+*> 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 DTPSV 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. The basic
+*> algorithm for A upper triangular is
+*>
+*> for j = 1, ..., n
+*> x(j) := ( b(j) - A[1:j-1,j]**T * x[1:j-1] ) / A(j,j)
+*> end
+*>
+*> We simultaneously compute two bounds
+*> G(j) = bound on ( b(i) - A[1:i-1,i]**T * 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 DTPSV if 1/M(n) and 1/G(n) are both greater
+*> than max(underflow, 1/overflow).
+*> \endverbatim
+*>
+* =====================================================================
SUBROUTINE DLATPS( UPLO, TRANS, DIAG, NORMIN, N, AP, X, SCALE,
$ CNORM, INFO )
*
-* -- LAPACK auxiliary routine (version 3.2) --
+* -- LAPACK auxiliary routine (version 3.4.2) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-* November 2006
+* September 2012
*
* .. Scalar Arguments ..
CHARACTER DIAG, NORMIN, TRANS, UPLO
@@ -15,152 +243,6 @@
DOUBLE PRECISION AP( * ), CNORM( * ), X( * )
* ..
*
-* Purpose
-* =======
-*
-* DLATPS solves one of the triangular systems
-*
-* A *x = s*b or A'*x = s*b
-*
-* with scaling to prevent overflow, where A is an upper or lower
-* triangular matrix stored in packed form. Here A' denotes the
-* 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
-* DTPSV 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'* x = s*b (Transpose)
-* = 'C': Solve A'* x = s*b (Conjugate transpose = 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.
-*
-* AP (input) DOUBLE PRECISION array, dimension (N*(N+1)/2)
-* The upper or lower triangular 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.
-*
-* X (input/output) DOUBLE PRECISION 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 or A'* 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, DTPSV
-* 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 DTPSV 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'*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 DTPSV if 1/M(n) and 1/G(n) are both greater
-* than max(underflow, 1/overflow).
-*
* =====================================================================
*
* .. Parameters ..
@@ -346,7 +428,7 @@
*
ELSE
*
-* Compute the growth in A' * x = b.
+* Compute the growth in A**T * x = b.
*
IF( UPPER ) THEN
JFIRST = 1
@@ -561,7 +643,7 @@
*
ELSE
*
-* Solve A' * x = b
+* Solve A**T * x = b
*
IP = JFIRST*( JFIRST+1 ) / 2
JLEN = 1
@@ -675,7 +757,7 @@
ELSE
*
* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and
-* scale = 0, and compute a solution to A'*x = 0.
+* scale = 0, and compute a solution to A**T*x = 0.
*
DO 140 I = 1, N
X( I ) = ZERO