Update LAPACK to 3.9.0

This commit is contained in:
Martin Kroeker 2019-12-30 13:13:42 +01:00 committed by GitHub
parent febeb9bcc1
commit d4859bb2cc
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
75 changed files with 1378 additions and 521 deletions

View File

@ -124,7 +124,7 @@
*> LDVT is INTEGER. *> LDVT is INTEGER.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] WORK *> \param[out] WORK
*> \verbatim *> \verbatim
*> WORK is COMPLEX*16 array. Workspace of size nb. *> WORK is COMPLEX*16 array. Workspace of size nb.
*> \endverbatim *> \endverbatim

View File

@ -19,7 +19,7 @@
* =========== * ===========
* *
* SUBROUTINE ZHECON_3( UPLO, N, A, LDA, E, IPIV, ANORM, RCOND, * SUBROUTINE ZHECON_3( UPLO, N, A, LDA, E, IPIV, ANORM, RCOND,
* WORK, IWORK, INFO ) * WORK, INFO )
* *
* .. Scalar Arguments .. * .. Scalar Arguments ..
* CHARACTER UPLO * CHARACTER UPLO
@ -27,7 +27,7 @@
* DOUBLE PRECISION ANORM, RCOND * DOUBLE PRECISION ANORM, RCOND
* .. * ..
* .. Array Arguments .. * .. Array Arguments ..
* INTEGER IPIV( * ), IWORK( * ) * INTEGER IPIV( * )
* COMPLEX*16 A( LDA, * ), E ( * ), WORK( * ) * COMPLEX*16 A( LDA, * ), E ( * ), WORK( * )
* .. * ..
* *
@ -129,11 +129,6 @@
*> WORK is COMPLEX*16 array, dimension (2*N) *> WORK is COMPLEX*16 array, dimension (2*N)
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] IWORK
*> \verbatim
*> IWORK is INTEGER array, dimension (N)
*> \endverbatim
*>
*> \param[out] INFO *> \param[out] INFO
*> \verbatim *> \verbatim
*> INFO is INTEGER *> INFO is INTEGER

View File

@ -210,7 +210,7 @@
*> eigenvalues are computed to high relative accuracy when *> eigenvalues are computed to high relative accuracy when
*> possible in future releases. The current code does not *> possible in future releases. The current code does not
*> make any guarantees about high relative accuracy, but *> make any guarantees about high relative accuracy, but
*> furutre releases will. See J. Barlow and J. Demmel, *> future releases will. See J. Barlow and J. Demmel,
*> "Computing Accurate Eigensystems of Scaled Diagonally *> "Computing Accurate Eigensystems of Scaled Diagonally
*> Dominant Matrices", LAPACK Working Note #7, for a discussion *> Dominant Matrices", LAPACK Working Note #7, for a discussion
*> of which matrices define their eigenvalues to high relative *> of which matrices define their eigenvalues to high relative

View File

@ -217,7 +217,7 @@
*> eigenvalues are computed to high relative accuracy when *> eigenvalues are computed to high relative accuracy when
*> possible in future releases. The current code does not *> possible in future releases. The current code does not
*> make any guarantees about high relative accuracy, but *> make any guarantees about high relative accuracy, but
*> furutre releases will. See J. Barlow and J. Demmel, *> future releases will. See J. Barlow and J. Demmel,
*> "Computing Accurate Eigensystems of Scaled Diagonally *> "Computing Accurate Eigensystems of Scaled Diagonally
*> Dominant Matrices", LAPACK Working Note #7, for a discussion *> Dominant Matrices", LAPACK Working Note #7, for a discussion
*> of which matrices define their eigenvalues to high relative *> of which matrices define their eigenvalues to high relative

View File

@ -97,6 +97,7 @@
*> B is COMPLEX*16 array, dimension (LDB,N) *> B is COMPLEX*16 array, dimension (LDB,N)
*> The triangular factor from the Cholesky factorization of B, *> The triangular factor from the Cholesky factorization of B,
*> as returned by ZPOTRF. *> as returned by ZPOTRF.
*> B is modified by the routine but restored on exit.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] LDB *> \param[in] LDB

View File

@ -97,6 +97,7 @@
*> B is COMPLEX*16 array, dimension (LDB,N) *> B is COMPLEX*16 array, dimension (LDB,N)
*> The triangular factor from the Cholesky factorization of B, *> The triangular factor from the Cholesky factorization of B,
*> as returned by ZPOTRF. *> as returned by ZPOTRF.
*> B is modified by the routine but restored on exit.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] LDB *> \param[in] LDB

View File

@ -102,7 +102,7 @@
*> \param[in] A *> \param[in] A
*> \verbatim *> \verbatim
*> A is COMPLEX*16 array, dimension (LDA,N) *> A is COMPLEX*16 array, dimension (LDA,N)
*> The symmetric matrix A. If UPLO = 'U', the leading N-by-N *> The Hermitian matrix A. If UPLO = 'U', the leading N-by-N
*> upper triangular part of A contains the upper triangular *> upper triangular part of A contains the upper triangular
*> part of the matrix A, and the strictly lower triangular *> part of the matrix A, and the strictly lower triangular
*> part of A is not referenced. If UPLO = 'L', the leading *> part of A is not referenced. If UPLO = 'L', the leading
@ -270,7 +270,7 @@
*> information as described below. There currently are up to three *> information as described below. There currently are up to three
*> pieces of information returned for each right-hand side. If *> pieces of information returned for each right-hand side. If
*> componentwise accuracy is not requested (PARAMS(3) = 0.0), then *> componentwise accuracy is not requested (PARAMS(3) = 0.0), then
*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most *> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS < 3, then at most
*> the first (:,N_ERR_BNDS) entries are returned. *> the first (:,N_ERR_BNDS) entries are returned.
*> *>
*> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith *> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
@ -306,14 +306,14 @@
*> \param[in] NPARAMS *> \param[in] NPARAMS
*> \verbatim *> \verbatim
*> NPARAMS is INTEGER *> NPARAMS is INTEGER
*> Specifies the number of parameters set in PARAMS. If .LE. 0, the *> Specifies the number of parameters set in PARAMS. If <= 0, the
*> PARAMS array is never referenced and default values are used. *> PARAMS array is never referenced and default values are used.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in,out] PARAMS *> \param[in,out] PARAMS
*> \verbatim *> \verbatim
*> PARAMS is DOUBLE PRECISION array, dimension NPARAMS *> PARAMS is DOUBLE PRECISION array, dimension NPARAMS
*> Specifies algorithm parameters. If an entry is .LT. 0.0, then *> Specifies algorithm parameters. If an entry is < 0.0, then
*> that entry will be filled with default value used for that *> that entry will be filled with default value used for that
*> parameter. Only positions up to NPARAMS are accessed; defaults *> parameter. Only positions up to NPARAMS are accessed; defaults
*> are used for higher-numbered parameters. *> are used for higher-numbered parameters.

View File

@ -42,7 +42,7 @@
*> matrices. *> matrices.
*> *>
*> Aasen's algorithm is used to factor A as *> Aasen's algorithm is used to factor A as
*> A = U * T * U**H, if UPLO = 'U', or *> A = U**H * T * U, if UPLO = 'U', or
*> A = L * T * L**H, if UPLO = 'L', *> A = L * T * L**H, if UPLO = 'L',
*> where U (or L) is a product of permutation and unit upper (lower) *> where U (or L) is a product of permutation and unit upper (lower)
*> triangular matrices, and T is Hermitian and tridiagonal. The factored form *> triangular matrices, and T is Hermitian and tridiagonal. The factored form
@ -86,7 +86,7 @@
*> *>
*> On exit, if INFO = 0, the tridiagonal matrix T and the *> On exit, if INFO = 0, the tridiagonal matrix T and the
*> multipliers used to obtain the factor U or L from the *> multipliers used to obtain the factor U or L from the
*> factorization A = U*T*U**H or A = L*T*L**H as computed by *> factorization A = U**H*T*U or A = L*T*L**H as computed by
*> ZHETRF_AA. *> ZHETRF_AA.
*> \endverbatim *> \endverbatim
*> *>
@ -230,7 +230,7 @@
RETURN RETURN
END IF END IF
* *
* Compute the factorization A = U*T*U**H or A = L*T*L**H. * Compute the factorization A = U**H*T*U or A = L*T*L**H.
* *
CALL ZHETRF_AA( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO ) CALL ZHETRF_AA( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO )
IF( INFO.EQ.0 ) THEN IF( INFO.EQ.0 ) THEN

View File

@ -44,7 +44,7 @@
*> matrices. *> matrices.
*> *>
*> Aasen's 2-stage algorithm is used to factor A as *> Aasen's 2-stage algorithm is used to factor A as
*> A = U * T * U**H, if UPLO = 'U', or *> A = U**H * T * U, if UPLO = 'U', or
*> A = L * T * L**H, if UPLO = 'L', *> A = L * T * L**H, if UPLO = 'L',
*> where U (or L) is a product of permutation and unit upper (lower) *> where U (or L) is a product of permutation and unit upper (lower)
*> triangular matrices, and T is Hermitian and band. The matrix T is *> triangular matrices, and T is Hermitian and band. The matrix T is
@ -211,9 +211,7 @@
* *
* .. Local Scalars .. * .. Local Scalars ..
LOGICAL UPPER, TQUERY, WQUERY LOGICAL UPPER, TQUERY, WQUERY
INTEGER I, J, K, I1, I2, TD INTEGER LWKOPT
INTEGER LDTB, LWKOPT, NB, KB, NT, IINFO
COMPLEX PIV
* .. * ..
* .. External Functions .. * .. External Functions ..
LOGICAL LSAME LOGICAL LSAME
@ -263,7 +261,7 @@
RETURN RETURN
END IF END IF
* *
* Compute the factorization A = U*T*U**H or A = L*T*L**H. * Compute the factorization A = U**H*T*U or A = L*T*L**H.
* *
CALL ZHETRF_AA_2STAGE( UPLO, N, A, LDA, TB, LTB, IPIV, IPIV2, CALL ZHETRF_AA_2STAGE( UPLO, N, A, LDA, TB, LTB, IPIV, IPIV2,
$ WORK, LWORK, INFO ) $ WORK, LWORK, INFO )

View File

@ -46,7 +46,7 @@
*> *>
*> ZHESVXX uses the diagonal pivoting factorization to compute the *> ZHESVXX uses the diagonal pivoting factorization to compute the
*> solution to a complex*16 system of linear equations A * X = B, where *> solution to a complex*16 system of linear equations A * X = B, where
*> A is an N-by-N symmetric matrix and X and B are N-by-NRHS *> A is an N-by-N Hermitian matrix and X and B are N-by-NRHS
*> matrices. *> matrices.
*> *>
*> If requested, both normwise and maximum componentwise error bounds *> If requested, both normwise and maximum componentwise error bounds
@ -88,7 +88,7 @@
*> A = L * D * L**T, if UPLO = 'L', *> A = L * D * L**T, if UPLO = 'L',
*> *>
*> where U (or L) is a product of permutation and unit upper (lower) *> where U (or L) is a product of permutation and unit upper (lower)
*> triangular matrices, and D is symmetric and block diagonal with *> triangular matrices, and D is Hermitian and block diagonal with
*> 1-by-1 and 2-by-2 diagonal blocks. *> 1-by-1 and 2-by-2 diagonal blocks.
*> *>
*> 3. If some D(i,i)=0, so that D is exactly singular, then the *> 3. If some D(i,i)=0, so that D is exactly singular, then the
@ -161,7 +161,7 @@
*> \param[in,out] A *> \param[in,out] A
*> \verbatim *> \verbatim
*> A is COMPLEX*16 array, dimension (LDA,N) *> A is COMPLEX*16 array, dimension (LDA,N)
*> The symmetric matrix A. If UPLO = 'U', the leading N-by-N *> The Hermitian matrix A. If UPLO = 'U', the leading N-by-N
*> upper triangular part of A contains the upper triangular *> upper triangular part of A contains the upper triangular
*> part of the matrix A, and the strictly lower triangular *> part of the matrix A, and the strictly lower triangular
*> part of A is not referenced. If UPLO = 'L', the leading *> part of A is not referenced. If UPLO = 'L', the leading
@ -378,7 +378,7 @@
*> information as described below. There currently are up to three *> information as described below. There currently are up to three
*> pieces of information returned for each right-hand side. If *> pieces of information returned for each right-hand side. If
*> componentwise accuracy is not requested (PARAMS(3) = 0.0), then *> componentwise accuracy is not requested (PARAMS(3) = 0.0), then
*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most *> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS < 3, then at most
*> the first (:,N_ERR_BNDS) entries are returned. *> the first (:,N_ERR_BNDS) entries are returned.
*> *>
*> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith *> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
@ -414,14 +414,14 @@
*> \param[in] NPARAMS *> \param[in] NPARAMS
*> \verbatim *> \verbatim
*> NPARAMS is INTEGER *> NPARAMS is INTEGER
*> Specifies the number of parameters set in PARAMS. If .LE. 0, the *> Specifies the number of parameters set in PARAMS. If <= 0, the
*> PARAMS array is never referenced and default values are used. *> PARAMS array is never referenced and default values are used.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in,out] PARAMS *> \param[in,out] PARAMS
*> \verbatim *> \verbatim
*> PARAMS is DOUBLE PRECISION array, dimension NPARAMS *> PARAMS is DOUBLE PRECISION array, dimension NPARAMS
*> Specifies algorithm parameters. If an entry is .LT. 0.0, then *> Specifies algorithm parameters. If an entry is < 0.0, then
*> that entry will be filled with default value used for that *> that entry will be filled with default value used for that
*> parameter. Only positions up to NPARAMS are accessed; defaults *> parameter. Only positions up to NPARAMS are accessed; defaults
*> are used for higher-numbered parameters. *> are used for higher-numbered parameters.

View File

@ -322,7 +322,7 @@
* *
* Factorize A as U*D*U**H using the upper triangle of A * Factorize A as U*D*U**H using the upper triangle of A
* *
* Initilize the first entry of array E, where superdiagonal * Initialize the first entry of array E, where superdiagonal
* elements of D are stored * elements of D are stored
* *
E( 1 ) = CZERO E( 1 ) = CZERO
@ -676,7 +676,7 @@
* *
* Factorize A as L*D*L**H using the lower triangle of A * Factorize A as L*D*L**H using the lower triangle of A
* *
* Initilize the unused last entry of the subdiagonal array E. * Initialize the unused last entry of the subdiagonal array E.
* *
E( N ) = CZERO E( N ) = CZERO
* *

View File

@ -123,23 +123,22 @@
*> *>
*> \param[out] HOUS2 *> \param[out] HOUS2
*> \verbatim *> \verbatim
*> HOUS2 is COMPLEX*16 array, dimension LHOUS2, that *> HOUS2 is COMPLEX*16 array, dimension (LHOUS2)
*> store the Householder representation of the stage2 *> Stores the Householder representation of the stage2
*> band to tridiagonal. *> band to tridiagonal.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] LHOUS2 *> \param[in] LHOUS2
*> \verbatim *> \verbatim
*> LHOUS2 is INTEGER *> LHOUS2 is INTEGER
*> The dimension of the array HOUS2. LHOUS2 = MAX(1, dimension) *> The dimension of the array HOUS2.
*> If LWORK = -1, or LHOUS2 = -1, *> If LWORK = -1, or LHOUS2 = -1,
*> then a query is assumed; the routine *> then a query is assumed; the routine
*> only calculates the optimal size of the HOUS2 array, returns *> only calculates the optimal size of the HOUS2 array, returns
*> this value as the first entry of the HOUS2 array, and no error *> this value as the first entry of the HOUS2 array, and no error
*> message related to LHOUS2 is issued by XERBLA. *> message related to LHOUS2 is issued by XERBLA.
*> LHOUS2 = MAX(1, dimension) where *> If VECT='N', LHOUS2 = max(1, 4*n);
*> dimension = 4*N if VECT='N' *> if VECT='V', option not yet available.
*> not available now if VECT='H'
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] WORK *> \param[out] WORK

View File

@ -50,9 +50,9 @@
* Arguments: * Arguments:
* ========== * ==========
* *
*> \param[in] STAGE *> \param[in] STAGE1
*> \verbatim *> \verbatim
*> STAGE is CHARACTER*1 *> STAGE1 is CHARACTER*1
*> = 'N': "No": to mention that the stage 1 of the reduction *> = 'N': "No": to mention that the stage 1 of the reduction
*> from dense to band using the zhetrd_he2hb routine *> from dense to band using the zhetrd_he2hb routine
*> was not called before this routine to reproduce AB. *> was not called before this routine to reproduce AB.
@ -512,8 +512,7 @@ C END IF
* *
* Call the kernel * Call the kernel
* *
#if defined(_OPENMP) && _OPENMP >= 201307 #if defined(_OPENMP)
IF( TTYPE.NE.1 ) THEN IF( TTYPE.NE.1 ) THEN
!$OMP TASK DEPEND(in:WORK(MYID+SHIFT-1)) !$OMP TASK DEPEND(in:WORK(MYID+SHIFT-1))
!$OMP$ DEPEND(in:WORK(MYID-1)) !$OMP$ DEPEND(in:WORK(MYID-1))

View File

@ -37,7 +37,7 @@
*> ZHETRF_AA computes the factorization of a complex hermitian matrix A *> ZHETRF_AA computes the factorization of a complex hermitian matrix A
*> using the Aasen's algorithm. The form of the factorization is *> using the Aasen's algorithm. The form of the factorization is
*> *>
*> A = U*T*U**H or A = L*T*L**H *> A = U**H*T*U or A = L*T*L**H
*> *>
*> where U (or L) is a product of permutation and unit upper (lower) *> where U (or L) is a product of permutation and unit upper (lower)
*> triangular matrices, and T is a hermitian tridiagonal matrix. *> triangular matrices, and T is a hermitian tridiagonal matrix.
@ -223,7 +223,7 @@
IF( UPPER ) THEN IF( UPPER ) THEN
* *
* ..................................................... * .....................................................
* Factorize A as L*D*L**H using the upper triangle of A * Factorize A as U**H*D*U using the upper triangle of A
* ..................................................... * .....................................................
* *
* copy first row A(1, 1:N) into H(1:n) (stored in WORK(1:N)) * copy first row A(1, 1:N) into H(1:n) (stored in WORK(1:N))
@ -256,7 +256,7 @@
$ A( MAX(1, J), J+1 ), LDA, $ A( MAX(1, J), J+1 ), LDA,
$ IPIV( J+1 ), WORK, N, WORK( N*NB+1 ) ) $ IPIV( J+1 ), WORK, N, WORK( N*NB+1 ) )
* *
* Ajust IPIV and apply it back (J-th step picks (J+1)-th pivot) * Adjust IPIV and apply it back (J-th step picks (J+1)-th pivot)
* *
DO J2 = J+2, MIN(N, J+JB+1) DO J2 = J+2, MIN(N, J+JB+1)
IPIV( J2 ) = IPIV( J2 ) + J IPIV( J2 ) = IPIV( J2 ) + J
@ -376,7 +376,7 @@
$ A( J+1, MAX(1, J) ), LDA, $ A( J+1, MAX(1, J) ), LDA,
$ IPIV( J+1 ), WORK, N, WORK( N*NB+1 ) ) $ IPIV( J+1 ), WORK, N, WORK( N*NB+1 ) )
* *
* Ajust IPIV and apply it back (J-th step picks (J+1)-th pivot) * Adjust IPIV and apply it back (J-th step picks (J+1)-th pivot)
* *
DO J2 = J+2, MIN(N, J+JB+1) DO J2 = J+2, MIN(N, J+JB+1)
IPIV( J2 ) = IPIV( J2 ) + J IPIV( J2 ) = IPIV( J2 ) + J

View File

