Update LAPACK to 3.9.0

This commit is contained in:
Martin Kroeker 2019-12-29 23:20:29 +01:00 committed by GitHub
parent a421ab9ce2
commit e3b07ca95d
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
30 changed files with 227 additions and 198 deletions

View File

@ -135,7 +135,7 @@
*> \param[in,out] S
*> \verbatim
*> S is DOUBLE PRECISION array, dimension (N)
*> The row scale factors for A. If EQUED = 'Y', A is multiplied on
*> The scale factors for A. If EQUED = 'Y', A is multiplied on
*> the left and right by diag(S). S is an input argument if FACT =
*> 'F'; otherwise, S is an output argument. If FACT = 'F' and EQUED
*> = 'Y', each element of S must be positive. If S is output, each
@ -263,7 +263,7 @@
*> information as described below. There currently are up to three
*> pieces of information returned for each right-hand side. If
*> 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 index in ERR_BNDS_COMP(i,:) corresponds to the ith
@ -299,14 +299,14 @@
*> \param[in] NPARAMS
*> \verbatim
*> 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.
*> \endverbatim
*>
*> \param[in,out] PARAMS
*> \verbatim
*> 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
*> parameter. Only positions up to NPARAMS are accessed; defaults
*> are used for higher-numbered parameters.

View File

@ -366,7 +366,7 @@
*> information as described below. There currently are up to three
*> pieces of information returned for each right-hand side. If
*> 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 index in ERR_BNDS_COMP(i,:) corresponds to the ith
@ -402,14 +402,14 @@
*> \param[in] NPARAMS
*> \verbatim
*> 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.
*> \endverbatim
*>
*> \param[in,out] PARAMS
*> \verbatim
*> 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
*> parameter. Only positions up to NPARAMS are accessed; defaults
*> are used for higher-numbered parameters.

View File

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

View File

@ -92,9 +92,9 @@
*> dimension (LDA,N)
*> On entry, the N-by-N coefficient matrix A.
*> On exit, if iterative refinement has been successfully used
*> (INFO.EQ.0 and ITER.GE.0, see description below), then A is
*> (INFO = 0 and ITER >= 0, see description below), then A is
*> unchanged, if double precision factorization has been used
*> (INFO.EQ.0 and ITER.LT.0, see description below), then the
*> (INFO = 0 and ITER < 0, see description below), then the
*> array A contains the factors L and U from the factorization
*> A = P*L*U; the unit diagonal elements of L are not stored.
*> \endverbatim
@ -111,8 +111,8 @@
*> The pivot indices that define the permutation matrix P;
*> row i of the matrix was interchanged with row IPIV(i).
*> Corresponds either to the single precision factorization
*> (if INFO.EQ.0 and ITER.GE.0) or the double precision
*> factorization (if INFO.EQ.0 and ITER.LT.0).
*> (if INFO = 0 and ITER >= 0) or the double precision
*> factorization (if INFO = 0 and ITER < 0).
*> \endverbatim
*>
*> \param[in] B
@ -406,7 +406,7 @@
30 CONTINUE
*
* If we are at this place of the code, this is because we have
* performed ITER=ITERMAX iterations and never satisified the
* performed ITER=ITERMAX iterations and never satisfied the
* stopping criterion, set up the ITER flag accordingly and follow up
* on double precision routine.
*

View File

@ -106,9 +106,9 @@
*> triangular part of the matrix A, and the strictly upper
*> triangular part of A is not referenced.
*> On exit, if iterative refinement has been successfully used
*> (INFO.EQ.0 and ITER.GE.0, see description below), then A is
*> (INFO = 0 and ITER >= 0, see description below), then A is
*> unchanged, if double precision factorization has been used
*> (INFO.EQ.0 and ITER.LT.0, see description below), then the
*> (INFO = 0 and ITER < 0, see description below), then the
*> array A contains the factor U or L from the Cholesky
*> factorization A = U**T*U or A = L*L**T.
*> \endverbatim
@ -413,7 +413,7 @@
30 CONTINUE
*
* If we are at this place of the code, this is because we have
* performed ITER=ITERMAX iterations and never satisified the
* performed ITER=ITERMAX iterations and never satisfied the
* stopping criterion, set up the ITER flag accordingly and follow
* up on double precision routine.
*