@ -38,7 +38,7 @@
*> ZHETRF_AA_2STAGE computes the factorization of a double hermitian matrix A *> ZHETRF_AA_2STAGE computes the factorization of a double hermitian matrix A
*> using the Aasen's algorithm. The form of the factorization is *> using the Aasen's algorithm. The form of the factorization is
*> *>
*> A = U*T*U**T or A = L*T*L**T *> A = U**H*T*U or A = L*T*L**H
*> *>
*> where U (or L) is a product of permutation and unit upper (lower) *> where U (or L) is a product of permutation and unit upper (lower)
*> triangular matrices, and T is a hermitian band matrix with the *> triangular matrices, and T is a hermitian band matrix with the
@ -66,7 +66,7 @@
*> *>
*> \param[in,out] A *> \param[in,out] A
*> \verbatim *> \verbatim
*> A is COMPLEX array, dimension (LDA,N) *> A is COMPLEX*16 array, dimension (LDA,N)
*> On entry, the hermitian matrix A. If UPLO = 'U', the leading *> On entry, the hermitian matrix A. If UPLO = 'U', the leading
*> N-by-N upper triangular part of A contains the upper *> N-by-N upper triangular part of A contains the upper
*> triangular part of the matrix A, and the strictly lower *> triangular part of the matrix A, and the strictly lower
@ -87,7 +87,7 @@
*> *>
*> \param[out] TB *> \param[out] TB
*> \verbatim *> \verbatim
*> TB is COMPLEX array, dimension (LTB) *> TB is COMPLEX*16 array, dimension (LTB)
*> On exit, details of the LU factorization of the band matrix. *> On exit, details of the LU factorization of the band matrix.
*> \endverbatim *> \endverbatim
*> *>
@ -121,7 +121,7 @@
*> *>
*> \param[out] WORK *> \param[out] WORK
*> \verbatim *> \verbatim
*> WORK is COMPLEX workspace of size LWORK *> WORK is COMPLEX*16 workspace of size LWORK
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] LWORK *> \param[in] LWORK
@ -276,7 +276,7 @@
IF( UPPER ) THEN IF( UPPER ) THEN
* *
* ..................................................... * .....................................................
* Factorize A as L*D*L**T using the upper triangle of A * Factorize A as U**H*D*U using the upper triangle of A
* ..................................................... * .....................................................
* *
DO J = 0, NT-1 DO J = 0, NT-1
@ -453,12 +453,15 @@ c END IF
CALL ZSWAP( K-1, A( (J+1)*NB+1, I1 ), 1, CALL ZSWAP( K-1, A( (J+1)*NB+1, I1 ), 1,
$ A( (J+1)*NB+1, I2 ), 1 ) $ A( (J+1)*NB+1, I2 ), 1 )
* > Swap A(I1+1:M, I1) with A(I2, I1+1:M) * > Swap A(I1+1:M, I1) with A(I2, I1+1:M)
IF( I2.GT.(I1+1) ) THEN
CALL ZSWAP( I2-I1-1, A( I1, I1+1 ), LDA, CALL ZSWAP( I2-I1-1, A( I1, I1+1 ), LDA,
$ A( I1+1, I2 ), 1 ) $ A( I1+1, I2 ), 1 )
CALL ZLACGV( I2-I1, A( I1, I1+1 ), LDA )
CALL ZLACGV( I2-I1-1, A( I1+1, I2 ), 1 ) CALL ZLACGV( I2-I1-1, A( I1+1, I2 ), 1 )
END IF
CALL ZLACGV( I2-I1, A( I1, I1+1 ), LDA )
* > Swap A(I2+1:M, I1) with A(I2+1:M, I2) * > Swap A(I2+1:M, I1) with A(I2+1:M, I2)
CALL ZSWAP( N-I2, A( I1, I2+1 ), LDA, IF( I2.LT.N )
$ CALL ZSWAP( N-I2, A( I1, I2+1 ), LDA,
$ A( I2, I2+1 ), LDA ) $ A( I2, I2+1 ), LDA )
* > Swap A(I1, I1) with A(I2, I2) * > Swap A(I1, I1) with A(I2, I2)
PIV = A( I1, I1 ) PIV = A( I1, I1 )
@ -476,7 +479,7 @@ c END IF
ELSE ELSE
* *
* ..................................................... * .....................................................
* Factorize A as L*D*L**T using the lower triangle of A * Factorize A as L*D*L**H using the lower triangle of A
* ..................................................... * .....................................................
* *
DO J = 0, NT-1 DO J = 0, NT-1
@ -630,12 +633,15 @@ c END IF
CALL ZSWAP( K-1, A( I1, (J+1)*NB+1 ), LDA, CALL ZSWAP( K-1, A( I1, (J+1)*NB+1 ), LDA,
$ A( I2, (J+1)*NB+1 ), LDA ) $ A( I2, (J+1)*NB+1 ), LDA )
* > Swap A(I1+1:M, I1) with A(I2, I1+1:M) * > Swap A(I1+1:M, I1) with A(I2, I1+1:M)
IF( I2.GT.(I1+1) ) THEN
CALL ZSWAP( I2-I1-1, A( I1+1, I1 ), 1, CALL ZSWAP( I2-I1-1, A( I1+1, I1 ), 1,
$ A( I2, I1+1 ), LDA ) $ A( I2, I1+1 ), LDA )
CALL ZLACGV( I2-I1, A( I1+1, I1 ), 1 )
CALL ZLACGV( I2-I1-1, A( I2, I1+1 ), LDA ) CALL ZLACGV( I2-I1-1, A( I2, I1+1 ), LDA )
END IF
CALL ZLACGV( I2-I1, A( I1+1, I1 ), 1 )
* > Swap A(I2+1:M, I1) with A(I2+1:M, I2) * > Swap A(I2+1:M, I1) with A(I2+1:M, I2)
CALL ZSWAP( N-I2, A( I2+1, I1 ), 1, IF( I2.LT.N )
$ CALL ZSWAP( N-I2, A( I2+1, I1 ), 1,
$ A( I2+1, I2 ), 1 ) $ A( I2+1, I2 ), 1 )
* > Swap A(I1, I1) with A(I2, I2) * > Swap A(I1, I1) with A(I2, I2)
PIV = A( I1, I1 ) PIV = A( I1, I1 )

View File

@ -62,7 +62,7 @@
*> \param[in,out] A *> \param[in,out] A
*> \verbatim *> \verbatim
*> A is COMPLEX*16 array, dimension (LDA,N) *> A is COMPLEX*16 array, dimension (LDA,N)
*> On entry, the NB diagonal matrix D and the multipliers *> On entry, the block diagonal matrix D and the multipliers
*> used to obtain the factor U or L as computed by ZHETRF. *> used to obtain the factor U or L as computed by ZHETRF.
*> *>
*> On exit, if INFO = 0, the (symmetric) inverse of the original *> On exit, if INFO = 0, the (symmetric) inverse of the original
@ -82,7 +82,7 @@
*> \param[in] IPIV *> \param[in] IPIV
*> \verbatim *> \verbatim
*> IPIV is INTEGER array, dimension (N) *> IPIV is INTEGER array, dimension (N)
*> Details of the interchanges and the NB structure of D *> Details of the interchanges and the block structure of D
*> as determined by ZHETRF. *> as determined by ZHETRF.
*> \endverbatim *> \endverbatim
*> *>

View File

@ -38,8 +38,8 @@
*> \verbatim *> \verbatim
*> *>
*> ZHETRS_AA solves a system of linear equations A*X = B with a complex *> ZHETRS_AA solves a system of linear equations A*X = B with a complex
*> hermitian matrix A using the factorization A = U*T*U**H or *> hermitian matrix A using the factorization A = U**H*T*U or
*> A = L*T*L**T computed by ZHETRF_AA. *> A = L*T*L**H computed by ZHETRF_AA.
*> \endverbatim *> \endverbatim
* *
* Arguments: * Arguments:
@ -50,7 +50,7 @@
*> UPLO is CHARACTER*1 *> UPLO is CHARACTER*1
*> Specifies whether the details of the factorization are stored *> Specifies whether the details of the factorization are stored
*> as an upper or lower triangular matrix. *> as an upper or lower triangular matrix.
*> = 'U': Upper triangular, form is A = U*T*U**H; *> = 'U': Upper triangular, form is A = U**H*T*U;
*> = 'L': Lower triangular, form is A = L*T*L**H. *> = 'L': Lower triangular, form is A = L*T*L**H.
*> \endverbatim *> \endverbatim
*> *>
@ -98,14 +98,16 @@
*> The leading dimension of the array B. LDB >= max(1,N). *> The leading dimension of the array B. LDB >= max(1,N).
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] WORK *> \param[out] WORK
*> \verbatim *> \verbatim
*> WORK is DOUBLE array, dimension (MAX(1,LWORK)) *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] LWORK *> \param[in] LWORK
*> \verbatim *> \verbatim
*> LWORK is INTEGER, LWORK >= MAX(1,3*N-2). *> LWORK is INTEGER
*> The dimension of the array WORK. LWORK >= max(1,3*N-2).
*> \endverbatim
*> *>
*> \param[out] INFO *> \param[out] INFO
*> \verbatim *> \verbatim
@ -199,9 +201,13 @@
* *
IF( UPPER ) THEN IF( UPPER ) THEN
* *
* Solve A*X = B, where A = U*T*U**T. * Solve A*X = B, where A = U**H*T*U.
* *
* Pivot, P**T * B * 1) Forward substitution with U**H
*
IF( N.GT.1 ) THEN
*
* Pivot, P**T * B -> B
* *
DO K = 1, N DO K = 1, N
KP = IPIV( K ) KP = IPIV( K )
@ -209,12 +215,15 @@
$ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
END DO END DO
* *
* Compute (U \P**T * B) -> B [ (U \P**T * B) ] * Compute U**H \ B -> B [ (U**H \P**T * B) ]
* *
CALL ZTRSM('L', 'U', 'C', 'U', N-1, NRHS, ONE, A( 1, 2 ), LDA, CALL ZTRSM( 'L', 'U', 'C', 'U', N-1, NRHS, ONE, A( 1, 2 ),
$ B( 2, 1 ), LDB) $ LDA, B( 2, 1 ), LDB )
END IF
* *
* Compute T \ B -> B [ T \ (U \P**T * B) ] * 2) Solve with triangular matrix T
*
* Compute T \ B -> B [ T \ (U**H \P**T * B) ]
* *
CALL ZLACPY( 'F', 1, N, A(1, 1), LDA+1, WORK(N), 1 ) CALL ZLACPY( 'F', 1, N, A(1, 1), LDA+1, WORK(N), 1 )
IF( N.GT.1 ) THEN IF( N.GT.1 ) THEN
@ -225,24 +234,33 @@
CALL ZGTSV( N, NRHS, WORK(1), WORK(N), WORK(2*N), B, LDB, CALL ZGTSV( N, NRHS, WORK(1), WORK(N), WORK(2*N), B, LDB,
$ INFO ) $ INFO )
* *
* Compute (U**T \ B) -> B [ U**T \ (T \ (U \P**T * B) ) ] * 3) Backward substitution with U
* *
CALL ZTRSM( 'L', 'U', 'N', 'U', N-1, NRHS, ONE, A( 1, 2 ), LDA, IF( N.GT.1 ) THEN
$ B(2, 1), LDB)
* *
* Pivot, P * B [ P * (U**T \ (T \ (U \P**T * B) )) ] * Compute U \ B -> B [ U \ (T \ (U**H \P**T * B) ) ]
*
CALL ZTRSM( 'L', 'U', 'N', 'U', N-1, NRHS, ONE, A( 1, 2 ),
$ LDA, B(2, 1), LDB)
*
* Pivot, P * B [ P * (U**H \ (T \ (U \P**T * B) )) ]
* *
DO K = N, 1, -1 DO K = N, 1, -1
KP = IPIV( K ) KP = IPIV( K )
IF( KP.NE.K ) IF( KP.NE.K )
$ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
END DO END DO
END IF
* *
ELSE ELSE
* *
* Solve A*X = B, where A = L*T*L**T. * Solve A*X = B, where A = L*T*L**H.
* *
* Pivot, P**T * B * 1) Forward substitution with L
*
IF( N.GT.1 ) THEN
*
* Pivot, P**T * B -> B
* *
DO K = 1, N DO K = 1, N
KP = IPIV( K ) KP = IPIV( K )
@ -250,10 +268,13 @@
$ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
END DO END DO
* *
* Compute (L \P**T * B) -> B [ (L \P**T * B) ] * Compute L \ B -> B [ (L \P**T * B) ]
* *
CALL ZTRSM( 'L', 'L', 'N', 'U', N-1, NRHS, ONE, A( 2, 1 ), LDA, CALL ZTRSM( 'L', 'L', 'N', 'U', N-1, NRHS, ONE, A( 2, 1 ),
$ B(2, 1), LDB) $ LDA, B(2, 1), LDB)
END IF
*
* 2) Solve with triangular matrix T
* *
* Compute T \ B -> B [ T \ (L \P**T * B) ] * Compute T \ B -> B [ T \ (L \P**T * B) ]
* *
@ -266,18 +287,23 @@
CALL ZGTSV(N, NRHS, WORK(1), WORK(N), WORK(2*N), B, LDB, CALL ZGTSV(N, NRHS, WORK(1), WORK(N), WORK(2*N), B, LDB,
$ INFO) $ INFO)
* *
* Compute (L**T \ B) -> B [ L**T \ (T \ (L \P**T * B) ) ] * 3) Backward substitution with L**H
* *
CALL ZTRSM( 'L', 'L', 'C', 'U', N-1, NRHS, ONE, A( 2, 1 ), LDA, IF( N.GT.1 ) THEN
$ B( 2, 1 ), LDB)
* *
* Pivot, P * B [ P * (L**T \ (T \ (L \P**T * B) )) ] * Compute L**H \ B -> B [ L**H \ (T \ (L \P**T * B) ) ]
*
CALL ZTRSM( 'L', 'L', 'C', 'U', N-1, NRHS, ONE, A( 2, 1 ),
$ LDA, B( 2, 1 ), LDB)
*
* Pivot, P * B [ P * (L**H \ (T \ (L \P**T * B) )) ]
* *
DO K = N, 1, -1 DO K = N, 1, -1
KP = IPIV( K ) KP = IPIV( K )
IF( KP.NE.K ) IF( KP.NE.K )
$ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
END DO END DO
END IF
* *
END IF END IF
* *

View File

@ -38,8 +38,8 @@
*> \verbatim *> \verbatim
*> *>
*> ZHETRS_AA_2STAGE solves a system of linear equations A*X = B with a *> ZHETRS_AA_2STAGE solves a system of linear equations A*X = B with a
*> hermitian matrix A using the factorization A = U*T*U**T or *> hermitian matrix A using the factorization A = U**H*T*U or
*> A = L*T*L**T computed by ZHETRF_AA_2STAGE. *> A = L*T*L**H computed by ZHETRF_AA_2STAGE.
*> \endverbatim *> \endverbatim
* *
* Arguments: * Arguments:
@ -50,8 +50,8 @@
*> UPLO is CHARACTER*1 *> UPLO is CHARACTER*1
*> Specifies whether the details of the factorization are stored *> Specifies whether the details of the factorization are stored
*> as an upper or lower triangular matrix. *> as an upper or lower triangular matrix.
*> = 'U': Upper triangular, form is A = U*T*U**T; *> = 'U': Upper triangular, form is A = U**H*T*U;
*> = 'L': Lower triangular, form is A = L*T*L**T. *> = 'L': Lower triangular, form is A = L*T*L**H.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] N *> \param[in] N
@ -210,33 +210,33 @@
* *
IF( UPPER ) THEN IF( UPPER ) THEN
* *
* Solve A*X = B, where A = U*T*U**T. * Solve A*X = B, where A = U**H*T*U.
* *
IF( N.GT.NB ) THEN IF( N.GT.NB ) THEN
* *
* Pivot, P**T * B * Pivot, P**T * B -> B
* *
CALL ZLASWP( NRHS, B, LDB, NB+1, N, IPIV, 1 ) CALL ZLASWP( NRHS, B, LDB, NB+1, N, IPIV, 1 )
* *
* Compute (U**T \P**T * B) -> B [ (U**T \P**T * B) ] * Compute (U**H \ B) -> B [ (U**H \P**T * B) ]
* *
CALL ZTRSM( 'L', 'U', 'C', 'U', N-NB, NRHS, ONE, A(1, NB+1), CALL ZTRSM( 'L', 'U', 'C', 'U', N-NB, NRHS, ONE, A(1, NB+1),
$ LDA, B(NB+1, 1), LDB) $ LDA, B(NB+1, 1), LDB)
* *
END IF END IF
* *
* Compute T \ B -> B [ T \ (U**T \P**T * B) ] * Compute T \ B -> B [ T \ (U**H \P**T * B) ]
* *
CALL ZGBTRS( 'N', N, NB, NB, NRHS, TB, LDTB, IPIV2, B, LDB, CALL ZGBTRS( 'N', N, NB, NB, NRHS, TB, LDTB, IPIV2, B, LDB,
$ INFO) $ INFO)
IF( N.GT.NB ) THEN IF( N.GT.NB ) THEN
* *
* Compute (U \ B) -> B [ U \ (T \ (U**T \P**T * B) ) ] * Compute (U \ B) -> B [ U \ (T \ (U**H \P**T * B) ) ]
* *
CALL ZTRSM( 'L', 'U', 'N', 'U', N-NB, NRHS, ONE, A(1, NB+1), CALL ZTRSM( 'L', 'U', 'N', 'U', N-NB, NRHS, ONE, A(1, NB+1),
$ LDA, B(NB+1, 1), LDB) $ LDA, B(NB+1, 1), LDB)
* *
* Pivot, P * B [ P * (U \ (T \ (U**T \P**T * B) )) ] * Pivot, P * B -> B [ P * (U \ (T \ (U**H \P**T * B) )) ]
* *
CALL ZLASWP( NRHS, B, LDB, NB+1, N, IPIV, -1 ) CALL ZLASWP( NRHS, B, LDB, NB+1, N, IPIV, -1 )
* *
@ -244,15 +244,15 @@
* *
ELSE ELSE
* *
* Solve A*X = B, where A = L*T*L**T. * Solve A*X = B, where A = L*T*L**H.
* *
IF( N.GT.NB ) THEN IF( N.GT.NB ) THEN
* *
* Pivot, P**T * B * Pivot, P**T * B -> B
* *
CALL ZLASWP( NRHS, B, LDB, NB+1, N, IPIV, 1 ) CALL ZLASWP( NRHS, B, LDB, NB+1, N, IPIV, 1 )
* *
* Compute (L \P**T * B) -> B [ (L \P**T * B) ] * Compute (L \ B) -> B [ (L \P**T * B) ]
* *
CALL ZTRSM( 'L', 'L', 'N', 'U', N-NB, NRHS, ONE, A(NB+1, 1), CALL ZTRSM( 'L', 'L', 'N', 'U', N-NB, NRHS, ONE, A(NB+1, 1),
$ LDA, B(NB+1, 1), LDB) $ LDA, B(NB+1, 1), LDB)
@ -265,12 +265,12 @@
$ INFO) $ INFO)
IF( N.GT.NB ) THEN IF( N.GT.NB ) THEN
* *
* Compute (L**T \ B) -> B [ L**T \ (T \ (L \P**T * B) ) ] * Compute (L**H \ B) -> B [ L**H \ (T \ (L \P**T * B) ) ]
* *
CALL ZTRSM( 'L', 'L', 'C', 'U', N-NB, NRHS, ONE, A(NB+1, 1), CALL ZTRSM( 'L', 'L', 'C', 'U', N-NB, NRHS, ONE, A(NB+1, 1),
$ LDA, B(NB+1, 1), LDB) $ LDA, B(NB+1, 1), LDB)
* *
* Pivot, P * B [ P * (L**T \ (T \ (L \P**T * B) )) ] * Pivot, P * B -> B [ P * (L**H \ (T \ (L \P**T * B) )) ]
* *
CALL ZLASWP( NRHS, B, LDB, NB+1, N, IPIV, -1 ) CALL ZLASWP( NRHS, B, LDB, NB+1, N, IPIV, -1 )
* *

View File