View File

@ -233,13 +233,13 @@
*> \param[in,out] TRYRAC
*> \verbatim
*> TRYRAC is LOGICAL
*> If TRYRAC.EQ..TRUE., indicates that the code should check whether
*> If TRYRAC = .TRUE., indicates that the code should check whether
*> the tridiagonal matrix defines its eigenvalues to high relative
*> accuracy. If so, the code uses relative-accuracy preserving
*> algorithms that might be (a bit) slower depending on the matrix.
*> If the matrix does not define its eigenvalues to high relative
*> accuracy, the code can uses possibly faster algorithms.
*> If TRYRAC.EQ..FALSE., the code is not required to guarantee
*> If TRYRAC = .FALSE., the code is not required to guarantee
*> relatively accurate eigenvalues and can use the fastest possible
*> techniques.
*> On exit, a .TRUE. TRYRAC will be set to .FALSE. if the matrix

View File

@ -291,7 +291,7 @@
*
* Convert PERMUTATIONS and IPIV
*
* Apply permutaions to submatrices of upper part of A
* Apply permutations to submatrices of upper part of A
* in factorization order where i decreases from N to 1
*
I = N
@ -344,7 +344,7 @@
*
* Revert PERMUTATIONS and IPIV
*
* Apply permutaions to submatrices of upper part of A
* Apply permutations to submatrices of upper part of A
* in reverse factorization order where i increases from 1 to N
*
I = 1
@ -435,7 +435,7 @@
*
* Convert PERMUTATIONS and IPIV
*
* Apply permutaions to submatrices of lower part of A
* Apply permutations to submatrices of lower part of A
* in factorization order where k increases from 1 to N
*
I = 1
@ -488,7 +488,7 @@
*
* Revert PERMUTATIONS and IPIV
*
* Apply permutaions to submatrices of lower part of A
* Apply permutations to submatrices of lower part of A
* in reverse factorization order where i decreases from N to 1
*
I = N

View File

@ -282,7 +282,7 @@
*
* Convert PERMUTATIONS
*
* Apply permutaions to submatrices of upper part of A
* Apply permutations to submatrices of upper part of A
* in factorization order where i decreases from N to 1
*
I = N
@ -333,7 +333,7 @@
*
* Revert PERMUTATIONS
*
* Apply permutaions to submatrices of upper part of A
* Apply permutations to submatrices of upper part of A
* in reverse factorization order where i increases from 1 to N
*
I = 1
@ -423,7 +423,7 @@
*
* Convert PERMUTATIONS
*
* Apply permutaions to submatrices of lower part of A
* Apply permutations to submatrices of lower part of A
* in factorization order where i increases from 1 to N
*
I = 1
@ -474,7 +474,7 @@
*
* Revert PERMUTATIONS
*
* Apply permutaions to submatrices of lower part of A
* Apply permutations to submatrices of lower part of A
* in reverse factorization order where i decreases from N to 1
*
I = N

View File