@ -69,7 +69,7 @@
*> \param[in] N *> \param[in] N
*> \verbatim *> \verbatim
*> N is INTEGER *> N is INTEGER
*> The order of the matrix H. N .GE. 0. *> The order of the matrix H. N >= 0.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] ILO *> \param[in] ILO
@ -86,7 +86,7 @@
*> set by a previous call to ZGEBAL, and then passed to ZGEHRD *> set by a previous call to ZGEBAL, and then passed to ZGEHRD
*> when the matrix output by ZGEBAL is reduced to Hessenberg *> when the matrix output by ZGEBAL is reduced to Hessenberg
*> form. Otherwise ILO and IHI should be set to 1 and N *> form. Otherwise ILO and IHI should be set to 1 and N
*> respectively. If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N. *> respectively. If N > 0, then 1 <= ILO <= IHI <= N.
*> If N = 0, then ILO = 1 and IHI = 0. *> If N = 0, then ILO = 1 and IHI = 0.
*> \endverbatim *> \endverbatim
*> *>
@ -98,17 +98,17 @@
*> triangular matrix T from the Schur decomposition (the *> triangular matrix T from the Schur decomposition (the
*> Schur form). If INFO = 0 and JOB = 'E', the contents of *> Schur form). If INFO = 0 and JOB = 'E', the contents of
*> H are unspecified on exit. (The output value of H when *> H are unspecified on exit. (The output value of H when
*> INFO.GT.0 is given under the description of INFO below.) *> INFO > 0 is given under the description of INFO below.)
*> *>
*> Unlike earlier versions of ZHSEQR, this subroutine may *> Unlike earlier versions of ZHSEQR, this subroutine may
*> explicitly H(i,j) = 0 for i.GT.j and j = 1, 2, ... ILO-1 *> explicitly H(i,j) = 0 for i > j and j = 1, 2, ... ILO-1
*> or j = IHI+1, IHI+2, ... N. *> or j = IHI+1, IHI+2, ... N.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] LDH *> \param[in] LDH
*> \verbatim *> \verbatim
*> LDH is INTEGER *> LDH is INTEGER
*> The leading dimension of the array H. LDH .GE. max(1,N). *> The leading dimension of the array H. LDH >= max(1,N).
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] W *> \param[out] W
@ -131,7 +131,7 @@
*> if INFO = 0, Z contains Q*Z. *> if INFO = 0, Z contains Q*Z.
*> Normally Q is the unitary matrix generated by ZUNGHR *> Normally Q is the unitary matrix generated by ZUNGHR
*> after the call to ZGEHRD which formed the Hessenberg matrix *> after the call to ZGEHRD which formed the Hessenberg matrix
*> H. (The output value of Z when INFO.GT.0 is given under *> H. (The output value of Z when INFO > 0 is given under
*> the description of INFO below.) *> the description of INFO below.)
*> \endverbatim *> \endverbatim
*> *>
@ -139,7 +139,7 @@
*> \verbatim *> \verbatim
*> LDZ is INTEGER *> LDZ is INTEGER
*> The leading dimension of the array Z. if COMPZ = 'I' or *> The leading dimension of the array Z. if COMPZ = 'I' or
*> COMPZ = 'V', then LDZ.GE.MAX(1,N). Otherwize, LDZ.GE.1. *> COMPZ = 'V', then LDZ >= MAX(1,N). Otherwise, LDZ >= 1.
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] WORK *> \param[out] WORK
@ -152,7 +152,7 @@
*> \param[in] LWORK *> \param[in] LWORK
*> \verbatim *> \verbatim
*> LWORK is INTEGER *> LWORK is INTEGER
*> The dimension of the array WORK. LWORK .GE. max(1,N) *> The dimension of the array WORK. LWORK >= max(1,N)
*> is sufficient and delivers very good and sometimes *> is sufficient and delivers very good and sometimes
*> optimal performance. However, LWORK as large as 11*N *> optimal performance. However, LWORK as large as 11*N
*> may be required for optimal performance. A workspace *> may be required for optimal performance. A workspace
@ -171,20 +171,20 @@
*> \verbatim *> \verbatim
*> INFO is INTEGER *> INFO is INTEGER
*> = 0: successful exit *> = 0: successful exit
*> .LT. 0: if INFO = -i, the i-th argument had an illegal *> < 0: if INFO = -i, the i-th argument had an illegal
*> value *> value
*> .GT. 0: if INFO = i, ZHSEQR failed to compute all of *> > 0: if INFO = i, ZHSEQR failed to compute all of
*> the eigenvalues. Elements 1:ilo-1 and i+1:n of WR *> the eigenvalues. Elements 1:ilo-1 and i+1:n of W
*> and WI contain those eigenvalues which have been *> contain those eigenvalues which have been
*> successfully computed. (Failures are rare.) *> successfully computed. (Failures are rare.)
*> *>
*> If INFO .GT. 0 and JOB = 'E', then on exit, the *> If INFO > 0 and JOB = 'E', then on exit, the
*> remaining unconverged eigenvalues are the eigen- *> remaining unconverged eigenvalues are the eigen-
*> values of the upper Hessenberg matrix rows and *> values of the upper Hessenberg matrix rows and
*> columns ILO through INFO of the final, output *> columns ILO through INFO of the final, output
*> value of H. *> value of H.
*> *>
*> If INFO .GT. 0 and JOB = 'S', then on exit *> If INFO > 0 and JOB = 'S', then on exit
*> *>
*> (*) (initial value of H)*U = U*(final value of H) *> (*) (initial value of H)*U = U*(final value of H)
*> *>
@ -192,19 +192,19 @@
*> value of H is upper Hessenberg and triangular in *> value of H is upper Hessenberg and triangular in
*> rows and columns INFO+1 through IHI. *> rows and columns INFO+1 through IHI.
*> *>
*> If INFO .GT. 0 and COMPZ = 'V', then on exit *> If INFO > 0 and COMPZ = 'V', then on exit
*> *>
*> (final value of Z) = (initial value of Z)*U *> (final value of Z) = (initial value of Z)*U
*> *>
*> where U is the unitary matrix in (*) (regard- *> where U is the unitary matrix in (*) (regard-
*> less of the value of JOB.) *> less of the value of JOB.)
*> *>
*> If INFO .GT. 0 and COMPZ = 'I', then on exit *> If INFO > 0 and COMPZ = 'I', then on exit
*> (final value of Z) = U *> (final value of Z) = U
*> where U is the unitary matrix in (*) (regard- *> where U is the unitary matrix in (*) (regard-
*> less of the value of JOB.) *> less of the value of JOB.)
*> *>
*> If INFO .GT. 0 and COMPZ = 'N', then Z is not *> If INFO > 0 and COMPZ = 'N', then Z is not
*> accessed. *> accessed.
*> \endverbatim *> \endverbatim
* *
@ -244,8 +244,8 @@
*> This depends on ILO, IHI and NS. NS is the *> This depends on ILO, IHI and NS. NS is the
*> number of simultaneous shifts returned *> number of simultaneous shifts returned
*> by ILAENV(ISPEC=15). (See ISPEC=15 below.) *> by ILAENV(ISPEC=15). (See ISPEC=15 below.)
*> The default for (IHI-ILO+1).LE.500 is NS. *> The default for (IHI-ILO+1) <= 500 is NS.
*> The default for (IHI-ILO+1).GT.500 is 3*NS/2. *> The default for (IHI-ILO+1) > 500 is 3*NS/2.
*> *>
*> ISPEC=14: Nibble crossover point. (See IPARMQ for *> ISPEC=14: Nibble crossover point. (See IPARMQ for
*> details.) Default: 14% of deflation window *> details.) Default: 14% of deflation window
@ -323,8 +323,8 @@
PARAMETER ( NTINY = 11 ) PARAMETER ( NTINY = 11 )
* *
* ==== NL allocates some local workspace to help small matrices * ==== NL allocates some local workspace to help small matrices
* . through a rare ZLAHQR failure. NL .GT. NTINY = 11 is * . through a rare ZLAHQR failure. NL > NTINY = 11 is
* . required and NL .LE. NMIN = ILAENV(ISPEC=12,...) is recom- * . required and NL <= NMIN = ILAENV(ISPEC=12,...) is recom-
* . mended. (The default value of NMIN is 75.) Using NL = 49 * . mended. (The default value of NMIN is 75.) Using NL = 49
* . allows up to six simultaneous shifts and a 16-by-16 * . allows up to six simultaneous shifts and a 16-by-16
* . deflation window. ==== * . deflation window. ====

View File

@ -133,13 +133,13 @@
*> i > 0: The ith argument is invalid. *> i > 0: The ith argument is invalid.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] WORK *> \param[out] WORK
*> \verbatim *> \verbatim
*> WORK is COMPLEX*16 array, dimension (2*N). *> WORK is COMPLEX*16 array, dimension (2*N).
*> Workspace. *> Workspace.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] RWORK *> \param[out] RWORK
*> \verbatim *> \verbatim
*> RWORK is DOUBLE PRECISION array, dimension (N). *> RWORK is DOUBLE PRECISION array, dimension (N).
*> Workspace. *> Workspace.

View File

@ -126,13 +126,13 @@
*> i > 0: The ith argument is invalid. *> i > 0: The ith argument is invalid.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] WORK *> \param[out] WORK
*> \verbatim *> \verbatim
*> WORK is COMPLEX*16 array, dimension (2*N). *> WORK is COMPLEX*16 array, dimension (2*N).
*> Workspace. *> Workspace.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] RWORK *> \param[out] RWORK
*> \verbatim *> \verbatim
*> RWORK is DOUBLE PRECISION array, dimension (N). *> RWORK is DOUBLE PRECISION array, dimension (N).
*> Workspace. *> Workspace.

View File

@ -65,19 +65,19 @@
*> \verbatim *> \verbatim
*> PREC_TYPE is INTEGER *> PREC_TYPE is INTEGER
*> Specifies the intermediate precision to be used in refinement. *> Specifies the intermediate precision to be used in refinement.
*> The value is defined by ILAPREC(P) where P is a CHARACTER and *> The value is defined by ILAPREC(P) where P is a CHARACTER and P
*> P = 'S': Single *> = 'S': Single
*> = 'D': Double *> = 'D': Double
*> = 'I': Indigenous *> = 'I': Indigenous
*> = 'X', 'E': Extra *> = 'X' or 'E': Extra
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] TRANS_TYPE *> \param[in] TRANS_TYPE
*> \verbatim *> \verbatim
*> TRANS_TYPE is INTEGER *> TRANS_TYPE is INTEGER
*> Specifies the transposition operation on A. *> Specifies the transposition operation on A.
*> The value is defined by ILATRANS(T) where T is a CHARACTER and *> The value is defined by ILATRANS(T) where T is a CHARACTER and T
*> T = 'N': No transpose *> = 'N': No transpose
*> = 'T': Transpose *> = 'T': Transpose
*> = 'C': Conjugate transpose *> = 'C': Conjugate transpose
*> \endverbatim *> \endverbatim
@ -269,7 +269,7 @@
*> information as described below. There currently are up to three *> information as described below. There currently are up to three
*> pieces of information returned for each right-hand side. If *> pieces of information returned for each right-hand side. If
*> componentwise accuracy is not requested (PARAMS(3) = 0.0), then *> componentwise accuracy is not requested (PARAMS(3) = 0.0), then
*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most *> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS < 3, then at most
*> the first (:,N_ERR_BNDS) entries are returned. *> the first (:,N_ERR_BNDS) entries are returned.
*> *>
*> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith *> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith

View File

@ -22,7 +22,7 @@
* LDAF, IPIV, C, CAPPLY, * LDAF, IPIV, C, CAPPLY,
* INFO, WORK, RWORK ) * INFO, WORK, RWORK )
* *
* .. Scalar Aguments .. * .. Scalar Arguments ..
* CHARACTER TRANS * CHARACTER TRANS
* LOGICAL CAPPLY * LOGICAL CAPPLY
* INTEGER N, LDA, LDAF, INFO * INTEGER N, LDA, LDAF, INFO
@ -114,13 +114,13 @@
*> i > 0: The ith argument is invalid. *> i > 0: The ith argument is invalid.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] WORK *> \param[out] WORK
*> \verbatim *> \verbatim
*> WORK is COMPLEX*16 array, dimension (2*N). *> WORK is COMPLEX*16 array, dimension (2*N).
*> Workspace. *> Workspace.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] RWORK *> \param[out] RWORK
*> \verbatim *> \verbatim
*> RWORK is DOUBLE PRECISION array, dimension (N). *> RWORK is DOUBLE PRECISION array, dimension (N).
*> Workspace. *> Workspace.
@ -148,7 +148,7 @@
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* December 2016 * December 2016
* *
* .. Scalar Aguments .. * .. Scalar Arguments ..
CHARACTER TRANS CHARACTER TRANS
LOGICAL CAPPLY LOGICAL CAPPLY
INTEGER N, LDA, LDAF, INFO INTEGER N, LDA, LDAF, INFO

View File

@ -107,13 +107,13 @@
*> i > 0: The ith argument is invalid. *> i > 0: The ith argument is invalid.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] WORK *> \param[out] WORK
*> \verbatim *> \verbatim
*> WORK is COMPLEX*16 array, dimension (2*N). *> WORK is COMPLEX*16 array, dimension (2*N).
*> Workspace. *> Workspace.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] RWORK *> \param[out] RWORK
*> \verbatim *> \verbatim
*> RWORK is DOUBLE PRECISION array, dimension (N). *> RWORK is DOUBLE PRECISION array, dimension (N).
*> Workspace. *> Workspace.

View File

@ -64,19 +64,19 @@
*> \verbatim *> \verbatim
*> PREC_TYPE is INTEGER *> PREC_TYPE is INTEGER
*> Specifies the intermediate precision to be used in refinement. *> Specifies the intermediate precision to be used in refinement.
*> The value is defined by ILAPREC(P) where P is a CHARACTER and *> The value is defined by ILAPREC(P) where P is a CHARACTER and P
*> P = 'S': Single *> = 'S': Single
*> = 'D': Double *> = 'D': Double
*> = 'I': Indigenous *> = 'I': Indigenous
*> = 'X', 'E': Extra *> = 'X' or 'E': Extra
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] TRANS_TYPE *> \param[in] TRANS_TYPE
*> \verbatim *> \verbatim
*> TRANS_TYPE is INTEGER *> TRANS_TYPE is INTEGER
*> Specifies the transposition operation on A. *> Specifies the transposition operation on A.
*> The value is defined by ILATRANS(T) where T is a CHARACTER and *> The value is defined by ILATRANS(T) where T is a CHARACTER and T
*> T = 'N': No transpose *> = 'N': No transpose
*> = 'T': Transpose *> = 'T': Transpose
*> = 'C': Conjugate transpose *> = 'C': Conjugate transpose
*> \endverbatim *> \endverbatim
@ -256,7 +256,7 @@
*> information as described below. There currently are up to three *> information as described below. There currently are up to three
*> pieces of information returned for each right-hand side. If *> pieces of information returned for each right-hand side. If
*> componentwise accuracy is not requested (PARAMS(3) = 0.0), then *> componentwise accuracy is not requested (PARAMS(3) = 0.0), then
*> ERRS_C is not accessed. If N_ERR_BNDS .LT. 3, then at most *> ERRS_C is not accessed. If N_ERR_BNDS < 3, then at most
*> the first (:,N_ERR_BNDS) entries are returned. *> the first (:,N_ERR_BNDS) entries are returned.
*> *>
*> The first index in ERRS_C(i,:) corresponds to the ith *> The first index in ERRS_C(i,:) corresponds to the ith

View File

@ -111,13 +111,13 @@
*> i > 0: The ith argument is invalid. *> i > 0: The ith argument is invalid.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] WORK *> \param[out] WORK
*> \verbatim *> \verbatim
*> WORK is COMPLEX*16 array, dimension (2*N). *> WORK is COMPLEX*16 array, dimension (2*N).
*> Workspace. *> Workspace.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] RWORK *> \param[out] RWORK
*> \verbatim *> \verbatim
*> RWORK is DOUBLE PRECISION array, dimension (N). *> RWORK is DOUBLE PRECISION array, dimension (N).
*> Workspace. *> Workspace.

View File

@ -104,13 +104,13 @@
*> i > 0: The ith argument is invalid. *> i > 0: The ith argument is invalid.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] WORK *> \param[out] WORK
*> \verbatim *> \verbatim
*> WORK is COMPLEX*16 array, dimension (2*N). *> WORK is COMPLEX*16 array, dimension (2*N).
*> Workspace. *> Workspace.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] RWORK *> \param[out] RWORK
*> \verbatim *> \verbatim
*> RWORK is DOUBLE PRECISION array, dimension (N). *> RWORK is DOUBLE PRECISION array, dimension (N).
*> Workspace. *> Workspace.

View File

@ -66,11 +66,11 @@
*> \verbatim *> \verbatim
*> PREC_TYPE is INTEGER *> PREC_TYPE is INTEGER
*> Specifies the intermediate precision to be used in refinement. *> Specifies the intermediate precision to be used in refinement.
*> The value is defined by ILAPREC(P) where P is a CHARACTER and *> The value is defined by ILAPREC(P) where P is a CHARACTER and P
*> P = 'S': Single *> = 'S': Single
*> = 'D': Double *> = 'D': Double
*> = 'I': Indigenous *> = 'I': Indigenous
*> = 'X', 'E': Extra *> = 'X' or 'E': Extra
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] UPLO *> \param[in] UPLO
@ -254,7 +254,7 @@
*> information as described below. There currently are up to three *> information as described below. There currently are up to three
*> pieces of information returned for each right-hand side. If *> pieces of information returned for each right-hand side. If
*> componentwise accuracy is not requested (PARAMS(3) = 0.0), then *> componentwise accuracy is not requested (PARAMS(3) = 0.0), then
*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most *> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS < 3, then at most
*> the first (:,N_ERR_BNDS) entries are returned. *> the first (:,N_ERR_BNDS) entries are returned.
*> *>
*> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith *> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith

View File

@ -102,7 +102,7 @@
*> as determined by ZHETRF. *> as determined by ZHETRF.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] WORK *> \param[out] WORK
*> \verbatim *> \verbatim
*> WORK is DOUBLE PRECISION array, dimension (2*N) *> WORK is DOUBLE PRECISION array, dimension (2*N)
*> \endverbatim *> \endverbatim

View File

@ -103,13 +103,13 @@
*> i > 0: The ith argument is invalid. *> i > 0: The ith argument is invalid.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] WORK *> \param[out] WORK
*> \verbatim *> \verbatim
*> WORK is COMPLEX*16 array, dimension (2*N). *> WORK is COMPLEX*16 array, dimension (2*N).
*> Workspace. *> Workspace.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] RWORK *> \param[out] RWORK
*> \verbatim *> \verbatim
*> RWORK is DOUBLE PRECISION array, dimension (N). *> RWORK is DOUBLE PRECISION array, dimension (N).
*> Workspace. *> Workspace.

View File

@ -96,13 +96,13 @@
*> i > 0: The ith argument is invalid. *> i > 0: The ith argument is invalid.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] WORK *> \param[out] WORK
*> \verbatim *> \verbatim
*> WORK is COMPLEX*16 array, dimension (2*N). *> WORK is COMPLEX*16 array, dimension (2*N).
*> Workspace. *> Workspace.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] RWORK *> \param[out] RWORK
*> \verbatim *> \verbatim
*> RWORK is DOUBLE PRECISION array, dimension (N). *> RWORK is DOUBLE PRECISION array, dimension (N).
*> Workspace. *> Workspace.

View File

@ -65,11 +65,11 @@
*> \verbatim *> \verbatim
*> PREC_TYPE is INTEGER *> PREC_TYPE is INTEGER
*> Specifies the intermediate precision to be used in refinement. *> Specifies the intermediate precision to be used in refinement.
*> The value is defined by ILAPREC(P) where P is a CHARACTER and *> The value is defined by ILAPREC(P) where P is a CHARACTER and P
*> P = 'S': Single *> = 'S': Single
*> = 'D': Double *> = 'D': Double
*> = 'I': Indigenous *> = 'I': Indigenous
*> = 'X', 'E': Extra *> = 'X' or 'E': Extra
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] UPLO *> \param[in] UPLO
@ -246,7 +246,7 @@
*> information as described below. There currently are up to three *> information as described below. There currently are up to three
*> pieces of information returned for each right-hand side. If *> pieces of information returned for each right-hand side. If
*> componentwise accuracy is not requested (PARAMS(3) = 0.0), then *> componentwise accuracy is not requested (PARAMS(3) = 0.0), then
*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most *> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS < 3, then at most
*> the first (:,N_ERR_BNDS) entries are returned. *> the first (:,N_ERR_BNDS) entries are returned.
*> *>
*> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith *> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith

View File

@ -86,7 +86,7 @@
*> The leading dimension of the array AF. LDAF >= max(1,N). *> The leading dimension of the array AF. LDAF >= max(1,N).
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] WORK *> \param[out] WORK
*> \verbatim *> \verbatim
*> WORK is DOUBLE PRECISION array, dimension (2*N) *> WORK is DOUBLE PRECISION array, dimension (2*N)
*> \endverbatim *> \endverbatim

View File

@ -111,13 +111,13 @@
*> i > 0: The ith argument is invalid. *> i > 0: The ith argument is invalid.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] WORK *> \param[out] WORK
*> \verbatim *> \verbatim
*> WORK is COMPLEX*16 array, dimension (2*N). *> WORK is COMPLEX*16 array, dimension (2*N).
*> Workspace. *> Workspace.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] RWORK *> \param[out] RWORK
*> \verbatim *> \verbatim
*> RWORK is DOUBLE PRECISION array, dimension (N). *> RWORK is DOUBLE PRECISION array, dimension (N).
*> Workspace. *> Workspace.

View File

@ -104,13 +104,13 @@
*> i > 0: The ith argument is invalid. *> i > 0: The ith argument is invalid.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] WORK *> \param[out] WORK
*> \verbatim *> \verbatim
*> WORK is COMPLEX*16 array, dimension (2*N). *> WORK is COMPLEX*16 array, dimension (2*N).
*> Workspace. *> Workspace.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] RWORK *> \param[out] RWORK
*> \verbatim *> \verbatim
*> RWORK is DOUBLE PRECISION array, dimension (N). *> RWORK is DOUBLE PRECISION array, dimension (N).
*> Workspace. *> Workspace.

View File

@ -66,11 +66,11 @@
*> \verbatim *> \verbatim
*> PREC_TYPE is INTEGER *> PREC_TYPE is INTEGER
*> Specifies the intermediate precision to be used in refinement. *> Specifies the intermediate precision to be used in refinement.
*> The value is defined by ILAPREC(P) where P is a CHARACTER and *> The value is defined by ILAPREC(P) where P is a CHARACTER and P
*> P = 'S': Single *> = 'S': Single
*> = 'D': Double *> = 'D': Double
*> = 'I': Indigenous *> = 'I': Indigenous
*> = 'X', 'E': Extra *> = 'X' or 'E': Extra
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] UPLO *> \param[in] UPLO
@ -254,7 +254,7 @@
*> information as described below. There currently are up to three *> information as described below. There currently are up to three
*> pieces of information returned for each right-hand side. If *> pieces of information returned for each right-hand side. If
*> componentwise accuracy is not requested (PARAMS(3) = 0.0), then *> componentwise accuracy is not requested (PARAMS(3) = 0.0), then
*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most *> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS < 3, then at most
*> the first (:,N_ERR_BNDS) entries are returned. *> the first (:,N_ERR_BNDS) entries are returned.
*> *>
*> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith *> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith

View File

@ -102,7 +102,7 @@
*> as determined by ZSYTRF. *> as determined by ZSYTRF.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] WORK *> \param[out] WORK
*> \verbatim *> \verbatim
*> WORK is DOUBLE PRECISION array, dimension (2*N) *> WORK is DOUBLE PRECISION array, dimension (2*N)
*> \endverbatim *> \endverbatim

View File

@ -36,7 +36,7 @@
*> ZLA_WWADDW adds a vector W into a doubled-single vector (X, Y). *> ZLA_WWADDW adds a vector W into a doubled-single vector (X, Y).
*> *>
*> This works for all extant IBM's hex and binary floating point *> This works for all extant IBM's hex and binary floating point
*> arithmetics, but not for decimal. *> arithmetic, but not for decimal.
*> \endverbatim *> \endverbatim
* *
* Arguments: * Arguments:

View File

@ -288,7 +288,8 @@
* *
* Swap A(I1, I2+1:N) with A(I2, I2+1:N) * Swap A(I1, I2+1:N) with A(I2, I2+1:N)
* *
CALL ZSWAP( M-I2, A( J1+I1-1, I2+1 ), LDA, IF( I2.LT.M )
$ CALL ZSWAP( M-I2, A( J1+I1-1, I2+1 ), LDA,
$ A( J1+I2-1, I2+1 ), LDA ) $ A( J1+I2-1, I2+1 ), LDA )
* *
* Swap A(I1, I1) with A(I2,I2) * Swap A(I1, I1) with A(I2,I2)
@ -329,6 +330,7 @@
* Compute L(J+2, J+1) = WORK( 3:N ) / T(J, J+1), * Compute L(J+2, J+1) = WORK( 3:N ) / T(J, J+1),
* where A(J, J+1) = T(J, J+1) and A(J+2:N, J) = L(J+2:N, J+1) * where A(J, J+1) = T(J, J+1) and A(J+2:N, J) = L(J+2:N, J+1)
* *
IF( J.LT.(M-1) ) THEN
IF( A( K, J+1 ).NE.ZERO ) THEN IF( A( K, J+1 ).NE.ZERO ) THEN
ALPHA = ONE / A( K, J+1 ) ALPHA = ONE / A( K, J+1 )
CALL ZCOPY( M-J-1, WORK( 3 ), 1, A( K, J+2 ), LDA ) CALL ZCOPY( M-J-1, WORK( 3 ), 1, A( K, J+2 ), LDA )
@ -338,6 +340,7 @@
$ A( K, J+2 ), LDA) $ A( K, J+2 ), LDA)
END IF END IF
END IF END IF
END IF
J = J + 1 J = J + 1
GO TO 10 GO TO 10
20 CONTINUE 20 CONTINUE
@ -440,7 +443,8 @@
* *
* Swap A(I2+1:N, I1) with A(I2+1:N, I2) * Swap A(I2+1:N, I1) with A(I2+1:N, I2)
* *
CALL ZSWAP( M-I2, A( I2+1, J1+I1-1 ), 1, IF( I2.LT.M )
$ CALL ZSWAP( M-I2, A( I2+1, J1+I1-1 ), 1,
$ A( I2+1, J1+I2-1 ), 1 ) $ A( I2+1, J1+I2-1 ), 1 )
* *
* Swap A(I1, I1) with A(I2, I2) * Swap A(I1, I1) with A(I2, I2)
@ -481,6 +485,7 @@
* Compute L(J+2, J+1) = WORK( 3:N ) / T(J, J+1), * Compute L(J+2, J+1) = WORK( 3:N ) / T(J, J+1),
* where A(J, J+1) = T(J, J+1) and A(J+2:N, J) = L(J+2:N, J+1) * where A(J, J+1) = T(J, J+1) and A(J+2:N, J) = L(J+2:N, J+1)
* *
IF( J.LT.(M-1) ) THEN
IF( A( J+1, K ).NE.ZERO ) THEN IF( A( J+1, K ).NE.ZERO ) THEN
ALPHA = ONE / A( J+1, K ) ALPHA = ONE / A( J+1, K )
CALL ZCOPY( M-J-1, WORK( 3 ), 1, A( J+2, K ), 1 ) CALL ZCOPY( M-J-1, WORK( 3 ), 1, A( J+2, K ), 1 )
@ -490,6 +495,7 @@
$ A( J+2, K ), LDA ) $ A( J+2, K ), LDA )
END IF END IF
END IF END IF
END IF
J = J + 1 J = J + 1
GO TO 30 GO TO 30
40 CONTINUE 40 CONTINUE

View File

@ -331,7 +331,7 @@
* Factorize the trailing columns of A using the upper triangle * Factorize the trailing columns of A using the upper triangle
* of A and working backwards, and compute the matrix W = U12*D * of A and working backwards, and compute the matrix W = U12*D
* for use in updating A11 (note that conjg(W) is actually stored) * for use in updating A11 (note that conjg(W) is actually stored)
* Initilize the first entry of array E, where superdiagonal * Initialize the first entry of array E, where superdiagonal
* elements of D are stored * elements of D are stored
* *
E( 1 ) = CZERO E( 1 ) = CZERO
@ -789,7 +789,7 @@
* of A and working forwards, and compute the matrix W = L21*D * of A and working forwards, and compute the matrix W = L21*D
* for use in updating A22 (note that conjg(W) is actually stored) * for use in updating A22 (note that conjg(W) is actually stored)
* *
* Initilize the unused last entry of the subdiagonal array E. * Initialize the unused last entry of the subdiagonal array E.
* *
E( N ) = CZERO E( N ) = CZERO
* *

View File

@ -139,25 +139,25 @@
*> \verbatim *> \verbatim
*> INFO is INTEGER *> INFO is INTEGER
*> = 0: successful exit *> = 0: successful exit
*> .GT. 0: if INFO = i, ZLAHQR failed to compute all the *> > 0: if INFO = i, ZLAHQR failed to compute all the
*> eigenvalues ILO to IHI in a total of 30 iterations *> eigenvalues ILO to IHI in a total of 30 iterations
*> per eigenvalue; elements i+1:ihi of W contain *> per eigenvalue; elements i+1:ihi of W contain
*> those eigenvalues which have been successfully *> those eigenvalues which have been successfully
*> computed. *> computed.
*> *>
*> If INFO .GT. 0 and WANTT is .FALSE., then on exit, *> If INFO > 0 and WANTT is .FALSE., then on exit,
*> the remaining unconverged eigenvalues are the *> the remaining unconverged eigenvalues are the
*> eigenvalues of the upper Hessenberg matrix *> eigenvalues of the upper Hessenberg matrix
*> rows and columns ILO thorugh INFO of the final, *> rows and columns ILO through INFO of the final,
*> output value of H. *> output value of H.
*> *>
*> If INFO .GT. 0 and WANTT is .TRUE., then on exit *> If INFO > 0 and WANTT is .TRUE., then on exit
*> (*) (initial value of H)*U = U*(final value of H) *> (*) (initial value of H)*U = U*(final value of H)
*> where U is an orthognal matrix. The final *> where U is an orthogonal matrix. The final
*> value of H is upper Hessenberg and triangular in *> value of H is upper Hessenberg and triangular in
*> rows and columns INFO+1 through IHI. *> rows and columns INFO+1 through IHI.
*> *>
*> If INFO .GT. 0 and WANTZ is .TRUE., then on exit *> If INFO > 0 and WANTZ is .TRUE., then on exit
*> (final value of Z) = (initial value of Z)*U *> (final value of Z) = (initial value of Z)*U
*> where U is the orthogonal matrix in (*) *> where U is the orthogonal matrix in (*)
*> (regardless of the value of WANTT.) *> (regardless of the value of WANTT.)

View File

@ -1,3 +1,4 @@
*> \brief \b ZLAMSWLQ
* *
* Definition: * Definition:
* =========== * ===========

View File

@ -1,3 +1,4 @@
*> \brief \b ZLAMTSQR
* *
* Definition: * Definition:
* =========== * ===========

View File

@ -130,6 +130,7 @@
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* December 2016 * December 2016
* *
IMPLICIT NONE
* .. Scalar Arguments .. * .. Scalar Arguments ..
CHARACTER NORM CHARACTER NORM
INTEGER KL, KU, LDAB, N INTEGER KL, KU, LDAB, N
@ -147,14 +148,17 @@
* .. * ..
* .. Local Scalars .. * .. Local Scalars ..
INTEGER I, J, K, L INTEGER I, J, K, L
DOUBLE PRECISION SCALE, SUM, VALUE, TEMP DOUBLE PRECISION SUM, VALUE, TEMP
* ..
* .. Local Arrays ..
DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 )
* .. * ..
* .. External Functions .. * .. External Functions ..
LOGICAL LSAME, DISNAN LOGICAL LSAME, DISNAN
EXTERNAL LSAME, DISNAN EXTERNAL LSAME, DISNAN
* .. * ..
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL ZLASSQ EXTERNAL ZLASSQ, DCOMBSSQ
* .. * ..
* .. Intrinsic Functions .. * .. Intrinsic Functions ..
INTRINSIC ABS, MAX, MIN, SQRT INTRINSIC ABS, MAX, MIN, SQRT
@ -207,15 +211,22 @@
ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN
* *
* Find normF(A). * Find normF(A).
* SSQ(1) is scale
* SSQ(2) is sum-of-squares
* For better accuracy, sum each column separately.
* *
SCALE = ZERO SSQ( 1 ) = ZERO
SUM = ONE SSQ( 2 ) = ONE
DO 90 J = 1, N DO 90 J = 1, N
L = MAX( 1, J-KU ) L = MAX( 1, J-KU )
K = KU + 1 - J + L K = KU + 1 - J + L
CALL ZLASSQ( MIN( N, J+KL )-L+1, AB( K, J ), 1, SCALE, SUM ) COLSSQ( 1 ) = ZERO
COLSSQ( 2 ) = ONE
CALL ZLASSQ( MIN( N, J+KL )-L+1, AB( K, J ), 1,
$ COLSSQ( 1 ), COLSSQ( 2 ) )
CALL DCOMBSSQ( SSQ, COLSSQ )
90 CONTINUE 90 CONTINUE
VALUE = SCALE*SQRT( SUM ) VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) )
END IF END IF
* *
ZLANGB = VALUE ZLANGB = VALUE

View File

@ -120,6 +120,7 @@
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* December 2016 * December 2016
* *
IMPLICIT NONE
* .. Scalar Arguments .. * .. Scalar Arguments ..
CHARACTER NORM CHARACTER NORM
INTEGER LDA, M, N INTEGER LDA, M, N
@ -137,14 +138,17 @@
* .. * ..
* .. Local Scalars .. * .. Local Scalars ..
INTEGER I, J INTEGER I, J
DOUBLE PRECISION SCALE, SUM, VALUE, TEMP DOUBLE PRECISION SUM, VALUE, TEMP
* ..
* .. Local Arrays ..
DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 )
* .. * ..
* .. External Functions .. * .. External Functions ..
LOGICAL LSAME, DISNAN LOGICAL LSAME, DISNAN
EXTERNAL LSAME, DISNAN EXTERNAL LSAME, DISNAN
* .. * ..
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL ZLASSQ EXTERNAL ZLASSQ, DCOMBSSQ
* .. * ..
* .. Intrinsic Functions .. * .. Intrinsic Functions ..
INTRINSIC ABS, MIN, SQRT INTRINSIC ABS, MIN, SQRT
@ -196,13 +200,19 @@
ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN
* *
* Find normF(A). * Find normF(A).
* SSQ(1) is scale
* SSQ(2) is sum-of-squares
* For better accuracy, sum each column separately.
* *
SCALE = ZERO SSQ( 1 ) = ZERO
SUM = ONE SSQ( 2 ) = ONE
DO 90 J = 1, N DO 90 J = 1, N
CALL ZLASSQ( M, A( 1, J ), 1, SCALE, SUM ) COLSSQ( 1 ) = ZERO
COLSSQ( 2 ) = ONE
CALL ZLASSQ( M, A( 1, J ), 1, COLSSQ( 1 ), COLSSQ( 2 ) )
CALL DCOMBSSQ( SSQ, COLSSQ )
90 CONTINUE 90 CONTINUE
VALUE = SCALE*SQRT( SUM ) VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) )
END IF END IF
* *
ZLANGE = VALUE ZLANGE = VALUE

View File

@ -137,6 +137,7 @@
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* December 2016 * December 2016
* *
IMPLICIT NONE
* .. Scalar Arguments .. * .. Scalar Arguments ..
CHARACTER NORM, UPLO CHARACTER NORM, UPLO
INTEGER K, LDAB, N INTEGER K, LDAB, N
@ -154,14 +155,17 @@
* .. * ..
* .. Local Scalars .. * .. Local Scalars ..
INTEGER I, J, L INTEGER I, J, L
DOUBLE PRECISION ABSA, SCALE, SUM, VALUE DOUBLE PRECISION ABSA, SUM, VALUE
* ..
* .. Local Arrays ..
DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 )
* .. * ..
* .. External Functions .. * .. External Functions ..
LOGICAL LSAME, DISNAN LOGICAL LSAME, DISNAN
EXTERNAL LSAME, DISNAN EXTERNAL LSAME, DISNAN
* .. * ..
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL ZLASSQ EXTERNAL ZLASSQ, DCOMBSSQ
* .. * ..
* .. Intrinsic Functions .. * .. Intrinsic Functions ..
INTRINSIC ABS, DBLE, MAX, MIN, SQRT INTRINSIC ABS, DBLE, MAX, MIN, SQRT
@ -233,39 +237,57 @@
ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN
* *
* Find normF(A). * Find normF(A).
* SSQ(1) is scale
* SSQ(2) is sum-of-squares
* For better accuracy, sum each column separately.
*
SSQ( 1 ) = ZERO
SSQ( 2 ) = ONE
*
* Sum off-diagonals
* *
SCALE = ZERO
SUM = ONE
IF( K.GT.0 ) THEN IF( K.GT.0 ) THEN
IF( LSAME( UPLO, 'U' ) ) THEN IF( LSAME( UPLO, 'U' ) ) THEN
DO 110 J = 2, N DO 110 J = 2, N
COLSSQ( 1 ) = ZERO
COLSSQ( 2 ) = ONE
CALL ZLASSQ( MIN( J-1, K ), AB( MAX( K+2-J, 1 ), J ), CALL ZLASSQ( MIN( J-1, K ), AB( MAX( K+2-J, 1 ), J ),
$ 1, SCALE, SUM ) $ 1, COLSSQ( 1 ), COLSSQ( 2 ) )
CALL DCOMBSSQ( SSQ, COLSSQ )
110 CONTINUE 110 CONTINUE
L = K + 1 L = K + 1
ELSE ELSE
DO 120 J = 1, N - 1 DO 120 J = 1, N - 1
CALL ZLASSQ( MIN( N-J, K ), AB( 2, J ), 1, SCALE, COLSSQ( 1 ) = ZERO
$ SUM ) COLSSQ( 2 ) = ONE
CALL ZLASSQ( MIN( N-J, K ), AB( 2, J ), 1,
$ COLSSQ( 1 ), COLSSQ( 2 ) )
CALL DCOMBSSQ( SSQ, COLSSQ )
120 CONTINUE 120 CONTINUE
L = 1 L = 1
END IF END IF
SUM = 2*SUM SSQ( 2 ) = 2*SSQ( 2 )
ELSE ELSE
L = 1 L = 1
END IF END IF
*
* Sum diagonal
*
COLSSQ( 1 ) = ZERO
COLSSQ( 2 ) = ONE
DO 130 J = 1, N DO 130 J = 1, N
IF( DBLE( AB( L, J ) ).NE.ZERO ) THEN IF( DBLE( AB( L, J ) ).NE.ZERO ) THEN
ABSA = ABS( DBLE( AB( L, J ) ) ) ABSA = ABS( DBLE( AB( L, J ) ) )
IF( SCALE.LT.ABSA ) THEN IF( COLSSQ( 1 ).LT.ABSA ) THEN
SUM = ONE + SUM*( SCALE / ABSA )**2 COLSSQ( 2 ) = ONE + COLSSQ(2)*( COLSSQ(1) / ABSA )**2
SCALE = ABSA COLSSQ( 1 ) = ABSA
ELSE ELSE
SUM = SUM + ( ABSA / SCALE )**2 COLSSQ( 2 ) = COLSSQ( 2 ) + ( ABSA / COLSSQ( 1 ) )**2
END IF END IF
END IF END IF
130 CONTINUE 130 CONTINUE
VALUE = SCALE*SQRT( SUM ) CALL DCOMBSSQ( SSQ, COLSSQ )
VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) )
END IF END IF
* *
ZLANHB = VALUE ZLANHB = VALUE

View File

@ -129,6 +129,7 @@
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* December 2016 * December 2016
* *
IMPLICIT NONE
* .. Scalar Arguments .. * .. Scalar Arguments ..
CHARACTER NORM, UPLO CHARACTER NORM, UPLO
INTEGER LDA, N INTEGER LDA, N
@ -146,14 +147,17 @@
* .. * ..
* .. Local Scalars .. * .. Local Scalars ..
INTEGER I, J INTEGER I, J
DOUBLE PRECISION ABSA, SCALE, SUM, VALUE DOUBLE PRECISION ABSA, SUM, VALUE
* ..
* .. Local Arrays ..
DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 )
* .. * ..
* .. External Functions .. * .. External Functions ..
LOGICAL LSAME, DISNAN LOGICAL LSAME, DISNAN
EXTERNAL LSAME, DISNAN EXTERNAL LSAME, DISNAN
* .. * ..
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL ZLASSQ EXTERNAL ZLASSQ, DCOMBSSQ
* .. * ..
* .. Intrinsic Functions .. * .. Intrinsic Functions ..
INTRINSIC ABS, DBLE, SQRT INTRINSIC ABS, DBLE, SQRT
@ -223,31 +227,48 @@
ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN
* *
* Find normF(A). * Find normF(A).
* SSQ(1) is scale
* SSQ(2) is sum-of-squares
* For better accuracy, sum each column separately.
*
SSQ( 1 ) = ZERO
SSQ( 2 ) = ONE
*
* Sum off-diagonals
* *
SCALE = ZERO
SUM = ONE
IF( LSAME( UPLO, 'U' ) ) THEN IF( LSAME( UPLO, 'U' ) ) THEN
DO 110 J = 2, N DO 110 J = 2, N
CALL ZLASSQ( J-1, A( 1, J ), 1, SCALE, SUM ) COLSSQ( 1 ) = ZERO
COLSSQ( 2 ) = ONE
CALL ZLASSQ( J-1, A( 1, J ), 1,
$ COLSSQ( 1 ), COLSSQ( 2 ) )
CALL DCOMBSSQ( SSQ, COLSSQ )
110 CONTINUE 110 CONTINUE
ELSE ELSE
DO 120 J = 1, N - 1 DO 120 J = 1, N - 1
CALL ZLASSQ( N-J, A( J+1, J ), 1, SCALE, SUM ) COLSSQ( 1 ) = ZERO
COLSSQ( 2 ) = ONE
CALL ZLASSQ( N-J, A( J+1, J ), 1,
$ COLSSQ( 1 ), COLSSQ( 2 ) )
CALL DCOMBSSQ( SSQ, COLSSQ )
120 CONTINUE 120 CONTINUE
END IF END IF
SUM = 2*SUM SSQ( 2 ) = 2*SSQ( 2 )
*
* Sum diagonal
*
DO 130 I = 1, N DO 130 I = 1, N
IF( DBLE( A( I, I ) ).NE.ZERO ) THEN IF( DBLE( A( I, I ) ).NE.ZERO ) THEN
ABSA = ABS( DBLE( A( I, I ) ) ) ABSA = ABS( DBLE( A( I, I ) ) )
IF( SCALE.LT.ABSA ) THEN IF( SSQ( 1 ).LT.ABSA ) THEN
SUM = ONE + SUM*( SCALE / ABSA )**2 SSQ( 2 ) = ONE + SSQ( 2 )*( SSQ( 1 ) / ABSA )**2
SCALE = ABSA SSQ( 1 ) = ABSA
ELSE ELSE
SUM = SUM + ( ABSA / SCALE )**2 SSQ( 2 ) = SSQ( 2 ) + ( ABSA / SSQ( 1 ) )**2
END IF END IF
END IF END IF
130 CONTINUE 130 CONTINUE
VALUE = SCALE*SQRT( SUM ) VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) )
END IF END IF
* *
ZLANHE = VALUE ZLANHE = VALUE

View File