@ -317,7 +317,7 @@
IF( .NOT.WANTZ ) THEN
CALL DSTERF( N, W, WORK( INDE ), INFO )
ELSE
* Not available in this release, and agrument checking should not
* Not available in this release, and argument checking should not
* let it getting here
RETURN
CALL DORGTR( UPLO, N, A, LDA, WORK( INDTAU ), WORK( INDWRK ),

View File

@ -385,7 +385,7 @@
IF( .NOT.WANTZ ) THEN
CALL DSTERF( N, W, WORK( INDE ), INFO )
ELSE
* Not available in this release, and agrument checking should not
* Not available in this release, and argument checking should not
* let it getting here
RETURN
CALL DSTEDC( 'I', N, W, WORK( INDE ), WORK( INDWRK ), N,

View File

@ -271,7 +271,7 @@
*> information as described below. There currently are up to three
*> pieces of information returned for each right-hand side. If
*> 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 index in ERR_BNDS_COMP(i,:) corresponds to the ith
@ -307,14 +307,14 @@
*> \param[in] NPARAMS
*> \verbatim
*> 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.
*> \endverbatim
*>
*> \param[in,out] PARAMS
*> \verbatim
*> 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
*> parameter. Only positions up to NPARAMS are accessed; defaults
*> are used for higher-numbered parameters.

View File

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

View File

@ -45,7 +45,7 @@
*> matrices.
*>
*> Aasen's 2-stage algorithm is used to factor A as
*> A = U * T * U**T, if UPLO = 'U', or
*> A = U**T * T * U, if UPLO = 'U', or
*> A = L * T * L**T, if UPLO = 'L',
*> where U (or L) is a product of permutation and unit upper (lower)
*> triangular matrices, and T is symmetric and band. The matrix T is
@ -259,7 +259,7 @@
END IF
*
*
* Compute the factorization A = U*T*U**T or A = L*T*L**T.
* Compute the factorization A = U**T*T*U or A = L*T*L**T.
*
CALL DSYTRF_AA_2STAGE( UPLO, N, A, LDA, TB, LTB, IPIV, IPIV2,
$ WORK, LWORK, INFO )

View File

@ -377,7 +377,7 @@
*> information as described below. There currently are up to three
*> pieces of information returned for each right-hand side. If
*> 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 index in ERR_BNDS_COMP(i,:) corresponds to the ith
@ -413,14 +413,14 @@
*> \param[in] NPARAMS
*> \verbatim
*> 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.
*> \endverbatim
*>
*> \param[in,out] PARAMS
*> \verbatim
*> 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
*> parameter. Only positions up to NPARAMS are accessed; defaults
*> are used for higher-numbered parameters.

View File

@ -312,7 +312,7 @@
*
* Factorize A as U*D*U**T 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
*
E( 1 ) = ZERO
@ -623,7 +623,7 @@
*
* Factorize A as L*D*L**T 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 ) = ZERO
*

View File

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

View File

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

View File

@ -39,7 +39,7 @@
*> the Bunch-Kaufman diagonal pivoting method. The form of the
*> factorization is
*>
*> A = U*D*U**T or A = L*D*L**T
*> A = U**T*D*U or A = L*D*L**T
*>
*> where U (or L) is a product of permutation and unit upper (lower)
*> triangular matrices, and D is symmetric and block diagonal with
@ -144,7 +144,7 @@
*>
*> \verbatim
*>
*> If UPLO = 'U', then A = U*D*U**T, where
*> If UPLO = 'U', then A = U**T*D*U, where
*> U = P(n)*U(n)* ... *P(k)U(k)* ...,
*> i.e., U is a product of terms P(k)*U(k), where k decreases from n to
*> 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
@ -262,7 +262,7 @@
*
IF( UPPER ) THEN
*
* Factorize A as U*D*U**T using the upper triangle of A
* Factorize A as U**T*D*U using the upper triangle of A
*
* K is the main loop index, decreasing from N to 1 in steps of
* KB, where KB is the number of columns factorized by DLASYF;

View File

@ -37,7 +37,7 @@
*> DSYTRF_AA computes the factorization of a real symmetric matrix A
*> using the Aasen's algorithm. The form of the factorization is
*>
*> A = U*T*U**T or A = L*T*L**T
*> A = U**T*T*U or A = L*T*L**T
*>
*> where U (or L) is a product of permutation and unit upper (lower)
*> triangular matrices, and T is a symmetric tridiagonal matrix.
@ -223,7 +223,7 @@
IF( UPPER ) THEN
*
* .....................................................
* Factorize A as L*D*L**T using the upper triangle of A
* Factorize A as U**T*D*U using the upper triangle of A
* .....................................................
*
* 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,
$ 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)
IPIV( J2 ) = IPIV( J2 ) + J
@ -375,7 +375,7 @@
$ A( J+1, MAX(1, J) ), LDA,
$ 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)
IPIV( J2 ) = IPIV( J2 ) + J