@ -122,6 +122,7 @@
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* December 2016 * December 2016
* *
IMPLICIT NONE
* .. Scalar Arguments .. * .. Scalar Arguments ..
CHARACTER NORM, UPLO CHARACTER NORM, UPLO
INTEGER N INTEGER N
@ -139,14 +140,17 @@
* .. * ..
* .. Local Scalars .. * .. Local Scalars ..
INTEGER I, J, K INTEGER I, J, K
DOUBLE PRECISION ABSA, SCALE, SUM, VALUE DOUBLE PRECISION ABSA, SUM, VALUE
* ..
* .. Local Arrays ..
DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 )
* .. * ..
* .. External Functions .. * .. External Functions ..
LOGICAL LSAME, DISNAN LOGICAL LSAME, DISNAN
EXTERNAL LSAME, DISNAN EXTERNAL LSAME, DISNAN
* .. * ..
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL ZLASSQ EXTERNAL ZLASSQ, DCOMBSSQ
* .. * ..
* .. Intrinsic Functions .. * .. Intrinsic Functions ..
INTRINSIC ABS, DBLE, SQRT INTRINSIC ABS, DBLE, SQRT
@ -225,31 +229,48 @@
ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN
* *
* Find normF(A). * Find normF(A).
* SSQ(1) is scale
* SSQ(2) is sum-of-squares
* For better accuracy, sum each column separately.
*
SSQ( 1 ) = ZERO
SSQ( 2 ) = ONE
*
* Sum off-diagonals
* *
SCALE = ZERO
SUM = ONE
K = 2 K = 2
IF( LSAME( UPLO, 'U' ) ) THEN IF( LSAME( UPLO, 'U' ) ) THEN
DO 110 J = 2, N DO 110 J = 2, N
CALL ZLASSQ( J-1, AP( K ), 1, SCALE, SUM ) COLSSQ( 1 ) = ZERO
COLSSQ( 2 ) = ONE
CALL ZLASSQ( J-1, AP( K ), 1, COLSSQ( 1 ), COLSSQ( 2 ) )
CALL DCOMBSSQ( SSQ, COLSSQ )
K = K + J K = K + J
110 CONTINUE 110 CONTINUE
ELSE ELSE
DO 120 J = 1, N - 1 DO 120 J = 1, N - 1
CALL ZLASSQ( N-J, AP( K ), 1, SCALE, SUM ) COLSSQ( 1 ) = ZERO
COLSSQ( 2 ) = ONE
CALL ZLASSQ( N-J, AP( K ), 1, COLSSQ( 1 ), COLSSQ( 2 ) )
CALL DCOMBSSQ( SSQ, COLSSQ )
K = K + N - J + 1 K = K + N - J + 1
120 CONTINUE 120 CONTINUE
END IF END IF
SUM = 2*SUM SSQ( 2 ) = 2*SSQ( 2 )
*
* Sum diagonal
*
K = 1 K = 1
COLSSQ( 1 ) = ZERO
COLSSQ( 2 ) = ONE
DO 130 I = 1, N DO 130 I = 1, N
IF( DBLE( AP( K ) ).NE.ZERO ) THEN IF( DBLE( AP( K ) ).NE.ZERO ) THEN
ABSA = ABS( DBLE( AP( K ) ) ) ABSA = ABS( DBLE( AP( K ) ) )
IF( SCALE.LT.ABSA ) THEN IF( COLSSQ( 1 ).LT.ABSA ) THEN
SUM = ONE + SUM*( SCALE / ABSA )**2 COLSSQ( 2 ) = ONE + COLSSQ(2)*( COLSSQ(1) / ABSA )**2
SCALE = ABSA COLSSQ( 1 ) = ABSA
ELSE ELSE
SUM = SUM + ( ABSA / SCALE )**2 COLSSQ( 2 ) = COLSSQ( 2 ) + ( ABSA / COLSSQ( 1 ) )**2
END IF END IF
END IF END IF
IF( LSAME( UPLO, 'U' ) ) THEN IF( LSAME( UPLO, 'U' ) ) THEN
@ -258,7 +279,8 @@
K = K + N - I + 1 K = K + N - I + 1
END IF END IF
130 CONTINUE 130 CONTINUE
VALUE = SCALE*SQRT( SUM ) CALL DCOMBSSQ( SSQ, COLSSQ )
VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) )
END IF END IF
* *
ZLANHP = VALUE ZLANHP = VALUE

View File

@ -114,6 +114,7 @@
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* December 2016 * December 2016
* *
IMPLICIT NONE
* .. Scalar Arguments .. * .. Scalar Arguments ..
CHARACTER NORM CHARACTER NORM
INTEGER LDA, N INTEGER LDA, N
@ -131,14 +132,17 @@
* .. * ..
* .. Local Scalars .. * .. Local Scalars ..
INTEGER I, J INTEGER I, J
DOUBLE PRECISION SCALE, SUM, VALUE DOUBLE PRECISION SUM, VALUE
* ..
* .. Local Arrays ..
DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 )
* .. * ..
* .. External Functions .. * .. External Functions ..
LOGICAL LSAME, DISNAN LOGICAL LSAME, DISNAN
EXTERNAL LSAME, DISNAN EXTERNAL LSAME, DISNAN
* .. * ..
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL ZLASSQ EXTERNAL ZLASSQ, DCOMBSSQ
* .. * ..
* .. Intrinsic Functions .. * .. Intrinsic Functions ..
INTRINSIC ABS, MIN, SQRT INTRINSIC ABS, MIN, SQRT
@ -190,13 +194,20 @@
ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN
* *
* Find normF(A). * Find normF(A).
* SSQ(1) is scale
* SSQ(2) is sum-of-squares
* For better accuracy, sum each column separately.
* *
SCALE = ZERO SSQ( 1 ) = ZERO
SUM = ONE SSQ( 2 ) = ONE
DO 90 J = 1, N DO 90 J = 1, N
CALL ZLASSQ( MIN( N, J+1 ), A( 1, J ), 1, SCALE, SUM ) COLSSQ( 1 ) = ZERO
COLSSQ( 2 ) = ONE
CALL ZLASSQ( MIN( N, J+1 ), A( 1, J ), 1,
$ COLSSQ( 1 ), COLSSQ( 2 ) )
CALL DCOMBSSQ( SSQ, COLSSQ )
90 CONTINUE 90 CONTINUE
VALUE = SCALE*SQRT( SUM ) VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) )
END IF END IF
* *
ZLANHS = VALUE ZLANHS = VALUE

View File

@ -135,6 +135,7 @@
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* December 2016 * December 2016
* *
IMPLICIT NONE
* .. Scalar Arguments .. * .. Scalar Arguments ..
CHARACTER NORM, UPLO CHARACTER NORM, UPLO
INTEGER K, LDAB, N INTEGER K, LDAB, N
@ -152,14 +153,17 @@
* .. * ..
* .. Local Scalars .. * .. Local Scalars ..
INTEGER I, J, L INTEGER I, J, L
DOUBLE PRECISION ABSA, SCALE, SUM, VALUE DOUBLE PRECISION ABSA, SUM, VALUE
* ..
* .. Local Arrays ..
DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 )
* .. * ..
* .. External Functions .. * .. External Functions ..
LOGICAL LSAME, DISNAN LOGICAL LSAME, DISNAN
EXTERNAL LSAME, DISNAN EXTERNAL LSAME, DISNAN
* .. * ..
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL ZLASSQ EXTERNAL ZLASSQ, DCOMBSSQ
* .. * ..
* .. Intrinsic Functions .. * .. Intrinsic Functions ..
INTRINSIC ABS, MAX, MIN, SQRT INTRINSIC ABS, MAX, MIN, SQRT
@ -227,29 +231,47 @@
ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN
* *
* Find normF(A). * Find normF(A).
* SSQ(1) is scale
* SSQ(2) is sum-of-squares
* For better accuracy, sum each column separately.
*
SSQ( 1 ) = ZERO
SSQ( 2 ) = ONE
*
* Sum off-diagonals
* *
SCALE = ZERO
SUM = ONE
IF( K.GT.0 ) THEN IF( K.GT.0 ) THEN
IF( LSAME( UPLO, 'U' ) ) THEN IF( LSAME( UPLO, 'U' ) ) THEN
DO 110 J = 2, N DO 110 J = 2, N
COLSSQ( 1 ) = ZERO
COLSSQ( 2 ) = ONE
CALL ZLASSQ( MIN( J-1, K ), AB( MAX( K+2-J, 1 ), J ), CALL ZLASSQ( MIN( J-1, K ), AB( MAX( K+2-J, 1 ), J ),
$ 1, SCALE, SUM ) $ 1, COLSSQ( 1 ), COLSSQ( 2 ) )
CALL DCOMBSSQ( SSQ, COLSSQ )
110 CONTINUE 110 CONTINUE
L = K + 1 L = K + 1
ELSE ELSE
DO 120 J = 1, N - 1 DO 120 J = 1, N - 1
CALL ZLASSQ( MIN( N-J, K ), AB( 2, J ), 1, SCALE, COLSSQ( 1 ) = ZERO
$ SUM ) COLSSQ( 2 ) = ONE
CALL ZLASSQ( MIN( N-J, K ), AB( 2, J ), 1,
$ COLSSQ( 1 ), COLSSQ( 2 ) )
CALL DCOMBSSQ( SSQ, COLSSQ )
120 CONTINUE 120 CONTINUE
L = 1 L = 1
END IF END IF
SUM = 2*SUM SSQ( 2 ) = 2*SSQ( 2 )
ELSE ELSE
L = 1 L = 1
END IF END IF
CALL ZLASSQ( N, AB( L, 1 ), LDAB, SCALE, SUM ) *
VALUE = SCALE*SQRT( SUM ) * Sum diagonal
*
COLSSQ( 1 ) = ZERO
COLSSQ( 2 ) = ONE
CALL ZLASSQ( N, AB( L, 1 ), LDAB, COLSSQ( 1 ), COLSSQ( 2 ) )
CALL DCOMBSSQ( SSQ, COLSSQ )
VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) )
END IF END IF
* *
ZLANSB = VALUE ZLANSB = VALUE

View File

@ -120,6 +120,7 @@
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* December 2016 * December 2016
* *
IMPLICIT NONE
* .. Scalar Arguments .. * .. Scalar Arguments ..
CHARACTER NORM, UPLO CHARACTER NORM, UPLO
INTEGER N INTEGER N
@ -137,14 +138,17 @@
* .. * ..
* .. Local Scalars .. * .. Local Scalars ..
INTEGER I, J, K INTEGER I, J, K
DOUBLE PRECISION ABSA, SCALE, SUM, VALUE DOUBLE PRECISION ABSA, SUM, VALUE
* ..
* .. Local Arrays ..
DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 )
* .. * ..
* .. External Functions .. * .. External Functions ..
LOGICAL LSAME, DISNAN LOGICAL LSAME, DISNAN
EXTERNAL LSAME, DISNAN EXTERNAL LSAME, DISNAN
* .. * ..
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL ZLASSQ EXTERNAL ZLASSQ, DCOMBSSQ
* .. * ..
* .. Intrinsic Functions .. * .. Intrinsic Functions ..
INTRINSIC ABS, DBLE, DIMAG, SQRT INTRINSIC ABS, DBLE, DIMAG, SQRT
@ -219,40 +223,57 @@
ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN
* *
* Find normF(A). * Find normF(A).
* SSQ(1) is scale
* SSQ(2) is sum-of-squares
* For better accuracy, sum each column separately.
*
SSQ( 1 ) = ZERO
SSQ( 2 ) = ONE
*
* Sum off-diagonals
* *
SCALE = ZERO
SUM = ONE
K = 2 K = 2
IF( LSAME( UPLO, 'U' ) ) THEN IF( LSAME( UPLO, 'U' ) ) THEN
DO 110 J = 2, N DO 110 J = 2, N
CALL ZLASSQ( J-1, AP( K ), 1, SCALE, SUM ) COLSSQ( 1 ) = ZERO
COLSSQ( 2 ) = ONE
CALL ZLASSQ( J-1, AP( K ), 1, COLSSQ( 1 ), COLSSQ( 2 ) )
CALL DCOMBSSQ( SSQ, COLSSQ )
K = K + J K = K + J
110 CONTINUE 110 CONTINUE
ELSE ELSE
DO 120 J = 1, N - 1 DO 120 J = 1, N - 1
CALL ZLASSQ( N-J, AP( K ), 1, SCALE, SUM ) COLSSQ( 1 ) = ZERO
COLSSQ( 2 ) = ONE
CALL ZLASSQ( N-J, AP( K ), 1, COLSSQ( 1 ), COLSSQ( 2 ) )
CALL DCOMBSSQ( SSQ, COLSSQ )
K = K + N - J + 1 K = K + N - J + 1
120 CONTINUE 120 CONTINUE
END IF END IF
SUM = 2*SUM SSQ( 2 ) = 2*SSQ( 2 )
*
* Sum diagonal
*
K = 1 K = 1
COLSSQ( 1 ) = ZERO
COLSSQ( 2 ) = ONE
DO 130 I = 1, N DO 130 I = 1, N
IF( DBLE( AP( K ) ).NE.ZERO ) THEN IF( DBLE( AP( K ) ).NE.ZERO ) THEN
ABSA = ABS( DBLE( AP( K ) ) ) ABSA = ABS( DBLE( AP( K ) ) )
IF( SCALE.LT.ABSA ) THEN IF( COLSSQ( 1 ).LT.ABSA ) THEN
SUM = ONE + SUM*( SCALE / ABSA )**2 COLSSQ( 2 ) = ONE + COLSSQ(2)*( COLSSQ(1) / ABSA )**2
SCALE = ABSA COLSSQ( 1 ) = ABSA
ELSE ELSE
SUM = SUM + ( ABSA / SCALE )**2 COLSSQ( 2 ) = COLSSQ( 2 ) + ( ABSA / COLSSQ( 1 ) )**2
END IF END IF
END IF END IF
IF( DIMAG( AP( K ) ).NE.ZERO ) THEN IF( DIMAG( AP( K ) ).NE.ZERO ) THEN
ABSA = ABS( DIMAG( AP( K ) ) ) ABSA = ABS( DIMAG( AP( K ) ) )
IF( SCALE.LT.ABSA ) THEN IF( COLSSQ( 1 ).LT.ABSA ) THEN
SUM = ONE + SUM*( SCALE / ABSA )**2 COLSSQ( 2 ) = ONE + COLSSQ(2)*( COLSSQ(1) / ABSA )**2
SCALE = ABSA COLSSQ( 1 ) = ABSA
ELSE ELSE
SUM = SUM + ( ABSA / SCALE )**2 COLSSQ( 2 ) = COLSSQ( 2 ) + ( ABSA / COLSSQ( 1 ) )**2
END IF END IF
END IF END IF
IF( LSAME( UPLO, 'U' ) ) THEN IF( LSAME( UPLO, 'U' ) ) THEN
@ -261,7 +282,8 @@
K = K + N - I + 1 K = K + N - I + 1
END IF END IF
130 CONTINUE 130 CONTINUE
VALUE = SCALE*SQRT( SUM ) CALL DCOMBSSQ( SSQ, COLSSQ )
VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) )
END IF END IF
* *
ZLANSP = VALUE ZLANSP = VALUE

View File

@ -128,6 +128,7 @@
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* December 2016 * December 2016
* *
IMPLICIT NONE
* .. Scalar Arguments .. * .. Scalar Arguments ..
CHARACTER NORM, UPLO CHARACTER NORM, UPLO
INTEGER LDA, N INTEGER LDA, N
@ -145,14 +146,17 @@
* .. * ..
* .. Local Scalars .. * .. Local Scalars ..
INTEGER I, J INTEGER I, J
DOUBLE PRECISION ABSA, SCALE, SUM, VALUE DOUBLE PRECISION ABSA, SUM, VALUE
* ..
* .. Local Arrays ..
DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 )
* .. * ..
* .. External Functions .. * .. External Functions ..
LOGICAL LSAME, DISNAN LOGICAL LSAME, DISNAN
EXTERNAL LSAME, DISNAN EXTERNAL LSAME, DISNAN
* .. * ..
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL ZLASSQ EXTERNAL ZLASSQ, DCOMBSSQ
* .. * ..
* .. Intrinsic Functions .. * .. Intrinsic Functions ..
INTRINSIC ABS, SQRT INTRINSIC ABS, SQRT
@ -218,21 +222,39 @@
ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN
* *
* Find normF(A). * Find normF(A).
* SSQ(1) is scale
* SSQ(2) is sum-of-squares
* For better accuracy, sum each column separately.
*
SSQ( 1 ) = ZERO
SSQ( 2 ) = ONE
*
* Sum off-diagonals
* *
SCALE = ZERO
SUM = ONE
IF( LSAME( UPLO, 'U' ) ) THEN IF( LSAME( UPLO, 'U' ) ) THEN
DO 110 J = 2, N DO 110 J = 2, N
CALL ZLASSQ( J-1, A( 1, J ), 1, SCALE, SUM ) COLSSQ( 1 ) = ZERO
COLSSQ( 2 ) = ONE
CALL ZLASSQ( J-1, A( 1, J ), 1, COLSSQ(1), COLSSQ(2) )
CALL DCOMBSSQ( SSQ, COLSSQ )
110 CONTINUE 110 CONTINUE
ELSE ELSE
DO 120 J = 1, N - 1 DO 120 J = 1, N - 1
CALL ZLASSQ( N-J, A( J+1, J ), 1, SCALE, SUM ) COLSSQ( 1 ) = ZERO
COLSSQ( 2 ) = ONE
CALL ZLASSQ( N-J, A( J+1, J ), 1, COLSSQ(1), COLSSQ(2) )
CALL DCOMBSSQ( SSQ, COLSSQ )
120 CONTINUE 120 CONTINUE
END IF END IF
SUM = 2*SUM SSQ( 2 ) = 2*SSQ( 2 )
CALL ZLASSQ( N, A, LDA+1, SCALE, SUM ) *
VALUE = SCALE*SQRT( SUM ) * Sum diagonal
*
COLSSQ( 1 ) = ZERO
COLSSQ( 2 ) = ONE
CALL ZLASSQ( N, A, LDA+1, COLSSQ( 1 ), COLSSQ( 2 ) )
CALL DCOMBSSQ( SSQ, COLSSQ )
VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) )
END IF END IF
* *
ZLANSY = VALUE ZLANSY = VALUE

View File

@ -146,6 +146,7 @@
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* December 2016 * December 2016
* *
IMPLICIT NONE
* .. Scalar Arguments .. * .. Scalar Arguments ..
CHARACTER DIAG, NORM, UPLO CHARACTER DIAG, NORM, UPLO
INTEGER K, LDAB, N INTEGER K, LDAB, N
@ -164,14 +165,17 @@
* .. Local Scalars .. * .. Local Scalars ..
LOGICAL UDIAG LOGICAL UDIAG
INTEGER I, J, L INTEGER I, J, L
DOUBLE PRECISION SCALE, SUM, VALUE DOUBLE PRECISION SUM, VALUE
* ..
* .. Local Arrays ..
DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 )
* .. * ..
* .. External Functions .. * .. External Functions ..
LOGICAL LSAME, DISNAN LOGICAL LSAME, DISNAN
EXTERNAL LSAME, DISNAN EXTERNAL LSAME, DISNAN
* .. * ..
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL ZLASSQ EXTERNAL ZLASSQ, DCOMBSSQ
* .. * ..
* .. Intrinsic Functions .. * .. Intrinsic Functions ..
INTRINSIC ABS, MAX, MIN, SQRT INTRINSIC ABS, MAX, MIN, SQRT
@ -313,46 +317,61 @@
ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN
* *
* Find normF(A). * Find normF(A).
* SSQ(1) is scale
* SSQ(2) is sum-of-squares
* For better accuracy, sum each column separately.
* *
IF( LSAME( UPLO, 'U' ) ) THEN IF( LSAME( UPLO, 'U' ) ) THEN
IF( LSAME( DIAG, 'U' ) ) THEN IF( LSAME( DIAG, 'U' ) ) THEN
SCALE = ONE SSQ( 1 ) = ONE
SUM = N SSQ( 2 ) = N
IF( K.GT.0 ) THEN IF( K.GT.0 ) THEN
DO 280 J = 2, N DO 280 J = 2, N
COLSSQ( 1 ) = ZERO
COLSSQ( 2 ) = ONE
CALL ZLASSQ( MIN( J-1, K ), CALL ZLASSQ( MIN( J-1, K ),
$ AB( MAX( K+2-J, 1 ), J ), 1, SCALE, $ AB( MAX( K+2-J, 1 ), J ), 1,
$ SUM ) $ COLSSQ( 1 ), COLSSQ( 2 ) )
CALL DCOMBSSQ( SSQ, COLSSQ )
280 CONTINUE 280 CONTINUE
END IF END IF
ELSE ELSE
SCALE = ZERO SSQ( 1 ) = ZERO
SUM = ONE SSQ( 2 ) = ONE
DO 290 J = 1, N DO 290 J = 1, N
COLSSQ( 1 ) = ZERO
COLSSQ( 2 ) = ONE
CALL ZLASSQ( MIN( J, K+1 ), AB( MAX( K+2-J, 1 ), J ), CALL ZLASSQ( MIN( J, K+1 ), AB( MAX( K+2-J, 1 ), J ),
$ 1, SCALE, SUM ) $ 1, COLSSQ( 1 ), COLSSQ( 2 ) )
CALL DCOMBSSQ( SSQ, COLSSQ )
290 CONTINUE 290 CONTINUE
END IF END IF
ELSE ELSE
IF( LSAME( DIAG, 'U' ) ) THEN IF( LSAME( DIAG, 'U' ) ) THEN
SCALE = ONE SSQ( 1 ) = ONE
SUM = N SSQ( 2 ) = N
IF( K.GT.0 ) THEN IF( K.GT.0 ) THEN
DO 300 J = 1, N - 1 DO 300 J = 1, N - 1
CALL ZLASSQ( MIN( N-J, K ), AB( 2, J ), 1, SCALE, COLSSQ( 1 ) = ZERO
$ SUM ) COLSSQ( 2 ) = ONE
CALL ZLASSQ( MIN( N-J, K ), AB( 2, J ), 1,
$ COLSSQ( 1 ), COLSSQ( 2 ) )
CALL DCOMBSSQ( SSQ, COLSSQ )
300 CONTINUE 300 CONTINUE
END IF END IF
ELSE ELSE
SCALE = ZERO SSQ( 1 ) = ZERO
SUM = ONE SSQ( 2 ) = ONE
DO 310 J = 1, N DO 310 J = 1, N
CALL ZLASSQ( MIN( N-J+1, K+1 ), AB( 1, J ), 1, SCALE, COLSSQ( 1 ) = ZERO
$ SUM ) COLSSQ( 2 ) = ONE
CALL ZLASSQ( MIN( N-J+1, K+1 ), AB( 1, J ), 1,
$ COLSSQ( 1 ), COLSSQ( 2 ) )
CALL DCOMBSSQ( SSQ, COLSSQ )
310 CONTINUE 310 CONTINUE
END IF END IF
END IF END IF
VALUE = SCALE*SQRT( SUM ) VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) )
END IF END IF
* *
ZLANTB = VALUE ZLANTB = VALUE

View File

@ -130,6 +130,7 @@
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* December 2016 * December 2016
* *
IMPLICIT NONE
* .. Scalar Arguments .. * .. Scalar Arguments ..
CHARACTER DIAG, NORM, UPLO CHARACTER DIAG, NORM, UPLO
INTEGER N INTEGER N
@ -148,14 +149,17 @@
* .. Local Scalars .. * .. Local Scalars ..
LOGICAL UDIAG LOGICAL UDIAG
INTEGER I, J, K INTEGER I, J, K
DOUBLE PRECISION SCALE, SUM, VALUE DOUBLE PRECISION SUM, VALUE
* ..
* .. Local Arrays ..
DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 )
* .. * ..
* .. External Functions .. * .. External Functions ..
LOGICAL LSAME, DISNAN LOGICAL LSAME, DISNAN
EXTERNAL LSAME, DISNAN EXTERNAL LSAME, DISNAN
* .. * ..
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL ZLASSQ EXTERNAL ZLASSQ, DCOMBSSQ
* .. * ..
* .. Intrinsic Functions .. * .. Intrinsic Functions ..
INTRINSIC ABS, SQRT INTRINSIC ABS, SQRT
@ -308,45 +312,64 @@
ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN
* *
* Find normF(A). * Find normF(A).
* SSQ(1) is scale
* SSQ(2) is sum-of-squares
* For better accuracy, sum each column separately.
* *
IF( LSAME( UPLO, 'U' ) ) THEN IF( LSAME( UPLO, 'U' ) ) THEN
IF( LSAME( DIAG, 'U' ) ) THEN IF( LSAME( DIAG, 'U' ) ) THEN
SCALE = ONE SSQ( 1 ) = ONE
SUM = N SSQ( 2 ) = N
K = 2 K = 2
DO 280 J = 2, N DO 280 J = 2, N
CALL ZLASSQ( J-1, AP( K ), 1, SCALE, SUM ) COLSSQ( 1 ) = ZERO
COLSSQ( 2 ) = ONE
CALL ZLASSQ( J-1, AP( K ), 1,
$ COLSSQ( 1 ), COLSSQ( 2 ) )
CALL DCOMBSSQ( SSQ, COLSSQ )
K = K + J K = K + J
280 CONTINUE 280 CONTINUE
ELSE ELSE
SCALE = ZERO SSQ( 1 ) = ZERO
SUM = ONE SSQ( 2 ) = ONE
K = 1 K = 1
DO 290 J = 1, N DO 290 J = 1, N
CALL ZLASSQ( J, AP( K ), 1, SCALE, SUM ) COLSSQ( 1 ) = ZERO
COLSSQ( 2 ) = ONE
CALL ZLASSQ( J, AP( K ), 1,
$ COLSSQ( 1 ), COLSSQ( 2 ) )
CALL DCOMBSSQ( SSQ, COLSSQ )
K = K + J K = K + J
290 CONTINUE 290 CONTINUE
END IF END IF
ELSE ELSE
IF( LSAME( DIAG, 'U' ) ) THEN IF( LSAME( DIAG, 'U' ) ) THEN
SCALE = ONE SSQ( 1 ) = ONE
SUM = N SSQ( 2 ) = N
K = 2 K = 2
DO 300 J = 1, N - 1 DO 300 J = 1, N - 1
CALL ZLASSQ( N-J, AP( K ), 1, SCALE, SUM ) COLSSQ( 1 ) = ZERO
COLSSQ( 2 ) = ONE
CALL ZLASSQ( N-J, AP( K ), 1,
$ COLSSQ( 1 ), COLSSQ( 2 ) )
CALL DCOMBSSQ( SSQ, COLSSQ )
K = K + N - J + 1 K = K + N - J + 1
300 CONTINUE 300 CONTINUE
ELSE ELSE
SCALE = ZERO SSQ( 1 ) = ZERO
SUM = ONE SSQ( 2 ) = ONE
K = 1 K = 1
DO 310 J = 1, N DO 310 J = 1, N
CALL ZLASSQ( N-J+1, AP( K ), 1, SCALE, SUM ) COLSSQ( 1 ) = ZERO
COLSSQ( 2 ) = ONE
CALL ZLASSQ( N-J+1, AP( K ), 1,
$ COLSSQ( 1 ), COLSSQ( 2 ) )
CALL DCOMBSSQ( SSQ, COLSSQ )
K = K + N - J + 1 K = K + N - J + 1
310 CONTINUE 310 CONTINUE
END IF END IF
END IF END IF
VALUE = SCALE*SQRT( SUM ) VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) )
END IF END IF
* *
ZLANTP = VALUE ZLANTP = VALUE

View File

@ -147,6 +147,7 @@
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* December 2016 * December 2016
* *
IMPLICIT NONE
* .. Scalar Arguments .. * .. Scalar Arguments ..
CHARACTER DIAG, NORM, UPLO CHARACTER DIAG, NORM, UPLO
INTEGER LDA, M, N INTEGER LDA, M, N
@ -165,14 +166,17 @@
* .. Local Scalars .. * .. Local Scalars ..
LOGICAL UDIAG LOGICAL UDIAG
INTEGER I, J INTEGER I, J
DOUBLE PRECISION SCALE, SUM, VALUE DOUBLE PRECISION SUM, VALUE
* ..
* .. Local Arrays ..
DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 )
* .. * ..
* .. External Functions .. * .. External Functions ..
LOGICAL LSAME, DISNAN LOGICAL LSAME, DISNAN
EXTERNAL LSAME, DISNAN EXTERNAL LSAME, DISNAN
* .. * ..
* .. External Subroutines .. * .. External Subroutines ..
EXTERNAL ZLASSQ EXTERNAL ZLASSQ, DCOMBSSQ
* .. * ..
* .. Intrinsic Functions .. * .. Intrinsic Functions ..
INTRINSIC ABS, MIN, SQRT INTRINSIC ABS, MIN, SQRT
@ -283,7 +287,7 @@
END IF END IF
ELSE ELSE
IF( LSAME( DIAG, 'U' ) ) THEN IF( LSAME( DIAG, 'U' ) ) THEN
DO 210 I = 1, N DO 210 I = 1, MIN( M, N )
WORK( I ) = ONE WORK( I ) = ONE
210 CONTINUE 210 CONTINUE
DO 220 I = N + 1, M DO 220 I = N + 1, M
@ -313,38 +317,56 @@
ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN
* *
* Find normF(A). * Find normF(A).
* SSQ(1) is scale
* SSQ(2) is sum-of-squares
* For better accuracy, sum each column separately.
* *
IF( LSAME( UPLO, 'U' ) ) THEN IF( LSAME( UPLO, 'U' ) ) THEN
IF( LSAME( DIAG, 'U' ) ) THEN IF( LSAME( DIAG, 'U' ) ) THEN
SCALE = ONE SSQ( 1 ) = ONE
SUM = MIN( M, N ) SSQ( 2 ) = MIN( M, N )
DO 290 J = 2, N DO 290 J = 2, N
CALL ZLASSQ( MIN( M, J-1 ), A( 1, J ), 1, SCALE, SUM ) COLSSQ( 1 ) = ZERO
COLSSQ( 2 ) = ONE
CALL ZLASSQ( MIN( M, J-1 ), A( 1, J ), 1,
$ COLSSQ( 1 ), COLSSQ( 2 ) )
CALL DCOMBSSQ( SSQ, COLSSQ )
290 CONTINUE 290 CONTINUE
ELSE ELSE
SCALE = ZERO SSQ( 1 ) = ZERO
SUM = ONE SSQ( 2 ) = ONE
DO 300 J = 1, N DO 300 J = 1, N
CALL ZLASSQ( MIN( M, J ), A( 1, J ), 1, SCALE, SUM ) COLSSQ( 1 ) = ZERO
COLSSQ( 2 ) = ONE
CALL ZLASSQ( MIN( M, J ), A( 1, J ), 1,
$ COLSSQ( 1 ), COLSSQ( 2 ) )
CALL DCOMBSSQ( SSQ, COLSSQ )
300 CONTINUE 300 CONTINUE
END IF END IF
ELSE ELSE
IF( LSAME( DIAG, 'U' ) ) THEN IF( LSAME( DIAG, 'U' ) ) THEN
SCALE = ONE SSQ( 1 ) = ONE
SUM = MIN( M, N ) SSQ( 2 ) = MIN( M, N )
DO 310 J = 1, N DO 310 J = 1, N
CALL ZLASSQ( M-J, A( MIN( M, J+1 ), J ), 1, SCALE, COLSSQ( 1 ) = ZERO
$ SUM ) COLSSQ( 2 ) = ONE
CALL ZLASSQ( M-J, A( MIN( M, J+1 ), J ), 1,
$ COLSSQ( 1 ), COLSSQ( 2 ) )
CALL DCOMBSSQ( SSQ, COLSSQ )
310 CONTINUE 310 CONTINUE
ELSE ELSE
SCALE = ZERO SSQ( 1 ) = ZERO
SUM = ONE SSQ( 2 ) = ONE
DO 320 J = 1, N DO 320 J = 1, N
CALL ZLASSQ( M-J+1, A( J, J ), 1, SCALE, SUM ) COLSSQ( 1 ) = ZERO
COLSSQ( 2 ) = ONE
CALL ZLASSQ( M-J+1, A( J, J ), 1,
$ COLSSQ( 1 ), COLSSQ( 2 ) )
CALL DCOMBSSQ( SSQ, COLSSQ )
320 CONTINUE 320 CONTINUE
END IF END IF
END IF END IF
VALUE = SCALE*SQRT( SUM ) VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) )
END IF END IF
* *
ZLANTR = VALUE ZLANTR = VALUE

View File

@ -127,7 +127,7 @@
*> \param[in,out] AUXV *> \param[in,out] AUXV
*> \verbatim *> \verbatim
*> AUXV is COMPLEX*16 array, dimension (NB) *> AUXV is COMPLEX*16 array, dimension (NB)
*> Auxiliar vector. *> Auxiliary vector.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in,out] F *> \param[in,out] F

View File

@ -66,7 +66,7 @@
*> \param[in] N *> \param[in] N
*> \verbatim *> \verbatim
*> N is INTEGER *> N is INTEGER
*> The order of the matrix H. N .GE. 0. *> The order of the matrix H. N >= 0.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] ILO *> \param[in] ILO
@ -79,12 +79,12 @@
*> IHI is INTEGER *> IHI is INTEGER
*> *>
*> It is assumed that H is already upper triangular in rows *> It is assumed that H is already upper triangular in rows
*> and columns 1:ILO-1 and IHI+1:N and, if ILO.GT.1, *> and columns 1:ILO-1 and IHI+1:N and, if ILO > 1,
*> H(ILO,ILO-1) is zero. ILO and IHI are normally set by a *> H(ILO,ILO-1) is zero. ILO and IHI are normally set by a
*> previous call to ZGEBAL, and then passed to ZGEHRD when the *> previous call to ZGEBAL, and then passed to ZGEHRD when the
*> matrix output by ZGEBAL is reduced to Hessenberg form. *> matrix output by ZGEBAL is reduced to Hessenberg form.
*> Otherwise, ILO and IHI should be set to 1 and N, *> Otherwise, ILO and IHI should be set to 1 and N,
*> respectively. If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N. *> respectively. If N > 0, then 1 <= ILO <= IHI <= N.
*> If N = 0, then ILO = 1 and IHI = 0. *> If N = 0, then ILO = 1 and IHI = 0.
*> \endverbatim *> \endverbatim
*> *>
@ -96,17 +96,17 @@
*> contains the upper triangular matrix T from the Schur *> contains the upper triangular matrix T from the Schur
*> decomposition (the Schur form). If INFO = 0 and WANT is *> decomposition (the Schur form). If INFO = 0 and WANT is
*> .FALSE., then the contents of H are unspecified on exit. *> .FALSE., then the contents of H are unspecified on exit.
*> (The output value of H when INFO.GT.0 is given under the *> (The output value of H when INFO > 0 is given under the
*> description of INFO below.) *> description of INFO below.)
*> *>
*> This subroutine may explicitly set H(i,j) = 0 for i.GT.j and *> This subroutine may explicitly set H(i,j) = 0 for i > j and
*> j = 1, 2, ... ILO-1 or j = IHI+1, IHI+2, ... N. *> j = 1, 2, ... ILO-1 or j = IHI+1, IHI+2, ... N.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] LDH *> \param[in] LDH
*> \verbatim *> \verbatim
*> LDH is INTEGER *> LDH is INTEGER
*> The leading dimension of the array H. LDH .GE. max(1,N). *> The leading dimension of the array H. LDH >= max(1,N).
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] W *> \param[out] W
@ -128,7 +128,7 @@
*> IHIZ is INTEGER *> IHIZ is INTEGER
*> Specify the rows of Z to which transformations must be *> Specify the rows of Z to which transformations must be
*> applied if WANTZ is .TRUE.. *> applied if WANTZ is .TRUE..
*> 1 .LE. ILOZ .LE. ILO; IHI .LE. IHIZ .LE. N. *> 1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in,out] Z *> \param[in,out] Z
@ -138,7 +138,7 @@
*> If WANTZ is .TRUE., then Z(ILO:IHI,ILOZ:IHIZ) is *> If WANTZ is .TRUE., then Z(ILO:IHI,ILOZ:IHIZ) is
*> replaced by Z(ILO:IHI,ILOZ:IHIZ)*U where U is the *> replaced by Z(ILO:IHI,ILOZ:IHIZ)*U where U is the
*> orthogonal Schur factor of H(ILO:IHI,ILO:IHI). *> orthogonal Schur factor of H(ILO:IHI,ILO:IHI).
*> (The output value of Z when INFO.GT.0 is given under *> (The output value of Z when INFO > 0 is given under
*> the description of INFO below.) *> the description of INFO below.)
*> \endverbatim *> \endverbatim
*> *>
@ -146,7 +146,7 @@
*> \verbatim *> \verbatim
*> LDZ is INTEGER *> LDZ is INTEGER
*> The leading dimension of the array Z. if WANTZ is .TRUE. *> The leading dimension of the array Z. if WANTZ is .TRUE.
*> then LDZ.GE.MAX(1,IHIZ). Otherwize, LDZ.GE.1. *> then LDZ >= MAX(1,IHIZ). Otherwise, LDZ >= 1.
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] WORK *> \param[out] WORK
@ -159,7 +159,7 @@
*> \param[in] LWORK *> \param[in] LWORK
*> \verbatim *> \verbatim
*> LWORK is INTEGER *> LWORK is INTEGER
*> The dimension of the array WORK. LWORK .GE. max(1,N) *> The dimension of the array WORK. LWORK >= max(1,N)
*> is sufficient, but LWORK typically as large as 6*N may *> is sufficient, but LWORK typically as large as 6*N may
*> be required for optimal performance. A workspace query *> be required for optimal performance. A workspace query
*> to determine the optimal workspace size is recommended. *> to determine the optimal workspace size is recommended.
@ -176,18 +176,18 @@
*> \verbatim *> \verbatim
*> INFO is INTEGER *> INFO is INTEGER
*> = 0: successful exit *> = 0: successful exit
*> .GT. 0: if INFO = i, ZLAQR0 failed to compute all of *> > 0: if INFO = i, ZLAQR0 failed to compute all of
*> the eigenvalues. Elements 1:ilo-1 and i+1:n of WR *> the eigenvalues. Elements 1:ilo-1 and i+1:n of WR
*> and WI contain those eigenvalues which have been *> and WI contain those eigenvalues which have been
*> successfully computed. (Failures are rare.) *> successfully computed. (Failures are rare.)
*> *>
*> If INFO .GT. 0 and WANT is .FALSE., then on exit, *> If INFO > 0 and WANT is .FALSE., then on exit,
*> the remaining unconverged eigenvalues are the eigen- *> the remaining unconverged eigenvalues are the eigen-
*> values of the upper Hessenberg matrix rows and *> values of the upper Hessenberg matrix rows and
*> columns ILO through INFO of the final, output *> columns ILO through INFO of the final, output
*> value of H. *> value of H.
*> *>
*> If INFO .GT. 0 and WANTT is .TRUE., then on exit *> If INFO > 0 and WANTT is .TRUE., then on exit
*> *>
*> (*) (initial value of H)*U = U*(final value of H) *> (*) (initial value of H)*U = U*(final value of H)
*> *>
@ -195,7 +195,7 @@
*> value of H is upper Hessenberg and triangular in *> value of H is upper Hessenberg and triangular in
*> rows and columns INFO+1 through IHI. *> rows and columns INFO+1 through IHI.
*> *>
*> If INFO .GT. 0 and WANTZ is .TRUE., then on exit *> If INFO > 0 and WANTZ is .TRUE., then on exit
*> *>
*> (final value of Z(ILO:IHI,ILOZ:IHIZ) *> (final value of Z(ILO:IHI,ILOZ:IHIZ)
*> = (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U *> = (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U
@ -203,7 +203,7 @@
*> where U is the unitary matrix in (*) (regard- *> where U is the unitary matrix in (*) (regard-
*> less of the value of WANTT.) *> less of the value of WANTT.)
*> *>
*> If INFO .GT. 0 and WANTZ is .FALSE., then Z is not *> If INFO > 0 and WANTZ is .FALSE., then Z is not
*> accessed. *> accessed.
*> \endverbatim *> \endverbatim
* *
@ -641,7 +641,7 @@
END IF END IF
END IF END IF
* *
* ==== Use up to NS of the the smallest magnatiude * ==== Use up to NS of the the smallest magnitude
* . shifts. If there aren't NS shifts available, * . shifts. If there aren't NS shifts available,
* . then use them all, possibly dropping one to * . then use them all, possibly dropping one to
* . make the number of shifts even. ==== * . make the number of shifts even. ====

View File

@ -64,7 +64,7 @@
*> \verbatim *> \verbatim
*> LDH is INTEGER *> LDH is INTEGER
*> The leading dimension of H as declared in *> The leading dimension of H as declared in
*> the calling procedure. LDH.GE.N *> the calling procedure. LDH >= N
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] S1 *> \param[in] S1

View File

@ -103,7 +103,7 @@
*> \param[in] NW *> \param[in] NW
*> \verbatim *> \verbatim
*> NW is INTEGER *> NW is INTEGER
*> Deflation window size. 1 .LE. NW .LE. (KBOT-KTOP+1). *> Deflation window size. 1 <= NW <= (KBOT-KTOP+1).
*> \endverbatim *> \endverbatim
*> *>
*> \param[in,out] H *> \param[in,out] H
@ -121,7 +121,7 @@
*> \verbatim *> \verbatim
*> LDH is INTEGER *> LDH is INTEGER
*> Leading dimension of H just as declared in the calling *> Leading dimension of H just as declared in the calling
*> subroutine. N .LE. LDH *> subroutine. N <= LDH
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] ILOZ *> \param[in] ILOZ
@ -133,7 +133,7 @@
*> \verbatim *> \verbatim
*> IHIZ is INTEGER *> IHIZ is INTEGER
*> Specify the rows of Z to which transformations must be *> Specify the rows of Z to which transformations must be
*> applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N. *> applied if WANTZ is .TRUE.. 1 <= ILOZ <= IHIZ <= N.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in,out] Z *> \param[in,out] Z
@ -149,7 +149,7 @@
*> \verbatim *> \verbatim
*> LDZ is INTEGER *> LDZ is INTEGER
*> The leading dimension of Z just as declared in the *> The leading dimension of Z just as declared in the
*> calling subroutine. 1 .LE. LDZ. *> calling subroutine. 1 <= LDZ.
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] NS *> \param[out] NS
@ -186,13 +186,13 @@
*> \verbatim *> \verbatim
*> LDV is INTEGER *> LDV is INTEGER
*> The leading dimension of V just as declared in the *> The leading dimension of V just as declared in the
*> calling subroutine. NW .LE. LDV *> calling subroutine. NW <= LDV
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] NH *> \param[in] NH
*> \verbatim *> \verbatim
*> NH is INTEGER *> NH is INTEGER
*> The number of columns of T. NH.GE.NW. *> The number of columns of T. NH >= NW.
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] T *> \param[out] T
@ -204,14 +204,14 @@
*> \verbatim *> \verbatim
*> LDT is INTEGER *> LDT is INTEGER
*> The leading dimension of T just as declared in the *> The leading dimension of T just as declared in the
*> calling subroutine. NW .LE. LDT *> calling subroutine. NW <= LDT
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] NV *> \param[in] NV
*> \verbatim *> \verbatim
*> NV is INTEGER *> NV is INTEGER
*> The number of rows of work array WV available for *> The number of rows of work array WV available for
*> workspace. NV.GE.NW. *> workspace. NV >= NW.
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] WV *> \param[out] WV
@ -223,7 +223,7 @@
*> \verbatim *> \verbatim
*> LDWV is INTEGER *> LDWV is INTEGER
*> The leading dimension of W just as declared in the *> The leading dimension of W just as declared in the
*> calling subroutine. NW .LE. LDV *> calling subroutine. NW <= LDV
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] WORK *> \param[out] WORK