View File

@ -38,7 +38,7 @@
*> DSYTRF_AA_2STAGE computes the factorization of a real symmetric matrix A
*> using the Aasen's algorithm. The form of the factorization is
*>
*> A = U*T*U**T or A = L*T*L**T
*> A = U**T*T*U or A = L*T*L**T
*>
*> where U (or L) is a product of permutation and unit upper (lower)
*> triangular matrices, and T is a symmetric band matrix with the
@ -103,6 +103,22 @@
*> no error message related to LTB is issued by XERBLA.
*> \endverbatim
*>
*> \param[out] IPIV
*> \verbatim
*> IPIV is INTEGER array, dimension (N)
*> On exit, it contains the details of the interchanges, i.e.,
*> the row and column k of A were interchanged with the
*> row and column IPIV(k).
*> \endverbatim
*>
*> \param[out] IPIV2
*> \verbatim
*> IPIV2 is INTEGER array, dimension (N)
*> On exit, it contains the details of the interchanges, i.e.,
*> the row and column k of T were interchanged with the
*> row and column IPIV2(k).
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*> WORK is DOUBLE PRECISION workspace of size LWORK
@ -120,22 +136,6 @@
*> no error message related to LWORK is issued by XERBLA.
*> \endverbatim
*>
*> \param[out] IPIV
*> \verbatim
*> IPIV is INTEGER array, dimension (N)
*> On exit, it contains the details of the interchanges, i.e.,
*> the row and column k of A were interchanged with the
*> row and column IPIV(k).
*> \endverbatim
*>
*> \param[out] IPIV2
*> \verbatim
*> IPIV2 is INTEGER array, dimension (N)
*> On exit, it contains the details of the interchanges, i.e.,
*> the row and column k of T were interchanged with the
*> row and column IPIV2(k).
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*> INFO is INTEGER
@ -275,7 +275,7 @@
IF( UPPER ) THEN
*
* .....................................................
* Factorize A as L*D*L**T using the upper triangle of A
* Factorize A as U**T*D*U using the upper triangle of A
* .....................................................
*
DO J = 0, NT-1
@ -443,10 +443,12 @@ c END IF
CALL DSWAP( K-1, A( (J+1)*NB+1, I1 ), 1,
$ A( (J+1)*NB+1, I2 ), 1 )
* > Swap A(I1+1:M, I1) with A(I2, I1+1:M)
CALL DSWAP( I2-I1-1, A( I1, I1+1 ), LDA,
IF( I2.GT.(I1+1) )
$ CALL DSWAP( I2-I1-1, A( I1, I1+1 ), LDA,
$ A( I1+1, I2 ), 1 )
* > Swap A(I2+1:M, I1) with A(I2+1:M, I2)
CALL DSWAP( N-I2, A( I1, I2+1 ), LDA,
IF( I2.LT.N )
$ CALL DSWAP( N-I2, A( I1, I2+1 ), LDA,
$ A( I2, I2+1 ), LDA )
* > Swap A(I1, I1) with A(I2, I2)
PIV = A( I1, I1 )
@ -616,10 +618,12 @@ c END IF
CALL DSWAP( K-1, A( I1, (J+1)*NB+1 ), LDA,
$ A( I2, (J+1)*NB+1 ), LDA )
* > Swap A(I1+1:M, I1) with A(I2, I1+1:M)
CALL DSWAP( I2-I1-1, A( I1+1, I1 ), 1,
IF( I2.GT.(I1+1) )
$ CALL DSWAP( I2-I1-1, A( I1+1, I1 ), 1,
$ A( I2, I1+1 ), LDA )
* > Swap A(I2+1:M, I1) with A(I2+1:M, I2)
CALL DSWAP( N-I2, A( I2+1, I1 ), 1,
IF( I2.LT.N )
$ CALL DSWAP( N-I2, A( I2+1, I1 ), 1,
$ A( I2+1, I2 ), 1 )
* > Swap A(I1, I1) with A(I2, I2)
PIV = A( I1, I1 )

View File

@ -62,7 +62,7 @@
*> \param[in,out] A
*> \verbatim
*> A is DOUBLE PRECISION 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 DSYTRF.
*>
*> On exit, if INFO = 0, the (symmetric) inverse of the original
@ -82,7 +82,7 @@
*> \param[in] IPIV
*> \verbatim
*> 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 DSYTRF.
*> \endverbatim
*>

View File

@ -37,7 +37,7 @@
*> \verbatim
*>
*> DSYTRS_AA solves a system of linear equations A*X = B with a real
*> symmetric matrix A using the factorization A = U*T*U**T or
*> symmetric matrix A using the factorization A = U**T*T*U or
*> A = L*T*L**T computed by DSYTRF_AA.
*> \endverbatim
*
@ -49,7 +49,7 @@
*> UPLO is CHARACTER*1
*> Specifies whether the details of the factorization are stored
*> as an upper or lower triangular matrix.
*> = 'U': Upper triangular, form is A = U*T*U**T;
*> = 'U': Upper triangular, form is A = U**T*T*U;
*> = 'L': Lower triangular, form is A = L*T*L**T.
*> \endverbatim
*>
@ -97,14 +97,16 @@
*> The leading dimension of the array B. LDB >= max(1,N).
*> \endverbatim
*>
*> \param[in] WORK
*> \param[out] WORK
*> \verbatim
*> WORK is DOUBLE array, dimension (MAX(1,LWORK))
*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
*> \endverbatim
*>
*> \param[in] LWORK
*> \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
*> \verbatim
@ -198,9 +200,13 @@
*
IF( UPPER ) THEN
*
* Solve A*X = B, where A = U*T*U**T.
* Solve A*X = B, where A = U**T*T*U.
*
* Pivot, P**T * B
* 1) Forward substitution with U**T
*
IF( N.GT.1 ) THEN
*
* Pivot, P**T * B -> B
*
DO K = 1, N
KP = IPIV( K )
@ -208,12 +214,15 @@
$ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
END DO
*
* Compute (U \P**T * B) -> B [ (U \P**T * B) ]
* Compute U**T \ B -> B [ (U**T \P**T * B) ]
*
CALL DTRSM('L', 'U', 'T', 'U', N-1, NRHS, ONE, A( 1, 2 ), LDA,
$ B( 2, 1 ), LDB)
CALL DTRSM('L', 'U', 'T', 'U', N-1, NRHS, ONE, A( 1, 2 ),
$ 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**T \P**T * B) ]
*
CALL DLACPY( 'F', 1, N, A( 1, 1 ), LDA+1, WORK( N ), 1)
IF( N.GT.1 ) THEN
@ -223,24 +232,33 @@
CALL DGTSV( N, NRHS, WORK( 1 ), WORK( N ), WORK( 2*N ), B, LDB,
$ INFO )
*
* Compute (U**T \ B) -> B [ U**T \ (T \ (U \P**T * B) ) ]
* 3) Backward substitution with U
*
CALL DTRSM( 'L', 'U', 'N', 'U', N-1, NRHS, ONE, A( 1, 2 ), LDA,
$ B( 2, 1 ), LDB)
IF( N.GT.1 ) THEN
*
* Pivot, P * B [ P * (U**T \ (T \ (U \P**T * B) )) ]
* Compute U \ B -> B [ U \ (T \ (U**T \P**T * B) ) ]
*
CALL DTRSM( 'L', 'U', 'N', 'U', N-1, NRHS, ONE, A( 1, 2 ),
$ LDA, B( 2, 1 ), LDB)
*
* Pivot, P * B -> B [ P * (U \ (T \ (U**T \P**T * B) )) ]
*
DO K = N, 1, -1
KP = IPIV( K )
IF( KP.NE.K )
$ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
END DO
END IF
*
ELSE
*
* Solve A*X = B, where A = L*T*L**T.
*
* Pivot, P**T * B
* 1) Forward substitution with L
*
IF( N.GT.1 ) THEN
*
* Pivot, P**T * B -> B
*
DO K = 1, N
KP = IPIV( K )
@ -248,10 +266,13 @@
$ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
END DO
*
* Compute (L \P**T * B) -> B [ (L \P**T * B) ]
* Compute L \ B -> B [ (L \P**T * B) ]
*
CALL DTRSM( 'L', 'L', 'N', 'U', N-1, NRHS, ONE, A( 2, 1 ), LDA,
$ B( 2, 1 ), LDB)
CALL DTRSM( 'L', 'L', 'N', 'U', N-1, NRHS, ONE, A( 2, 1 ),
$ LDA, B( 2, 1 ), LDB)
END IF
*
* 2) Solve with triangular matrix T
*
* Compute T \ B -> B [ T \ (L \P**T * B) ]
*
@ -263,18 +284,23 @@
CALL DGTSV( N, NRHS, WORK( 1 ), WORK(N), WORK( 2*N ), B, LDB,
$ INFO)
*
* 3) Backward substitution with L**T
*
IF( N.GT.1 ) THEN
*
* Compute (L**T \ B) -> B [ L**T \ (T \ (L \P**T * B) ) ]
*
CALL DTRSM( 'L', 'L', 'T', 'U', N-1, NRHS, ONE, A( 2, 1 ), LDA,
$ B( 2, 1 ), LDB)
CALL DTRSM( 'L', 'L', 'T', 'U', N-1, NRHS, ONE, A( 2, 1 ),
$ LDA, B( 2, 1 ), LDB)
*
* Pivot, P * B [ P * (L**T \ (T \ (L \P**T * B) )) ]
* Pivot, P * B -> B [ P * (L**T \ (T \ (L \P**T * B) )) ]
*
DO K = N, 1, -1
KP = IPIV( K )
IF( KP.NE.K )
$ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
END DO
END IF
*
END IF
*