View File

@ -100,7 +100,7 @@
*> \param[in] NW *> \param[in] NW
*> \verbatim *> \verbatim
*> NW is INTEGER *> NW is INTEGER
*> Deflation window size. 1 .LE. NW .LE. (KBOT-KTOP+1). *> Deflation window size. 1 <= NW <= (KBOT-KTOP+1).
*> \endverbatim *> \endverbatim
*> *>
*> \param[in,out] H *> \param[in,out] H
@ -118,7 +118,7 @@
*> \verbatim *> \verbatim
*> LDH is INTEGER *> LDH is INTEGER
*> Leading dimension of H just as declared in the calling *> Leading dimension of H just as declared in the calling
*> subroutine. N .LE. LDH *> subroutine. N <= LDH
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] ILOZ *> \param[in] ILOZ
@ -130,7 +130,7 @@
*> \verbatim *> \verbatim
*> IHIZ is INTEGER *> IHIZ is INTEGER
*> Specify the rows of Z to which transformations must be *> Specify the rows of Z to which transformations must be
*> applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N. *> applied if WANTZ is .TRUE.. 1 <= ILOZ <= IHIZ <= N.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in,out] Z *> \param[in,out] Z
@ -146,7 +146,7 @@
*> \verbatim *> \verbatim
*> LDZ is INTEGER *> LDZ is INTEGER
*> The leading dimension of Z just as declared in the *> The leading dimension of Z just as declared in the
*> calling subroutine. 1 .LE. LDZ. *> calling subroutine. 1 <= LDZ.
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] NS *> \param[out] NS
@ -183,13 +183,13 @@
*> \verbatim *> \verbatim
*> LDV is INTEGER *> LDV is INTEGER
*> The leading dimension of V just as declared in the *> The leading dimension of V just as declared in the
*> calling subroutine. NW .LE. LDV *> calling subroutine. NW <= LDV
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] NH *> \param[in] NH
*> \verbatim *> \verbatim
*> NH is INTEGER *> NH is INTEGER
*> The number of columns of T. NH.GE.NW. *> The number of columns of T. NH >= NW.
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] T *> \param[out] T
@ -201,14 +201,14 @@
*> \verbatim *> \verbatim
*> LDT is INTEGER *> LDT is INTEGER
*> The leading dimension of T just as declared in the *> The leading dimension of T just as declared in the
*> calling subroutine. NW .LE. LDT *> calling subroutine. NW <= LDT
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] NV *> \param[in] NV
*> \verbatim *> \verbatim
*> NV is INTEGER *> NV is INTEGER
*> The number of rows of work array WV available for *> The number of rows of work array WV available for
*> workspace. NV.GE.NW. *> workspace. NV >= NW.
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] WV *> \param[out] WV
@ -220,7 +220,7 @@
*> \verbatim *> \verbatim
*> LDWV is INTEGER *> LDWV is INTEGER
*> The leading dimension of W just as declared in the *> The leading dimension of W just as declared in the
*> calling subroutine. NW .LE. LDV *> calling subroutine. NW <= LDV
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] WORK *> \param[out] WORK

View File

@ -73,7 +73,7 @@
*> \param[in] N *> \param[in] N
*> \verbatim *> \verbatim
*> N is INTEGER *> N is INTEGER
*> The order of the matrix H. N .GE. 0. *> The order of the matrix H. N >= 0.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] ILO *> \param[in] ILO
@ -85,12 +85,12 @@
*> \verbatim *> \verbatim
*> IHI is INTEGER *> IHI is INTEGER
*> It is assumed that H is already upper triangular in rows *> It is assumed that H is already upper triangular in rows
*> and columns 1:ILO-1 and IHI+1:N and, if ILO.GT.1, *> and columns 1:ILO-1 and IHI+1:N and, if ILO > 1,
*> H(ILO,ILO-1) is zero. ILO and IHI are normally set by a *> H(ILO,ILO-1) is zero. ILO and IHI are normally set by a
*> previous call to ZGEBAL, and then passed to ZGEHRD when the *> previous call to ZGEBAL, and then passed to ZGEHRD when the
*> matrix output by ZGEBAL is reduced to Hessenberg form. *> matrix output by ZGEBAL is reduced to Hessenberg form.
*> Otherwise, ILO and IHI should be set to 1 and N, *> Otherwise, ILO and IHI should be set to 1 and N,
*> respectively. If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N. *> respectively. If N > 0, then 1 <= ILO <= IHI <= N.
*> If N = 0, then ILO = 1 and IHI = 0. *> If N = 0, then ILO = 1 and IHI = 0.
*> \endverbatim *> \endverbatim
*> *>
@ -102,17 +102,17 @@
*> contains the upper triangular matrix T from the Schur *> contains the upper triangular matrix T from the Schur
*> decomposition (the Schur form). If INFO = 0 and WANT is *> decomposition (the Schur form). If INFO = 0 and WANT is
*> .FALSE., then the contents of H are unspecified on exit. *> .FALSE., then the contents of H are unspecified on exit.
*> (The output value of H when INFO.GT.0 is given under the *> (The output value of H when INFO > 0 is given under the
*> description of INFO below.) *> description of INFO below.)
*> *>
*> This subroutine may explicitly set H(i,j) = 0 for i.GT.j and *> This subroutine may explicitly set H(i,j) = 0 for i > j and
*> j = 1, 2, ... ILO-1 or j = IHI+1, IHI+2, ... N. *> j = 1, 2, ... ILO-1 or j = IHI+1, IHI+2, ... N.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] LDH *> \param[in] LDH
*> \verbatim *> \verbatim
*> LDH is INTEGER *> LDH is INTEGER
*> The leading dimension of the array H. LDH .GE. max(1,N). *> The leading dimension of the array H. LDH >= max(1,N).
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] W *> \param[out] W
@ -134,7 +134,7 @@
*> IHIZ is INTEGER *> IHIZ is INTEGER
*> Specify the rows of Z to which transformations must be *> Specify the rows of Z to which transformations must be
*> applied if WANTZ is .TRUE.. *> applied if WANTZ is .TRUE..
*> 1 .LE. ILOZ .LE. ILO; IHI .LE. IHIZ .LE. N. *> 1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in,out] Z *> \param[in,out] Z
@ -144,7 +144,7 @@
*> If WANTZ is .TRUE., then Z(ILO:IHI,ILOZ:IHIZ) is *> If WANTZ is .TRUE., then Z(ILO:IHI,ILOZ:IHIZ) is
*> replaced by Z(ILO:IHI,ILOZ:IHIZ)*U where U is the *> replaced by Z(ILO:IHI,ILOZ:IHIZ)*U where U is the
*> orthogonal Schur factor of H(ILO:IHI,ILO:IHI). *> orthogonal Schur factor of H(ILO:IHI,ILO:IHI).
*> (The output value of Z when INFO.GT.0 is given under *> (The output value of Z when INFO > 0 is given under
*> the description of INFO below.) *> the description of INFO below.)
*> \endverbatim *> \endverbatim
*> *>
@ -152,7 +152,7 @@
*> \verbatim *> \verbatim
*> LDZ is INTEGER *> LDZ is INTEGER
*> The leading dimension of the array Z. if WANTZ is .TRUE. *> The leading dimension of the array Z. if WANTZ is .TRUE.
*> then LDZ.GE.MAX(1,IHIZ). Otherwize, LDZ.GE.1. *> then LDZ >= MAX(1,IHIZ). Otherwise, LDZ >= 1.
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] WORK *> \param[out] WORK
@ -165,7 +165,7 @@
*> \param[in] LWORK *> \param[in] LWORK
*> \verbatim *> \verbatim
*> LWORK is INTEGER *> LWORK is INTEGER
*> The dimension of the array WORK. LWORK .GE. max(1,N) *> The dimension of the array WORK. LWORK >= max(1,N)
*> is sufficient, but LWORK typically as large as 6*N may *> is sufficient, but LWORK typically as large as 6*N may
*> be required for optimal performance. A workspace query *> be required for optimal performance. A workspace query
*> to determine the optimal workspace size is recommended. *> to determine the optimal workspace size is recommended.
@ -182,18 +182,18 @@
*> \verbatim *> \verbatim
*> INFO is INTEGER *> INFO is INTEGER
*> = 0: successful exit *> = 0: successful exit
*> .GT. 0: if INFO = i, ZLAQR4 failed to compute all of *> > 0: if INFO = i, ZLAQR4 failed to compute all of
*> the eigenvalues. Elements 1:ilo-1 and i+1:n of WR *> the eigenvalues. Elements 1:ilo-1 and i+1:n of WR
*> and WI contain those eigenvalues which have been *> and WI contain those eigenvalues which have been
*> successfully computed. (Failures are rare.) *> successfully computed. (Failures are rare.)
*> *>
*> If INFO .GT. 0 and WANT is .FALSE., then on exit, *> If INFO > 0 and WANT is .FALSE., then on exit,
*> the remaining unconverged eigenvalues are the eigen- *> the remaining unconverged eigenvalues are the eigen-
*> values of the upper Hessenberg matrix rows and *> values of the upper Hessenberg matrix rows and
*> columns ILO through INFO of the final, output *> columns ILO through INFO of the final, output
*> value of H. *> value of H.
*> *>
*> If INFO .GT. 0 and WANTT is .TRUE., then on exit *> If INFO > 0 and WANTT is .TRUE., then on exit
*> *>
*> (*) (initial value of H)*U = U*(final value of H) *> (*) (initial value of H)*U = U*(final value of H)
*> *>
@ -201,7 +201,7 @@
*> value of H is upper Hessenberg and triangular in *> value of H is upper Hessenberg and triangular in
*> rows and columns INFO+1 through IHI. *> rows and columns INFO+1 through IHI.
*> *>
*> If INFO .GT. 0 and WANTZ is .TRUE., then on exit *> If INFO > 0 and WANTZ is .TRUE., then on exit
*> *>
*> (final value of Z(ILO:IHI,ILOZ:IHIZ) *> (final value of Z(ILO:IHI,ILOZ:IHIZ)
*> = (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U *> = (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U
@ -209,7 +209,7 @@
*> where U is the unitary matrix in (*) (regard- *> where U is the unitary matrix in (*) (regard-
*> less of the value of WANTT.) *> less of the value of WANTT.)
*> *>
*> If INFO .GT. 0 and WANTZ is .FALSE., then Z is not *> If INFO > 0 and WANTZ is .FALSE., then Z is not
*> accessed. *> accessed.
*> \endverbatim *> \endverbatim
* *
@ -641,7 +641,7 @@
END IF END IF
END IF END IF
* *
* ==== Use up to NS of the the smallest magnatiude * ==== Use up to NS of the the smallest magnitude
* . shifts. If there aren't NS shifts available, * . shifts. If there aren't NS shifts available,
* . then use them all, possibly dropping one to * . then use them all, possibly dropping one to
* . make the number of shifts even. ==== * . make the number of shifts even. ====

View File

@ -125,7 +125,7 @@
*> \verbatim *> \verbatim
*> LDH is INTEGER *> LDH is INTEGER
*> LDH is the leading dimension of H just as declared in the *> LDH is the leading dimension of H just as declared in the
*> calling procedure. LDH.GE.MAX(1,N). *> calling procedure. LDH >= MAX(1,N).
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] ILOZ *> \param[in] ILOZ
@ -137,7 +137,7 @@
*> \verbatim *> \verbatim
*> IHIZ is INTEGER *> IHIZ is INTEGER
*> Specify the rows of Z to which transformations must be *> Specify the rows of Z to which transformations must be
*> applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N *> applied if WANTZ is .TRUE.. 1 <= ILOZ <= IHIZ <= N
*> \endverbatim *> \endverbatim
*> *>
*> \param[in,out] Z *> \param[in,out] Z
@ -153,7 +153,7 @@
*> \verbatim *> \verbatim
*> LDZ is INTEGER *> LDZ is INTEGER
*> LDA is the leading dimension of Z just as declared in *> LDA is the leading dimension of Z just as declared in
*> the calling procedure. LDZ.GE.N. *> the calling procedure. LDZ >= N.
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] V *> \param[out] V
@ -165,7 +165,7 @@
*> \verbatim *> \verbatim
*> LDV is INTEGER *> LDV is INTEGER
*> LDV is the leading dimension of V as declared in the *> LDV is the leading dimension of V as declared in the
*> calling procedure. LDV.GE.3. *> calling procedure. LDV >= 3.
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] U *> \param[out] U
@ -177,33 +177,14 @@
*> \verbatim *> \verbatim
*> LDU is INTEGER *> LDU is INTEGER
*> LDU is the leading dimension of U just as declared in the *> LDU is the leading dimension of U just as declared in the
*> in the calling subroutine. LDU.GE.3*NSHFTS-3. *> in the calling subroutine. LDU >= 3*NSHFTS-3.
*> \endverbatim
*>
*> \param[in] NH
*> \verbatim
*> NH is INTEGER
*> NH is the number of columns in array WH available for
*> workspace. NH.GE.1.
*> \endverbatim
*>
*> \param[out] WH
*> \verbatim
*> WH is COMPLEX*16 array, dimension (LDWH,NH)
*> \endverbatim
*>
*> \param[in] LDWH
*> \verbatim
*> LDWH is INTEGER
*> Leading dimension of WH just as declared in the
*> calling procedure. LDWH.GE.3*NSHFTS-3.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] NV *> \param[in] NV
*> \verbatim *> \verbatim
*> NV is INTEGER *> NV is INTEGER
*> NV is the number of rows in WV agailable for workspace. *> NV is the number of rows in WV agailable for workspace.
*> NV.GE.1. *> NV >= 1.
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] WV *> \param[out] WV
@ -215,9 +196,28 @@
*> \verbatim *> \verbatim
*> LDWV is INTEGER *> LDWV is INTEGER
*> LDWV is the leading dimension of WV as declared in the *> LDWV is the leading dimension of WV as declared in the
*> in the calling subroutine. LDWV.GE.NV. *> in the calling subroutine. LDWV >= NV.
*> \endverbatim *> \endverbatim
* *
*> \param[in] NH
*> \verbatim
*> NH is INTEGER
*> NH is the number of columns in array WH available for
*> workspace. NH >= 1.
*> \endverbatim
*>
*> \param[out] WH
*> \verbatim
*> WH is COMPLEX*16 array, dimension (LDWH,NH)
*> \endverbatim
*>
*> \param[in] LDWH
*> \verbatim
*> LDWH is INTEGER
*> Leading dimension of WH just as declared in the
*> calling procedure. LDWH >= 3*NSHFTS-3.
*> \endverbatim
*>
* Authors: * Authors:
* ======== * ========
* *

View File

@ -92,6 +92,8 @@
*> K is INTEGER *> K is INTEGER
*> The order of the matrix T (= the number of elementary *> The order of the matrix T (= the number of elementary
*> reflectors whose product defines the block reflector). *> reflectors whose product defines the block reflector).
*> If SIDE = 'L', M >= K >= 0;
*> if SIDE = 'R', N >= K >= 0.
*> \endverbatim *> \endverbatim
*> *>
*> \param[in] V *> \param[in] V

View File

@ -94,7 +94,7 @@
*> \param[in] LDC *> \param[in] LDC
*> \verbatim *> \verbatim
*> LDC is INTEGER *> LDC is INTEGER
*> The leading dimension of the array C. LDA >= max(1,M). *> The leading dimension of the array C. LDC >= max(1,M).
*> \endverbatim *> \endverbatim
*> *>
*> \param[out] WORK *> \param[out] WORK

View File

@ -103,7 +103,7 @@
* *
*> \date December 2016 *> \date December 2016
* *
*> \ingroup complex16_eig *> \ingroup complex16OTHERauxiliary
* *
* ===================================================================== * =====================================================================
SUBROUTINE ZLARFY( UPLO, N, V, INCV, TAU, C, LDC, WORK ) SUBROUTINE ZLARFY( UPLO, N, V, INCV, TAU, C, LDC, WORK )

View File

@ -143,7 +143,7 @@
*> RTOL2 is DOUBLE PRECISION *> RTOL2 is DOUBLE PRECISION
*> Parameters for bisection. *> Parameters for bisection.
*> An interval [LEFT,RIGHT] has converged if *> An interval [LEFT,RIGHT] has converged if
*> RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) ) *> RIGHT-LEFT < MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
*> \endverbatim *> \endverbatim
*> *>
*> \param[in,out] W *> \param[in,out] W

View File

@ -41,7 +41,7 @@
*> where x( i ) = abs( X( 1 + ( i - 1 )*INCX ) ). The value of sumsq is *> where x( i ) = abs( X( 1 + ( i - 1 )*INCX ) ). The value of sumsq is
*> assumed to be at least unity and the value of ssq will then satisfy *> assumed to be at least unity and the value of ssq will then satisfy
*> *>
*> 1.0 .le. ssq .le. ( sumsq + 2*n ). *> 1.0 <= ssq <= ( sumsq + 2*n ).
*> *>
*> scale is assumed to be non-negative and scl returns the value *> scale is assumed to be non-negative and scl returns the value
*> *>
@ -65,7 +65,7 @@
*> *>
*> \param[in] X *> \param[in] X
*> \verbatim *> \verbatim
*> X is COMPLEX*16 array, dimension (N) *> X is COMPLEX*16 array, dimension (1+(N-1)*INCX)
*> The vector x as described above. *> The vector x as described above.
*> x( i ) = X( 1 + ( i - 1 )*INCX ), 1 <= i <= n. *> x( i ) = X( 1 + ( i - 1 )*INCX ), 1 <= i <= n.
*> \endverbatim *> \endverbatim

View File

@ -1,3 +1,4 @@
*> \brief \b ZLASWLQ
* *
* Definition: * Definition:
* =========== * ===========
@ -18,9 +19,20 @@
*> *>
*> \verbatim *> \verbatim
*> *>
*> ZLASWLQ computes a blocked Short-Wide LQ factorization of a *> ZLASWLQ computes a blocked Tall-Skinny LQ factorization of
*> M-by-N matrix A, where N >= M: *> a complexx M-by-N matrix A for M <= N:
*> A = L * Q *>
*> A = ( L 0 ) * Q,
*>
*> where:
*>
*> Q is a n-by-N orthogonal matrix, stored on exit in an implicit
*> form in the elements above the digonal of the array A and in
*> the elemenst of the array T;
*> L is an lower-triangular M-by-M matrix stored on exit in
*> the elements on and below the diagonal of the array A.
*> 0 is a M-by-(N-M) zero matrix, if M < N, and is not stored.
*>
*> \endverbatim *> \endverbatim
* *
* Arguments: * Arguments:
@ -150,7 +162,7 @@
SUBROUTINE ZLASWLQ( M, N, MB, NB, A, LDA, T, LDT, WORK, LWORK, SUBROUTINE ZLASWLQ( M, N, MB, NB, A, LDA, T, LDT, WORK, LWORK,
$ INFO) $ INFO)
* *
* -- LAPACK computational routine (version 3.7.1) -- * -- LAPACK computational routine (version 3.9.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd. -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd. --
* June 2017 * June 2017

View File

@ -284,7 +284,8 @@
* *
* Swap A(I1, I2+1:M) with A(I2, I2+1:M) * Swap A(I1, I2+1:M) with A(I2, I2+1:M)
* *
CALL ZSWAP( M-I2, A( J1+I1-1, I2+1 ), LDA, IF( I2.LT.M )
$ CALL ZSWAP( M-I2, A( J1+I1-1, I2+1 ), LDA,
$ A( J1+I2-1, I2+1 ), LDA ) $ A( J1+I2-1, I2+1 ), LDA )
* *
* Swap A(I1, I1) with A(I2,I2) * Swap A(I1, I1) with A(I2,I2)
@ -325,6 +326,7 @@
* Compute L(J+2, J+1) = WORK( 3:M ) / T(J, J+1), * Compute L(J+2, J+1) = WORK( 3:M ) / T(J, J+1),
* where A(J, J+1) = T(J, J+1) and A(J+2:M, J) = L(J+2:M, J+1) * where A(J, J+1) = T(J, J+1) and A(J+2:M, J) = L(J+2:M, J+1)
* *
IF( J.LT.(M-1) ) THEN
IF( A( K, J+1 ).NE.ZERO ) THEN IF( A( K, J+1 ).NE.ZERO ) THEN
ALPHA = ONE / A( K, J+1 ) ALPHA = ONE / A( K, J+1 )
CALL ZCOPY( M-J-1, WORK( 3 ), 1, A( K, J+2 ), LDA ) CALL ZCOPY( M-J-1, WORK( 3 ), 1, A( K, J+2 ), LDA )
@ -334,6 +336,7 @@
$ A( K, J+2 ), LDA) $ A( K, J+2 ), LDA)
END IF END IF
END IF END IF
END IF
J = J + 1 J = J + 1
GO TO 10 GO TO 10
20 CONTINUE 20 CONTINUE
@ -432,7 +435,8 @@
* *
* Swap A(I2+1:M, I1) with A(I2+1:M, I2) * Swap A(I2+1:M, I1) with A(I2+1:M, I2)
* *
CALL ZSWAP( M-I2, A( I2+1, J1+I1-1 ), 1, IF( I2.LT.M )
$ CALL ZSWAP( M-I2, A( I2+1, J1+I1-1 ), 1,
$ A( I2+1, J1+I2-1 ), 1 ) $ A( I2+1, J1+I2-1 ), 1 )
* *
* Swap A(I1, I1) with A(I2, I2) * Swap A(I1, I1) with A(I2, I2)
@ -473,6 +477,7 @@
* Compute L(J+2, J+1) = WORK( 3:M ) / T(J, J+1), * Compute L(J+2, J+1) = WORK( 3:M ) / T(J, J+1),
* where A(J, J+1) = T(J, J+1) and A(J+2:M, J) = L(J+2:M, J+1) * where A(J, J+1) = T(J, J+1) and A(J+2:M, J) = L(J+2:M, J+1)
* *
IF( J.LT.(M-1) ) THEN
IF( A( J+1, K ).NE.ZERO ) THEN IF( A( J+1, K ).NE.ZERO ) THEN
ALPHA = ONE / A( J+1, K ) ALPHA = ONE / A( J+1, K )
CALL ZCOPY( M-J-1, WORK( 3 ), 1, A( J+2, K ), 1 ) CALL ZCOPY( M-J-1, WORK( 3 ), 1, A( J+2, K ), 1 )
@ -482,6 +487,7 @@
$ A( J+2, K ), LDA ) $ A( J+2, K ), LDA )
END IF END IF
END IF END IF
END IF
J = J + 1 J = J + 1
GO TO 30 GO TO 30
40 CONTINUE 40 CONTINUE