View File

@ -36,7 +36,7 @@
*> \verbatim
*>
*> DSYTRS_AA_2STAGE solves a system of linear equations A*X = B with a real
*> symmetric matrix A using the factorization A = U*T*U**T or
*> symmetric matrix A using the factorization A = U**T*T*U or
*> A = L*T*L**T computed by DSYTRF_AA_2STAGE.
*> \endverbatim
*
@ -48,7 +48,7 @@
*> UPLO is CHARACTER*1
*> Specifies whether the details of the factorization are stored
*> as an upper or lower triangular matrix.
*> = 'U': Upper triangular, form is A = U*T*U**T;
*> = 'U': Upper triangular, form is A = U**T*T*U;
*> = 'L': Lower triangular, form is A = L*T*L**T.
*> \endverbatim
*>
@ -208,15 +208,15 @@
*
IF( UPPER ) THEN
*
* Solve A*X = B, where A = U*T*U**T.
* Solve A*X = B, where A = U**T*T*U.
*
IF( N.GT.NB ) THEN
*
* Pivot, P**T * B
* Pivot, P**T * B -> B
*
CALL DLASWP( NRHS, B, LDB, NB+1, N, IPIV, 1 )
*
* Compute (U**T \P**T * B) -> B [ (U**T \P**T * B) ]
* Compute (U**T \ B) -> B [ (U**T \P**T * B) ]
*
CALL DTRSM( 'L', 'U', 'T', 'U', N-NB, NRHS, ONE, A(1, NB+1),
$ LDA, B(NB+1, 1), LDB)
@ -234,7 +234,7 @@
CALL DTRSM( 'L', 'U', 'N', 'U', N-NB, NRHS, ONE, A(1, NB+1),
$ 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**T \P**T * B) )) ]
*
CALL DLASWP( NRHS, B, LDB, NB+1, N, IPIV, -1 )
*
@ -246,11 +246,11 @@
*
IF( N.GT.NB ) THEN
*
* Pivot, P**T * B
* Pivot, P**T * B -> B
*
CALL DLASWP( 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 DTRSM( 'L', 'L', 'N', 'U', N-NB, NRHS, ONE, A(NB+1, 1),
$ LDA, B(NB+1, 1), LDB)
@ -268,7 +268,7 @@
CALL DTRSM( 'L', 'L', 'T', 'U', N-NB, NRHS, ONE, A(NB+1, 1),
$ LDA, B(NB+1, 1), LDB)
*
* Pivot, P * B [ P * (L**T \ (T \ (L \P**T * B) )) ]
* Pivot, P * B -> B [ P * (L**T \ (T \ (L \P**T * B) )) ]
*
CALL DLASWP( NRHS, B, LDB, NB+1, N, IPIV, -1 )
*

View File

@ -71,7 +71,7 @@
*> R * B**T + L * E**T = scale * -F
*>
*> This case is used to compute an estimate of Dif[(A, D), (B, E)] =
*> sigma_min(Z) using reverse communicaton with DLACON.
*> sigma_min(Z) using reverse communication with DLACON.
*>
*> DTGSY2 also (IJOB >= 1) contributes to the computation in DTGSYL
*> of an upper bound on the separation between to matrix pairs. Then
@ -85,7 +85,7 @@
*> \param[in] TRANS
*> \verbatim
*> TRANS is CHARACTER*1
*> = 'N', solve the generalized Sylvester equation (1).
*> = 'N': solve the generalized Sylvester equation (1).
*> = 'T': solve the 'transposed' system (3).
*> \endverbatim
*>

View File

@ -88,8 +88,8 @@
*> \param[in] TRANS
*> \verbatim
*> TRANS is CHARACTER*1
*> = 'N', solve the generalized Sylvester equation (1).
*> = 'T', solve the 'transposed' system (3).
*> = 'N': solve the generalized Sylvester equation (1).
*> = 'T': solve the 'transposed' system (3).
*> \endverbatim
*>
*> \param[in] IJOB

View File

@ -94,7 +94,7 @@
*>
*> \param[in] V
*> \verbatim
*> V is DOUBLE PRECISION array, dimension (LDA,K)
*> V is DOUBLE PRECISION array, dimension (LDV,K)
*> The i-th row must contain the vector which defines the
*> elementary reflector H(i), for i = 1,2,...,k, as returned by
*> DTPLQT in B. See Further Details.

View File

@ -94,7 +94,7 @@
*>
*> \param[in] V
*> \verbatim
*> V is DOUBLE PRECISION array, dimension (LDA,K)
*> V is DOUBLE PRECISION array, dimension (LDV,K)
*> The i-th column must contain the vector which defines the
*> elementary reflector H(i), for i = 1,2,...,k, as returned by
*> CTPQRT in B. See Further Details.

View File

@ -152,8 +152,8 @@
*> \verbatim
*> LDA is INTEGER
*> The leading dimension of the array A.
*> If SIDE = 'L', LDC >= max(1,K);
*> If SIDE = 'R', LDC >= max(1,M).
*> If SIDE = 'L', LDA >= max(1,K);
*> If SIDE = 'R', LDA >= max(1,M).
*> \endverbatim
*>
*> \param[in,out] B