View File

@ -330,7 +330,7 @@
* of A and working backwards, and compute the matrix W = U12*D * of A and working backwards, and compute the matrix W = U12*D
* for use in updating A11 * for use in updating A11
* *
* Initilize the first entry of array E, where superdiagonal * Initialize the first entry of array E, where superdiagonal
* elements of D are stored * elements of D are stored
* *
E( 1 ) = CZERO E( 1 ) = CZERO
@ -658,7 +658,7 @@
* of A and working forwards, and compute the matrix W = L21*D * of A and working forwards, and compute the matrix W = L21*D
* for use in updating A22 * for use in updating A22
* *
* Initilize the unused last entry of the subdiagonal array E. * Initialize the unused last entry of the subdiagonal array E.
* *
E( N ) = CZERO E( N ) = CZERO
* *

View File

@ -261,7 +261,7 @@
* *
* Solve for U- part, lockahead for RHS(N) = +-1. This is not done * Solve for U- part, lockahead for RHS(N) = +-1. This is not done
* In BSOLVE and will hopefully give us a better estimate because * In BSOLVE and will hopefully give us a better estimate because
* any ill-conditioning of the original matrix is transfered to U * any ill-conditioning of the original matrix is transferred to U
* and not to L. U(N, N) is an approximation to sigma_min(LU). * and not to L. U(N, N) is an approximation to sigma_min(LU).
* *
CALL ZCOPY( N-1, RHS, 1, WORK, 1 ) CALL ZCOPY( N-1, RHS, 1, WORK, 1 )

View File

@ -1,3 +1,4 @@
*> \brief \b ZLATSQR
* *
* Definition: * Definition:
* =========== * ===========
@ -18,9 +19,23 @@
*> *>
*> \verbatim *> \verbatim
*> *>
*> SLATSQR computes a blocked Tall-Skinny QR factorization of *> ZLATSQR computes a blocked Tall-Skinny QR factorization of
*> an M-by-N matrix A, where M >= N: *> a complex M-by-N matrix A for M >= N:
*> A = Q * R . *>
*> A = Q * ( R ),
*> ( 0 )
*>
*> where:
*>
*> Q is a M-by-M orthogonal matrix, stored on exit in an implicit
*> form in the elements below the digonal of the array A and in
*> the elemenst of the array T;
*>
*> R is an upper-triangular N-by-N matrix, stored on exit in
*> the elements on and above the diagonal of the array A.
*>
*> 0 is a (M-N)-by-N zero matrix, and is not stored.
*>
*> \endverbatim *> \endverbatim
* *
* Arguments: * Arguments:
@ -149,10 +164,10 @@
SUBROUTINE ZLATSQR( M, N, MB, NB, A, LDA, T, LDT, WORK, SUBROUTINE ZLATSQR( M, N, MB, NB, A, LDA, T, LDT, WORK,
$ LWORK, INFO) $ LWORK, INFO)
* *
* -- LAPACK computational routine (version 3.7.0) -- * -- LAPACK computational routine (version 3.9.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd. -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd. --
* December 2016 * November 2019
* *
* .. Scalar Arguments .. * .. Scalar Arguments ..
INTEGER INFO, LDA, M, N, MB, NB, LDT, LWORK INTEGER INFO, LDA, M, N, MB, NB, LDT, LWORK

View File

@ -0,0 +1,248 @@
*> \brief \b ZLAUNHR_COL_GETRFNP
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download ZLAUNHR_COL_GETRFNP + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlaunhr_col_getrfnp.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlaunhr_col_getrfnp.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlaunhr_col_getrfnp.f">
*> [TXT]</a>
*> \endhtmlonly
*
* Definition:
* ===========
*
* SUBROUTINE ZLAUNHR_COL_GETRFNP( M, N, A, LDA, D, INFO )
*
* .. Scalar Arguments ..
* INTEGER INFO, LDA, M, N
* ..
* .. Array Arguments ..
* COMPLEX*16 A( LDA, * ), D( * )
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> ZLAUNHR_COL_GETRFNP computes the modified LU factorization without
*> pivoting of a complex general M-by-N matrix A. The factorization has
*> the form:
*>
*> A - S = L * U,
*>
*> where:
*> S is a m-by-n diagonal sign matrix with the diagonal D, so that
*> D(i) = S(i,i), 1 <= i <= min(M,N). The diagonal D is constructed
*> as D(i)=-SIGN(A(i,i)), where A(i,i) is the value after performing
*> i-1 steps of Gaussian elimination. This means that the diagonal
*> element at each step of "modified" Gaussian elimination is
*> at least one in absolute value (so that division-by-zero not
*> not possible during the division by the diagonal element);
*>
*> L is a M-by-N lower triangular matrix with unit diagonal elements
*> (lower trapezoidal if M > N);
*>
*> and U is a M-by-N upper triangular matrix
*> (upper trapezoidal if M < N).
*>
*> This routine is an auxiliary routine used in the Householder
*> reconstruction routine ZUNHR_COL. In ZUNHR_COL, this routine is
*> applied to an M-by-N matrix A with orthonormal columns, where each
*> element is bounded by one in absolute value. With the choice of
*> the matrix S above, one can show that the diagonal element at each
*> step of Gaussian elimination is the largest (in absolute value) in
*> the column on or below the diagonal, so that no pivoting is required
*> for numerical stability [1].
*>
*> For more details on the Householder reconstruction algorithm,
*> including the modified LU factorization, see [1].
*>
*> This is the blocked right-looking version of the algorithm,
*> calling Level 3 BLAS to update the submatrix. To factorize a block,
*> this routine calls the recursive routine ZLAUNHR_COL_GETRFNP2.
*>
*> [1] "Reconstructing Householder vectors from tall-skinny QR",
*> G. Ballard, J. Demmel, L. Grigori, M. Jacquelin, H.D. Nguyen,
*> E. Solomonik, J. Parallel Distrib. Comput.,
*> vol. 85, pp. 3-31, 2015.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] M
*> \verbatim
*> M is INTEGER
*> The number of rows of the matrix A. M >= 0.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> The number of columns of the matrix A. N >= 0.
*> \endverbatim
*>
*> \param[in,out] A
*> \verbatim
*> A is COMPLEX*16 array, dimension (LDA,N)
*> On entry, the M-by-N matrix to be factored.
*> On exit, the factors L and U from the factorization
*> A-S=L*U; the unit diagonal elements of L are not stored.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*> LDA is INTEGER
*> The leading dimension of the array A. LDA >= max(1,M).
*> \endverbatim
*>
*> \param[out] D
*> \verbatim
*> D is COMPLEX*16 array, dimension min(M,N)
*> The diagonal elements of the diagonal M-by-N sign matrix S,
*> D(i) = S(i,i), where 1 <= i <= min(M,N). The elements can be
*> only ( +1.0, 0.0 ) or (-1.0, 0.0 ).
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*> INFO is INTEGER
*> = 0: successful exit
*> < 0: if INFO = -i, the i-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 November 2019
*
*> \ingroup complex16GEcomputational
*
*> \par Contributors:
* ==================
*>
*> \verbatim
*>
*> November 2019, Igor Kozachenko,
*> Computer Science Division,
*> University of California, Berkeley
*>
*> \endverbatim
*
* =====================================================================
SUBROUTINE ZLAUNHR_COL_GETRFNP( M, N, A, LDA, D, INFO )
IMPLICIT NONE
*
* -- LAPACK computational routine (version 3.9.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* November 2019
*
* .. Scalar Arguments ..
INTEGER INFO, LDA, M, N
* ..
* .. Array Arguments ..
COMPLEX*16 A( LDA, * ), D( * )
* ..
*
* =====================================================================
*
* .. Parameters ..
COMPLEX*16 CONE
PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) )
* ..
* .. Local Scalars ..
INTEGER IINFO, J, JB, NB
* ..
* .. External Subroutines ..
EXTERNAL ZGEMM, ZLAUNHR_COL_GETRFNP2, ZTRSM, XERBLA
* ..
* .. External Functions ..
INTEGER ILAENV
EXTERNAL ILAENV
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX, MIN
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
INFO = 0
IF( M.LT.0 ) THEN
INFO = -1
ELSE IF( N.LT.0 ) THEN
INFO = -2
ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
INFO = -4
END IF
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'ZLAUNHR_COL_GETRFNP', -INFO )
RETURN
END IF
*
* Quick return if possible
*
IF( MIN( M, N ).EQ.0 )
$ RETURN
*
* Determine the block size for this environment.
*
NB = ILAENV( 1, 'ZLAUNHR_COL_GETRFNP', ' ', M, N, -1, -1 )
IF( NB.LE.1 .OR. NB.GE.MIN( M, N ) ) THEN
*
* Use unblocked code.
*
CALL ZLAUNHR_COL_GETRFNP2( M, N, A, LDA, D, INFO )
ELSE
*
* Use blocked code.
*
DO J = 1, MIN( M, N ), NB
JB = MIN( MIN( M, N )-J+1, NB )
*
* Factor diagonal and subdiagonal blocks.
*
CALL ZLAUNHR_COL_GETRFNP2( M-J+1, JB, A( J, J ), LDA,
$ D( J ), IINFO )
*
IF( J+JB.LE.N ) THEN
*
* Compute block row of U.
*
CALL ZTRSM( 'Left', 'Lower', 'No transpose', 'Unit', JB,
$ N-J-JB+1, CONE, A( J, J ), LDA, A( J, J+JB ),
$ LDA )
IF( J+JB.LE.M ) THEN
*
* Update trailing submatrix.
*
CALL ZGEMM( 'No transpose', 'No transpose', M-J-JB+1,
$ N-J-JB+1, JB, -CONE, A( J+JB, J ), LDA,
$ A( J, J+JB ), LDA, CONE, A( J+JB, J+JB ),
$ LDA )
END IF
END IF
END DO
END IF
RETURN
*
* End of ZLAUNHR_COL_GETRFNP
*
END

View File

@ -0,0 +1,314 @@
*> \brief \b ZLAUNHR_COL_GETRFNP2
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download ZLAUNHR_COL_GETRFNP2 + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlaunhr_col_getrfnp2.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlaunhr_col_getrfnp2.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlaunhr_col_getrfnp2.f">
*> [TXT]</a>
*> \endhtmlonly
*
* Definition:
* ===========
*
* RECURSIVE SUBROUTINE ZLAUNHR_COL_GETRFNP2( M, N, A, LDA, D, INFO )
*
* .. Scalar Arguments ..
* INTEGER INFO, LDA, M, N
* ..
* .. Array Arguments ..
* COMPLEX*16 A( LDA, * ), D( * )
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> ZLAUNHR_COL_GETRFNP2 computes the modified LU factorization without
*> pivoting of a complex general M-by-N matrix A. The factorization has
*> the form:
*>
*> A - S = L * U,
*>
*> where:
*> S is a m-by-n diagonal sign matrix with the diagonal D, so that
*> D(i) = S(i,i), 1 <= i <= min(M,N). The diagonal D is constructed
*> as D(i)=-SIGN(A(i,i)), where A(i,i) is the value after performing
*> i-1 steps of Gaussian elimination. This means that the diagonal
*> element at each step of "modified" Gaussian elimination is at
*> least one in absolute value (so that division-by-zero not
*> possible during the division by the diagonal element);
*>
*> L is a M-by-N lower triangular matrix with unit diagonal elements
*> (lower trapezoidal if M > N);
*>
*> and U is a M-by-N upper triangular matrix
*> (upper trapezoidal if M < N).
*>
*> This routine is an auxiliary routine used in the Householder
*> reconstruction routine ZUNHR_COL. In ZUNHR_COL, this routine is
*> applied to an M-by-N matrix A with orthonormal columns, where each
*> element is bounded by one in absolute value. With the choice of
*> the matrix S above, one can show that the diagonal element at each
*> step of Gaussian elimination is the largest (in absolute value) in
*> the column on or below the diagonal, so that no pivoting is required
*> for numerical stability [1].
*>
*> For more details on the Householder reconstruction algorithm,
*> including the modified LU factorization, see [1].
*>
*> This is the recursive version of the LU factorization algorithm.
*> Denote A - S by B. The algorithm divides the matrix B into four
*> submatrices:
*>
*> [ B11 | B12 ] where B11 is n1 by n1,
*> B = [ -----|----- ] B21 is (m-n1) by n1,
*> [ B21 | B22 ] B12 is n1 by n2,
*> B22 is (m-n1) by n2,
*> with n1 = min(m,n)/2, n2 = n-n1.
*>
*>
*> The subroutine calls itself to factor B11, solves for B21,
*> solves for B12, updates B22, then calls itself to factor B22.
*>
*> For more details on the recursive LU algorithm, see [2].
*>
*> ZLAUNHR_COL_GETRFNP2 is called to factorize a block by the blocked
*> routine ZLAUNHR_COL_GETRFNP, which uses blocked code calling
*. Level 3 BLAS to update the submatrix. However, ZLAUNHR_COL_GETRFNP2
*> is self-sufficient and can be used without ZLAUNHR_COL_GETRFNP.
*>
*> [1] "Reconstructing Householder vectors from tall-skinny QR",
*> G. Ballard, J. Demmel, L. Grigori, M. Jacquelin, H.D. Nguyen,
*> E. Solomonik, J. Parallel Distrib. Comput.,
*> vol. 85, pp. 3-31, 2015.
*>
*> [2] "Recursion leads to automatic variable blocking for dense linear
*> algebra algorithms", F. Gustavson, IBM J. of Res. and Dev.,
*> vol. 41, no. 6, pp. 737-755, 1997.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] M
*> \verbatim
*> M is INTEGER
*> The number of rows of the matrix A. M >= 0.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> The number of columns of the matrix A. N >= 0.
*> \endverbatim
*>
*> \param[in,out] A
*> \verbatim
*> A is COMPLEX*16 array, dimension (LDA,N)
*> On entry, the M-by-N matrix to be factored.
*> On exit, the factors L and U from the factorization
*> A-S=L*U; the unit diagonal elements of L are not stored.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*> LDA is INTEGER
*> The leading dimension of the array A. LDA >= max(1,M).
*> \endverbatim
*>
*> \param[out] D
*> \verbatim
*> D is COMPLEX*16 array, dimension min(M,N)
*> The diagonal elements of the diagonal M-by-N sign matrix S,
*> D(i) = S(i,i), where 1 <= i <= min(M,N). The elements can be
*> only ( +1.0, 0.0 ) or (-1.0, 0.0 ).
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*> INFO is INTEGER
*> = 0: successful exit
*> < 0: if INFO = -i, the i-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 November 2019
*
*> \ingroup complex16GEcomputational
*
*> \par Contributors:
* ==================
*>
*> \verbatim
*>
*> November 2019, Igor Kozachenko,
*> Computer Science Division,
*> University of California, Berkeley
*>
*> \endverbatim
*
* =====================================================================
RECURSIVE SUBROUTINE ZLAUNHR_COL_GETRFNP2( M, N, A, LDA, D, INFO )
IMPLICIT NONE
*
* -- LAPACK computational routine (version 3.9.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* November 2019
*
* .. Scalar Arguments ..
INTEGER INFO, LDA, M, N
* ..
* .. Array Arguments ..
COMPLEX*16 A( LDA, * ), D( * )
* ..
*
* =====================================================================
*
* .. Parameters ..
DOUBLE PRECISION ONE
PARAMETER ( ONE = 1.0D+0 )
COMPLEX*16 CONE
PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) )
* ..
* .. Local Scalars ..
DOUBLE PRECISION SFMIN
INTEGER I, IINFO, N1, N2
COMPLEX*16 Z
* ..
* .. External Functions ..
DOUBLE PRECISION DLAMCH
EXTERNAL DLAMCH
* ..
* .. External Subroutines ..
EXTERNAL ZGEMM, ZSCAL, ZTRSM, XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, DBLE, DCMPLX, DIMAG, DSIGN, MAX, MIN
* ..
* .. Statement Functions ..
DOUBLE PRECISION CABS1
* ..
* .. Statement Function definitions ..
CABS1( Z ) = ABS( DBLE( Z ) ) + ABS( DIMAG( Z ) )
* ..
* .. Executable Statements ..
*
* Test the input parameters
*
INFO = 0
IF( M.LT.0 ) THEN
INFO = -1
ELSE IF( N.LT.0 ) THEN
INFO = -2
ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
INFO = -4
END IF
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'ZLAUNHR_COL_GETRFNP2', -INFO )
RETURN
END IF
*
* Quick return if possible
*
IF( MIN( M, N ).EQ.0 )
$ RETURN
IF ( M.EQ.1 ) THEN
*
* One row case, (also recursion termination case),
* use unblocked code
*
* Transfer the sign
*
D( 1 ) = DCMPLX( -DSIGN( ONE, DBLE( A( 1, 1 ) ) ) )
*
* Construct the row of U
*
A( 1, 1 ) = A( 1, 1 ) - D( 1 )
*
ELSE IF( N.EQ.1 ) THEN
*
* One column case, (also recursion termination case),
* use unblocked code
*
* Transfer the sign
*
D( 1 ) = DCMPLX( -DSIGN( ONE, DBLE( A( 1, 1 ) ) ) )
*
* Construct the row of U
*
A( 1, 1 ) = A( 1, 1 ) - D( 1 )
*
* Scale the elements 2:M of the column
*
* Determine machine safe minimum
*
SFMIN = DLAMCH('S')
*
* Construct the subdiagonal elements of L
*
IF( CABS1( A( 1, 1 ) ) .GE. SFMIN ) THEN
CALL ZSCAL( M-1, CONE / A( 1, 1 ), A( 2, 1 ), 1 )
ELSE
DO I = 2, M
A( I, 1 ) = A( I, 1 ) / A( 1, 1 )
END DO
END IF
*
ELSE
*
* Divide the matrix B into four submatrices
*
N1 = MIN( M, N ) / 2
N2 = N-N1
*
* Factor B11, recursive call
*
CALL ZLAUNHR_COL_GETRFNP2( N1, N1, A, LDA, D, IINFO )
*
* Solve for B21
*
CALL ZTRSM( 'R', 'U', 'N', 'N', M-N1, N1, CONE, A, LDA,
$ A( N1+1, 1 ), LDA )
*
* Solve for B12
*
CALL ZTRSM( 'L', 'L', 'N', 'U', N1, N2, CONE, A, LDA,
$ A( 1, N1+1 ), LDA )
*
* Update B22, i.e. compute the Schur complement
* B22 := B22 - B21*B12
*
CALL ZGEMM( 'N', 'N', M-N1, N2, N1, -CONE, A( N1+1, 1 ), LDA,
$ A( 1, N1+1 ), LDA, CONE, A( N1+1, N1+1 ), LDA )
*
* Factor B22, recursive call
*
CALL ZLAUNHR_COL_GETRFNP2( M-N1, N2, A( N1+1, N1+1 ), LDA,
$ D( N1+1 ), IINFO )
*
END IF
RETURN
*
* End of ZLAUNHR_COL_GETRFNP2
*
END