diff options
| author | julielangou <julie@cs.utk.edu> | 2016-11-23 16:27:46 -0800 |
|---|---|---|
| committer | GitHub <noreply@github.com> | 2016-11-23 16:27:46 -0800 |
| commit | 9564fbbc99be911d182a616d8535f03ecc2259e3 (patch) | |
| tree | 76014a5a8d834be819dafadc0efa1d75b38e8668 | |
| parent | 14f49ebfde6908a959f7bcefbdcb2a95ab68c1f3 (diff) | |
| parent | 6714497b881128cd61207ee9f715477c1faaefd5 (diff) | |
Merge pull request #89 from iyamazaki/lapack-aasen
Lapack aasen , Add Aasen's for complex symmetric matrix, and a few cleanups/fixes in testers.
52 files changed, 6166 insertions, 483 deletions
diff --git a/SRC/Makefile b/SRC/Makefile index dfb62ae7..9bb40582 100644 --- a/SRC/Makefile +++ b/SRC/Makefile @@ -229,7 +229,7 @@ CLASRC = \ clarf.o clarfb.o clarfg.o clarft.o clarfgp.o \ clarfx.o clarfy.o clargv.o clarnv.o clarrv.o clartg.o clartv.o \ clarz.o clarzb.o clarzt.o clascl.o claset.o clasr.o classq.o \ - claswp.o clasyf.o clasyf_rook.o clasyf_rk.o \ + claswp.o clasyf.o clasyf_rook.o clasyf_rk.o clasyf_aa.o \ clatbs.o clatdf.o clatps.o clatrd.o clatrs.o clatrz.o \ clauu2.o clauum.o cpbcon.o cpbequ.o cpbrfs.o cpbstf.o cpbsv.o \ cpbsvx.o cpbtf2.o cpbtrf.o cpbtrs.o cpocon.o cpoequ.o cporfs.o \ @@ -245,8 +245,8 @@ CLASRC = \ csyconv.o csyconvf.o csyconvf_rook.o \ csytf2_rook.o csytrf_rook.o csytrs_rook.o \ csytri_rook.o csycon_rook.o csysv_rook.o \ - csytf2_rk.o csytrf_rk.o csytrs_3.o \ - csytri_3.o csytri_3x.o csycon_3.o csysv_rk.o \ + csytf2_rk.o csytrf_rk.o csytrf_aa.o csytrs_3.o csytrs_aa.o \ + csytri_3.o csytri_3x.o csycon_3.o csysv_rk.o csysv_aa.o \ ctbcon.o ctbrfs.o ctbtrs.o ctgevc.o ctgex2.o \ ctgexc.o ctgsen.o ctgsja.o ctgsna.o ctgsy2.o ctgsyl.o ctpcon.o \ ctprfs.o ctptri.o \ @@ -424,7 +424,7 @@ ZLASRC = \ zlarfg.o zlarft.o zlarfgp.o \ zlarfx.o zlarfy.o zlargv.o zlarnv.o zlarrv.o zlartg.o zlartv.o \ zlarz.o zlarzb.o zlarzt.o zlascl.o zlaset.o zlasr.o \ - zlassq.o zlaswp.o zlasyf.o zlasyf_rook.o zlasyf_rk.o \ + zlassq.o zlaswp.o zlasyf.o zlasyf_rook.o zlasyf_rk.o zlasyf_aa.o \ zlatbs.o zlatdf.o zlatps.o zlatrd.o zlatrs.o zlatrz.o zlauu2.o \ zlauum.o zpbcon.o zpbequ.o zpbrfs.o zpbstf.o zpbsv.o \ zpbsvx.o zpbtf2.o zpbtrf.o zpbtrs.o zpocon.o zpoequ.o zporfs.o \ @@ -438,10 +438,10 @@ ZLASRC = \ zsyr.o zsyrfs.o zsysv.o zsysvx.o zsytf2.o zsytrf.o zsytri.o zsytri2.o zsytri2x.o \ zsyswapr.o zsytrs.o zsytrs2.o \ zsyconv.o zsyconvf.o zsyconvf_rook.o \ - zsytf2_rook.o zsytrf_rook.o zsytrs_rook.o \ + zsytf2_rook.o zsytrf_rook.o zsytrs_rook.o zsytrs_aa.o \ zsytri_rook.o zsycon_rook.o zsysv_rook.o \ - zsytf2_rk.o zsytrf_rk.o zsytrs_3.o \ - zsytri_3.o zsytri_3x.o zsycon_3.o zsysv_rk.o \ + zsytf2_rk.o zsytrf_rk.o zsytrf_aa.o zsytrs_3.o \ + zsytri_3.o zsytri_3x.o zsycon_3.o zsysv_rk.o zsysv_aa.o \ ztbcon.o ztbrfs.o ztbtrs.o ztgevc.o ztgex2.o \ ztgexc.o ztgsen.o ztgsja.o ztgsna.o ztgsy2.o ztgsyl.o ztpcon.o \ ztprfs.o ztptri.o \ diff --git a/SRC/chesv_aa.f b/SRC/chesv_aa.f index aae45e60..642c9932 100644 --- a/SRC/chesv_aa.f +++ b/SRC/chesv_aa.f @@ -19,7 +19,7 @@ * =========== * * SUBROUTINE CHESV_AA( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK, -* LWORK, INFO ) +* LWORK, INFO ) * * .. Scalar Arguments .. * CHARACTER UPLO @@ -27,7 +27,7 @@ * .. * .. Array Arguments .. * INTEGER IPIV( * ) -* COMPLEX A( LDA, * ), B( LDB, * ), WORK( * ) +* COMPLEX A( LDA, * ), B( LDB, * ), WORK( * ) * .. * * @@ -126,9 +126,9 @@ *> \param[in] LWORK *> \verbatim *> LWORK is INTEGER -*> The length of WORK. LWORK >= 1, and for best performance -*> LWORK >= max(1,N*NB), where NB is the optimal blocksize for -*> CHETRF. +*> The length of WORK. LWORK >= MAX(1,2*N,3*N-2), and for best +*> performance LWORK >= MAX(1,N*NB), where NB is the optimal +*> blocksize for CHETRF. *> for LWORK < N, TRS will be done with Level BLAS 2 *> for LWORK >= N, TRS will be done with Level BLAS 3 *> @@ -162,7 +162,7 @@ * * ===================================================================== SUBROUTINE CHESV_AA( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK, - $ LWORK, INFO ) + $ LWORK, INFO ) * * -- LAPACK driver routine (version 3.7.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -175,7 +175,7 @@ * .. * .. Array Arguments .. INTEGER IPIV( * ) - COMPLEX A( LDA, * ), B( LDB, * ), WORK( * ) + COMPLEX A( LDA, * ), B( LDB, * ), WORK( * ) * .. * * ===================================================================== diff --git a/SRC/chetrf_aa.f b/SRC/chetrf_aa.f index 883e9f00..3dd7c036 100644 --- a/SRC/chetrf_aa.f +++ b/SRC/chetrf_aa.f @@ -37,7 +37,7 @@ *> CHETRF_AA computes the factorization of a complex hermitian 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*U**H or A = L*T*L**H *> *> where U (or L) is a product of permutation and unit upper (lower) *> triangular matrices, and T is a hermitian tridiagonal matrix. @@ -230,7 +230,7 @@ IF( UPPER ) THEN * * ..................................................... -* Factorize A as L*D*L**T using the upper triangle of A +* Factorize A as L*D*L**H using the upper triangle of A * ..................................................... * * copy first row A(1, 1:N) into H(1:n) (stored in WORK(1:N)) @@ -353,7 +353,7 @@ 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 * ..................................................... * * copy first column A(1:N, 1) into H(1:N, 1) diff --git a/SRC/chetrs_aa.f b/SRC/chetrs_aa.f index 47d7825d..ef55c8f0 100644 --- a/SRC/chetrs_aa.f +++ b/SRC/chetrs_aa.f @@ -27,7 +27,7 @@ * .. * .. Array Arguments .. * INTEGER IPIV( * ) -* COMPLEX A( LDA, * ), B( LDB, * ), WORK( * ) +* COMPLEX A( LDA, * ), B( LDB, * ), WORK( * ) * .. * * @@ -36,9 +36,9 @@ *> *> \verbatim *> -*> CHETRS_AA solves a system of linear equations A*X = B with a real -*> hermitian matrix A using the factorization A = U*T*U**T or -*> A = L*T*L**T computed by CHETRF_AA. +*> CHETRS_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 +*> A = L*T*L**H computed by CHETRF_AA. *> \endverbatim * * Arguments: @@ -49,8 +49,8 @@ *> 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; -*> = 'L': Lower triangular, form is A = L*T*L**T. +*> = 'U': Upper triangular, form is A = U*T*U**H; +*> = 'L': Lower triangular, form is A = L*T*L**H. *> \endverbatim *> *> \param[in] N @@ -104,7 +104,7 @@ *> *> \param[in] LWORK *> \verbatim -*> LWORK is INTEGER, LWORK >= 3*N-2. +*> LWORK is INTEGER, LWORK >= MAX(1,3*N-2). *> *> \param[out] INFO *> \verbatim @@ -142,12 +142,12 @@ * .. * .. Array Arguments .. INTEGER IPIV( * ) - COMPLEX A( LDA, * ), B( LDB, * ), WORK( * ) + COMPLEX A( LDA, * ), B( LDB, * ), WORK( * ) * .. * * ===================================================================== * - COMPLEX ONE + COMPLEX ONE PARAMETER ( ONE = 1.0E+0 ) * .. * .. Local Scalars .. @@ -179,7 +179,7 @@ INFO = -5 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN INFO = -8 - ELSE IF( LWORK.LT.(3*N-2) .AND. .NOT.LQUERY ) THEN + ELSE IF( LWORK.LT.MAX( 1, 3*N-2 ) .AND. .NOT.LQUERY ) THEN INFO = -10 END IF IF( INFO.NE.0 ) THEN diff --git a/SRC/clahef_aa.f b/SRC/clahef_aa.f index 01e8f98c..81ca9023 100644 --- a/SRC/clahef_aa.f +++ b/SRC/clahef_aa.f @@ -19,7 +19,7 @@ * =========== * * SUBROUTINE CLAHEF_AA( UPLO, J1, M, NB, A, LDA, IPIV, -* H, LDH, WORK, INFO ) +* H, LDH, WORK, INFO ) * * .. Scalar Arguments .. * CHARACTER UPLO @@ -27,7 +27,7 @@ * .. * .. Array Arguments .. * INTEGER IPIV( * ) -* COMPLEX A( LDA, * ), H( LDH, * ), WORK( * ) +* COMPLEX A( LDA, * ), H( LDH, * ), WORK( * ) * .. * * @@ -36,7 +36,7 @@ *> *> \verbatim *> -*> DLATRF_AA factorizes a panel of a real hermitian matrix A using +*> CLAHEF_AA factorizes a panel of a complex hermitian matrix A using *> the Aasen's algorithm. The panel consists of a set of NB rows of A *> when UPLO is U, or a set of NB columns when UPLO is L. *> @@ -46,7 +46,7 @@ *> which is used to factorize the first panel. *> *> The resulting J-th row of U, or J-th column of L, is stored in the -*> (J-1)-th row, or column, of A (without the unit diatonals), while +*> (J-1)-th row, or column, of A (without the unit diagonals), while *> the diagonal and subdiagonal of A are overwritten by those of T. *> *> \endverbatim @@ -152,7 +152,7 @@ * * ===================================================================== SUBROUTINE CLAHEF_AA( UPLO, J1, M, NB, A, LDA, IPIV, - $ H, LDH, WORK, INFO ) + $ H, LDH, WORK, INFO ) * * -- LAPACK computational routine (version 3.7.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -167,17 +167,17 @@ * .. * .. Array Arguments .. INTEGER IPIV( * ) - COMPLEX A( LDA, * ), H( LDH, * ), WORK( * ) + COMPLEX A( LDA, * ), H( LDH, * ), WORK( * ) * .. * * ===================================================================== * .. Parameters .. - COMPLEX ZERO, ONE + COMPLEX ZERO, ONE PARAMETER ( ZERO = (0.0E+0, 0.0E+0), ONE = (1.0E+0, 0.0E+0) ) * * .. Local Scalars .. INTEGER J, K, K1, I1, I2 - COMPLEX PIV, ALPHA + COMPLEX PIV, ALPHA * .. * .. External Functions .. LOGICAL LSAME diff --git a/SRC/clasyf_aa.f b/SRC/clasyf_aa.f new file mode 100644 index 00000000..b69cc547 --- /dev/null +++ b/SRC/clasyf_aa.f @@ -0,0 +1,506 @@ +*> \brief \b CLASYF_AA +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download CLASYF_AA + dependencies +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clasyf_aa.f"> +*> [TGZ]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clasyf_aa.f"> +*> [ZIP]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clasyf_aa.f"> +*> [TXT]</a> +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE CLASYF_AA( UPLO, J1, M, NB, A, LDA, IPIV, +* H, LDH, WORK, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER UPLO +* INTEGER J1, M, NB, LDA, LDH, INFO +* .. +* .. Array Arguments .. +* INTEGER IPIV( * ) +* COMPLEX A( LDA, * ), H( LDH, * ), WORK( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DLATRF_AA factorizes a panel of a complex symmetric matrix A using +*> the Aasen's algorithm. The panel consists of a set of NB rows of A +*> when UPLO is U, or a set of NB columns when UPLO is L. +*> +*> In order to factorize the panel, the Aasen's algorithm requires the +*> last row, or column, of the previous panel. The first row, or column, +*> of A is set to be the first row, or column, of an identity matrix, +*> which is used to factorize the first panel. +*> +*> The resulting J-th row of U, or J-th column of L, is stored in the +*> (J-1)-th row, or column, of A (without the unit diagonals), while +*> the diagonal and subdiagonal of A are overwritten by those of T. +*> +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] UPLO +*> \verbatim +*> UPLO is CHARACTER*1 +*> = 'U': Upper triangle of A is stored; +*> = 'L': Lower triangle of A is stored. +*> \endverbatim +*> +*> \param[in] J1 +*> \verbatim +*> J1 is INTEGER +*> The location of the first row, or column, of the panel +*> within the submatrix of A, passed to this routine, e.g., +*> when called by CSYTRF_AA, for the first panel, J1 is 1, +*> while for the remaining panels, J1 is 2. +*> \endverbatim +*> +*> \param[in] M +*> \verbatim +*> M is INTEGER +*> The dimension of the submatrix. M >= 0. +*> \endverbatim +*> +*> \param[in] NB +*> \verbatim +*> NB is INTEGER +*> The dimension of the panel to be facotorized. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is REAL array, dimension (LDA,M) for +*> the first panel, while dimension (LDA,M+1) for the +*> remaining panels. +*> +*> On entry, A contains the last row, or column, of +*> the previous panel, and the trailing submatrix of A +*> to be factorized, except for the first panel, only +*> the panel is passed. +*> +*> On exit, the leading panel is factorized. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N). +*> \endverbatim +*> +*> \param[out] IPIV +*> \verbatim +*> IPIV is INTEGER array, dimension (N) +*> Details of the row and column interchanges, +*> the row and column k were interchanged with the row and +*> column IPIV(k). +*> \endverbatim +*> +*> \param[in,out] H +*> \verbatim +*> H is REAL workspace, dimension (LDH,NB). +*> +*> \endverbatim +*> +*> \param[in] LDH +*> \verbatim +*> LDH is INTEGER +*> The leading dimension of the workspace H. LDH >= max(1,M). +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is REAL workspace, dimension (M). +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit +*> < 0: if INFO = -i, the i-th argument had an illegal value +*> > 0: if INFO = i, D(i,i) is exactly zero. The factorization +*> has been completed, but the block diagonal matrix D is +*> exactly singular, and division by zero will occur if it +*> is used to solve a system of equations. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2016 +* +*> \ingroup complexSYcomputational +* +* ===================================================================== + SUBROUTINE CLASYF_AA( UPLO, J1, M, NB, A, LDA, IPIV, + $ H, LDH, WORK, INFO ) +* +* -- LAPACK computational routine (version 3.7.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2016 +* + IMPLICIT NONE +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER M, NB, J1, LDA, LDH, INFO +* .. +* .. Array Arguments .. + INTEGER IPIV( * ) + COMPLEX A( LDA, * ), H( LDH, * ), WORK( * ) +* .. +* +* ===================================================================== +* .. Parameters .. + COMPLEX ZERO, ONE + PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) +* +* .. Local Scalars .. + INTEGER J, K, K1, I1, I2 + COMPLEX PIV, ALPHA +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER ICAMAX, ILAENV + EXTERNAL LSAME, ILAENV, ICAMAX +* .. +* .. External Subroutines .. + EXTERNAL XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX +* .. +* .. Executable Statements .. +* + INFO = 0 + J = 1 +* +* K1 is the first column of the panel to be factorized +* i.e., K1 is 2 for the first block column, and 1 for the rest of the blocks +* + K1 = (2-J1)+1 +* + IF( LSAME( UPLO, 'U' ) ) THEN +* +* ..................................................... +* Factorize A as U**T*D*U using the upper triangle of A +* ..................................................... +* + 10 CONTINUE + IF ( J.GT.MIN(M, NB) ) + $ GO TO 20 +* +* K is the column to be factorized +* when being called from CSYTRF_AA, +* > for the first block column, J1 is 1, hence J1+J-1 is J, +* > for the rest of the columns, J1 is 2, and J1+J-1 is J+1, +* + K = J1+J-1 +* +* H(J:N, J) := A(J, J:N) - H(J:N, 1:(J-1)) * L(J1:(J-1), J), +* where H(J:N, J) has been initialized to be A(J, J:N) +* + IF( K.GT.2 ) THEN +* +* K is the column to be factorized +* > for the first block column, K is J, skipping the first two +* columns +* > for the rest of the columns, K is J+1, skipping only the +* first column +* + CALL CGEMV( 'No transpose', M-J+1, J-K1, + $ -ONE, H( J, K1 ), LDH, + $ A( 1, J ), 1, + $ ONE, H( J, J ), 1 ) + END IF +* +* Copy H(i:n, i) into WORK +* + CALL CCOPY( M-J+1, H( J, J ), 1, WORK( 1 ), 1 ) +* + IF( J.GT.K1 ) THEN +* +* Compute WORK := WORK - L(J-1, J:N) * T(J-1,J), +* where A(J-1, J) stores T(J-1, J) and A(J-2, J:N) stores U(J-1, J:N) +* + ALPHA = -A( K-1, J ) + CALL CAXPY( M-J+1, ALPHA, A( K-2, J ), LDA, WORK( 1 ), 1 ) + END IF +* +* Set A(J, J) = T(J, J) +* + A( K, J ) = WORK( 1 ) +* + IF( J.LT.M ) THEN +* +* Compute WORK(2:N) = T(J, J) L(J, (J+1):N) +* where A(J, J) stores T(J, J) and A(J-1, (J+1):N) stores U(J, (J+1):N) +* + IF( K.GT.1 ) THEN + ALPHA = -A( K, J ) + CALL CAXPY( M-J, ALPHA, A( K-1, J+1 ), LDA, + $ WORK( 2 ), 1 ) + ENDIF +* +* Find max(|WORK(2:n)|) +* + I2 = ICAMAX( M-J, WORK( 2 ), 1 ) + 1 + PIV = WORK( I2 ) +* +* Apply symmetric pivot +* + IF( (I2.NE.2) .AND. (PIV.NE.0) ) THEN +* +* Swap WORK(I1) and WORK(I2) +* + I1 = 2 + WORK( I2 ) = WORK( I1 ) + WORK( I1 ) = PIV +* +* Swap A(I1, I1+1:N) with A(I1+1:N, I2) +* + I1 = I1+J-1 + I2 = I2+J-1 + CALL CSWAP( I2-I1-1, A( J1+I1-1, I1+1 ), LDA, + $ A( J1+I1, I2 ), 1 ) +* +* Swap A(I1, I2+1:N) with A(I2, I2+1:N) +* + CALL CSWAP( M-I2, A( J1+I1-1, I2+1 ), LDA, + $ A( J1+I2-1, I2+1 ), LDA ) +* +* Swap A(I1, I1) with A(I2,I2) +* + PIV = A( I1+J1-1, I1 ) + A( J1+I1-1, I1 ) = A( J1+I2-1, I2 ) + A( J1+I2-1, I2 ) = PIV +* +* Swap H(I1, 1:J1) with H(I2, 1:J1) +* + CALL CSWAP( I1-1, H( I1, 1 ), LDH, H( I2, 1 ), LDH ) + IPIV( I1 ) = I2 +* + IF( I1.GT.(K1-1) ) THEN +* +* Swap L(1:I1-1, I1) with L(1:I1-1, I2), +* skipping the first column +* + CALL CSWAP( I1-K1+1, A( 1, I1 ), 1, + $ A( 1, I2 ), 1 ) + END IF + ELSE + IPIV( J+1 ) = J+1 + ENDIF +* +* Set A(J, J+1) = T(J, J+1) +* + A( K, J+1 ) = WORK( 2 ) + IF( (A( K, J ).EQ.ZERO ) .AND. + $ ( (J.EQ.M) .OR. (A( K, J+1 ).EQ.ZERO))) THEN + IF(INFO .EQ. 0) THEN + INFO = J + ENDIF + END IF +* + IF( J.LT.NB ) THEN +* +* Copy A(J+1:N, J+1) into H(J:N, J), +* + CALL CCOPY( M-J, A( K+1, J+1 ), LDA, + $ H( J+1, J+1 ), 1 ) + END IF +* +* 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) +* + IF( A( K, J+1 ).NE.ZERO ) THEN + ALPHA = ONE / A( K, J+1 ) + CALL CCOPY( M-J-1, WORK( 3 ), 1, A( K, J+2 ), LDA ) + CALL CSCAL( M-J-1, ALPHA, A( K, J+2 ), LDA ) + ELSE + CALL CLASET( 'Full', 1, M-J-1, ZERO, ZERO, + $ A( K, J+2 ), LDA) + END IF + ELSE + IF( (A( K, J ).EQ.ZERO) .AND. (INFO.EQ.0) ) THEN + INFO = J + END IF + END IF + J = J + 1 + GO TO 10 + 20 CONTINUE +* + ELSE +* +* ..................................................... +* Factorize A as L*D*L**T using the lower triangle of A +* ..................................................... +* + 30 CONTINUE + IF( J.GT.MIN( M, NB ) ) + $ GO TO 40 +* +* K is the column to be factorized +* when being called from CSYTRF_AA, +* > for the first block column, J1 is 1, hence J1+J-1 is J, +* > for the rest of the columns, J1 is 2, and J1+J-1 is J+1, +* + K = J1+J-1 +* +* H(J:N, J) := A(J:N, J) - H(J:N, 1:(J-1)) * L(J, J1:(J-1))^T, +* where H(J:N, J) has been initialized to be A(J:N, J) +* + IF( K.GT.2 ) THEN +* +* K is the column to be factorized +* > for the first block column, K is J, skipping the first two +* columns +* > for the rest of the columns, K is J+1, skipping only the +* first column +* + CALL CGEMV( 'No transpose', M-J+1, J-K1, + $ -ONE, H( J, K1 ), LDH, + $ A( J, 1 ), LDA, + $ ONE, H( J, J ), 1 ) + END IF +* +* Copy H(J:N, J) into WORK +* + CALL CCOPY( M-J+1, H( J, J ), 1, WORK( 1 ), 1 ) +* + IF( J.GT.K1 ) THEN +* +* Compute WORK := WORK - L(J:N, J-1) * T(J-1,J), +* where A(J-1, J) = T(J-1, J) and A(J, J-2) = L(J, J-1) +* + ALPHA = -A( J, K-1 ) + CALL CAXPY( M-J+1, ALPHA, A( J, K-2 ), 1, WORK( 1 ), 1 ) + END IF +* +* Set A(J, J) = T(J, J) +* + A( J, K ) = WORK( 1 ) +* + IF( J.LT.M ) THEN +* +* Compute WORK(2:N) = T(J, J) L((J+1):N, J) +* where A(J, J) = T(J, J) and A((J+1):N, J-1) = L((J+1):N, J) +* + IF( K.GT.1 ) THEN + ALPHA = -A( J, K ) + CALL CAXPY( M-J, ALPHA, A( J+1, K-1 ), 1, + $ WORK( 2 ), 1 ) + ENDIF +* +* Find max(|WORK(2:n)|) +* + I2 = ICAMAX( M-J, WORK( 2 ), 1 ) + 1 + PIV = WORK( I2 ) +* +* Apply symmetric pivot +* + IF( (I2.NE.2) .AND. (PIV.NE.0) ) THEN +* +* Swap WORK(I1) and WORK(I2) +* + I1 = 2 + WORK( I2 ) = WORK( I1 ) + WORK( I1 ) = PIV +* +* Swap A(I1+1:N, I1) with A(I2, I1+1:N) +* + I1 = I1+J-1 + I2 = I2+J-1 + CALL CSWAP( I2-I1-1, A( I1+1, J1+I1-1 ), 1, + $ A( I2, J1+I1 ), LDA ) +* +* Swap A(I2+1:N, I1) with A(I2+1:N, I2) +* + CALL CSWAP( M-I2, A( I2+1, J1+I1-1 ), 1, + $ A( I2+1, J1+I2-1 ), 1 ) +* +* Swap A(I1, I1) with A(I2, I2) +* + PIV = A( I1, J1+I1-1 ) + A( I1, J1+I1-1 ) = A( I2, J1+I2-1 ) + A( I2, J1+I2-1 ) = PIV +* +* Swap H(I1, I1:J1) with H(I2, I2:J1) +* + CALL CSWAP( I1-1, H( I1, 1 ), LDH, H( I2, 1 ), LDH ) + IPIV( I1 ) = I2 +* + IF( I1.GT.(K1-1) ) THEN +* +* Swap L(1:I1-1, I1) with L(1:I1-1, I2), +* skipping the first column +* + CALL CSWAP( I1-K1+1, A( I1, 1 ), LDA, + $ A( I2, 1 ), LDA ) + END IF + ELSE + IPIV( J+1 ) = J+1 + ENDIF +* +* Set A(J+1, J) = T(J+1, J) +* + A( J+1, K ) = WORK( 2 ) + IF( (A( J, K ).EQ.ZERO) .AND. + $ ( (J.EQ.M) .OR. (A( J+1, K ).EQ.ZERO)) ) THEN + IF (INFO .EQ. 0) + $ INFO = J + END IF +* + IF( J.LT.NB ) THEN +* +* Copy A(J+1:N, J+1) into H(J+1:N, J), +* + CALL CCOPY( M-J, A( J+1, K+1 ), 1, + $ H( J+1, J+1 ), 1 ) + END IF +* +* 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) +* + IF( A( J+1, K ).NE.ZERO ) THEN + ALPHA = ONE / A( J+1, K ) + CALL CCOPY( M-J-1, WORK( 3 ), 1, A( J+2, K ), 1 ) + CALL CSCAL( M-J-1, ALPHA, A( J+2, K ), 1 ) + ELSE + CALL CLASET( 'Full', M-J-1, 1, ZERO, ZERO, + $ A( J+2, K ), LDA ) + END IF + ELSE + IF( (A( J, K ).EQ.ZERO) .AND. (INFO.EQ.0) ) THEN + INFO = J + END IF + END IF + J = J + 1 + GO TO 30 + 40 CONTINUE + END IF + RETURN +* +* End of CLASYF_AA +* + END diff --git a/SRC/csysv_aa.f b/SRC/csysv_aa.f new file mode 100644 index 00000000..7c82a400 --- /dev/null +++ b/SRC/csysv_aa.f @@ -0,0 +1,254 @@ +*> \brief <b> CSYSV_AA computes the solution to system of linear equations A * X = B for SY matrices</b> +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download CSYSV_AA + dependencies +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/csysv_aa.f"> +*> [TGZ]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/csysv_aa.f"> +*> [ZIP]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/csysv_aa.f"> +*> [TXT]</a> +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE CSYSV_AA( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK, +* LWORK, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER UPLO +* INTEGER N, NRHS, LDA, LDB, LWORK, INFO +* .. +* .. Array Arguments .. +* INTEGER IPIV( * ) +* COMPLEX A( LDA, * ), B( LDB, * ), WORK( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> CSYSV computes the solution to a complex system of linear equations +*> A * X = B, +*> where A is an N-by-N symmetric matrix and X and B are N-by-NRHS +*> matrices. +*> +*> Aasen's algorithm is used to factor A as +*> A = U * T * U**T, 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 +*> form of A is then used to solve the system of equations A * X = B. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] UPLO +*> \verbatim +*> UPLO is CHARACTER*1 +*> = 'U': Upper triangle of A is stored; +*> = 'L': Lower triangle of A is stored. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of linear equations, i.e., the order of the +*> matrix A. N >= 0. +*> \endverbatim +*> +*> \param[in] NRHS +*> \verbatim +*> NRHS is INTEGER +*> The number of right hand sides, i.e., the number of columns +*> of the matrix B. NRHS >= 0. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is REAL array, dimension (LDA,N) +*> On entry, the symmetric matrix A. If UPLO = 'U', the leading +*> N-by-N upper triangular part of A contains the upper +*> triangular part of the matrix A, and the strictly lower +*> triangular part of A is not referenced. If UPLO = 'L', the +*> leading N-by-N lower triangular part of A contains the lower +*> triangular part of the matrix A, and the strictly upper +*> triangular part of A is not referenced. +*> +*> 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 +*> CSYTRF. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N). +*> \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[in,out] B +*> \verbatim +*> B is REAL array, dimension (LDB,NRHS) +*> On entry, the N-by-NRHS right hand side matrix B. +*> On exit, if INFO = 0, the N-by-NRHS solution matrix X. +*> \endverbatim +*> +*> \param[in] LDB +*> \verbatim +*> LDB is INTEGER +*> The leading dimension of the array B. LDB >= max(1,N). +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is REAL array, dimension (MAX(1,LWORK)) +*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. +*> \endverbatim +*> +*> \param[in] LWORK +*> \verbatim +*> LWORK is INTEGER +*> The length of WORK. LWORK >= MAX(2*N, 3*N-2), and for +*> the best performance, LWORK >= max(1,N*NB), where NB is +*> the optimal blocksize for CSYTRF_AA. +*> +*> If LWORK = -1, then a workspace query is assumed; the routine +*> only calculates the optimal size of the WORK array, returns +*> this value as the first entry of the WORK array, and no error +*> message related to LWORK is issued by XERBLA. +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit +*> < 0: if INFO = -i, the i-th argument had an illegal value +*> > 0: if INFO = i, D(i,i) is exactly zero. The factorization +*> has been completed, but the block diagonal matrix D is +*> exactly singular, so the solution could not be computed. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2016 +* +*> \ingroup complexSYsolve +* +* ===================================================================== + SUBROUTINE CSYSV_AA( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK, + $ LWORK, INFO ) +* +* -- LAPACK driver routine (version 3.7.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2016 +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER INFO, LDA, LDB, LWORK, N, NRHS +* .. +* .. Array Arguments .. + INTEGER IPIV( * ) + COMPLEX A( LDA, * ), B( LDB, * ), WORK( * ) +* .. +* +* ===================================================================== +* +* .. Local Scalars .. + LOGICAL LQUERY + INTEGER LWKOPT, LWKOPT_SYTRF, LWKOPT_SYTRS +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER ILAENV + EXTERNAL ILAENV, LSAME +* .. +* .. External Subroutines .. + EXTERNAL XERBLA, CSYTRF, CSYTRS, CSYTRS2 +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX +* .. +* .. Executable Statements .. +* +* Test the input parameters. +* + INFO = 0 + LQUERY = ( LWORK.EQ.-1 ) + IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( NRHS.LT.0 ) THEN + INFO = -3 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -5 + ELSE IF( LDB.LT.MAX( 1, N ) ) THEN + INFO = -8 + ELSE IF( LWORK.LT.MAX(2*N, 3*N-2) .AND. .NOT.LQUERY ) THEN + INFO = -10 + END IF +* + IF( INFO.EQ.0 ) THEN + CALL CSYTRF_AA( UPLO, N, A, LDA, IPIV, WORK, -1, INFO ) + LWKOPT_SYTRF = INT( WORK(1) ) + CALL CSYTRS_AA( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK, + $ -1, INFO ) + LWKOPT_SYTRS = INT( WORK(1) ) + LWKOPT = MAX( LWKOPT_SYTRF, LWKOPT_SYTRS ) + WORK( 1 ) = LWKOPT + IF( LWORK.LT.LWKOPT .AND. .NOT.LQUERY ) THEN + INFO = -10 + END IF + END IF +* + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'CSYSV_AA ', -INFO ) + RETURN + ELSE IF( LQUERY ) THEN + RETURN + END IF +* +* Compute the factorization A = U*T*U**T or A = L*T*L**T. +* + CALL CSYTRF_AA( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO ) + IF( INFO.EQ.0 ) THEN +* +* Solve the system A*X = B, overwriting B with X. +* + CALL CSYTRS_AA( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK, + $ LWORK, INFO ) +* + END IF +* + WORK( 1 ) = LWKOPT +* + RETURN +* +* End of CSYSV_AA +* + END diff --git a/SRC/csytrf_aa.f b/SRC/csytrf_aa.f new file mode 100644 index 00000000..c6b76137 --- /dev/null +++ b/SRC/csytrf_aa.f @@ -0,0 +1,480 @@ +*> \brief \b CSYTRF_AA +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download CSYTRF_AA + dependencies +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/csytrf_aa.f"> +*> [TGZ]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/csytrf_aa.f"> +*> [ZIP]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/csytrf_aa.f"> +*> [TXT]</a> +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE CSYTRF_AA( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER UPLO +* INTEGER N, LDA, LWORK, INFO +* .. +* .. Array Arguments .. +* INTEGER IPIV( * ) +* COMPLEX A( LDA, * ), WORK( * ) +* .. +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> CSYTRF_AA computes the factorization of a complex 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 +*> +*> where U (or L) is a product of permutation and unit upper (lower) +*> triangular matrices, and T is a complex symmetric tridiagonal matrix. +*> +*> This is the blocked version of the algorithm, calling Level 3 BLAS. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] UPLO +*> \verbatim +*> UPLO is CHARACTER*1 +*> = 'U': Upper triangle of A is stored; +*> = 'L': Lower triangle of A is stored. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the matrix A. N >= 0. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is REAL array, dimension (LDA,N) +*> On entry, the symmetric matrix A. If UPLO = 'U', the leading +*> N-by-N upper triangular part of A contains the upper +*> triangular part of the matrix A, and the strictly lower +*> triangular part of A is not referenced. If UPLO = 'L', the +*> leading N-by-N lower triangular part of A contains the lower +*> triangular part of the matrix A, and the strictly upper +*> triangular part of A is not referenced. +*> +*> On exit, the tridiagonal matrix is stored in the diagonals +*> and the subdiagonals of A just below (or above) the diagonals, +*> and L is stored below (or above) the subdiaonals, when UPLO +*> is 'L' (or 'U'). +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N). +*> \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] WORK +*> \verbatim +*> WORK is REAL array, dimension (MAX(1,LWORK)) +*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. +*> \endverbatim +*> +*> \param[in] LWORK +*> \verbatim +*> LWORK is INTEGER +*> The length of WORK. LWORK >= MAX(1,2*N). For optimum performance +*> LWORK >= N*(1+NB), where NB is the optimal blocksize. +*> +*> If LWORK = -1, then a workspace query is assumed; the routine +*> only calculates the optimal size of the WORK array, returns +*> this value as the first entry of the WORK array, and no error +*> message related to LWORK is issued by XERBLA. +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit +*> < 0: if INFO = -i, the i-th argument had an illegal value +*> > 0: if INFO = i, D(i,i) is exactly zero. The factorization +*> has been completed, but the block diagonal matrix D is +*> exactly singular, and division by zero will occur if it +*> is used to solve a system of equations. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2016 +* +*> \ingroup complexSYcomputational +* +* ===================================================================== + SUBROUTINE CSYTRF_AA( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO) +* +* -- LAPACK computational routine (version 3.7.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2016 +* + IMPLICIT NONE +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER N, LDA, LWORK, INFO +* .. +* .. Array Arguments .. + INTEGER IPIV( * ) + COMPLEX A( LDA, * ), WORK( * ) +* .. +* +* ===================================================================== +* .. Parameters .. + COMPLEX ZERO, ONE + PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) +* +* .. Local Scalars .. + LOGICAL LQUERY, UPPER + INTEGER J, LWKOPT, IINFO + INTEGER NB, MJ, NJ, K1, K2, J1, J2, J3, JB + COMPLEX ALPHA +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER ILAENV + EXTERNAL LSAME, ILAENV +* .. +* .. External Subroutines .. + EXTERNAL XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX +* .. +* .. Executable Statements .. +* +* Determine the block size +* + NB = ILAENV( 1, 'CSYTRF', UPLO, N, -1, -1, -1 ) +* +* Test the input parameters. +* + INFO = 0 + UPPER = LSAME( UPLO, 'U' ) + LQUERY = ( LWORK.EQ.-1 ) + IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -4 + ELSE IF( LWORK.LT.MAX( 1, 2*N ) .AND. .NOT.LQUERY ) THEN + INFO = -7 + END IF +* + IF( INFO.EQ.0 ) THEN + LWKOPT = (NB+1)*N + WORK( 1 ) = LWKOPT + END IF +* + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'CSYTRF_AA', -INFO ) + RETURN + ELSE IF( LQUERY ) THEN + RETURN + END IF +* +* Quick return +* + IF ( N.EQ.0 ) THEN + RETURN + ENDIF + IPIV( 1 ) = 1 + IF ( N.EQ.1 ) THEN + IF ( A( 1, 1 ).EQ.ZERO ) THEN + INFO = 1 + END IF + RETURN + END IF +* +* Adjubst block size based on the workspace size +* + IF( LWORK.LT.((1+NB)*N) ) THEN + NB = ( LWORK-N ) / N + END IF +* + IF( UPPER ) THEN +* +* ..................................................... +* Factorize A as L*D*L**T using the upper triangle of A +* ..................................................... +* +* Copy first row A(1, 1:N) into H(1:n) (stored in WORK(1:N)) +* + CALL CCOPY( N, A( 1, 1 ), LDA, WORK( 1 ), 1 ) +* +* J is the main loop index, increasing from 1 to N in steps of +* JB, where JB is the number of columns factorized by CLASYF; +* JB is either NB, or N-J+1 for the last block +* + J = 0 + 10 CONTINUE + IF( J.GE.N ) + $ GO TO 20 +* +* each step of the main loop +* J is the last column of the previous panel +* J1 is the first column of the current panel +* K1 identifies if the previous column of the panel has been +* explicitly stored, e.g., K1=1 for the first panel, and +* K1=0 for the rest +* + J1 = J + 1 + JB = MIN( N-J1+1, NB ) + K1 = MAX(1, J)-J +* +* Panel factorization +* + CALL CLASYF_AA( UPLO, 2-K1, N-J, JB, + $ A( MAX(1, J), J+1 ), LDA, + $ IPIV( J+1 ), WORK, N, WORK( N*NB+1 ), + $ IINFO ) + IF( (IINFO.GT.0) .AND. (INFO.EQ.0) ) THEN + INFO = IINFO+J + ENDIF +* +* Ajust 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 + IF( (J2.NE.IPIV(J2)) .AND. ((J1-K1).GT.2) ) THEN + CALL CSWAP( J1-K1-2, A( 1, J2 ), 1, + $ A( 1, IPIV(J2) ), 1 ) + END IF + END DO + J = J + JB +* +* Trailing submatrix update, where +* the row A(J1-1, J2-1:N) stores U(J1, J2+1:N) and +* WORK stores the current block of the auxiriarly matrix H +* + IF( J.LT.N ) THEN +* +* If first panel and JB=1 (NB=1), then nothing to do +* + IF( J1.GT.1 .OR. JB.GT.1 ) THEN +* +* Merge rank-1 update with BLAS-3 update +* + ALPHA = A( J, J+1 ) + A( J, J+1 ) = ONE + CALL CCOPY( N-J, A( J-1, J+1 ), LDA, + $ WORK( (J+1-J1+1)+JB*N ), 1 ) + CALL CSCAL( N-J, ALPHA, WORK( (J+1-J1+1)+JB*N ), 1 ) +* +* K1 identifies if the previous column of the panel has been +* explicitly stored, e.g., K1=1 and K2= 0 for the first panel, +* while K1=0 and K2=1 for the rest +* + IF( J1.GT.1 ) THEN +* +* Not first panel +* + K2 = 1 + ELSE +* +* First panel +* + K2 = 0 +* +* First update skips the first column +* + JB = JB - 1 + END IF +* + DO J2 = J+1, N, NB + NJ = MIN( NB, N-J2+1 ) +* +* Update (J2, J2) diagonal block with CGEMV +* + J3 = J2 + DO MJ = NJ-1, 1, -1 + CALL CGEMV( 'No transpose', MJ, JB+1, + $ -ONE, WORK( J3-J1+1+K1*N ), N, + $ A( J1-K2, J3 ), 1, + $ ONE, A( J3, J3 ), LDA ) + J3 = J3 + 1 + END DO +* +* Update off-diagonal block of J2-th block row with CGEMM +* + CALL CGEMM( 'Transpose', 'Transpose', + $ NJ, N-J3+1, JB+1, + $ -ONE, A( J1-K2, J2 ), LDA, + $ WORK( J3-J1+1+K1*N ), N, + $ ONE, A( J2, J3 ), LDA ) + END DO +* +* Recover T( J, J+1 ) +* + A( J, J+1 ) = ALPHA + END IF +* +* WORK(J+1, 1) stores H(J+1, 1) +* + CALL CCOPY( N-J, A( J+1, J+1 ), LDA, WORK( 1 ), 1 ) + END IF + GO TO 10 + ELSE +* +* ..................................................... +* Factorize A as L*D*L**T using the lower triangle of A +* ..................................................... +* +* copy first column A(1:N, 1) into H(1:N, 1) +* (stored in WORK(1:N)) +* + CALL CCOPY( N, A( 1, 1 ), 1, WORK( 1 ), 1 ) +* +* J is the main loop index, increasing from 1 to N in steps of +* JB, where JB is the number of columns factorized by CLASYF; +* JB is either NB, or N-J+1 for the last block +* + J = 0 + 11 CONTINUE + IF( J.GE.N ) + $ GO TO 20 +* +* each step of the main loop +* J is the last column of the previous panel +* J1 is the first column of the current panel +* K1 identifies if the previous column of the panel has been +* explicitly stored, e.g., K1=1 for the first panel, and +* K1=0 for the rest +* + J1 = J+1 + JB = MIN( N-J1+1, NB ) + K1 = MAX(1, J)-J +* +* Panel factorization +* + CALL CLASYF_AA( UPLO, 2-K1, N-J, JB, + $ A( J+1, MAX(1, J) ), LDA, + $ IPIV( J+1 ), WORK, N, WORK( N*NB+1 ), IINFO) + IF( (IINFO.GT.0) .AND. (INFO.EQ.0) ) THEN + INFO = IINFO+J + ENDIF +* +* Ajust 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 + IF( (J2.NE.IPIV(J2)) .AND. ((J1-K1).GT.2) ) THEN + CALL CSWAP( J1-K1-2, A( J2, 1 ), LDA, + $ A( IPIV(J2), 1 ), LDA ) + END IF + END DO + J = J + JB +* +* Trailing submatrix update, where +* A(J2+1, J1-1) stores L(J2+1, J1) and +* WORK(J2+1, 1) stores H(J2+1, 1) +* + IF( J.LT.N ) THEN +* +* if first panel and JB=1 (NB=1), then nothing to do +* + IF( J1.GT.1 .OR. JB.GT.1 ) THEN +* +* Merge rank-1 update with BLAS-3 update +* + ALPHA = A( J+1, J ) + A( J+1, J ) = ONE + CALL CCOPY( N-J, A( J+1, J-1 ), 1, + $ WORK( (J+1-J1+1)+JB*N ), 1 ) + CALL CSCAL( N-J, ALPHA, WORK( (J+1-J1+1)+JB*N ), 1 ) +* +* K1 identifies if the previous column of the panel has been +* explicitly stored, e.g., K1=1 and K2= 0 for the first panel, +* while K1=0 and K2=1 for the rest +* + IF( J1.GT.1 ) THEN +* +* Not first panel +* + K2 = 1 + ELSE +* +* First panel +* + K2 = 0 +* +* First update skips the first column +* + JB = JB - 1 + END IF +* + DO J2 = J+1, N, NB + NJ = MIN( NB, N-J2+1 ) +* +* Update (J2, J2) diagonal block with CGEMV +* + J3 = J2 + DO MJ = NJ-1, 1, -1 + CALL CGEMV( 'No transpose', MJ, JB+1, + $ -ONE, WORK( J3-J1+1+K1*N ), N, + $ A( J3, J1-K2 ), LDA, + $ ONE, A( J3, J3 ), 1 ) + J3 = J3 + 1 + END DO +* +* Update off-diagonal block in J2-th block column with CGEMM +* + CALL CGEMM( 'No transpose', 'Transpose', + $ N-J3+1, NJ, JB+1, + $ -ONE, WORK( J3-J1+1+K1*N ), N, + $ A( J2, J1-K2 ), LDA, + $ ONE, A( J3, J2 ), LDA ) + END DO +* +* Recover T( J+1, J ) +* + A( J+1, J ) = ALPHA + END IF +* +* WORK(J+1, 1) stores H(J+1, 1) +* + CALL CCOPY( N-J, A( J+1, J+1 ), 1, WORK( 1 ), 1 ) + END IF + GO TO 11 + END IF +* + 20 CONTINUE + RETURN +* +* End of CSYTRF_AA +* + END diff --git a/SRC/csytrs_aa.f b/SRC/csytrs_aa.f new file mode 100644 index 00000000..fd75ba4c --- /dev/null +++ b/SRC/csytrs_aa.f @@ -0,0 +1,285 @@ +*> \brief \b CSYTRS_AA +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download CSYTRS_AA + dependencies +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/csytrs_aa.f"> +*> [TGZ]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/csytrs_aa.f"> +*> [ZIP]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/csytrs_aa.f"> +*> [TXT]</a> +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE CSYTRS_AA( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, +* WORK, LWORK, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER UPLO +* INTEGER N, NRHS, LDA, LDB, LWORK, INFO +* .. +* .. Array Arguments .. +* INTEGER IPIV( * ) +* COMPLEX A( LDA, * ), B( LDB, * ), WORK( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> CSYTRS_AA solves a system of linear equations A*X = B with a complex +*> symmetric matrix A using the factorization A = U*T*U**T or +*> A = L*T*L**T computed by CSYTRF_AA. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] UPLO +*> \verbatim +*> 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; +*> = 'L': Lower triangular, form is A = L*T*L**T. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the matrix A. N >= 0. +*> \endverbatim +*> +*> \param[in] NRHS +*> \verbatim +*> NRHS is INTEGER +*> The number of right hand sides, i.e., the number of columns +*> of the matrix B. NRHS >= 0. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is REAL array, dimension (LDA,N) +*> Details of factors computed by CSYTRF_AA. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N). +*> \endverbatim +*> +*> \param[in] IPIV +*> \verbatim +*> IPIV is INTEGER array, dimension (N) +*> Details of the interchanges as computed by CSYTRF_AA. +*> \endverbatim +*> +*> \param[in,out] B +*> \verbatim +*> B is REAL array, dimension (LDB,NRHS) +*> On entry, the right hand side matrix B. +*> On exit, the solution matrix X. +*> \endverbatim +*> +*> \param[in] LDB +*> \verbatim +*> LDB is INTEGER +*> The leading dimension of the array B. LDB >= max(1,N). +*> \endverbatim +*> +*> \param[in] WORK +*> \verbatim +*> WORK is DOUBLE array, dimension (MAX(1,LWORK)) +*> \endverbatim +*> +*> \param[in] LWORK +*> \verbatim +*> LWORK is INTEGER, LWORK >= MAX(1,3*N-2). +*> +*> \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 2016 +* +*> \ingroup complexSYcomputational +* +* ===================================================================== + SUBROUTINE CSYTRS_AA( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, + $ WORK, LWORK, INFO ) +* +* -- LAPACK computational routine (version 3.7.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2016 +* + IMPLICIT NONE +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER N, NRHS, LDA, LDB, LWORK, INFO +* .. +* .. Array Arguments .. + INTEGER IPIV( * ) + COMPLEX A( LDA, * ), B( LDB, * ), WORK( * ) +* .. +* +* ===================================================================== +* + COMPLEX ONE + PARAMETER ( ONE = 1.0E+0 ) +* .. +* .. Local Scalars .. + LOGICAL LQUERY, UPPER + INTEGER K, KP, LWKOPT +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. External Subroutines .. + EXTERNAL CGTSV, CSWAP, CTRSM, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX +* .. +* .. Executable Statements .. +* + INFO = 0 + UPPER = LSAME( UPLO, 'U' ) + LQUERY = ( LWORK.EQ.-1 ) + IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( NRHS.LT.0 ) THEN + INFO = -3 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -5 + ELSE IF( LDB.LT.MAX( 1, N ) ) THEN + INFO = -8 + ELSE IF( LWORK.LT.MAX( 1, 3*N-2 ) .AND. .NOT.LQUERY ) THEN + INFO = -10 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'CSYTRS_AA', -INFO ) + RETURN + ELSE IF( LQUERY ) THEN + LWKOPT = (3*N-2) + WORK( 1 ) = LWKOPT + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 .OR. NRHS.EQ.0 ) + $ RETURN +* + IF( UPPER ) THEN +* +* Solve A*X = B, where A = U*T*U**T. +* +* Pivot, P**T * B +* + DO K = 1, N + KP = IPIV( K ) + IF( KP.NE.K ) + $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) + END DO +* +* Compute (U \P**T * B) -> B [ (U \P**T * B) ] +* + CALL CTRSM('L', 'U', 'T', 'U', N-1, NRHS, ONE, A( 1, 2 ), LDA, + $ B( 2, 1 ), LDB) +* +* Compute T \ B -> B [ T \ (U \P**T * B) ] +* + CALL CLACPY( 'F', 1, N, A( 1, 1 ), LDA+1, WORK( N ), 1) + IF( N.GT.1 ) THEN + CALL CLACPY( 'F', 1, N-1, A( 1, 2 ), LDA+1, WORK( 1 ), 1 ) + CALL CLACPY( 'F', 1, N-1, A( 1, 2 ), LDA+1, WORK( 2*N ), 1 ) + END IF + CALL CGTSV( N, NRHS, WORK( 1 ), WORK( N ), WORK( 2*N ), B, LDB, + $ INFO ) +* +* Compute (U**T \ B) -> B [ U**T \ (T \ (U \P**T * B) ) ] +* + CALL CTRSM( 'L', 'U', 'N', 'U', N-1, NRHS, ONE, A( 1, 2 ), LDA, + $ B( 2, 1 ), LDB) +* +* Pivot, P * B [ P * (U**T \ (T \ (U \P**T * B) )) ] +* + DO K = N, 1, -1 + KP = IPIV( K ) + IF( KP.NE.K ) + $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) + END DO +* + ELSE +* +* Solve A*X = B, where A = L*T*L**T. +* +* Pivot, P**T * B +* + DO K = 1, N + KP = IPIV( K ) + IF( KP.NE.K ) + $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) + END DO +* +* Compute (L \P**T * B) -> B [ (L \P**T * B) ] +* + CALL CTRSM( 'L', 'L', 'N', 'U', N-1, NRHS, ONE, A( 2, 1 ), LDA, + $ B( 2, 1 ), LDB) +* +* Compute T \ B -> B [ T \ (L \P**T * B) ] +* + CALL CLACPY( 'F', 1, N, A(1, 1), LDA+1, WORK(N), 1) + IF( N.GT.1 ) THEN + CALL CLACPY( 'F', 1, N-1, A( 2, 1 ), LDA+1, WORK( 1 ), 1 ) + CALL CLACPY( 'F', 1, N-1, A( 2, 1 ), LDA+1, WORK( 2*N ), 1 ) + END IF + CALL CGTSV( N, NRHS, WORK( 1 ), WORK(N), WORK( 2*N ), B, LDB, + $ INFO) +* +* Compute (L**T \ B) -> B [ L**T \ (T \ (L \P**T * B) ) ] +* + CALL CTRSM( '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) )) ] +* + DO K = N, 1, -1 + KP = IPIV( K ) + IF( KP.NE.K ) + $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) + END DO +* + END IF +* + RETURN +* +* End of CSYTRS_AA +* + END diff --git a/SRC/dlasyf_aa.f b/SRC/dlasyf_aa.f index cc0b80f2..393b50a8 100644 --- a/SRC/dlasyf_aa.f +++ b/SRC/dlasyf_aa.f @@ -19,7 +19,7 @@ * =========== * * SUBROUTINE DLASYF_AA( UPLO, J1, M, NB, A, LDA, IPIV, -* H, LDH, WORK, INFO ) +* H, LDH, WORK, INFO ) * * .. Scalar Arguments .. * CHARACTER UPLO @@ -46,7 +46,7 @@ *> which is used to factorize the first panel. *> *> The resulting J-th row of U, or J-th column of L, is stored in the -*> (J-1)-th row, or column, of A (without the unit diatonals), while +*> (J-1)-th row, or column, of A (without the unit diagonals), while *> the diagonal and subdiagonal of A are overwritten by those of T. *> *> \endverbatim @@ -152,7 +152,7 @@ * * ===================================================================== SUBROUTINE DLASYF_AA( UPLO, J1, M, NB, A, LDA, IPIV, - $ H, LDH, WORK, INFO ) + $ H, LDH, WORK, INFO ) * * -- LAPACK computational routine (version 3.7.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- diff --git a/SRC/dsysv_aa.f b/SRC/dsysv_aa.f index 9f9969fb..055097fb 100644 --- a/SRC/dsysv_aa.f +++ b/SRC/dsysv_aa.f @@ -19,7 +19,7 @@ * =========== * * SUBROUTINE DSYSV_AA( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK, -* LWORK, INFO ) +* LWORK, INFO ) * * .. Scalar Arguments .. * CHARACTER UPLO @@ -126,8 +126,8 @@ *> \param[in] LWORK *> \verbatim *> LWORK is INTEGER -*> The length of WORK. LWORK >= MAX(2*N, 3*N-2), and for -*> the best performance, LWORK >= max(1,N*NB), where NB is +*> The length of WORK. LWORK >= MAX(1,2*N,3*N-2), and for +*> the best performance, LWORK >= MAX(1,N*NB), where NB is *> the optimal blocksize for DSYTRF_AA. *> *> If LWORK = -1, then a workspace query is assumed; the routine @@ -156,11 +156,13 @@ * *> \date November 2016 * +* @precisions fortran d -> z c +* *> \ingroup doubleSYsolve * * ===================================================================== SUBROUTINE DSYSV_AA( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK, - $ LWORK, INFO ) + $ LWORK, INFO ) * * -- LAPACK driver routine (version 3.7.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- diff --git a/SRC/dsytrf_aa.f b/SRC/dsytrf_aa.f index 07919a2c..0e038806 100644 --- a/SRC/dsytrf_aa.f +++ b/SRC/dsytrf_aa.f @@ -101,7 +101,7 @@ *> \param[in] LWORK *> \verbatim *> LWORK is INTEGER -*> The length of WORK. LWORK >=2*N. For optimum performance +*> The length of WORK. LWORK >= MAX(1,2*N). For optimum performance *> LWORK >= N*(1+NB), where NB is the optimal blocksize. *> *> If LWORK = -1, then a workspace query is assumed; the routine @@ -191,7 +191,7 @@ INFO = -2 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN INFO = -4 - ELSE IF( LWORK.LT.( 2*N ) .AND. .NOT.LQUERY ) THEN + ELSE IF( LWORK.LT.MAX( 1, 2*N ) .AND. .NOT.LQUERY ) THEN INFO = -7 END IF * diff --git a/SRC/dsytrs_aa.f b/SRC/dsytrs_aa.f index ddb9d3fc..6c56b919 100644 --- a/SRC/dsytrs_aa.f +++ b/SRC/dsytrs_aa.f @@ -104,7 +104,7 @@ *> *> \param[in] LWORK *> \verbatim -*> LWORK is INTEGER, LWORK >= 3*N-2. +*> LWORK is INTEGER, LWORK >= MAX(1,3*N-2). *> *> \param[out] INFO *> \verbatim @@ -179,7 +179,7 @@ INFO = -5 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN INFO = -8 - ELSE IF( LWORK.LT.(3*N-2) .AND. .NOT.LQUERY ) THEN + ELSE IF( LWORK.LT.MAX( 1, 3*N-2 ) .AND. .NOT.LQUERY ) THEN INFO = -10 END IF IF( INFO.NE.0 ) THEN diff --git a/SRC/slasyf_aa.f b/SRC/slasyf_aa.f index 953d574d..3869f2a6 100644 --- a/SRC/slasyf_aa.f +++ b/SRC/slasyf_aa.f @@ -19,7 +19,7 @@ * =========== * * SUBROUTINE SLASYF_AA( UPLO, J1, M, NB, A, LDA, IPIV, -* H, LDH, WORK, INFO ) +* H, LDH, WORK, INFO ) * * .. Scalar Arguments .. * CHARACTER UPLO @@ -27,7 +27,7 @@ * .. * .. Array Arguments .. * INTEGER IPIV( * ) -* REAL A( LDA, * ), H( LDH, * ), WORK( * ) +* REAL A( LDA, * ), H( LDH, * ), WORK( * ) * .. * * @@ -46,7 +46,7 @@ *> which is used to factorize the first panel. *> *> The resulting J-th row of U, or J-th column of L, is stored in the -*> (J-1)-th row, or column, of A (without the unit diatonals), while +*> (J-1)-th row, or column, of A (without the unit diagonals), while *> the diagonal and subdiagonal of A are overwritten by those of T. *> *> \endverbatim @@ -152,7 +152,7 @@ * * ===================================================================== SUBROUTINE SLASYF_AA( UPLO, J1, M, NB, A, LDA, IPIV, - $ H, LDH, WORK, INFO ) + $ H, LDH, WORK, INFO ) * * -- LAPACK computational routine (version 3.7.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -167,7 +167,7 @@ * .. * .. Array Arguments .. INTEGER IPIV( * ) - REAL A( LDA, * ), H( LDH, * ), WORK( * ) + REAL A( LDA, * ), H( LDH, * ), WORK( * ) * .. * * ===================================================================== @@ -177,7 +177,7 @@ * * .. Local Scalars .. INTEGER J, K, K1, I1, I2 - REAL PIV, ALPHA + REAL PIV, ALPHA * .. * .. External Functions .. LOGICAL LSAME diff --git a/SRC/ssysv_aa.f b/SRC/ssysv_aa.f index 52f6eb5f..d61a346c 100644 --- a/SRC/ssysv_aa.f +++ b/SRC/ssysv_aa.f @@ -19,7 +19,7 @@ * =========== * * SUBROUTINE SSYSV_AA( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK, -* LWORK, INFO ) +* LWORK, INFO ) * * .. Scalar Arguments .. * CHARACTER UPLO @@ -27,7 +27,7 @@ * .. * .. Array Arguments .. * INTEGER IPIV( * ) -* REAL A( LDA, * ), B( LDB, * ), WORK( * ) +* REAL A( LDA, * ), B( LDB, * ), WORK( * ) * .. * * @@ -126,8 +126,8 @@ *> \param[in] LWORK *> \verbatim *> LWORK is INTEGER -*> The length of WORK. LWORK >= MAX(2*N, 3*N-2), and for -*> the best performance, LWORK >= max(1,N*NB), where NB is +*> The length of WORK. LWORK >= MAX(1,2*N,3*N-2), and for +*> the best performance, LWORK >= MAX(1,N*NB), where NB is *> the optimal blocksize for SSYTRF_AA. *> *> If LWORK = -1, then a workspace query is assumed; the routine @@ -160,7 +160,7 @@ * * ===================================================================== SUBROUTINE SSYSV_AA( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK, - $ LWORK, INFO ) + $ LWORK, INFO ) * * -- LAPACK driver routine (version 3.7.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -173,7 +173,7 @@ * .. * .. Array Arguments .. INTEGER IPIV( * ) - REAL A( LDA, * ), B( LDB, * ), WORK( * ) + REAL A( LDA, * ), B( LDB, * ), WORK( * ) * .. * * ===================================================================== diff --git a/SRC/ssytrf_aa.f b/SRC/ssytrf_aa.f index 13498c9b..a22ff05d 100644 --- a/SRC/ssytrf_aa.f +++ b/SRC/ssytrf_aa.f @@ -101,7 +101,7 @@ *> \param[in] LWORK *> \verbatim *> LWORK is INTEGER -*> The length of WORK. LWORK >=2*N. For optimum performance +*> The length of WORK. LWORK >= MAX(1,2*N). For optimum performance *> LWORK >= N*(1+NB), where NB is the optimal blocksize. *> *> If LWORK = -1, then a workspace query is assumed; the routine @@ -191,7 +191,7 @@ INFO = -2 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN INFO = -4 - ELSE IF( LWORK.LT.( 2*N ) .AND. .NOT.LQUERY ) THEN + ELSE IF( LWORK.LT.MAX( 1, 2*N ) .AND. .NOT.LQUERY ) THEN INFO = -7 END IF * diff --git a/SRC/ssytrs_aa.f b/SRC/ssytrs_aa.f index 06c793ae..911016e0 100644 --- a/SRC/ssytrs_aa.f +++ b/SRC/ssytrs_aa.f @@ -104,7 +104,7 @@ *> *> \param[in] LWORK *> \verbatim -*> LWORK is INTEGER, LWORK >= 3*N-2. +*> LWORK is INTEGER, LWORK >= MAX(1,3*N-2). *> *> \param[out] INFO *> \verbatim @@ -179,7 +179,7 @@ INFO = -5 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN INFO = -8 - ELSE IF( LWORK.LT.(3*N-2) .AND. .NOT.LQUERY ) THEN + ELSE IF( LWORK.LT.MAX( 1, 3*N-2 ) .AND. .NOT.LQUERY ) THEN INFO = -10 END IF IF( INFO.NE.0 ) THEN diff --git a/SRC/zhesv_aa.f b/SRC/zhesv_aa.f index 50314547..dffa4754 100644 --- a/SRC/zhesv_aa.f +++ b/SRC/zhesv_aa.f @@ -19,7 +19,7 @@ * =========== * * SUBROUTINE ZHESV_AA( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK, -* LWORK, INFO ) +* LWORK, INFO ) * * .. Scalar Arguments .. * CHARACTER UPLO @@ -126,9 +126,9 @@ *> \param[in] LWORK *> \verbatim *> LWORK is INTEGER -*> The length of WORK. LWORK >= 1, and for best performance -*> LWORK >= max(1,N*NB), where NB is the optimal blocksize for -*> ZHETRF. +*> The length of WORK. LWORK >= MAX(1,2*N,3*N-2), and for best +*> performance LWORK >= max(1,N*NB), where NB is the optimal +*> blocksize for ZHETRF. *> for LWORK < N, TRS will be done with Level BLAS 2 *> for LWORK >= N, TRS will be done with Level BLAS 3 *> @@ -162,7 +162,7 @@ * * ===================================================================== SUBROUTINE ZHESV_AA( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK, - $ LWORK, INFO ) + $ LWORK, INFO ) * * -- LAPACK driver routine (version 3.7.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- diff --git a/SRC/zhetrf_aa.f b/SRC/zhetrf_aa.f index 73a8383a..a3dd0950 100644 --- a/SRC/zhetrf_aa.f +++ b/SRC/zhetrf_aa.f @@ -37,7 +37,7 @@ *> ZHETRF_AA computes the factorization of a complex hermitian 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*U**H or A = L*T*L**H *> *> where U (or L) is a product of permutation and unit upper (lower) *> triangular matrices, and T is a hermitian tridiagonal matrix. @@ -101,7 +101,7 @@ *> \param[in] LWORK *> \verbatim *> LWORK is INTEGER -*> The length of WORK. LWORK >= 2*N. For optimum performance +*> The length of WORK. LWORK >= MAX(1,2*N). For optimum performance *> LWORK >= N*(1+NB), where NB is the optimal blocksize. *> *> If LWORK = -1, then a workspace query is assumed; the routine @@ -191,7 +191,7 @@ INFO = -2 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN INFO = -4 - ELSE IF( LWORK.LT.( 2*N ) .AND. .NOT.LQUERY ) THEN + ELSE IF( LWORK.LT.MAX( 1, 2*N ) .AND. .NOT.LQUERY ) THEN INFO = -7 END IF * @@ -230,7 +230,7 @@ IF( UPPER ) THEN * * ..................................................... -* Factorize A as L*D*L**T using the upper triangle of A +* Factorize A as L*D*L**H using the upper triangle of A * ..................................................... * * copy first row A(1, 1:N) into H(1:n) (stored in WORK(1:N)) @@ -353,7 +353,7 @@ 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 * ..................................................... * * copy first column A(1:N, 1) into H(1:N, 1) diff --git a/SRC/zhetrs_aa.f b/SRC/zhetrs_aa.f index fd819d5a..0a02b8a6 100644 --- a/SRC/zhetrs_aa.f +++ b/SRC/zhetrs_aa.f @@ -27,7 +27,7 @@ * .. * .. Array Arguments .. * INTEGER IPIV( * ) -* COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) +* COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) * .. * * @@ -37,8 +37,8 @@ *> *> \verbatim *> -*> ZHETRS_AA solves a system of linear equations A*X = B with a real -*> hermitian matrix A using the factorization A = U*T*U**T or +*> 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 *> A = L*T*L**T computed by ZHETRF_AA. *> \endverbatim * @@ -50,8 +50,8 @@ *> 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; -*> = 'L': Lower triangular, form is A = L*T*L**T. +*> = 'U': Upper triangular, form is A = U*T*U**H; +*> = 'L': Lower triangular, form is A = L*T*L**H. *> \endverbatim *> *> \param[in] N @@ -105,7 +105,7 @@ *> *> \param[in] LWORK *> \verbatim -*> LWORK is INTEGER, LWORK >= 3*N-2. +*> LWORK is INTEGER, LWORK >= MAX(1,3*N-2). *> *> \param[out] INFO *> \verbatim @@ -143,12 +143,12 @@ * .. * .. Array Arguments .. INTEGER IPIV( * ) - COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) + COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) * .. * * ===================================================================== * - COMPLEX*16 ONE + COMPLEX*16 ONE PARAMETER ( ONE = 1.0D+0 ) * .. * .. Local Scalars .. @@ -180,7 +180,7 @@ INFO = -5 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN INFO = -8 - ELSE IF( LWORK.LT.(3*N-2) .AND. .NOT.LQUERY ) THEN + ELSE IF( LWORK.LT.MAX( 1, 3*N-2 ) .AND. .NOT.LQUERY ) THEN INFO = -10 END IF IF( INFO.NE.0 ) THEN diff --git a/SRC/zlahef_aa.f b/SRC/zlahef_aa.f index 45d1b67c..ef42f752 100644 --- a/SRC/zlahef_aa.f +++ b/SRC/zlahef_aa.f @@ -19,7 +19,7 @@ * =========== * * SUBROUTINE ZLAHEF_AA( UPLO, J1, M, NB, A, LDA, IPIV, -* H, LDH, WORK, INFO ) +* H, LDH, WORK, INFO ) * * .. Scalar Arguments .. * CHARACTER UPLO @@ -36,7 +36,7 @@ *> *> \verbatim *> -*> DLATRF_AA factorizes a panel of a real hermitian matrix A using +*> DLAHEF_AA factorizes a panel of a complex hermitian matrix A using *> the Aasen's algorithm. The panel consists of a set of NB rows of A *> when UPLO is U, or a set of NB columns when UPLO is L. *> @@ -46,7 +46,7 @@ *> which is used to factorize the first panel. *> *> The resulting J-th row of U, or J-th column of L, is stored in the -*> (J-1)-th row, or column, of A (without the unit diatonals), while +*> (J-1)-th row, or column, of A (without the unit diagonals), while *> the diagonal and subdiagonal of A are overwritten by those of T. *> *> \endverbatim @@ -152,7 +152,7 @@ * * ===================================================================== SUBROUTINE ZLAHEF_AA( UPLO, J1, M, NB, A, LDA, IPIV, - $ H, LDH, WORK, INFO ) + $ H, LDH, WORK, INFO ) * * -- LAPACK computational routine (version 3.7.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- diff --git a/SRC/zlasyf_aa.f b/SRC/zlasyf_aa.f new file mode 100644 index 00000000..fb914662 --- /dev/null +++ b/SRC/zlasyf_aa.f @@ -0,0 +1,506 @@ +*> \brief \b ZLASYF_AA +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download ZLASYF_AA + dependencies +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlasyf_aa.f"> +*> [TGZ]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlasyf_aa.f"> +*> [ZIP]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlasyf_aa.f"> +*> [TXT]</a> +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE ZLASYF_AA( UPLO, J1, M, NB, A, LDA, IPIV, +* H, LDH, WORK, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER UPLO +* INTEGER J1, M, NB, LDA, LDH, INFO +* .. +* .. Array Arguments .. +* INTEGER IPIV( * ) +* COMPLEX*16 A( LDA, * ), H( LDH, * ), WORK( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DLATRF_AA factorizes a panel of a complex symmetric matrix A using +*> the Aasen's algorithm. The panel consists of a set of NB rows of A +*> when UPLO is U, or a set of NB columns when UPLO is L. +*> +*> In order to factorize the panel, the Aasen's algorithm requires the +*> last row, or column, of the previous panel. The first row, or column, +*> of A is set to be the first row, or column, of an identity matrix, +*> which is used to factorize the first panel. +*> +*> The resulting J-th row of U, or J-th column of L, is stored in the +*> (J-1)-th row, or column, of A (without the unit diagonals), while +*> the diagonal and subdiagonal of A are overwritten by those of T. +*> +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] UPLO +*> \verbatim +*> UPLO is CHARACTER*1 +*> = 'U': Upper triangle of A is stored; +*> = 'L': Lower triangle of A is stored. +*> \endverbatim +*> +*> \param[in] J1 +*> \verbatim +*> J1 is INTEGER +*> The location of the first row, or column, of the panel +*> within the submatrix of A, passed to this routine, e.g., +*> when called by ZSYTRF_AA, for the first panel, J1 is 1, +*> while for the remaining panels, J1 is 2. +*> \endverbatim +*> +*> \param[in] M +*> \verbatim +*> M is INTEGER +*> The dimension of the submatrix. M >= 0. +*> \endverbatim +*> +*> \param[in] NB +*> \verbatim +*> NB is INTEGER +*> The dimension of the panel to be facotorized. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is COMPLEX*16 array, dimension (LDA,M) for +*> the first panel, while dimension (LDA,M+1) for the +*> remaining panels. +*> +*> On entry, A contains the last row, or column, of +*> the previous panel, and the trailing submatrix of A +*> to be factorized, except for the first panel, only +*> the panel is passed. +*> +*> On exit, the leading panel is factorized. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N). +*> \endverbatim +*> +*> \param[out] IPIV +*> \verbatim +*> IPIV is INTEGER array, dimension (N) +*> Details of the row and column interchanges, +*> the row and column k were interchanged with the row and +*> column IPIV(k). +*> \endverbatim +*> +*> \param[in,out] H +*> \verbatim +*> H is COMPLEX*16 workspace, dimension (LDH,NB). +*> +*> \endverbatim +*> +*> \param[in] LDH +*> \verbatim +*> LDH is INTEGER +*> The leading dimension of the workspace H. LDH >= max(1,M). +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is COMPLEX*16 workspace, dimension (M). +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit +*> < 0: if INFO = -i, the i-th argument had an illegal value +*> > 0: if INFO = i, D(i,i) is exactly zero. The factorization +*> has been completed, but the block diagonal matrix D is +*> exactly singular, and division by zero will occur if it +*> is used to solve a system of equations. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2016 +* +*> \ingroup complex16SYcomputational +* +* ===================================================================== + SUBROUTINE ZLASYF_AA( UPLO, J1, M, NB, A, LDA, IPIV, + $ H, LDH, WORK, INFO ) +* +* -- LAPACK computational routine (version 3.7.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2016 +* + IMPLICIT NONE +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER M, NB, J1, LDA, LDH, INFO +* .. +* .. Array Arguments .. + INTEGER IPIV( * ) + COMPLEX*16 A( LDA, * ), H( LDH, * ), WORK( * ) +* .. +* +* ===================================================================== +* .. Parameters .. + COMPLEX*16 ZERO, ONE + PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) +* +* .. Local Scalars .. + INTEGER J, K, K1, I1, I2 + COMPLEX*16 PIV, ALPHA +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER IZAMAX, ILAENV + EXTERNAL LSAME, ILAENV, IZAMAX +* .. +* .. External Subroutines .. + EXTERNAL XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX +* .. +* .. Executable Statements .. +* + INFO = 0 + J = 1 +* +* K1 is the first column of the panel to be factorized +* i.e., K1 is 2 for the first block column, and 1 for the rest of the blocks +* + K1 = (2-J1)+1 +* + IF( LSAME( UPLO, 'U' ) ) THEN +* +* ..................................................... +* Factorize A as U**T*D*U using the upper triangle of A +* ..................................................... +* + 10 CONTINUE + IF ( J.GT.MIN(M, NB) ) + $ GO TO 20 +* +* K is the column to be factorized +* when being called from ZSYTRF_AA, +* > for the first block column, J1 is 1, hence J1+J-1 is J, +* > for the rest of the columns, J1 is 2, and J1+J-1 is J+1, +* + K = J1+J-1 +* +* H(J:N, J) := A(J, J:N) - H(J:N, 1:(J-1)) * L(J1:(J-1), J), +* where H(J:N, J) has been initialized to be A(J, J:N) +* + IF( K.GT.2 ) THEN +* +* K is the column to be factorized +* > for the first block column, K is J, skipping the first two +* columns +* > for the rest of the columns, K is J+1, skipping only the +* first column +* + CALL ZGEMV( 'No transpose', M-J+1, J-K1, + $ -ONE, H( J, K1 ), LDH, + $ A( 1, J ), 1, + $ ONE, H( J, J ), 1 ) + END IF +* +* Copy H(i:n, i) into WORK +* + CALL ZCOPY( M-J+1, H( J, J ), 1, WORK( 1 ), 1 ) +* + IF( J.GT.K1 ) THEN +* +* Compute WORK := WORK - L(J-1, J:N) * T(J-1,J), +* where A(J-1, J) stores T(J-1, J) and A(J-2, J:N) stores U(J-1, J:N) +* + ALPHA = -A( K-1, J ) + CALL ZAXPY( M-J+1, ALPHA, A( K-2, J ), LDA, WORK( 1 ), 1 ) + END IF +* +* Set A(J, J) = T(J, J) +* + A( K, J ) = WORK( 1 ) +* + IF( J.LT.M ) THEN +* +* Compute WORK(2:N) = T(J, J) L(J, (J+1):N) +* where A(J, J) stores T(J, J) and A(J-1, (J+1):N) stores U(J, (J+1):N) +* + IF( K.GT.1 ) THEN + ALPHA = -A( K, J ) + CALL ZAXPY( M-J, ALPHA, A( K-1, J+1 ), LDA, + $ WORK( 2 ), 1 ) + ENDIF +* +* Find max(|WORK(2:n)|) +* + I2 = IZAMAX( M-J, WORK( 2 ), 1 ) + 1 + PIV = WORK( I2 ) +* +* Apply symmetric pivot +* + IF( (I2.NE.2) .AND. (PIV.NE.0) ) THEN +* +* Swap WORK(I1) and WORK(I2) +* + I1 = 2 + WORK( I2 ) = WORK( I1 ) + WORK( I1 ) = PIV +* +* Swap A(I1, I1+1:N) with A(I1+1:N, I2) +* + I1 = I1+J-1 + I2 = I2+J-1 + CALL ZSWAP( I2-I1-1, A( J1+I1-1, I1+1 ), LDA, + $ A( J1+I1, I2 ), 1 ) +* +* Swap A(I1, I2+1:N) with A(I2, I2+1:N) +* + CALL ZSWAP( M-I2, A( J1+I1-1, I2+1 ), LDA, + $ A( J1+I2-1, I2+1 ), LDA ) +* +* Swap A(I1, I1) with A(I2,I2) +* + PIV = A( I1+J1-1, I1 ) + A( J1+I1-1, I1 ) = A( J1+I2-1, I2 ) + A( J1+I2-1, I2 ) = PIV +* +* Swap H(I1, 1:J1) with H(I2, 1:J1) +* + CALL ZSWAP( I1-1, H( I1, 1 ), LDH, H( I2, 1 ), LDH ) + IPIV( I1 ) = I2 +* + IF( I1.GT.(K1-1) ) THEN +* +* Swap L(1:I1-1, I1) with L(1:I1-1, I2), +* skipping the first column +* + CALL ZSWAP( I1-K1+1, A( 1, I1 ), 1, + $ A( 1, I2 ), 1 ) + END IF + ELSE + IPIV( J+1 ) = J+1 + ENDIF +* +* Set A(J, J+1) = T(J, J+1) +* + A( K, J+1 ) = WORK( 2 ) + IF( (A( K, J ).EQ.ZERO ) .AND. + $ ( (J.EQ.M) .OR. (A( K, J+1 ).EQ.ZERO))) THEN + IF(INFO .EQ. 0) THEN + INFO = J + ENDIF + END IF +* + IF( J.LT.NB ) THEN +* +* Copy A(J+1:N, J+1) into H(J:N, J), +* + CALL ZCOPY( M-J, A( K+1, J+1 ), LDA, + $ H( J+1, J+1 ), 1 ) + END IF +* +* 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) +* + IF( A( K, J+1 ).NE.ZERO ) THEN + ALPHA = ONE / A( K, J+1 ) + CALL ZCOPY( M-J-1, WORK( 3 ), 1, A( K, J+2 ), LDA ) + CALL ZSCAL( M-J-1, ALPHA, A( K, J+2 ), LDA ) + ELSE + CALL ZLASET( 'Full', 1, M-J-1, ZERO, ZERO, + $ A( K, J+2 ), LDA) + END IF + ELSE + IF( (A( K, J ).EQ.ZERO) .AND. (INFO.EQ.0) ) THEN + INFO = J + END IF + END IF + J = J + 1 + GO TO 10 + 20 CONTINUE +* + ELSE +* +* ..................................................... +* Factorize A as L*D*L**T using the lower triangle of A +* ..................................................... +* + 30 CONTINUE + IF( J.GT.MIN( M, NB ) ) + $ GO TO 40 +* +* K is the column to be factorized +* when being called from ZSYTRF_AA, +* > for the first block column, J1 is 1, hence J1+J-1 is J, +* > for the rest of the columns, J1 is 2, and J1+J-1 is J+1, +* + K = J1+J-1 +* +* H(J:N, J) := A(J:N, J) - H(J:N, 1:(J-1)) * L(J, J1:(J-1))^T, +* where H(J:N, J) has been initialized to be A(J:N, J) +* + IF( K.GT.2 ) THEN +* +* K is the column to be factorized +* > for the first block column, K is J, skipping the first two +* columns +* > for the rest of the columns, K is J+1, skipping only the +* first column +* + CALL ZGEMV( 'No transpose', M-J+1, J-K1, + $ -ONE, H( J, K1 ), LDH, + $ A( J, 1 ), LDA, + $ ONE, H( J, J ), 1 ) + END IF +* +* Copy H(J:N, J) into WORK +* + CALL ZCOPY( M-J+1, H( J, J ), 1, WORK( 1 ), 1 ) +* + IF( J.GT.K1 ) THEN +* +* Compute WORK := WORK - L(J:N, J-1) * T(J-1,J), +* where A(J-1, J) = T(J-1, J) and A(J, J-2) = L(J, J-1) +* + ALPHA = -A( J, K-1 ) + CALL ZAXPY( M-J+1, ALPHA, A( J, K-2 ), 1, WORK( 1 ), 1 ) + END IF +* +* Set A(J, J) = T(J, J) +* + A( J, K ) = WORK( 1 ) +* + IF( J.LT.M ) THEN +* +* Compute WORK(2:N) = T(J, J) L((J+1):N, J) +* where A(J, J) = T(J, J) and A((J+1):N, J-1) = L((J+1):N, J) +* + IF( K.GT.1 ) THEN + ALPHA = -A( J, K ) + CALL ZAXPY( M-J, ALPHA, A( J+1, K-1 ), 1, + $ WORK( 2 ), 1 ) + ENDIF +* +* Find max(|WORK(2:n)|) +* + I2 = IZAMAX( M-J, WORK( 2 ), 1 ) + 1 + PIV = WORK( I2 ) +* +* Apply symmetric pivot +* + IF( (I2.NE.2) .AND. (PIV.NE.0) ) THEN +* +* Swap WORK(I1) and WORK(I2) +* + I1 = 2 + WORK( I2 ) = WORK( I1 ) + WORK( I1 ) = PIV +* +* Swap A(I1+1:N, I1) with A(I2, I1+1:N) +* + I1 = I1+J-1 + I2 = I2+J-1 + CALL ZSWAP( I2-I1-1, A( I1+1, J1+I1-1 ), 1, + $ A( I2, J1+I1 ), LDA ) +* +* Swap A(I2+1:N, I1) with A(I2+1:N, I2) +* + CALL ZSWAP( M-I2, A( I2+1, J1+I1-1 ), 1, + $ A( I2+1, J1+I2-1 ), 1 ) +* +* Swap A(I1, I1) with A(I2, I2) +* + PIV = A( I1, J1+I1-1 ) + A( I1, J1+I1-1 ) = A( I2, J1+I2-1 ) + A( I2, J1+I2-1 ) = PIV +* +* Swap H(I1, I1:J1) with H(I2, I2:J1) +* + CALL ZSWAP( I1-1, H( I1, 1 ), LDH, H( I2, 1 ), LDH ) + IPIV( I1 ) = I2 +* + IF( I1.GT.(K1-1) ) THEN +* +* Swap L(1:I1-1, I1) with L(1:I1-1, I2), +* skipping the first column +* + CALL ZSWAP( I1-K1+1, A( I1, 1 ), LDA, + $ A( I2, 1 ), LDA ) + END IF + ELSE + IPIV( J+1 ) = J+1 + ENDIF +* +* Set A(J+1, J) = T(J+1, J) +* + A( J+1, K ) = WORK( 2 ) + IF( (A( J, K ).EQ.ZERO) .AND. + $ ( (J.EQ.M) .OR. (A( J+1, K ).EQ.ZERO)) ) THEN + IF (INFO .EQ. 0) + $ INFO = J + END IF +* + IF( J.LT.NB ) THEN +* +* Copy A(J+1:N, J+1) into H(J+1:N, J), +* + CALL ZCOPY( M-J, A( J+1, K+1 ), 1, + $ H( J+1, J+1 ), 1 ) + END IF +* +* 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) +* + IF( A( J+1, K ).NE.ZERO ) THEN + ALPHA = ONE / A( J+1, K ) + CALL ZCOPY( M-J-1, WORK( 3 ), 1, A( J+2, K ), 1 ) + CALL ZSCAL( M-J-1, ALPHA, A( J+2, K ), 1 ) + ELSE + CALL ZLASET( 'Full', M-J-1, 1, ZERO, ZERO, + $ A( J+2, K ), LDA ) + END IF + ELSE + IF( (A( J, K ).EQ.ZERO) .AND. (INFO.EQ.0) ) THEN + INFO = J + END IF + END IF + J = J + 1 + GO TO 30 + 40 CONTINUE + END IF + RETURN +* +* End of ZLASYF_AA +* + END diff --git a/SRC/zsysv_aa.f b/SRC/zsysv_aa.f new file mode 100644 index 00000000..6c767148 --- /dev/null +++ b/SRC/zsysv_aa.f @@ -0,0 +1,254 @@ +*> \brief <b> ZSYSV_AA computes the solution to system of linear equations A * X = B for SY matrices</b> +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download ZSYSV_AA + dependencies +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zsysv_aa.f"> +*> [TGZ]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zsysv_aa.f"> +*> [ZIP]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zsysv_aa.f"> +*> [TXT]</a> +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE ZSYSV_AA( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK, +* LWORK, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER UPLO +* INTEGER N, NRHS, LDA, LDB, LWORK, INFO +* .. +* .. Array Arguments .. +* INTEGER IPIV( * ) +* COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> ZSYSV computes the solution to a complex system of linear equations +*> A * X = B, +*> where A is an N-by-N symmetric matrix and X and B are N-by-NRHS +*> matrices. +*> +*> Aasen's algorithm is used to factor A as +*> A = U * T * U**T, 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 +*> form of A is then used to solve the system of equations A * X = B. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] UPLO +*> \verbatim +*> UPLO is CHARACTER*1 +*> = 'U': Upper triangle of A is stored; +*> = 'L': Lower triangle of A is stored. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of linear equations, i.e., the order of the +*> matrix A. N >= 0. +*> \endverbatim +*> +*> \param[in] NRHS +*> \verbatim +*> NRHS is INTEGER +*> The number of right hand sides, i.e., the number of columns +*> of the matrix B. NRHS >= 0. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is COMPLEX*16 array, dimension (LDA,N) +*> On entry, the symmetric matrix A. If UPLO = 'U', the leading +*> N-by-N upper triangular part of A contains the upper +*> triangular part of the matrix A, and the strictly lower +*> triangular part of A is not referenced. If UPLO = 'L', the +*> leading N-by-N lower triangular part of A contains the lower +*> triangular part of the matrix A, and the strictly upper +*> triangular part of A is not referenced. +*> +*> 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 +*> ZSYTRF. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N). +*> \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[in,out] B +*> \verbatim +*> B is COMPLEX*16 array, dimension (LDB,NRHS) +*> On entry, the N-by-NRHS right hand side matrix B. +*> On exit, if INFO = 0, the N-by-NRHS solution matrix X. +*> \endverbatim +*> +*> \param[in] LDB +*> \verbatim +*> LDB is INTEGER +*> The leading dimension of the array B. LDB >= max(1,N). +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK)) +*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. +*> \endverbatim +*> +*> \param[in] LWORK +*> \verbatim +*> LWORK is INTEGER +*> The length of WORK. LWORK >= MAX(1,2*N,3*N-2), and for +*> the best performance, LWORK >= MAX(1,N*NB), where NB is +*> the optimal blocksize for ZSYTRF_AA. +*> +*> If LWORK = -1, then a workspace query is assumed; the routine +*> only calculates the optimal size of the WORK array, returns +*> this value as the first entry of the WORK array, and no error +*> message related to LWORK is issued by XERBLA. +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit +*> < 0: if INFO = -i, the i-th argument had an illegal value +*> > 0: if INFO = i, D(i,i) is exactly zero. The factorization +*> has been completed, but the block diagonal matrix D is +*> exactly singular, so the solution could not be computed. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2016 +* +*> \ingroup complex16SYsolve +* +* ===================================================================== + SUBROUTINE ZSYSV_AA( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK, + $ LWORK, INFO ) +* +* -- LAPACK driver routine (version 3.7.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2016 +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER INFO, LDA, LDB, LWORK, N, NRHS +* .. +* .. Array Arguments .. + INTEGER IPIV( * ) + COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) +* .. +* +* ===================================================================== +* +* .. Local Scalars .. + LOGICAL LQUERY + INTEGER LWKOPT, LWKOPT_SYTRF, LWKOPT_SYTRS +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER ILAENV + EXTERNAL ILAENV, LSAME +* .. +* .. External Subroutines .. + EXTERNAL XERBLA, ZSYTRF, ZSYTRS, ZSYTRS2 +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX +* .. +* .. Executable Statements .. +* +* Test the input parameters. +* + INFO = 0 + LQUERY = ( LWORK.EQ.-1 ) + IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( NRHS.LT.0 ) THEN + INFO = -3 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -5 + ELSE IF( LDB.LT.MAX( 1, N ) ) THEN + INFO = -8 + ELSE IF( LWORK.LT.MAX(2*N, 3*N-2) .AND. .NOT.LQUERY ) THEN + INFO = -10 + END IF +* + IF( INFO.EQ.0 ) THEN + CALL ZSYTRF_AA( UPLO, N, A, LDA, IPIV, WORK, -1, INFO ) + LWKOPT_SYTRF = INT( WORK(1) ) + CALL ZSYTRS_AA( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK, + $ -1, INFO ) + LWKOPT_SYTRS = INT( WORK(1) ) + LWKOPT = MAX( LWKOPT_SYTRF, LWKOPT_SYTRS ) + WORK( 1 ) = LWKOPT + IF( LWORK.LT.LWKOPT .AND. .NOT.LQUERY ) THEN + INFO = -10 + END IF + END IF +* + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'ZSYSV_AA ', -INFO ) + RETURN + ELSE IF( LQUERY ) THEN + RETURN + END IF +* +* Compute the factorization A = U*T*U**T or A = L*T*L**T. +* + CALL ZSYTRF_AA( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO ) + IF( INFO.EQ.0 ) THEN +* +* Solve the system A*X = B, overwriting B with X. +* + CALL ZSYTRS_AA( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK, + $ LWORK, INFO ) +* + END IF +* + WORK( 1 ) = LWKOPT +* + RETURN +* +* End of ZSYSV_AA +* + END diff --git a/SRC/zsytrf_aa.f b/SRC/zsytrf_aa.f new file mode 100644 index 00000000..f82e5139 --- /dev/null +++ b/SRC/zsytrf_aa.f @@ -0,0 +1,480 @@ +*> \brief \b ZSYTRF_AA +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download ZSYTRF_AA + dependencies +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zsytrf_aa.f"> +*> [TGZ]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zsytrf_aa.f"> +*> [ZIP]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zsytrf_aa.f"> +*> [TXT]</a> +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE ZSYTRF_AA( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER UPLO +* INTEGER N, LDA, LWORK, INFO +* .. +* .. Array Arguments .. +* INTEGER IPIV( * ) +* COMPLEX*16 A( LDA, * ), WORK( * ) +* .. +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> ZSYTRF_AA computes the factorization of a complex 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 +*> +*> where U (or L) is a product of permutation and unit upper (lower) +*> triangular matrices, and T is a complex symmetric tridiagonal matrix. +*> +*> This is the blocked version of the algorithm, calling Level 3 BLAS. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] UPLO +*> \verbatim +*> UPLO is CHARACTER*1 +*> = 'U': Upper triangle of A is stored; +*> = 'L': Lower triangle of A is stored. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the matrix A. N >= 0. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is COMPLEX*16 array, dimension (LDA,N) +*> On entry, the symmetric matrix A. If UPLO = 'U', the leading +*> N-by-N upper triangular part of A contains the upper +*> triangular part of the matrix A, and the strictly lower +*> triangular part of A is not referenced. If UPLO = 'L', the +*> leading N-by-N lower triangular part of A contains the lower +*> triangular part of the matrix A, and the strictly upper +*> triangular part of A is not referenced. +*> +*> On exit, the tridiagonal matrix is stored in the diagonals +*> and the subdiagonals of A just below (or above) the diagonals, +*> and L is stored below (or above) the subdiaonals, when UPLO +*> is 'L' (or 'U'). +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N). +*> \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] WORK +*> \verbatim +*> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK)) +*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. +*> \endverbatim +*> +*> \param[in] LWORK +*> \verbatim +*> LWORK is INTEGER +*> The length of WORK. LWORK >=MAX(1,2*N). For optimum performance +*> LWORK >= N*(1+NB), where NB is the optimal blocksize. +*> +*> If LWORK = -1, then a workspace query is assumed; the routine +*> only calculates the optimal size of the WORK array, returns +*> this value as the first entry of the WORK array, and no error +*> message related to LWORK is issued by XERBLA. +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit +*> < 0: if INFO = -i, the i-th argument had an illegal value +*> > 0: if INFO = i, D(i,i) is exactly zero. The factorization +*> has been completed, but the block diagonal matrix D is +*> exactly singular, and division by zero will occur if it +*> is used to solve a system of equations. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2016 +* +*> \ingroup complex16SYcomputational +* +* ===================================================================== + SUBROUTINE ZSYTRF_AA( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO) +* +* -- LAPACK computational routine (version 3.7.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2016 +* + IMPLICIT NONE +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER N, LDA, LWORK, INFO +* .. +* .. Array Arguments .. + INTEGER IPIV( * ) + COMPLEX*16 A( LDA, * ), WORK( * ) +* .. +* +* ===================================================================== +* .. Parameters .. + COMPLEX*16 ZERO, ONE + PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) +* +* .. Local Scalars .. + LOGICAL LQUERY, UPPER + INTEGER J, LWKOPT, IINFO + INTEGER NB, MJ, NJ, K1, K2, J1, J2, J3, JB + COMPLEX*16 ALPHA +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER ILAENV + EXTERNAL LSAME, ILAENV +* .. +* .. External Subroutines .. + EXTERNAL XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX +* .. +* .. Executable Statements .. +* +* Determine the block size +* + NB = ILAENV( 1, 'ZSYTRF', UPLO, N, -1, -1, -1 ) +* +* Test the input parameters. +* + INFO = 0 + UPPER = LSAME( UPLO, 'U' ) + LQUERY = ( LWORK.EQ.-1 ) + IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -4 + ELSE IF( LWORK.LT.MAX( 1, 2*N ) .AND. .NOT.LQUERY ) THEN + INFO = -7 + END IF +* + IF( INFO.EQ.0 ) THEN + LWKOPT = (NB+1)*N + WORK( 1 ) = LWKOPT + END IF +* + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'ZSYTRF_AA', -INFO ) + RETURN + ELSE IF( LQUERY ) THEN + RETURN + END IF +* +* Quick return +* + IF ( N.EQ.0 ) THEN + RETURN + ENDIF + IPIV( 1 ) = 1 + IF ( N.EQ.1 ) THEN + IF ( A( 1, 1 ).EQ.ZERO ) THEN + INFO = 1 + END IF + RETURN + END IF +* +* Adjubst block size based on the workspace size +* + IF( LWORK.LT.((1+NB)*N) ) THEN + NB = ( LWORK-N ) / N + END IF +* + IF( UPPER ) THEN +* +* ..................................................... +* Factorize A as L*D*L**T using the upper triangle of A +* ..................................................... +* +* Copy first row A(1, 1:N) into H(1:n) (stored in WORK(1:N)) +* + CALL ZCOPY( N, A( 1, 1 ), LDA, WORK( 1 ), 1 ) +* +* J is the main loop index, increasing from 1 to N in steps of +* JB, where JB is the number of columns factorized by ZLASYF; +* JB is either NB, or N-J+1 for the last block +* + J = 0 + 10 CONTINUE + IF( J.GE.N ) + $ GO TO 20 +* +* each step of the main loop +* J is the last column of the previous panel +* J1 is the first column of the current panel +* K1 identifies if the previous column of the panel has been +* explicitly stored, e.g., K1=1 for the first panel, and +* K1=0 for the rest +* + J1 = J + 1 + JB = MIN( N-J1+1, NB ) + K1 = MAX(1, J)-J +* +* Panel factorization +* + CALL ZLASYF_AA( UPLO, 2-K1, N-J, JB, + $ A( MAX(1, J), J+1 ), LDA, + $ IPIV( J+1 ), WORK, N, WORK( N*NB+1 ), + $ IINFO ) + IF( (IINFO.GT.0) .AND. (INFO.EQ.0) ) THEN + INFO = IINFO+J + ENDIF +* +* Ajust 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 + IF( (J2.NE.IPIV(J2)) .AND. ((J1-K1).GT.2) ) THEN + CALL ZSWAP( J1-K1-2, A( 1, J2 ), 1, + $ A( 1, IPIV(J2) ), 1 ) + END IF + END DO + J = J + JB +* +* Trailing submatrix update, where +* the row A(J1-1, J2-1:N) stores U(J1, J2+1:N) and +* WORK stores the current block of the auxiriarly matrix H +* + IF( J.LT.N ) THEN +* +* If first panel and JB=1 (NB=1), then nothing to do +* + IF( J1.GT.1 .OR. JB.GT.1 ) THEN +* +* Merge rank-1 update with BLAS-3 update +* + ALPHA = A( J, J+1 ) + A( J, J+1 ) = ONE + CALL ZCOPY( N-J, A( J-1, J+1 ), LDA, + $ WORK( (J+1-J1+1)+JB*N ), 1 ) + CALL ZSCAL( N-J, ALPHA, WORK( (J+1-J1+1)+JB*N ), 1 ) +* +* K1 identifies if the previous column of the panel has been +* explicitly stored, e.g., K1=1 and K2= 0 for the first panel, +* while K1=0 and K2=1 for the rest +* + IF( J1.GT.1 ) THEN +* +* Not first panel +* + K2 = 1 + ELSE +* +* First panel +* + K2 = 0 +* +* First update skips the first column +* + JB = JB - 1 + END IF +* + DO J2 = J+1, N, NB + NJ = MIN( NB, N-J2+1 ) +* +* Update (J2, J2) diagonal block with ZGEMV +* + J3 = J2 + DO MJ = NJ-1, 1, -1 + CALL ZGEMV( 'No transpose', MJ, JB+1, + $ -ONE, WORK( J3-J1+1+K1*N ), N, + $ A( J1-K2, J3 ), 1, + $ ONE, A( J3, J3 ), LDA ) + J3 = J3 + 1 + END DO +* +* Update off-diagonal block of J2-th block row with ZGEMM +* + CALL ZGEMM( 'Transpose', 'Transpose', + $ NJ, N-J3+1, JB+1, + $ -ONE, A( J1-K2, J2 ), LDA, + $ WORK( J3-J1+1+K1*N ), N, + $ ONE, A( J2, J3 ), LDA ) + END DO +* +* Recover T( J, J+1 ) +* + A( J, J+1 ) = ALPHA + END IF +* +* WORK(J+1, 1) stores H(J+1, 1) +* + CALL ZCOPY( N-J, A( J+1, J+1 ), LDA, WORK( 1 ), 1 ) + END IF + GO TO 10 + ELSE +* +* ..................................................... +* Factorize A as L*D*L**T using the lower triangle of A +* ..................................................... +* +* copy first column A(1:N, 1) into H(1:N, 1) +* (stored in WORK(1:N)) +* + CALL ZCOPY( N, A( 1, 1 ), 1, WORK( 1 ), 1 ) +* +* J is the main loop index, increasing from 1 to N in steps of +* JB, where JB is the number of columns factorized by ZLASYF; +* JB is either NB, or N-J+1 for the last block +* + J = 0 + 11 CONTINUE + IF( J.GE.N ) + $ GO TO 20 +* +* each step of the main loop +* J is the last column of the previous panel +* J1 is the first column of the current panel +* K1 identifies if the previous column of the panel has been +* explicitly stored, e.g., K1=1 for the first panel, and +* K1=0 for the rest +* + J1 = J+1 + JB = MIN( N-J1+1, NB ) + K1 = MAX(1, J)-J +* +* Panel factorization +* + CALL ZLASYF_AA( UPLO, 2-K1, N-J, JB, + $ A( J+1, MAX(1, J) ), LDA, + $ IPIV( J+1 ), WORK, N, WORK( N*NB+1 ), IINFO) + IF( (IINFO.GT.0) .AND. (INFO.EQ.0) ) THEN + INFO = IINFO+J + ENDIF +* +* Ajust 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 + IF( (J2.NE.IPIV(J2)) .AND. ((J1-K1).GT.2) ) THEN + CALL ZSWAP( J1-K1-2, A( J2, 1 ), LDA, + $ A( IPIV(J2), 1 ), LDA ) + END IF + END DO + J = J + JB +* +* Trailing submatrix update, where +* A(J2+1, J1-1) stores L(J2+1, J1) and +* WORK(J2+1, 1) stores H(J2+1, 1) +* + IF( J.LT.N ) THEN +* +* if first panel and JB=1 (NB=1), then nothing to do +* + IF( J1.GT.1 .OR. JB.GT.1 ) THEN +* +* Merge rank-1 update with BLAS-3 update +* + ALPHA = A( J+1, J ) + A( J+1, J ) = ONE + CALL ZCOPY( N-J, A( J+1, J-1 ), 1, + $ WORK( (J+1-J1+1)+JB*N ), 1 ) + CALL ZSCAL( N-J, ALPHA, WORK( (J+1-J1+1)+JB*N ), 1 ) +* +* K1 identifies if the previous column of the panel has been +* explicitly stored, e.g., K1=1 and K2= 0 for the first panel, +* while K1=0 and K2=1 for the rest +* + IF( J1.GT.1 ) THEN +* +* Not first panel +* + K2 = 1 + ELSE +* +* First panel +* + K2 = 0 +* +* First update skips the first column +* + JB = JB - 1 + END IF +* + DO J2 = J+1, N, NB + NJ = MIN( NB, N-J2+1 ) +* +* Update (J2, J2) diagonal block with ZGEMV +* + J3 = J2 + DO MJ = NJ-1, 1, -1 + CALL ZGEMV( 'No transpose', MJ, JB+1, + $ -ONE, WORK( J3-J1+1+K1*N ), N, + $ A( J3, J1-K2 ), LDA, + $ ONE, A( J3, J3 ), 1 ) + J3 = J3 + 1 + END DO +* +* Update off-diagonal block in J2-th block column with ZGEMM +* + CALL ZGEMM( 'No transpose', 'Transpose', + $ N-J3+1, NJ, JB+1, + $ -ONE, WORK( J3-J1+1+K1*N ), N, + $ A( J2, J1-K2 ), LDA, + $ ONE, A( J3, J2 ), LDA ) + END DO +* +* Recover T( J+1, J ) +* + A( J+1, J ) = ALPHA + END IF +* +* WORK(J+1, 1) stores H(J+1, 1) +* + CALL ZCOPY( N-J, A( J+1, J+1 ), 1, WORK( 1 ), 1 ) + END IF + GO TO 11 + END IF +* + 20 CONTINUE + RETURN +* +* End of ZSYTRF_AA +* + END diff --git a/SRC/zsytrs_aa.f b/SRC/zsytrs_aa.f new file mode 100644 index 00000000..cae83a76 --- /dev/null +++ b/SRC/zsytrs_aa.f @@ -0,0 +1,285 @@ +*> \brief \b ZSYTRS_AA +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download ZSYTRS_AA + dependencies +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zsytrs_aa.f"> +*> [TGZ]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zsytrs_aa.f"> +*> [ZIP]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zsytrs_aa.f"> +*> [TXT]</a> +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE ZSYTRS_AA( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, +* WORK, LWORK, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER UPLO +* INTEGER N, NRHS, LDA, LDB, LWORK, INFO +* .. +* .. Array Arguments .. +* INTEGER IPIV( * ) +* COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> ZSYTRS_AA solves a system of linear equations A*X = B with a complex +*> symmetric matrix A using the factorization A = U*T*U**T or +*> A = L*T*L**T computed by ZSYTRF_AA. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] UPLO +*> \verbatim +*> 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; +*> = 'L': Lower triangular, form is A = L*T*L**T. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the matrix A. N >= 0. +*> \endverbatim +*> +*> \param[in] NRHS +*> \verbatim +*> NRHS is INTEGER +*> The number of right hand sides, i.e., the number of columns +*> of the matrix B. NRHS >= 0. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is COMPLEX*16 array, dimension (LDA,N) +*> Details of factors computed by ZSYTRF_AA. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N). +*> \endverbatim +*> +*> \param[in] IPIV +*> \verbatim +*> IPIV is INTEGER array, dimension (N) +*> Details of the interchanges as computed by ZSYTRF_AA. +*> \endverbatim +*> +*> \param[in,out] B +*> \verbatim +*> B is COMPLEX*16 array, dimension (LDB,NRHS) +*> On entry, the right hand side matrix B. +*> On exit, the solution matrix X. +*> \endverbatim +*> +*> \param[in] LDB +*> \verbatim +*> LDB is INTEGER +*> The leading dimension of the array B. LDB >= max(1,N). +*> \endverbatim +*> +*> \param[in] WORK +*> \verbatim +*> WORK is DOUBLE array, dimension (MAX(1,LWORK)) +*> \endverbatim +*> +*> \param[in] LWORK +*> \verbatim +*> LWORK is INTEGER, LWORK >= MAX(1,3*N-2). +*> +*> \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 2016 +* +*> \ingroup complex16SYcomputational +* +* ===================================================================== + SUBROUTINE ZSYTRS_AA( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, + $ WORK, LWORK, INFO ) +* +* -- LAPACK computational routine (version 3.7.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2016 +* + IMPLICIT NONE +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER N, NRHS, LDA, LDB, LWORK, INFO +* .. +* .. Array Arguments .. + INTEGER IPIV( * ) + COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) +* .. +* +* ===================================================================== +* + COMPLEX*16 ONE + PARAMETER ( ONE = 1.0D+0 ) +* .. +* .. Local Scalars .. + LOGICAL LQUERY, UPPER + INTEGER K, KP, LWKOPT +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. External Subroutines .. + EXTERNAL ZGTSV, ZSWAP, ZTRSM, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX +* .. +* .. Executable Statements .. +* + INFO = 0 + UPPER = LSAME( UPLO, 'U' ) + LQUERY = ( LWORK.EQ.-1 ) + IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( NRHS.LT.0 ) THEN + INFO = -3 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -5 + ELSE IF( LDB.LT.MAX( 1, N ) ) THEN + INFO = -8 + ELSE IF( LWORK.LT.MAX( 1, 3*N-2 ) .AND. .NOT.LQUERY ) THEN + INFO = -10 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'ZSYTRS_AA', -INFO ) + RETURN + ELSE IF( LQUERY ) THEN + LWKOPT = (3*N-2) + WORK( 1 ) = LWKOPT + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 .OR. NRHS.EQ.0 ) + $ RETURN +* + IF( UPPER ) THEN +* +* Solve A*X = B, where A = U*T*U**T. +* +* Pivot, P**T * B +* + DO K = 1, N + KP = IPIV( K ) + IF( KP.NE.K ) + $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) + END DO +* +* Compute (U \P**T * B) -> B [ (U \P**T * B) ] +* + CALL ZTRSM('L', 'U', 'T', 'U', N-1, NRHS, ONE, A( 1, 2 ), LDA, + $ B( 2, 1 ), LDB) +* +* Compute T \ B -> B [ T \ (U \P**T * B) ] +* + CALL ZLACPY( 'F', 1, N, A( 1, 1 ), LDA+1, WORK( N ), 1) + IF( N.GT.1 ) THEN + CALL ZLACPY( 'F', 1, N-1, A( 1, 2 ), LDA+1, WORK( 1 ), 1 ) + CALL ZLACPY( 'F', 1, N-1, A( 1, 2 ), LDA+1, WORK( 2*N ), 1 ) + END IF + CALL ZGTSV( N, NRHS, WORK( 1 ), WORK( N ), WORK( 2*N ), B, LDB, + $ INFO ) +* +* Compute (U**T \ B) -> B [ U**T \ (T \ (U \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**T \ (T \ (U \P**T * B) )) ] +* + DO K = N, 1, -1 + KP = IPIV( K ) + IF( KP.NE.K ) + $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) + END DO +* + ELSE +* +* Solve A*X = B, where A = L*T*L**T. +* +* Pivot, P**T * B +* + DO K = 1, N + KP = IPIV( K ) + IF( KP.NE.K ) + $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) + END DO +* +* Compute (L \P**T * B) -> B [ (L \P**T * B) ] +* + CALL ZTRSM( 'L', 'L', 'N', 'U', N-1, NRHS, ONE, A( 2, 1 ), LDA, + $ B( 2, 1 ), LDB) +* +* Compute T \ B -> B [ T \ (L \P**T * B) ] +* + CALL ZLACPY( 'F', 1, N, A(1, 1), LDA+1, WORK(N), 1) + IF( N.GT.1 ) THEN + CALL ZLACPY( 'F', 1, N-1, A( 2, 1 ), LDA+1, WORK( 1 ), 1 ) + CALL ZLACPY( 'F', 1, N-1, A( 2, 1 ), LDA+1, WORK( 2*N ), 1 ) + END IF + CALL ZGTSV( N, NRHS, WORK( 1 ), WORK(N), WORK( 2*N ), B, LDB, + $ INFO) +* +* Compute (L**T \ B) -> B [ L**T \ (T \ (L \P**T * B) ) ] +* + CALL ZTRSM( '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) )) ] +* + DO K = N, 1, -1 + KP = IPIV( K ) + IF( KP.NE.K ) + $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) + END DO +* + END IF +* + RETURN +* +* End of ZSYTRS_AA +* + END diff --git a/TESTING/LIN/Makefile b/TESTING/LIN/Makefile index 15d5e94f..879cee2e 100644 --- a/TESTING/LIN/Makefile +++ b/TESTING/LIN/Makefile @@ -90,11 +90,11 @@ CLINTST = cchkaa.o \ cchkeq.o cchkgb.o cchkge.o cchkgt.o \ cchkhe.o cchkhe_rook.o cchkhe_rk.o cchkhe_aa.o cchkhp.o cchklq.o cchkpb.o \ cchkpo.o cchkps.o cchkpp.o cchkpt.o cchkq3.o cchkql.o \ - cchkqr.o cchkrq.o cchksp.o cchksy.o cchksy_rook.o cchksy_rk.o cchktb.o \ + cchkqr.o cchkrq.o cchksp.o cchksy.o cchksy_rook.o cchksy_rk.o cchksy_aa.o cchktb.o \ cchktp.o cchktr.o cchktz.o \ cdrvgt.o cdrvhe_rook.o cdrvhe_rk.o cdrvhe_aa.o cdrvhp.o \ cdrvls.o cdrvpb.o cdrvpp.o cdrvpt.o \ - cdrvsp.o cdrvsy_rook.o cdrvsy_rk.o \ + cdrvsp.o cdrvsy_rook.o cdrvsy_rk.o cdrvsy_aa.o \ cerrgt.o cerrlq.o \ cerrls.o cerrps.o cerrql.o cerrqp.o \ cerrqr.o cerrrq.o cerrtr.o cerrtz.o \ @@ -113,7 +113,7 @@ CLINTST = cchkaa.o \ cqrt12.o cqrt13.o cqrt14.o cqrt15.o cqrt16.o \ cqrt17.o crqt01.o crqt02.o crqt03.o crzt01.o crzt02.o \ csbmv.o cspt01.o \ - cspt02.o cspt03.o csyt01.o csyt01_rook.o csyt01_3.o csyt02.o csyt03.o \ + cspt02.o cspt03.o csyt01.o csyt01_rook.o csyt01_3.o csyt01_aa.o csyt02.o csyt03.o \ ctbt02.o ctbt03.o ctbt05.o ctbt06.o ctpt01.o \ ctpt02.o ctpt03.o ctpt05.o ctpt06.o ctrt01.o \ ctrt02.o ctrt03.o ctrt05.o ctrt06.o \ @@ -174,11 +174,11 @@ ZLINTST = zchkaa.o \ zchkeq.o zchkgb.o zchkge.o zchkgt.o \ zchkhe.o zchkhe_rook.o zchkhe_rk.o zchkhe_aa.o zchkhp.o zchklq.o zchkpb.o \ zchkpo.o zchkps.o zchkpp.o zchkpt.o zchkq3.o zchkql.o \ - zchkqr.o zchkrq.o zchksp.o zchksy.o zchksy_rook.o zchksy_rk.o zchktb.o \ + zchkqr.o zchkrq.o zchksp.o zchksy.o zchksy_rook.o zchksy_rk.o zchksy_aa.o zchktb.o \ zchktp.o zchktr.o zchktz.o \ zdrvgt.o zdrvhe_rook.o zdrvhe_rk.o zdrvhe_aa.o zdrvhp.o \ zdrvls.o zdrvpb.o zdrvpp.o zdrvpt.o \ - zdrvsp.o zdrvsy_rook.o zdrvsy_rk.o \ + zdrvsp.o zdrvsy_rook.o zdrvsy_rk.o zdrvsy_aa.o \ zerrgt.o zerrlq.o \ zerrls.o zerrps.o zerrql.o zerrqp.o \ zerrqr.o zerrrq.o zerrtr.o zerrtz.o \ @@ -197,7 +197,7 @@ ZLINTST = zchkaa.o \ zqrt12.o zqrt13.o zqrt14.o zqrt15.o zqrt16.o \ zqrt17.o zrqt01.o zrqt02.o zrqt03.o zrzt01.o zrzt02.o \ zsbmv.o zspt01.o \ - zspt02.o zspt03.o zsyt01.o zsyt01_rook.o zsyt01_3.o zsyt02.o zsyt03.o \ + zspt02.o zspt03.o zsyt01.o zsyt01_rook.o zsyt01_3.o zsyt01_aa.o zsyt02.o zsyt03.o \ ztbt02.o ztbt03.o ztbt05.o ztbt06.o ztpt01.o \ ztpt02.o ztpt03.o ztpt05.o ztpt06.o ztrt01.o \ ztrt02.o ztrt03.o ztrt05.o ztrt06.o \ diff --git a/TESTING/LIN/cchkaa.f b/TESTING/LIN/cchkaa.f index cf04e78d..f2ef59f1 100644 --- a/TESTING/LIN/cchkaa.f +++ b/TESTING/LIN/cchkaa.f @@ -165,12 +165,12 @@ $ CCHKHE_ROOK, CCHKHE_RK, CCHKHE_AA, CCHKLQ, $ CCHKPB,CCHKPO, CCHKPS, CCHKPP, CCHKPT, CCHKQ3, $ CCHKQL, CCHKQR, CCHKRQ, CCHKSP, CCHKSY, - $ CCHKSY_ROOK, CCHKSY_RK, CCHKTB, CCHKTP, - $ CCHKTR, CCHKTZ, CDRVGB, CDRVGE, CDRVGT, CDRVHE, - $ CDRVHE_ROOK, CDRVHE_RK, CDRVHE_AA, CDRVHP, - $ CDRVLS, CDRVPB, CDRVPO, CDRVPP, CDRVPT, CDRVSP, - $ CDRVSY, CDRVSY_ROOK, CDRVSY_RK, ILAVER, CCHKQRT, - $ CCHKQRTP + $ CCHKSY_ROOK, CCHKSY_RK, CCHKSY_AA, CCHKTB, + $ CCHKTP, CCHKTR, CCHKTZ, CDRVGB, CDRVGE, CDRVGT, + $ CDRVHE, CDRVHE_ROOK, CDRVHE_RK, CDRVHE_AA, + $ CDRVHP, CDRVLS, CDRVPB, CDRVPO, CDRVPP, CDRVPT, + $ CDRVSP, CDRVSY, CDRVSY_ROOK, CDRVSY_RK, + $ CDRVSY_AA, ILAVER, CCHKQRT, CCHKQRTP * .. * .. Scalars in Common .. LOGICAL LERR, OK @@ -830,6 +830,31 @@ WRITE( NOUT, FMT = 9988 )PATH END IF * + ELSE IF( LSAMEN( 2, C2, 'SA' ) ) THEN +* +* SA: symmetric indefinite matrices with Aasen's algorithm, +* + NTYPES = 11 + CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) +* + IF( TSTCHK ) THEN + CALL CCHKSY_AA( DOTYPE, NN, NVAL, NNB2, NBVAL2, NNS, NSVAL, + $ THRESH, TSTERR, LDA, A( 1, 1 ), A( 1, 2 ), + $ A( 1, 3 ), B( 1, 1 ), B( 1, 2 ), + $ B( 1, 3 ), WORK, RWORK, IWORK, NOUT ) + ELSE + WRITE( NOUT, FMT = 9989 )PATH + END IF +* + IF( TSTDRV ) THEN + CALL CDRVSY_AA( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, + $ LDA, A( 1, 1 ), A( 1, 2 ), A( 1, 3 ), + $ B( 1, 1 ), B( 1, 2 ), B( 1, 3 ), WORK, + $ RWORK, IWORK, NOUT ) + ELSE + WRITE( NOUT, FMT = 9988 )PATH + END IF +* ELSE IF( LSAMEN( 2, C2, 'SP' ) ) THEN * * SP: symmetric indefinite packed matrices, diff --git a/TESTING/LIN/cchkhe_aa.f b/TESTING/LIN/cchkhe_aa.f index 702677bd..ca23c0f2 100644 --- a/TESTING/LIN/cchkhe_aa.f +++ b/TESTING/LIN/cchkhe_aa.f @@ -205,13 +205,13 @@ PARAMETER ( NTESTS = 9 ) * .. * .. Local Scalars .. - LOGICAL TRFCON, ZEROT + LOGICAL ZEROT CHARACTER DIST, TYPE, UPLO, XTYPE CHARACTER*3 PATH, MATPATH INTEGER I, I1, I2, IMAT, IN, INB, INFO, IOFF, IRHS, $ IUPLO, IZERO, J, K, KL, KU, LDA, LWORK, MODE, $ N, NB, NERRS, NFAIL, NIMAT, NRHS, NRUN, NT - REAL ANORM, CNDNUM, RCOND, RCONDC + REAL ANORM, CNDNUM * .. * .. Local Arrays .. CHARACTER UPLOS( 2 ) @@ -224,7 +224,7 @@ * .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, XLAENV, CERRHE, CGET04, - $ ZHECON, CHERFS, CHET01, CHETRF_AA, ZHETRI2, + $ ZHECON, CHERFS, CHET01_AA, CHETRF_AA, ZHETRI2, $ CHETRS_AA, CLACPY, CLAIPD, CLARHS, CLATB4, $ CLATMS, CPOT02, ZPOT03, ZPOT05 * .. @@ -431,10 +431,10 @@ * the block structure of D. AINV is a work array for * block factorization, LWORK is the length of AINV. * - LWORK = ( NB+1 )*LDA + LWORK = MAX( 1, ( NB+1 )*LDA ) SRNAMT = 'CHETRF_AA' CALL CHETRF_AA( UPLO, N, AFAC, LDA, IWORK, AINV, - $ LWORK, INFO ) + $ LWORK, INFO ) * * Adjust the expected value of INFO to account for * pivoting. @@ -464,19 +464,11 @@ $ NOUT ) END IF * -* Set the condition estimate flag if the INFO is not 0. -* - IF( INFO.NE.0 ) THEN - TRFCON = .TRUE. - ELSE - TRFCON = .FALSE. - END IF -* *+ TEST 1 * Reconstruct matrix from factors and compute residual. * CALL CHET01_AA( UPLO, N, A, LDA, AFAC, LDA, IWORK, - $ AINV, LDA, RWORK, RESULT( 1 ) ) + $ AINV, LDA, RWORK, RESULT( 1 ) ) NT = 1 * * @@ -494,10 +486,9 @@ 110 CONTINUE NRUN = NRUN + NT * -* Do only the condition estimate if INFO is not 0. +* Skip solver test if INFO is not 0. * - IF( TRFCON ) THEN - RCONDC = ZERO + IF( INFO.NE.0 ) THEN GO TO 140 END IF * @@ -506,7 +497,7 @@ DO 130 IRHS = 1, NNS NRHS = NSVAL( IRHS ) * -*+ TEST 3 (Using TRS) +*+ TEST 2 (Using TRS) * Solve and compute residual for A * X = B. * * Choose a set of NRHS random solution vectors @@ -519,9 +510,9 @@ CALL CLACPY( 'Full', N, NRHS, B, LDA, X, LDA ) * SRNAMT = 'CHETRS_AA' - LWORK = 3*N-2 + LWORK = MAX( 1, 3*N-2 ) CALL CHETRS_AA( UPLO, N, NRHS, AFAC, LDA, IWORK, - $ X, LDA, WORK, LWORK, INFO ) + $ X, LDA, WORK, LWORK, INFO ) * * Check error code from CHETRS and handle error. * diff --git a/TESTING/LIN/cchksy_aa.f b/TESTING/LIN/cchksy_aa.f new file mode 100644 index 00000000..534be92f --- /dev/null +++ b/TESTING/LIN/cchksy_aa.f @@ -0,0 +1,572 @@ +*> \brief \b CCHKSY_AA +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE CCHKSY_AA( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL, +* THRESH, TSTERR, NMAX, A, AFAC, AINV, B, X, +* XACT, WORK, RWORK, IWORK, NOUT ) +* +* .. Scalar Arguments .. +* LOGICAL TSTERR +* INTEGER NMAX, NN, NNB, NNS, NOUT +* REAL THRESH +* .. +* .. Array Arguments .. +* LOGICAL DOTYPE( * ) +* INTEGER IWORK( * ), NBVAL( * ), NSVAL( * ), NVAL( * ) +* COMPLEX A( * ), AFAC( * ), AINV( * ), B( * ), +* $ RWORK( * ), WORK( * ), X( * ), XACT( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> CCHKSY_AA tests CSYTRF_AA, -TRS_AA. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] DOTYPE +*> \verbatim +*> DOTYPE is LOGICAL array, dimension (NTYPES) +*> The matrix types to be used for testing. Matrices of type j +*> (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) = +*> .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used. +*> \endverbatim +*> +*> \param[in] NN +*> \verbatim +*> NN is INTEGER +*> The number of values of N contained in the vector NVAL. +*> \endverbatim +*> +*> \param[in] NVAL +*> \verbatim +*> NVAL is INTEGER array, dimension (NN) +*> The values of the matrix dimension N. +*> \endverbatim +*> +*> \param[in] NNB +*> \verbatim +*> NNB is INTEGER +*> The number of values of NB contained in the vector NBVAL. +*> \endverbatim +*> +*> \param[in] NBVAL +*> \verbatim +*> NBVAL is INTEGER array, dimension (NBVAL) +*> The values of the blocksize NB. +*> \endverbatim +*> +*> \param[in] NNS +*> \verbatim +*> NNS is INTEGER +*> The number of values of NRHS contained in the vector NSVAL. +*> \endverbatim +*> +*> \param[in] NSVAL +*> \verbatim +*> NSVAL is INTEGER array, dimension (NNS) +*> The values of the number of right hand sides NRHS. +*> \endverbatim +*> +*> \param[in] THRESH +*> \verbatim +*> THRESH is REAL +*> The threshold value for the test ratios. A result is +*> included in the output file if RESULT >= THRESH. To have +*> every test ratio printed, use THRESH = 0. +*> \endverbatim +*> +*> \param[in] TSTERR +*> \verbatim +*> TSTERR is LOGICAL +*> Flag that indicates whether error exits are to be tested. +*> \endverbatim +*> +*> \param[in] NMAX +*> \verbatim +*> NMAX is INTEGER +*> The maximum value permitted for N, used in dimensioning the +*> work arrays. +*> \endverbatim +*> +*> \param[out] A +*> \verbatim +*> A is REAL array, dimension (NMAX*NMAX) +*> \endverbatim +*> +*> \param[out] AFAC +*> \verbatim +*> AFAC is REAL array, dimension (NMAX*NMAX) +*> \endverbatim +*> +*> \param[out] AINV +*> \verbatim +*> AINV is REAL array, dimension (NMAX*NMAX) +*> \endverbatim +*> +*> \param[out] B +*> \verbatim +*> B is REAL array, dimension (NMAX*NSMAX) +*> where NSMAX is the largest entry in NSVAL. +*> \endverbatim +*> +*> \param[out] X +*> \verbatim +*> X is REAL array, dimension (NMAX*NSMAX) +*> \endverbatim +*> +*> \param[out] XACT +*> \verbatim +*> XACT is REAL array, dimension (NMAX*NSMAX) +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is REAL array, dimension (NMAX*max(3,NSMAX)) +*> \endverbatim +*> +*> \param[out] RWORK +*> \verbatim +*> RWORK is REAL array, dimension (max(NMAX,2*NSMAX)) +*> \endverbatim +*> +*> \param[out] IWORK +*> \verbatim +*> IWORK is INTEGER array, dimension (2*NMAX) +*> \endverbatim +*> +*> \param[in] NOUT +*> \verbatim +*> NOUT is INTEGER +*> The unit number for output. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2016 +* +* @generated from LIN/dchksy_aa.f, fortran d -> c, Wed Nov 16 21:34:18 2016 +* +*> \ingroup complex_lin +* +* ===================================================================== + SUBROUTINE CCHKSY_AA( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL, + $ THRESH, TSTERR, NMAX, A, AFAC, AINV, B, + $ X, XACT, WORK, RWORK, IWORK, NOUT ) +* +* -- LAPACK test routine (version 3.7.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2016 +* + IMPLICIT NONE +* +* .. Scalar Arguments .. + LOGICAL TSTERR + INTEGER NN, NNB, NNS, NMAX, NOUT + REAL THRESH +* .. +* .. Array Arguments .. + LOGICAL DOTYPE( * ) + INTEGER IWORK( * ), NBVAL( * ), NSVAL( * ), NVAL( * ) + REAL RWORK( * ) + COMPLEX A( * ), AFAC( * ), AINV( * ), B( * ), + $ WORK( * ), X( * ), XACT( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + REAL ZERO + PARAMETER ( ZERO = 0.0D+0 ) + COMPLEX CZERO + PARAMETER ( CZERO = 0.0E+0 ) + INTEGER NTYPES + PARAMETER ( NTYPES = 10 ) + INTEGER NTESTS + PARAMETER ( NTESTS = 9 ) +* .. +* .. Local Scalars .. + LOGICAL ZEROT + CHARACTER DIST, TYPE, UPLO, XTYPE + CHARACTER*3 PATH, MATPATH + INTEGER I, I1, I2, IMAT, IN, INB, INFO, IOFF, IRHS, + $ IUPLO, IZERO, J, K, KL, KU, LDA, LWORK, MODE, + $ N, NB, NERRS, NFAIL, NIMAT, NRHS, NRUN, NT + REAL ANORM, CNDNUM +* .. +* .. Local Arrays .. + CHARACTER UPLOS( 2 ) + INTEGER ISEED( 4 ), ISEEDY( 4 ) + REAL RESULT( NTESTS ) +* .. +* .. External Functions .. + REAL DGET06, CLANSY + EXTERNAL DGET06, CLANSY +* .. +* .. External Subroutines .. + EXTERNAL ALAERH, ALAHD, ALASUM, CERRSY, CGET04, CLACPY, + $ CLARHS, CLATB4, CLATMS, CSYT02, DSYT03, DSYT05, + $ DSYCON, CSYRFS, CSYT01_AA, CSYTRF_AA, + $ DSYTRI2, CSYTRS_AA, XLAENV +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. Scalars in Common .. + LOGICAL LERR, OK + CHARACTER*32 SRNAMT + INTEGER INFOT, NUNIT +* .. +* .. Common blocks .. + COMMON / INFOC / INFOT, NUNIT, OK, LERR + COMMON / SRNAMC / SRNAMT +* .. +* .. Data statements .. + DATA ISEEDY / 1988, 1989, 1990, 1991 / + DATA UPLOS / 'U', 'L' / +* .. +* .. Executable Statements .. +* +* Initialize constants and the random number seed. +* +* Test path +* + PATH( 1: 1 ) = 'Complex precision' + PATH( 2: 3 ) = 'SA' +* +* Path to generate matrices +* + MATPATH( 1: 1 ) = 'Complex precision' + MATPATH( 2: 3 ) = 'SY' + NRUN = 0 + NFAIL = 0 + NERRS = 0 + DO 10 I = 1, 4 + ISEED( I ) = ISEEDY( I ) + 10 CONTINUE +* +* Test the error exits +* + IF( TSTERR ) + $ CALL CERRSY( PATH, NOUT ) + INFOT = 0 +* +* Set the minimum block size for which the block routine should +* be used, which will be later returned by ILAENV +* + CALL XLAENV( 2, 2 ) +* +* Do for each value of N in NVAL +* + DO 180 IN = 1, NN + N = NVAL( IN ) + IF( N .GT. NMAX ) THEN + NFAIL = NFAIL + 1 + WRITE(NOUT, 9995) 'M ', N, NMAX + GO TO 180 + END IF + LDA = MAX( N, 1 ) + XTYPE = 'N' + NIMAT = NTYPES + IF( N.LE.0 ) + $ NIMAT = 1 +* + IZERO = 0 +* +* Do for each value of matrix type IMAT +* + DO 170 IMAT = 1, NIMAT +* +* Do the tests only if DOTYPE( IMAT ) is true. +* + IF( .NOT.DOTYPE( IMAT ) ) + $ GO TO 170 +* +* Skip types 3, 4, 5, or 6 if the matrix size is too small. +* + ZEROT = IMAT.GE.3 .AND. IMAT.LE.6 + IF( ZEROT .AND. N.LT.IMAT-2 ) + $ GO TO 170 +* +* Do first for UPLO = 'U', then for UPLO = 'L' +* + DO 160 IUPLO = 1, 2 + UPLO = UPLOS( IUPLO ) +* +* Begin generate the test matrix A. +* +* +* Set up parameters with CLATB4 for the matrix generator +* based on the type of matrix to be generated. +* + CALL CLATB4( MATPATH, IMAT, N, N, TYPE, KL, KU, + $ ANORM, MODE, CNDNUM, DIST ) +* +* Generate a matrix with CLATMS. +* + SRNAMT = 'CLATMS' + CALL CLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE, + $ CNDNUM, ANORM, KL, KU, UPLO, A, LDA, WORK, + $ INFO ) +* +* Check error code from CLATMS and handle error. +* + IF( INFO.NE.0 ) THEN + CALL ALAERH( PATH, 'CLATMS', INFO, 0, UPLO, N, N, -1, + $ -1, -1, IMAT, NFAIL, NERRS, NOUT ) +* +* Skip all tests for this generated matrix +* + GO TO 160 + END IF +* +* For matrix types 3-6, zero one or more rows and +* columns of the matrix to test that INFO is returned +* correctly. +* + IF( ZEROT ) THEN + IF( IMAT.EQ.3 ) THEN + IZERO = 1 + ELSE IF( IMAT.EQ.4 ) THEN + IZERO = N + ELSE + IZERO = N / 2 + 1 + END IF +* + IF( IMAT.LT.6 ) THEN +* +* Set row and column IZERO to zero. +* + IF( IUPLO.EQ.1 ) THEN + IOFF = ( IZERO-1 )*LDA + DO 20 I = 1, IZERO - 1 + A( IOFF+I ) = CZERO + 20 CONTINUE + IOFF = IOFF + IZERO + DO 30 I = IZERO, N + A( IOFF ) = CZERO + IOFF = IOFF + LDA + 30 CONTINUE + ELSE + IOFF = IZERO + DO 40 I = 1, IZERO - 1 + A( IOFF ) = CZERO + IOFF = IOFF + LDA + 40 CONTINUE + IOFF = IOFF - IZERO + DO 50 I = IZERO, N + A( IOFF+I ) = CZERO + 50 CONTINUE + END IF + ELSE + IF( IUPLO.EQ.1 ) THEN +* +* Set the first IZERO rows and columns to zero. +* + IOFF = 0 + DO 70 J = 1, N + I2 = MIN( J, IZERO ) + DO 60 I = 1, I2 + A( IOFF+I ) = CZERO + 60 CONTINUE + IOFF = IOFF + LDA + 70 CONTINUE + IZERO = 1 + ELSE +* +* Set the last IZERO rows and columns to zero. +* + IOFF = 0 + DO 90 J = 1, N + I1 = MAX( J, IZERO ) + DO 80 I = I1, N + A( IOFF+I ) = CZERO + 80 CONTINUE + IOFF = IOFF + LDA + 90 CONTINUE + END IF + END IF + ELSE + IZERO = 0 + END IF +* +* End generate the test matrix A. +* +* Do for each value of NB in NBVAL +* + DO 150 INB = 1, NNB +* +* Set the optimal blocksize, which will be later +* returned by ILAENV. +* + NB = NBVAL( INB ) + CALL XLAENV( 1, NB ) +* +* Copy the test matrix A into matrix AFAC which +* will be factorized in place. This is needed to +* preserve the test matrix A for subsequent tests. +* + CALL CLACPY( UPLO, N, N, A, LDA, AFAC, LDA ) +* +* Compute the L*D*L**T or U*D*U**T factorization of the +* matrix. IWORK stores details of the interchanges and +* the block structure of D. AINV is a work array for +* block factorization, LWORK is the length of AINV. +* + SRNAMT = 'CSYTRF_AA' + LWORK = MAX( 1, N*NB + N ) + CALL CSYTRF_AA( UPLO, N, AFAC, LDA, IWORK, AINV, + $ LWORK, INFO ) +* +* Adjust the expected value of INFO to account for +* pivoting. +* + IF( IZERO.GT.0 ) THEN + J = 1 + K = IZERO + 100 CONTINUE + IF( J.EQ.K ) THEN + K = IWORK( J ) + ELSE IF( IWORK( J ).EQ.K ) THEN + K = J + END IF + IF( J.LT.K ) THEN + J = J + 1 + GO TO 100 + END IF + ELSE + K = 0 + END IF +* +* Check error code from CSYTRF and handle error. +* + IF( INFO.NE.K ) THEN + CALL ALAERH( PATH, 'CSYTRF_AA', INFO, K, UPLO, + $ N, N, -1, -1, NB, IMAT, NFAIL, NERRS, + $ NOUT ) + END IF +* +*+ TEST 1 +* Reconstruct matrix from factors and compute residual. +* + CALL CSYT01_AA( UPLO, N, A, LDA, AFAC, LDA, IWORK, + $ AINV, LDA, RWORK, RESULT( 1 ) ) + NT = 1 +* +* +* Print information about the tests that did not pass +* the threshold. +* + DO 110 K = 1, NT + IF( RESULT( K ).GE.THRESH ) THEN + IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) + $ CALL ALAHD( NOUT, PATH ) + WRITE( NOUT, FMT = 9999 )UPLO, N, NB, IMAT, K, + $ RESULT( K ) + NFAIL = NFAIL + 1 + END IF + 110 CONTINUE + NRUN = NRUN + NT +* +* Skip solver test if INFO is not 0. +* + IF( INFO.NE.0 ) THEN + GO TO 140 + END IF +* +* Do for each value of NRHS in NSVAL. +* + DO 130 IRHS = 1, NNS + NRHS = NSVAL( IRHS ) +* +*+ TEST 2 (Using TRS) +* Solve and compute residual for A * X = B. +* +* Choose a set of NRHS random solution vectors +* stored in XACT and set up the right hand side B +* + SRNAMT = 'CLARHS' + CALL CLARHS( MATPATH, XTYPE, UPLO, ' ', N, N, + $ KL, KU, NRHS, A, LDA, XACT, LDA, + $ B, LDA, ISEED, INFO ) + CALL CLACPY( 'Full', N, NRHS, B, LDA, X, LDA ) +* + SRNAMT = 'CSYTRS_AA' + LWORK = MAX( 1, 3*N-2 ) + CALL CSYTRS_AA( UPLO, N, NRHS, AFAC, LDA, + $ IWORK, X, LDA, WORK, LWORK, + $ INFO ) +* +* Check error code from CSYTRS and handle error. +* + IF( INFO.NE.0 ) THEN + CALL ALAERH( PATH, 'CSYTRS_AA', INFO, 0, + $ UPLO, N, N, -1, -1, NRHS, IMAT, + $ NFAIL, NERRS, NOUT ) + END IF +* + CALL CLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA ) +* +* Compute the residual for the solution +* + CALL CSYT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK, + $ LDA, RWORK, RESULT( 2 ) ) +* +* +* Print information about the tests that did not pass +* the threshold. +* + DO 120 K = 2, 2 + IF( RESULT( K ).GE.THRESH ) THEN + IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) + $ CALL ALAHD( NOUT, PATH ) + WRITE( NOUT, FMT = 9998 )UPLO, N, NRHS, + $ IMAT, K, RESULT( K ) + NFAIL = NFAIL + 1 + END IF + 120 CONTINUE + NRUN = NRUN + 1 +* +* End do for each value of NRHS in NSVAL. +* + 130 CONTINUE + 140 CONTINUE + 150 CONTINUE + 160 CONTINUE + 170 CONTINUE + 180 CONTINUE +* +* Print a summary of the results. +* + CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS ) +* + 9999 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', NB =', I4, ', type ', + $ I2, ', test ', I2, ', ratio =', G12.5 ) + 9998 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', NRHS=', I3, ', type ', + $ I2, ', test(', I2, ') =', G12.5 ) + 9995 FORMAT( ' Invalid input value: ', A4, '=', I6, '; must be <=', + $ I6 ) + RETURN +* +* End of CCHKSY_AA +* + END diff --git a/TESTING/LIN/cdrvhe_aa.f b/TESTING/LIN/cdrvhe_aa.f index 38ebca59..4e4f73bb 100644 --- a/TESTING/LIN/cdrvhe_aa.f +++ b/TESTING/LIN/cdrvhe_aa.f @@ -9,8 +9,8 @@ * =========== * * SUBROUTINE CDRVHE_AA( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, NMAX, -* A, AFAC, AINV, B, X, XACT, WORK, RWORK, IWORK, -* NOUT ) +* A, AFAC, AINV, B, X, XACT, WORK, RWORK, IWORK, +* NOUT ) * * .. Scalar Arguments .. * LOGICAL TSTERR @@ -186,9 +186,9 @@ CHARACTER DIST, FACT, TYPE, UPLO, XTYPE CHARACTER*3 MATPATH, PATH INTEGER I, I1, I2, IFACT, IMAT, IN, INFO, IOFF, IUPLO, - $ IZERO, J, K, K1, KL, KU, LDA, LWORK, MODE, N, + $ IZERO, J, K, KL, KU, LDA, LWORK, MODE, N, $ NB, NBMIN, NERRS, NFAIL, NIMAT, NRUN, NT - REAL AINVNM, ANORM, CNDNUM, RCOND, RCONDC + REAL ANORM, CNDNUM * .. * .. Local Arrays .. CHARACTER FACTS( NFACT ), UPLOS( 2 ) @@ -385,45 +385,6 @@ * FACT = FACTS( IFACT ) * -* Compute the condition number for comparison with -* the value returned by CHESVX. -* - IF( ZEROT ) THEN - IF( IFACT.EQ.1 ) - $ GO TO 150 - RCONDC = ZERO -* - ELSE IF( IFACT.EQ.1 ) THEN -* -* Compute the 1-norm of A. -* - ANORM = CLANHE( '1', UPLO, N, A, LDA, RWORK ) -* -* Factor the matrix A. -* -c CALL CLACPY( UPLO, N, N, A, LDA, AFAC, LDA ) -c SRNAMT = 'CHETRF_AA' -c CALL CHETRF_AA( UPLO, N, AFAC, LDA, IWORK, -c $ WORK, LWORK, INFO ) -* -* Compute inv(A) and take its norm. -* -c CALL CLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA ) -c LWORK = (N+NB+1)*(NB+3) -c SRNAMT = 'CHETRI2' -c CALL CHETRI2( UPLO, N, AINV, LDA, IWORK, WORK, -c $ LWORK, INFO ) -c AINVNM = CLANHE( '1', UPLO, N, AINV, LDA, RWORK ) -* -* Compute the 1-norm condition number of A. -* -c IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN -c RCONDC = ONE -c ELSE -c RCONDC = ( ONE / ANORM ) / AINVNM -c END IF - END IF -* * Form an exact solution and set the right hand side. * SRNAMT = 'CLARHS' @@ -487,12 +448,7 @@ c END IF CALL CLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA ) CALL CPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK, $ LDA, RWORK, RESULT( 2 ) ) -* -* Check solution from generated exact solution. -* - CALL CGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, - $ RESULT( 3 ) ) - NT = 3 + NT = 2 * * Print information about the tests that did not pass * the threshold. diff --git a/TESTING/LIN/cdrvsy_aa.f b/TESTING/LIN/cdrvsy_aa.f new file mode 100644 index 00000000..69a4e556 --- /dev/null +++ b/TESTING/LIN/cdrvsy_aa.f @@ -0,0 +1,480 @@ +*> \brief \b CDRVSY_AA +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE CDRVSY_AA( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, NMAX, +* A, AFAC, AINV, B, X, XACT, WORK, RWORK, IWORK, +* NOUT ) +* +* .. Scalar Arguments .. +* LOGICAL TSTERR +* INTEGER NMAX, NN, NOUT, NRHS +* REAL THRESH +* .. +* .. Array Arguments .. +* LOGICAL DOTYPE( * ) +* INTEGER IWORK( * ), NVAL( * ) +* REAL RWORK( * ) +* COMPLEX A( * ), AFAC( * ), AINV( * ), B( * ), +* $ WORK( * ), X( * ), XACT( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> CDRVSY_AA tests the driver routine CSYSV_AA. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] DOTYPE +*> \verbatim +*> DOTYPE is LOGICAL array, dimension (NTYPES) +*> The matrix types to be used for testing. Matrices of type j +*> (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) = +*> .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used. +*> \endverbatim +*> +*> \param[in] NN +*> \verbatim +*> NN is INTEGER +*> The number of values of N contained in the vector NVAL. +*> \endverbatim +*> +*> \param[in] NVAL +*> \verbatim +*> NVAL is INTEGER array, dimension (NN) +*> The values of the matrix dimension N. +*> \endverbatim +*> +*> \param[in] NRHS +*> \verbatim +*> NRHS is INTEGER +*> The number of right hand side vectors to be generated for +*> each linear system. +*> \endverbatim +*> +*> \param[in] THRESH +*> \verbatim +*> THRESH is REAL +*> The threshold value for the test ratios. A result is +*> included in the output file if RESULT >= THRESH. To have +*> every test ratio printed, use THRESH = 0. +*> \endverbatim +*> +*> \param[in] TSTERR +*> \verbatim +*> TSTERR is LOGICAL +*> Flag that indicates whether error exits are to be tested. +*> \endverbatim +*> +*> \param[in] NMAX +*> \verbatim +*> NMAX is INTEGER +*> The maximum value permitted for N, used in dimensioning the +*> work arrays. +*> \endverbatim +*> +*> \param[out] A +*> \verbatim +*> A is REAL array, dimension (NMAX*NMAX) +*> \endverbatim +*> +*> \param[out] AFAC +*> \verbatim +*> AFAC is REAL array, dimension (NMAX*NMAX) +*> \endverbatim +*> +*> \param[out] AINV +*> \verbatim +*> AINV is REAL array, dimension (NMAX*NMAX) +*> \endverbatim +*> +*> \param[out] B +*> \verbatim +*> B is REAL array, dimension (NMAX*NRHS) +*> \endverbatim +*> +*> \param[out] X +*> \verbatim +*> X is REAL array, dimension (NMAX*NRHS) +*> \endverbatim +*> +*> \param[out] XACT +*> \verbatim +*> XACT is REAL array, dimension (NMAX*NRHS) +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is REAL array, dimension (NMAX*max(2,NRHS)) +*> \endverbatim +*> +*> \param[out] RWORK +*> \verbatim +*> RWORK is REAL array, dimension (NMAX+2*NRHS) +*> \endverbatim +*> +*> \param[out] IWORK +*> \verbatim +*> IWORK is INTEGER array, dimension (2*NMAX) +*> \endverbatim +*> +*> \param[in] NOUT +*> \verbatim +*> NOUT is INTEGER +*> The unit number for output. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2016 +* +* @generated from LIN/ddrvsy_aa.f, fortran d -> c, Thu Nov 17 12:14:51 2016 +* +*> \ingroup complex_lin +* +* ===================================================================== + SUBROUTINE CDRVSY_AA( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, + $ NMAX, A, AFAC, AINV, B, X, XACT, WORK, + $ RWORK, IWORK, NOUT ) +* +* -- LAPACK test routine (version 3.7.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2016 +* +* .. Scalar Arguments .. + LOGICAL TSTERR + INTEGER NMAX, NN, NOUT, NRHS + REAL THRESH +* .. +* .. Array Arguments .. + LOGICAL DOTYPE( * ) + INTEGER IWORK( * ), NVAL( * ) + REAL RWORK( * ) + COMPLEX A( * ), AFAC( * ), AINV( * ), B( * ), + $ WORK( * ), X( * ), XACT( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + REAL ZERO + PARAMETER ( ZERO = 0.0D+0 ) + COMPLEX CZERO + PARAMETER ( CZERO = 0.0E+0 ) + INTEGER NTYPES, NTESTS + PARAMETER ( NTYPES = 10, NTESTS = 3 ) + INTEGER NFACT + PARAMETER ( NFACT = 2 ) +* .. +* .. Local Scalars .. + LOGICAL ZEROT + CHARACTER DIST, FACT, TYPE, UPLO, XTYPE + CHARACTER*3 MATPATH, PATH + INTEGER I, I1, I2, IFACT, IMAT, IN, INFO, IOFF, IUPLO, + $ IZERO, J, K, KL, KU, LDA, LWORK, MODE, N, + $ NB, NBMIN, NERRS, NFAIL, NIMAT, NRUN, NT + REAL ANORM, CNDNUM +* .. +* .. Local Arrays .. + CHARACTER FACTS( NFACT ), UPLOS( 2 ) + INTEGER ISEED( 4 ), ISEEDY( 4 ) + REAL RESULT( NTESTS ) +* .. +* .. External Functions .. + REAL DGET06, CLANSY + EXTERNAL DGET06, CLANSY +* .. +* .. External Subroutines .. + EXTERNAL ALADHD, ALAERH, ALASVM, DERRVX, CGET04, CLACPY, + $ CLARHS, CLASET, CLATB4, CLATMS, CSYT02, DSYT05, + $ CSYSV_AA, CSYT01_AA, CSYTRF_AA, XLAENV +* .. +* .. Scalars in Common .. + LOGICAL LERR, OK + CHARACTER*32 SRNAMT + INTEGER INFOT, NUNIT +* .. +* .. Common blocks .. + COMMON / INFOC / INFOT, NUNIT, OK, LERR + COMMON / SRNAMC / SRNAMT +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. Data statements .. + DATA ISEEDY / 1988, 1989, 1990, 1991 / + DATA UPLOS / 'U', 'L' / , FACTS / 'F', 'N' / +* .. +* .. Executable Statements .. +* +* Initialize constants and the random number seed. +* +* Test path +* + PATH( 1: 1 ) = 'Complex precision' + PATH( 2: 3 ) = 'SA' +* +* Path to generate matrices +* + MATPATH( 1: 1 ) = 'Complex precision' + MATPATH( 2: 3 ) = 'SY' +* + NRUN = 0 + NFAIL = 0 + NERRS = 0 + DO 10 I = 1, 4 + ISEED( I ) = ISEEDY( I ) + 10 CONTINUE + LWORK = MAX( 2*NMAX, NMAX*NRHS ) +* +* Test the error exits +* + IF( TSTERR ) + $ CALL CERRVX( PATH, NOUT ) + INFOT = 0 +* +* Set the block size and minimum block size for testing. +* + NB = 1 + NBMIN = 2 + CALL XLAENV( 1, NB ) + CALL XLAENV( 2, NBMIN ) +* +* Do for each value of N in NVAL +* + DO 180 IN = 1, NN + N = NVAL( IN ) + LDA = MAX( N, 1 ) + XTYPE = 'N' + NIMAT = NTYPES + IF( N.LE.0 ) + $ NIMAT = 1 +* + DO 170 IMAT = 1, NIMAT +* +* Do the tests only if DOTYPE( IMAT ) is true. +* + IF( .NOT.DOTYPE( IMAT ) ) + $ GO TO 170 +* +* Skip types 3, 4, 5, or 6 if the matrix size is too small. +* + ZEROT = IMAT.GE.3 .AND. IMAT.LE.6 + IF( ZEROT .AND. N.LT.IMAT-2 ) + $ GO TO 170 +* +* Do first for UPLO = 'U', then for UPLO = 'L' +* + DO 160 IUPLO = 1, 2 + UPLO = UPLOS( IUPLO ) +* +* Set up parameters with CLATB4 and generate a test matrix +* with CLATMS. +* + CALL CLATB4( MATPATH, IMAT, N, N, TYPE, KL, KU, ANORM, + $ MODE, CNDNUM, DIST ) +* + SRNAMT = 'CLATMS' + CALL CLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE, + $ CNDNUM, ANORM, KL, KU, UPLO, A, LDA, WORK, + $ INFO ) +* +* Check error code from CLATMS. +* + IF( INFO.NE.0 ) THEN + CALL ALAERH( PATH, 'CLATMS', INFO, 0, UPLO, N, N, -1, + $ -1, -1, IMAT, NFAIL, NERRS, NOUT ) + GO TO 160 + END IF +* +* For types 3-6, zero one or more rows and columns of the +* matrix to test that INFO is returned correctly. +* + IF( ZEROT ) THEN + IF( IMAT.EQ.3 ) THEN + IZERO = 1 + ELSE IF( IMAT.EQ.4 ) THEN + IZERO = N + ELSE + IZERO = N / 2 + 1 + END IF +* + IF( IMAT.LT.6 ) THEN +* +* Set row and column IZERO to zero. +* + IF( IUPLO.EQ.1 ) THEN + IOFF = ( IZERO-1 )*LDA + DO 20 I = 1, IZERO - 1 + A( IOFF+I ) = CZERO + 20 CONTINUE + IOFF = IOFF + IZERO + DO 30 I = IZERO, N + A( IOFF ) = CZERO + IOFF = IOFF + LDA + 30 CONTINUE + ELSE + IOFF = IZERO + DO 40 I = 1, IZERO - 1 + A( IOFF ) = CZERO + IOFF = IOFF + LDA + 40 CONTINUE + IOFF = IOFF - IZERO + DO 50 I = IZERO, N + A( IOFF+I ) = CZERO + 50 CONTINUE + END IF + ELSE + IOFF = 0 + IF( IUPLO.EQ.1 ) THEN +* +* Set the first IZERO rows and columns to zero. +* + DO 70 J = 1, N + I2 = MIN( J, IZERO ) + DO 60 I = 1, I2 + A( IOFF+I ) = CZERO + 60 CONTINUE + IOFF = IOFF + LDA + 70 CONTINUE + IZERO = 1 + ELSE +* +* Set the last IZERO rows and columns to zero. +* + DO 90 J = 1, N + I1 = MAX( J, IZERO ) + DO 80 I = I1, N + A( IOFF+I ) = CZERO + 80 CONTINUE + IOFF = IOFF + LDA + 90 CONTINUE + END IF + END IF + ELSE + IZERO = 0 + END IF +* + DO 150 IFACT = 1, NFACT +* +* Do first for FACT = 'F', then for other values. +* + FACT = FACTS( IFACT ) +* +* Form an exact solution and set the right hand side. +* + SRNAMT = 'CLARHS' + CALL CLARHS( MATPATH, XTYPE, UPLO, ' ', N, N, KL, KU, + $ NRHS, A, LDA, XACT, LDA, B, LDA, ISEED, + $ INFO ) + XTYPE = 'C' +* +* --- Test CSYSV_AA --- +* + IF( IFACT.EQ.2 ) THEN + CALL CLACPY( UPLO, N, N, A, LDA, AFAC, LDA ) + CALL CLACPY( 'Full', N, NRHS, B, LDA, X, LDA ) +* +* Factor the matrix and solve the system using CSYSV_AA. +* + SRNAMT = 'CSYSV_AA' + CALL CSYSV_AA( UPLO, N, NRHS, AFAC, LDA, IWORK, + $ X, LDA, WORK, LWORK, INFO ) +* +* Adjust the expected value of INFO to account for +* pivoting. +* + IF( IZERO.GT.0 ) THEN + J = 1 + K = IZERO + 100 CONTINUE + IF( J.EQ.K ) THEN + K = IWORK( J ) + ELSE IF( IWORK( J ).EQ.K ) THEN + K = J + END IF + IF( J.LT.K ) THEN + J = J + 1 + GO TO 100 + END IF + ELSE + K = 0 + END IF +* +* Check error code from CSYSV_AA . +* + IF( INFO.NE.K ) THEN + CALL ALAERH( PATH, 'CSYSV_AA ', INFO, K, + $ UPLO, N, N, -1, -1, NRHS, + $ IMAT, NFAIL, NERRS, NOUT ) + GO TO 120 + ELSE IF( INFO.NE.0 ) THEN + GO TO 120 + END IF +* +* Reconstruct matrix from factors and compute +* residual. +* + CALL CSYT01_AA( UPLO, N, A, LDA, AFAC, LDA, + $ IWORK, AINV, LDA, RWORK, + $ RESULT( 1 ) ) +* +* Compute residual of the computed solution. +* + CALL CLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA ) + CALL CSYT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK, + $ LDA, RWORK, RESULT( 2 ) ) + NT = 2 +* +* Print information about the tests that did not pass +* the threshold. +* + DO 110 K = 1, NT + IF( RESULT( K ).GE.THRESH ) THEN + IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) + $ CALL ALADHD( NOUT, PATH ) + WRITE( NOUT, FMT = 9999 )'CSYSV_AA ', + $ UPLO, N, IMAT, K, RESULT( K ) + NFAIL = NFAIL + 1 + END IF + 110 CONTINUE + NRUN = NRUN + NT + 120 CONTINUE + END IF +* + 150 CONTINUE +* + 160 CONTINUE + 170 CONTINUE + 180 CONTINUE +* +* Print a summary of the results. +* + CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS ) +* + 9999 FORMAT( 1X, A, ', UPLO=''', A1, ''', N =', I5, ', type ', I2, + $ ', test ', I2, ', ratio =', G12.5 ) + RETURN +* +* End of CDRVSY_AA +* + END diff --git a/TESTING/LIN/cerrhe.f b/TESTING/LIN/cerrhe.f index 3711b8e3..2bc50c0d 100644 --- a/TESTING/LIN/cerrhe.f +++ b/TESTING/LIN/cerrhe.f @@ -488,6 +488,12 @@ INFOT = 4 CALL CHETRF_AA( 'U', 2, A, 1, IP, W, 4, INFO ) CALL CHKXER( 'CHETRF_AA', INFOT, NOUT, LERR, OK ) + INFOT = 7 + CALL CHETRF_AA( 'U', 0, A, 1, IP, W, 0, INFO ) + CALL CHKXER( 'CHETRF_AA', INFOT, NOUT, LERR, OK ) + INFOT = 7 + CALL CHETRF_AA( 'U', 0, A, 1, IP, W, -2, INFO ) + CALL CHKXER( 'CHETRF_AA', INFOT, NOUT, LERR, OK ) * * CHETRS_AA * @@ -507,6 +513,12 @@ INFOT = 8 CALL CHETRS_AA( 'U', 2, 1, A, 2, IP, B, 1, W, 1, INFO ) CALL CHKXER( 'CHETRS_AA', INFOT, NOUT, LERR, OK ) + INFOT = 10 + CALL CHETRS_AA( 'U', 0, 1, A, 1, IP, B, 1, W, 0, INFO ) + CALL CHKXER( 'CHETRS_AA', INFOT, NOUT, LERR, OK ) + INFOT = 10 + CALL CHETRS_AA( 'U', 0, 1, A, 1, IP, B, 1, W, -2, INFO ) + CALL CHKXER( 'CHETRS_AA', INFOT, NOUT, LERR, OK ) * * Test error exits of the routines that use factorization * of a Hermitian indefinite packed matrix with patrial diff --git a/TESTING/LIN/cerrsy.f b/TESTING/LIN/cerrsy.f index c7613bd6..e4bdc1dd 100644 --- a/TESTING/LIN/cerrsy.f +++ b/TESTING/LIN/cerrsy.f @@ -547,6 +547,56 @@ INFOT = 5 CALL CSPCON( 'U', 1, A, IP, -ANRM, RCOND, W, INFO ) CALL CHKXER( 'CSPCON', INFOT, NOUT, LERR, OK ) +* + ELSE IF( LSAMEN( 2, C2, 'SA' ) ) THEN +* +* Test error exits of the routines that use factorization +* of a symmetric indefinite matrix with Aasen's algorithm +* +* CSYTRF_AA +* + SRNAMT = 'CSYTRF_AA' + INFOT = 1 + CALL CSYTRF_AA( '/', 0, A, 1, IP, W, 1, INFO ) + CALL CHKXER( 'CSYTRF_AA', INFOT, NOUT, LERR, OK ) + INFOT = 2 + CALL CSYTRF_AA( 'U', -1, A, 1, IP, W, 1, INFO ) + CALL CHKXER( 'CSYTRF_AA', INFOT, NOUT, LERR, OK ) + INFOT = 4 + CALL CSYTRF_AA( 'U', 2, A, 1, IP, W, 4, INFO ) + CALL CHKXER( 'CSYTRF_AA', INFOT, NOUT, LERR, OK ) + INFOT = 7 + CALL CSYTRF_AA( 'U', 0, A, 1, IP, W, 0, INFO ) + CALL CHKXER( 'CSYTRF_AA', INFOT, NOUT, LERR, OK ) + INFOT = 7 + CALL CSYTRF_AA( 'U', 0, A, 1, IP, W, -2, INFO ) + CALL CHKXER( 'CSYTRF_AA', INFOT, NOUT, LERR, OK ) +* +* CSYTRS_AA +* + SRNAMT = 'CSYTRS_AA' + INFOT = 1 + CALL CSYTRS_AA( '/', 0, 0, A, 1, IP, B, 1, W, 1, INFO ) + CALL CHKXER( 'CSYTRS_AA', INFOT, NOUT, LERR, OK ) + INFOT = 2 + CALL CSYTRS_AA( 'U', -1, 0, A, 1, IP, B, 1, W, 1, INFO ) + CALL CHKXER( 'CSYTRS_AA', INFOT, NOUT, LERR, OK ) + INFOT = 3 + CALL CSYTRS_AA( 'U', 0, -1, A, 1, IP, B, 1, W, 1, INFO ) + CALL CHKXER( 'CSYTRS_AA', INFOT, NOUT, LERR, OK ) + INFOT = 5 + CALL CSYTRS_AA( 'U', 2, 1, A, 1, IP, B, 2, W, 1, INFO ) + CALL CHKXER( 'CSYTRS_AA', INFOT, NOUT, LERR, OK ) + INFOT = 8 + CALL CSYTRS_AA( 'U', 2, 1, A, 2, IP, B, 1, W, 1, INFO ) + CALL CHKXER( 'CSYTRS_AA', INFOT, NOUT, LERR, OK ) + INFOT = 10 + CALL CSYTRS_AA( 'U', 0, 1, A, 1, IP, B, 1, W, 0, INFO ) + CALL CHKXER( 'CSYTRS_AA', INFOT, NOUT, LERR, OK ) + INFOT = 10 + CALL CSYTRS_AA( 'U', 0, 1, A, 1, IP, B, 1, W, -2, INFO ) + CALL CHKXER( 'CSYTRS_AA', INFOT, NOUT, LERR, OK ) +* END IF * * Print a summary line. diff --git a/TESTING/LIN/chet01_aa.f b/TESTING/LIN/chet01_aa.f index 8f797f11..31b504d2 100644 --- a/TESTING/LIN/chet01_aa.f +++ b/TESTING/LIN/chet01_aa.f @@ -9,17 +9,17 @@ * =========== * * SUBROUTINE CHET01_AA( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, -* C, LDC, RWORK, RESID ) +* C, LDC, RWORK, RESID ) * * .. Scalar Arguments .. * CHARACTER UPLO * INTEGER LDA, LDAFAC, LDC, N -* COMPLEX RESID +* REAL RESID * .. * .. Array Arguments .. * INTEGER IPIV( * ) -* COMPLEX A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ), -* $ RWORK( * ) +* REAL RWORK( * ) +* COMPLEX A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ) * .. * * @@ -123,7 +123,7 @@ * * ===================================================================== SUBROUTINE CHET01_AA( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, C, - $ LDC, RWORK, RESID ) + $ LDC, RWORK, RESID ) * * -- LAPACK test routine (version 3.7.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -137,8 +137,8 @@ * .. * .. Array Arguments .. INTEGER IPIV( * ) - COMPLEX A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ), - $ RWORK( * ) + REAL RWORK( * ) + COMPLEX A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ) * .. * * ===================================================================== @@ -197,27 +197,29 @@ $ LDC+1 ) CALL CLACGV( N-1, C( 1, 2 ), LDC+1 ) ENDIF - ENDIF * -* Call CTRMM to form the product U' * D (or L * D ). +* Call CTRMM to form the product U' * D (or L * D ). * - IF( LSAME( UPLO, 'U' ) ) THEN - CALL CTRMM( 'Left', UPLO, 'Conjugate transpose', 'Unit', N-1, - $ N, CONE, AFAC( 1, 2 ), LDAFAC, C( 2, 1 ), LDC ) - ELSE - CALL CTRMM( 'Left', UPLO, 'No transpose', 'Unit', N-1, N, - $ CONE, AFAC( 2, 1 ), LDAFAC, C( 2, 1 ), LDC ) - END IF + IF( LSAME( UPLO, 'U' ) ) THEN + CALL CTRMM( 'Left', UPLO, 'Conjugate transpose', 'Unit', + $ N-1, N, CONE, AFAC( 1, 2 ), LDAFAC, C( 2, 1 ), + $ LDC ) + ELSE + CALL CTRMM( 'Left', UPLO, 'No transpose', 'Unit', N-1, N, + $ CONE, AFAC( 2, 1 ), LDAFAC, C( 2, 1 ), LDC ) + END IF * -* Call CTRMM again to multiply by U (or L ). +* Call CTRMM again to multiply by U (or L ). * - IF( LSAME( UPLO, 'U' ) ) THEN - CALL CTRMM( 'Right', UPLO, 'No transpose', 'Unit', N, N-1, - $ CONE, AFAC( 1, 2 ), LDAFAC, C( 1, 2 ), LDC ) - ELSE - CALL CTRMM( 'Right', UPLO, 'Conjugate transpose', 'Unit', N, - $ N-1, CONE, AFAC( 2, 1 ), LDAFAC, C( 1, 2 ), LDC ) - END IF + IF( LSAME( UPLO, 'U' ) ) THEN + CALL CTRMM( 'Right', UPLO, 'No transpose', 'Unit', N, N-1, + $ CONE, AFAC( 1, 2 ), LDAFAC, C( 1, 2 ), LDC ) + ELSE + CALL CTRMM( 'Right', UPLO, 'Conjugate transpose', 'Unit', N, + $ N-1, CONE, AFAC( 2, 1 ), LDAFAC, C( 1, 2 ), + $ LDC ) + END IF + ENDIF * * Apply hermitian pivots * diff --git a/TESTING/LIN/csyt01_aa.f b/TESTING/LIN/csyt01_aa.f new file mode 100644 index 00000000..7c7382a3 --- /dev/null +++ b/TESTING/LIN/csyt01_aa.f @@ -0,0 +1,265 @@ +*> \brief \b CSYT01 +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE CSYT01_AA( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, C, LDC, +* RWORK, RESID ) +* +* .. Scalar Arguments .. +* CHARACTER UPLO +* INTEGER LDA, LDAFAC, LDC, N +* REAL RESID +* .. +* .. Array Arguments .. +* INTEGER IPIV( * ) +* REAL RWORK( * ) +* COMPLEX A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> CSYT01 reconstructs a hermitian indefinite matrix A from its +*> block L*D*L' or U*D*U' factorization and computes the residual +*> norm( C - A ) / ( N * norm(A) * EPS ), +*> where C is the reconstructed matrix and EPS is the machine epsilon. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] UPLO +*> \verbatim +*> UPLO is CHARACTER*1 +*> Specifies whether the upper or lower triangular part of the +*> hermitian matrix A is stored: +*> = 'U': Upper triangular +*> = 'L': Lower triangular +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of rows and columns of the matrix A. N >= 0. +*> \endverbatim +*> +*> \param[in] A +*> \verbatim +*> A is REAL array, dimension (LDA,N) +*> The original hermitian matrix A. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N) +*> \endverbatim +*> +*> \param[in] AFAC +*> \verbatim +*> AFAC is REAL array, dimension (LDAFAC,N) +*> The factored form of the matrix A. AFAC contains the block +*> diagonal matrix D and the multipliers used to obtain the +*> factor L or U from the block L*D*L' or U*D*U' factorization +*> as computed by CSYTRF. +*> \endverbatim +*> +*> \param[in] LDAFAC +*> \verbatim +*> LDAFAC is INTEGER +*> The leading dimension of the array AFAC. LDAFAC >= max(1,N). +*> \endverbatim +*> +*> \param[in] IPIV +*> \verbatim +*> IPIV is INTEGER array, dimension (N) +*> The pivot indices from CSYTRF. +*> \endverbatim +*> +*> \param[out] C +*> \verbatim +*> C is REAL array, dimension (LDC,N) +*> \endverbatim +*> +*> \param[in] LDC +*> \verbatim +*> LDC is INTEGER +*> The leading dimension of the array C. LDC >= max(1,N). +*> \endverbatim +*> +*> \param[out] RWORK +*> \verbatim +*> RWORK is REAL array, dimension (N) +*> \endverbatim +*> +*> \param[out] RESID +*> \verbatim +*> RESID is REAL +*> If UPLO = 'L', norm(L*D*L' - A) / ( N * norm(A) * EPS ) +*> If UPLO = 'U', norm(U*D*U' - A) / ( N * norm(A) * EPS ) +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2016 +* +* @generated from LIN/dsyt01_aa.f, fortran d -> c, Thu Nov 17 13:01:50 2016 +* +*> \ingroup complex_lin +* +* ===================================================================== + SUBROUTINE CSYT01_AA( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, C, + $ LDC, RWORK, RESID ) +* +* -- LAPACK test routine (version 3.5.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2016 +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER LDA, LDAFAC, LDC, N + REAL RESID +* .. +* .. Array Arguments .. + INTEGER IPIV( * ) + COMPLEX A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ) + REAL RWORK( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + REAL ZERO, ONE + PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) + COMPLEX CZERO, CONE + PARAMETER ( CZERO = 0.0E+0, CONE = 1.0E+0 ) +* .. +* .. Local Scalars .. + INTEGER I, J + REAL ANORM, EPS +* .. +* .. External Functions .. + LOGICAL LSAME + REAL SLAMCH, CLANSY + EXTERNAL LSAME, SLAMCH, CLANSY +* .. +* .. External Subroutines .. + EXTERNAL CLASET, CLAVSY +* .. +* .. Intrinsic Functions .. + INTRINSIC DBLE +* .. +* .. Executable Statements .. +* +* Quick exit if N = 0. +* + IF( N.LE.0 ) THEN + RESID = ZERO + RETURN + END IF +* +* Determine EPS and the norm of A. +* + EPS = SLAMCH( 'Epsilon' ) + ANORM = CLANSY( '1', UPLO, N, A, LDA, RWORK ) +* +* Initialize C to the tridiagonal matrix T. +* + CALL CLASET( 'Full', N, N, CZERO, CZERO, C, LDC ) + CALL CLACPY( 'F', 1, N, AFAC( 1, 1 ), LDAFAC+1, C( 1, 1 ), LDC+1 ) + IF( N.GT.1 ) THEN + IF( LSAME( UPLO, 'U' ) ) THEN + CALL CLACPY( 'F', 1, N-1, AFAC( 1, 2 ), LDAFAC+1, C( 1, 2 ), + $ LDC+1 ) + CALL CLACPY( 'F', 1, N-1, AFAC( 1, 2 ), LDAFAC+1, C( 2, 1 ), + $ LDC+1 ) + ELSE + CALL CLACPY( 'F', 1, N-1, AFAC( 2, 1 ), LDAFAC+1, C( 1, 2 ), + $ LDC+1 ) + CALL CLACPY( 'F', 1, N-1, AFAC( 2, 1 ), LDAFAC+1, C( 2, 1 ), + $ LDC+1 ) + ENDIF +* +* Call CTRMM to form the product U' * D (or L * D ). +* + IF( LSAME( UPLO, 'U' ) ) THEN + CALL CTRMM( 'Left', UPLO, 'Transpose', 'Unit', N-1, N, + $ CONE, AFAC( 1, 2 ), LDAFAC, C( 2, 1 ), LDC ) + ELSE + CALL CTRMM( 'Left', UPLO, 'No transpose', 'Unit', N-1, N, + $ CONE, AFAC( 2, 1 ), LDAFAC, C( 2, 1 ), LDC ) + END IF +* +* Call CTRMM again to multiply by U (or L ). +* + IF( LSAME( UPLO, 'U' ) ) THEN + CALL CTRMM( 'Right', UPLO, 'No transpose', 'Unit', N, N-1, + $ CONE, AFAC( 1, 2 ), LDAFAC, C( 1, 2 ), LDC ) + ELSE + CALL CTRMM( 'Right', UPLO, 'Transpose', 'Unit', N, N-1, + $ CONE, AFAC( 2, 1 ), LDAFAC, C( 1, 2 ), LDC ) + END IF + ENDIF +* +* Apply symmetric pivots +* + DO J = N, 1, -1 + I = IPIV( J ) + IF( I.NE.J ) + $ CALL CSWAP( N, C( J, 1 ), LDC, C( I, 1 ), LDC ) + END DO + DO J = N, 1, -1 + I = IPIV( J ) + IF( I.NE.J ) + $ CALL CSWAP( N, C( 1, J ), 1, C( 1, I ), 1 ) + END DO +* +* +* Compute the difference C - A . +* + IF( LSAME( UPLO, 'U' ) ) THEN + DO J = 1, N + DO I = 1, J + C( I, J ) = C( I, J ) - A( I, J ) + END DO + END DO + ELSE + DO J = 1, N + DO I = J, N + C( I, J ) = C( I, J ) - A( I, J ) + END DO + END DO + END IF +* +* Compute norm( C - A ) / ( N * norm(A) * EPS ) +* + RESID = CLANSY( '1', UPLO, N, C, LDC, RWORK ) +* + IF( ANORM.LE.ZERO ) THEN + IF( RESID.NE.ZERO ) + $ RESID = ONE / EPS + ELSE + RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS + END IF +* + RETURN +* +* End of CSYT01 +* + END diff --git a/TESTING/LIN/dchksy_aa.f b/TESTING/LIN/dchksy_aa.f index a596bcb9..a01f4fa0 100644 --- a/TESTING/LIN/dchksy_aa.f +++ b/TESTING/LIN/dchksy_aa.f @@ -163,6 +163,7 @@ * *> \date November 2016 * +* @precisions fortran d -> z c * *> \ingroup double_lin * @@ -201,13 +202,13 @@ PARAMETER ( NTESTS = 9 ) * .. * .. Local Scalars .. - LOGICAL TRFCON, ZEROT + LOGICAL ZEROT CHARACTER DIST, TYPE, UPLO, XTYPE CHARACTER*3 PATH, MATPATH INTEGER I, I1, I2, IMAT, IN, INB, INFO, IOFF, IRHS, $ IUPLO, IZERO, J, K, KL, KU, LDA, LWORK, MODE, $ N, NB, NERRS, NFAIL, NIMAT, NRHS, NRUN, NT - DOUBLE PRECISION ANORM, CNDNUM, RCONDC + DOUBLE PRECISION ANORM, CNDNUM * .. * .. Local Arrays .. CHARACTER UPLOS( 2 ) @@ -221,7 +222,7 @@ * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, DERRSY, DGET04, DLACPY, $ DLARHS, DLATB4, DLATMS, DPOT02, DPOT03, DPOT05, - $ DSYCON, DSYRFS, DSYT01, DSYTRF_AA, + $ DSYCON, DSYRFS, DSYT01_AA, DSYTRF_AA, $ DSYTRI2, DSYTRS_AA, XLAENV * .. * .. Intrinsic Functions .. @@ -429,9 +430,9 @@ * block factorization, LWORK is the length of AINV. * SRNAMT = 'DSYTRF_AA' - LWORK = N*NB + N + LWORK = MAX( 1, N*NB + N ) CALL DSYTRF_AA( UPLO, N, AFAC, LDA, IWORK, AINV, - $ LWORK, INFO ) + $ LWORK, INFO ) * * Adjust the expected value of INFO to account for * pivoting. @@ -461,19 +462,11 @@ $ NOUT ) END IF * -* Set the condition estimate flag if the INFO is not 0. -* - IF( INFO.NE.0 ) THEN - TRFCON = .TRUE. - ELSE - TRFCON = .FALSE. - END IF -* *+ TEST 1 * Reconstruct matrix from factors and compute residual. * CALL DSYT01_AA( UPLO, N, A, LDA, AFAC, LDA, IWORK, - $ AINV, LDA, RWORK, RESULT( 1 ) ) + $ AINV, LDA, RWORK, RESULT( 1 ) ) NT = 1 * * @@ -491,10 +484,9 @@ 110 CONTINUE NRUN = NRUN + NT * -* Do only the condition estimate if INFO is not 0. +* Skip solver test if INFO is not 0. * - IF( TRFCON ) THEN - RCONDC = ZERO + IF( INFO.NE.0 ) THEN GO TO 140 END IF * @@ -503,7 +495,7 @@ DO 130 IRHS = 1, NNS NRHS = NSVAL( IRHS ) * -*+ TEST 3 ( Using TRS) +*+ TEST 2 (Using TRS) * Solve and compute residual for A * X = B. * * Choose a set of NRHS random solution vectors @@ -516,10 +508,10 @@ CALL DLACPY( 'Full', N, NRHS, B, LDA, X, LDA ) * SRNAMT = 'DSYTRS_AA' - LWORK = 3*N-2 + LWORK = MAX( 1, 3*N-2 ) CALL DSYTRS_AA( UPLO, N, NRHS, AFAC, LDA, - $ IWORK, X, LDA, WORK, LWORK, - $ INFO ) + $ IWORK, X, LDA, WORK, LWORK, + $ INFO ) * * Check error code from DSYTRS and handle error. * diff --git a/TESTING/LIN/ddrvsy_aa.f b/TESTING/LIN/ddrvsy_aa.f index be5d6eb3..25a58292 100644 --- a/TESTING/LIN/ddrvsy_aa.f +++ b/TESTING/LIN/ddrvsy_aa.f @@ -9,8 +9,8 @@ * =========== * * SUBROUTINE DDRVSY_AA( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, NMAX, -* A, AFAC, AINV, B, X, XACT, WORK, RWORK, IWORK, -* NOUT ) +* A, AFAC, AINV, B, X, XACT, WORK, RWORK, IWORK, +* NOUT ) * * .. Scalar Arguments .. * LOGICAL TSTERR @@ -145,12 +145,14 @@ * *> \date November 2016 * +* @precisions fortran d -> z c +* *> \ingroup double_lin * * ===================================================================== SUBROUTINE DDRVSY_AA( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, - $ NMAX, A, AFAC, AINV, B, X, XACT, WORK, - $ RWORK, IWORK, NOUT ) + $ NMAX, A, AFAC, AINV, B, X, XACT, WORK, + $ RWORK, IWORK, NOUT ) * * -- LAPACK test routine (version 3.7.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -184,9 +186,9 @@ CHARACTER DIST, FACT, TYPE, UPLO, XTYPE CHARACTER*3 MATPATH, PATH INTEGER I, I1, I2, IFACT, IMAT, IN, INFO, IOFF, IUPLO, - $ IZERO, J, K, K1, KL, KU, LDA, LWORK, MODE, N, + $ IZERO, J, K, KL, KU, LDA, LWORK, MODE, N, $ NB, NBMIN, NERRS, NFAIL, NIMAT, NRUN, NT - DOUBLE PRECISION AINVNM, ANORM, CNDNUM, RCOND, RCONDC + DOUBLE PRECISION ANORM, CNDNUM * .. * .. Local Arrays .. CHARACTER FACTS( NFACT ), UPLOS( 2 ) @@ -374,44 +376,6 @@ * FACT = FACTS( IFACT ) * -* Compute the condition number for comparison with -* the value returned by DSYSVX. -* - IF( ZEROT ) THEN - IF( IFACT.EQ.1 ) - $ GO TO 150 - RCONDC = ZERO -* - ELSE IF( IFACT.EQ.1 ) THEN -* -* Compute the 1-norm of A. -* - ANORM = DLANSY( '1', UPLO, N, A, LDA, RWORK ) -* -* Factor the matrix A. -* -c CALL DLACPY( UPLO, N, N, A, LDA, AFAC, LDA ) -c CALL DSYTRF( UPLO, N, AFAC, LDA, IWORK, WORK, -c $ LWORK, INFO ) -* -* Compute inv(A) and take its norm. -* -c CALL DLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA ) -c LWORK = (N+NB+1)*(NB+3) -c SRNAMT = 'DSYTRI2' -c CALL DSYTRI2( UPLO, N, AINV, LDA, IWORK, WORK, -c $ LWORK, INFO ) -c AINVNM = DLANSY( '1', UPLO, N, AINV, LDA, RWORK ) -* -* Compute the 1-norm condition number of A. -* -c IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN -c RCONDC = ONE -c ELSE -c RCONDC = ( ONE / ANORM ) / AINVNM -c END IF - END IF -* * Form an exact solution and set the right hand side. * SRNAMT = 'DLARHS' @@ -475,12 +439,7 @@ c END IF CALL DLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA ) CALL DPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK, $ LDA, RWORK, RESULT( 2 ) ) -* -* Check solution from generated exact solution. -* - CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, - $ RESULT( 3 ) ) - NT = 3 + NT = 2 * * Print information about the tests that did not pass * the threshold. diff --git a/TESTING/LIN/derrsy.f b/TESTING/LIN/derrsy.f index 056e931b..7fe74de3 100644 --- a/TESTING/LIN/derrsy.f +++ b/TESTING/LIN/derrsy.f @@ -493,6 +493,12 @@ INFOT = 4 CALL DSYTRF_AA( 'U', 2, A, 1, IP, W, 4, INFO ) CALL CHKXER( 'DSYTRF_AA', INFOT, NOUT, LERR, OK ) + INFOT = 7 + CALL DSYTRF_AA( 'U', 0, A, 1, IP, W, 0, INFO ) + CALL CHKXER( 'DSYTRF_AA', INFOT, NOUT, LERR, OK ) + INFOT = 7 + CALL DSYTRF_AA( 'U', 0, A, 1, IP, W, -2, INFO ) + CALL CHKXER( 'DSYTRF_AA', INFOT, NOUT, LERR, OK ) * * DSYTRS_AA * @@ -512,6 +518,12 @@ INFOT = 8 CALL DSYTRS_AA( 'U', 2, 1, A, 2, IP, B, 1, W, 1, INFO ) CALL CHKXER( 'DSYTRS_AA', INFOT, NOUT, LERR, OK ) + INFOT = 10 + CALL DSYTRS_AA( 'U', 0, 1, A, 2, IP, B, 1, W, 0, INFO ) + CALL CHKXER( 'DSYTRS_AA', INFOT, NOUT, LERR, OK ) + INFOT = 10 + CALL DSYTRS_AA( 'U', 0, 1, A, 2, IP, B, 1, W, -2, INFO ) + CALL CHKXER( 'DSYTRS_AA', INFOT, NOUT, LERR, OK ) * ELSE IF( LSAMEN( 2, C2, 'SP' ) ) THEN * diff --git a/TESTING/LIN/dsyt01_aa.f b/TESTING/LIN/dsyt01_aa.f index bc30df38..3a704de3 100644 --- a/TESTING/LIN/dsyt01_aa.f +++ b/TESTING/LIN/dsyt01_aa.f @@ -8,8 +8,8 @@ * Definition: * =========== * -* SUBROUTINE DSYT01( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, C, LDC, -* RWORK, RESID ) +* SUBROUTINE DSYT01_AA( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, C, LDC, +* RWORK, RESID ) * * .. Scalar Arguments .. * CHARACTER UPLO @@ -118,7 +118,7 @@ * *> \date November 2016 * -* @precisions fortran d -> s +* @precisions fortran d -> z c * *> \ingroup double_lin * @@ -193,27 +193,27 @@ CALL DLACPY( 'F', 1, N-1, AFAC( 2, 1 ), LDAFAC+1, C( 2, 1 ), $ LDC+1 ) ENDIF - ENDIF * -* Call DTRMM to form the product U' * D (or L * D ). +* Call DTRMM to form the product U' * D (or L * D ). * - IF( LSAME( UPLO, 'U' ) ) THEN - CALL DTRMM( 'Left', UPLO, 'Transpose', 'Unit', N-1, N, - $ ONE, AFAC( 1, 2 ), LDAFAC, C( 2, 1 ), LDC ) - ELSE - CALL DTRMM( 'Left', UPLO, 'No transpose', 'Unit', N-1, N, - $ ONE, AFAC( 2, 1 ), LDAFAC, C( 2, 1 ), LDC ) - END IF + IF( LSAME( UPLO, 'U' ) ) THEN + CALL DTRMM( 'Left', UPLO, 'Transpose', 'Unit', N-1, N, + $ ONE, AFAC( 1, 2 ), LDAFAC, C( 2, 1 ), LDC ) + ELSE + CALL DTRMM( 'Left', UPLO, 'No transpose', 'Unit', N-1, N, + $ ONE, AFAC( 2, 1 ), LDAFAC, C( 2, 1 ), LDC ) + END IF * -* Call DTRMM again to multiply by U (or L ). +* Call DTRMM again to multiply by U (or L ). * - IF( LSAME( UPLO, 'U' ) ) THEN - CALL DTRMM( 'Right', UPLO, 'No transpose', 'Unit', N, N-1, - $ ONE, AFAC( 1, 2 ), LDAFAC, C( 1, 2 ), LDC ) - ELSE - CALL DTRMM( 'Right', UPLO, 'Transpose', 'Unit', N, N-1, - $ ONE, AFAC( 2, 1 ), LDAFAC, C( 1, 2 ), LDC ) - END IF + IF( LSAME( UPLO, 'U' ) ) THEN + CALL DTRMM( 'Right', UPLO, 'No transpose', 'Unit', N, N-1, + $ ONE, AFAC( 1, 2 ), LDAFAC, C( 1, 2 ), LDC ) + ELSE + CALL DTRMM( 'Right', UPLO, 'Transpose', 'Unit', N, N-1, + $ ONE, AFAC( 2, 1 ), LDAFAC, C( 1, 2 ), LDC ) + END IF + ENDIF * * Apply symmetric pivots * diff --git a/TESTING/LIN/schksy_aa.f b/TESTING/LIN/schksy_aa.f index 0f668723..ee00fdb0 100644 --- a/TESTING/LIN/schksy_aa.f +++ b/TESTING/LIN/schksy_aa.f @@ -163,6 +163,7 @@ * *> \date November 2016 * +* @precisions fortran d -> z c * *> \ingroup real_lin * @@ -201,13 +202,13 @@ PARAMETER ( NTESTS = 9 ) * .. * .. Local Scalars .. - LOGICAL TRFCON, ZEROT + LOGICAL ZEROT CHARACTER DIST, TYPE, UPLO, XTYPE CHARACTER*3 PATH, MATPATH INTEGER I, I1, I2, IMAT, IN, INB, INFO, IOFF, IRHS, $ IUPLO, IZERO, J, K, KL, KU, LDA, LWORK, MODE, $ N, NB, NERRS, NFAIL, NIMAT, NRHS, NRUN, NT - REAL ANORM, CNDNUM, RCONDC + REAL ANORM, CNDNUM * .. * .. Local Arrays .. CHARACTER UPLOS( 2 ) @@ -430,9 +431,9 @@ * block factorization, LWORK is the length of AINV. * SRNAMT = 'SSYTRF_AA' - LWORK = N*NB + N + LWORK = MAX( 1, N*NB + N ) CALL SSYTRF_AA( UPLO, N, AFAC, LDA, IWORK, AINV, - $ LWORK, INFO ) + $ LWORK, INFO ) * * Adjust the expected value of INFO to account for * pivoting. @@ -462,19 +463,11 @@ $ NOUT ) END IF * -* Set the condition estimate flag if the INFO is not 0. -* - IF( INFO.NE.0 ) THEN - TRFCON = .TRUE. - ELSE - TRFCON = .FALSE. - END IF -* *+ TEST 1 * Reconstruct matrix from factors and compute residual. * CALL SSYT01_AA( UPLO, N, A, LDA, AFAC, LDA, IWORK, - $ AINV, LDA, RWORK, RESULT( 1 ) ) + $ AINV, LDA, RWORK, RESULT( 1 ) ) NT = 1 * * @@ -492,10 +485,9 @@ 110 CONTINUE NRUN = NRUN + NT * -* Do only the condition estimate if INFO is not 0. +* Skip solver test if INFO is not 0. * - IF( TRFCON ) THEN - RCONDC = ZERO + IF( INFO.NE.0 ) THEN GO TO 140 END IF * @@ -504,7 +496,7 @@ DO 130 IRHS = 1, NNS NRHS = NSVAL( IRHS ) * -*+ TEST 3 ( Using TRS) +*+ TEST 2 (Using TRS) * Solve and compute residual for A * X = B. * * Choose a set of NRHS random solution vectors @@ -517,10 +509,10 @@ CALL SLACPY( 'Full', N, NRHS, B, LDA, X, LDA ) * SRNAMT = 'SSYTRS_AA' - LWORK = 3*N-2 + LWORK = MAX( 1, 3*N-2 ) CALL SSYTRS_AA( UPLO, N, NRHS, AFAC, LDA, - $ IWORK, X, LDA, WORK, LWORK, - $ INFO ) + $ IWORK, X, LDA, WORK, LWORK, + $ INFO ) * * Check error code from SSYTRS and handle error. * diff --git a/TESTING/LIN/sdrvsy_aa.f b/TESTING/LIN/sdrvsy_aa.f index 3fef3c70..da5cf8a2 100644 --- a/TESTING/LIN/sdrvsy_aa.f +++ b/TESTING/LIN/sdrvsy_aa.f @@ -9,8 +9,8 @@ * =========== * * SUBROUTINE SDRVSY_AA( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, NMAX, -* A, AFAC, AINV, B, X, XACT, WORK, RWORK, IWORK, -* NOUT ) +* A, AFAC, AINV, B, X, XACT, WORK, RWORK, IWORK, +* NOUT ) * * .. Scalar Arguments .. * LOGICAL TSTERR @@ -149,8 +149,8 @@ * * ===================================================================== SUBROUTINE SDRVSY_AA( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, - $ NMAX, A, AFAC, AINV, B, X, XACT, WORK, - $ RWORK, IWORK, NOUT ) + $ NMAX, A, AFAC, AINV, B, X, XACT, WORK, + $ RWORK, IWORK, NOUT ) * * -- LAPACK test routine (version 3.7.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -184,9 +184,9 @@ CHARACTER DIST, FACT, TYPE, UPLO, XTYPE CHARACTER*3 MATPATH, PATH INTEGER I, I1, I2, IFACT, IMAT, IN, INFO, IOFF, IUPLO, - $ IZERO, J, K, K1, KL, KU, LDA, LWORK, MODE, N, + $ IZERO, J, K, KL, KU, LDA, LWORK, MODE, N, $ NB, NBMIN, NERRS, NFAIL, NIMAT, NRUN, NT - REAL AINVNM, ANORM, CNDNUM, RCOND, RCONDC + REAL ANORM, CNDNUM * .. * .. Local Arrays .. CHARACTER FACTS( NFACT ), UPLOS( 2 ) @@ -374,44 +374,6 @@ * FACT = FACTS( IFACT ) * -* Compute the condition number for comparison with -* the value returned by SSYSVX. -* - IF( ZEROT ) THEN - IF( IFACT.EQ.1 ) - $ GO TO 150 - RCONDC = ZERO -* - ELSE IF( IFACT.EQ.1 ) THEN -* -* Compute the 1-norm of A. -* - ANORM = SLANSY( '1', UPLO, N, A, LDA, RWORK ) -* -* Factor the matrix A. -* -c CALL SLACPY( UPLO, N, N, A, LDA, AFAC, LDA ) -c CALL SSYTRF( UPLO, N, AFAC, LDA, IWORK, WORK, -c $ LWORK, INFO ) -* -* Compute inv(A) and take its norm. -* -c CALL SLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA ) -c LWORK = (N+NB+1)*(NB+3) -c SRNAMT = 'DSYTRI2' -c CALL DSYTRI2( UPLO, N, AINV, LDA, IWORK, WORK, -c $ LWORK, INFO ) -c AINVNM = SLANSY( '1', UPLO, N, AINV, LDA, RWORK ) -* -* Compute the 1-norm condition number of A. -* -c IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN -c RCONDC = ONE -c ELSE -c RCONDC = ( ONE / ANORM ) / AINVNM -c END IF - END IF -* * Form an exact solution and set the right hand side. * SRNAMT = 'SLARHS' @@ -475,12 +437,7 @@ c END IF CALL SLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA ) CALL SPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK, $ LDA, RWORK, RESULT( 2 ) ) -* -* Check solution from generated exact solution. -* - CALL SGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, - $ RESULT( 3 ) ) - NT = 3 + NT = 2 * * Print information about the tests that did not pass * the threshold. diff --git a/TESTING/LIN/serrsy.f b/TESTING/LIN/serrsy.f index bf69893a..25309611 100644 --- a/TESTING/LIN/serrsy.f +++ b/TESTING/LIN/serrsy.f @@ -92,7 +92,8 @@ $ SSYTF2_RK, SSYTF2_ROOK, SSYTRF, SSYTRF_RK, $ SSYTRF_ROOK, SSYTRI, SSYTF2, SSYTRI_3, $ SSYTRI_3X, SSYTRI_ROOK, SSYTRF_AA, SSYTRI2, - $ SYTRI2X, SSYTRS, SSYTRS_3, SSYTRS_ROOK, SSYTRS_AA + $ SYTRI2X, SSYTRS, SSYTRS_3, SSYTRS_ROOK, + $ SSYTRS_AA * .. * .. Scalars in Common .. LOGICAL LERR, OK @@ -492,6 +493,12 @@ INFOT = 4 CALL SSYTRF_AA( 'U', 2, A, 1, IP, W, 4, INFO ) CALL CHKXER( 'SSYTRF_AA', INFOT, NOUT, LERR, OK ) + INFOT = 7 + CALL SSYTRF_AA( 'U', 0, A, 1, IP, W, 0, INFO ) + CALL CHKXER( 'SSYTRF_AA', INFOT, NOUT, LERR, OK ) + INFOT = 7 + CALL SSYTRF_AA( 'U', 0, A, 1, IP, W, -2, INFO ) + CALL CHKXER( 'SSYTRF_AA', INFOT, NOUT, LERR, OK ) * * SSYTRS_AA * @@ -511,6 +518,12 @@ INFOT = 8 CALL SSYTRS_AA( 'U', 2, 1, A, 2, IP, B, 1, W, 1, INFO ) CALL CHKXER( 'SSYTRS_AA', INFOT, NOUT, LERR, OK ) + INFOT = 10 + CALL SSYTRS_AA( 'U', 0, 1, A, 2, IP, B, 1, W, 0, INFO ) + CALL CHKXER( 'SSYTRS_AA', INFOT, NOUT, LERR, OK ) + INFOT = 10 + CALL SSYTRS_AA( 'U', 0, 1, A, 2, IP, B, 1, W, -2, INFO ) + CALL CHKXER( 'SSYTRS_AA', INFOT, NOUT, LERR, OK ) * ELSE IF( LSAMEN( 2, C2, 'SP' ) ) THEN * diff --git a/TESTING/LIN/ssyt01_aa.f b/TESTING/LIN/ssyt01_aa.f index 0e72fa71..5855ba22 100644 --- a/TESTING/LIN/ssyt01_aa.f +++ b/TESTING/LIN/ssyt01_aa.f @@ -9,7 +9,7 @@ * =========== * * SUBROUTINE SSYT01_AA( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, -* C, LDC, RWORK, RESID ) +* C, LDC, RWORK, RESID ) * * .. Scalar Arguments .. * CHARACTER UPLO @@ -123,7 +123,7 @@ * * ===================================================================== SUBROUTINE SSYT01_AA( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, C, - $ LDC, RWORK, RESID ) + $ LDC, RWORK, RESID ) * * -- LAPACK test routine (version 3.7.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -144,7 +144,7 @@ * ===================================================================== * * .. Parameters .. - REAL ZERO, ONE + REAL ZERO, ONE PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) * .. * .. Local Scalars .. @@ -192,27 +192,27 @@ CALL SLACPY( 'F', 1, N-1, AFAC( 2, 1 ), LDAFAC+1, C( 2, 1 ), $ LDC+1 ) ENDIF - ENDIF * -* Call STRMM to form the product U' * D (or L * D ). +* Call STRMM to form the product U' * D (or L * D ). * - IF( LSAME( UPLO, 'U' ) ) THEN - CALL STRMM( 'Left', UPLO, 'Transpose', 'Unit', N-1, N, - $ ONE, AFAC( 1, 2 ), LDAFAC, C( 2, 1 ), LDC ) - ELSE - CALL STRMM( 'Left', UPLO, 'No transpose', 'Unit', N-1, N, - $ ONE, AFAC( 2, 1 ), LDAFAC, C( 2, 1 ), LDC ) - END IF + IF( LSAME( UPLO, 'U' ) ) THEN + CALL STRMM( 'Left', UPLO, 'Transpose', 'Unit', N-1, N, + $ ONE, AFAC( 1, 2 ), LDAFAC, C( 2, 1 ), LDC ) + ELSE + CALL STRMM( 'Left', UPLO, 'No transpose', 'Unit', N-1, N, + $ ONE, AFAC( 2, 1 ), LDAFAC, C( 2, 1 ), LDC ) + END IF * -* Call STRMM again to multiply by U (or L ). +* Call STRMM again to multiply by U (or L ). * - IF( LSAME( UPLO, 'U' ) ) THEN - CALL STRMM( 'Right', UPLO, 'No transpose', 'Unit', N, N-1, - $ ONE, AFAC( 1, 2 ), LDAFAC, C( 1, 2 ), LDC ) - ELSE - CALL STRMM( 'Right', UPLO, 'Transpose', 'Unit', N, N-1, - $ ONE, AFAC( 2, 1 ), LDAFAC, C( 1, 2 ), LDC ) - END IF + IF( LSAME( UPLO, 'U' ) ) THEN + CALL STRMM( 'Right', UPLO, 'No transpose', 'Unit', N, N-1, + $ ONE, AFAC( 1, 2 ), LDAFAC, C( 1, 2 ), LDC ) + ELSE + CALL STRMM( 'Right', UPLO, 'Transpose', 'Unit', N, N-1, + $ ONE, AFAC( 2, 1 ), LDAFAC, C( 1, 2 ), LDC ) + END IF + ENDIF * * Apply symmetric pivots * diff --git a/TESTING/LIN/zchkaa.f b/TESTING/LIN/zchkaa.f index f9be8451..818f1e63 100644 --- a/TESTING/LIN/zchkaa.f +++ b/TESTING/LIN/zchkaa.f @@ -165,11 +165,12 @@ $ ZCHKHE_ROOK, ZCHKHE_RK, ZCHKHE_AA, ZCHKHP, $ ZCHKLQ, ZCHKPB, ZCHKPO, ZCHKPS, ZCHKPP, ZCHKPT, $ ZCHKQ3, ZCHKQL, ZCHKQR, ZCHKRQ, ZCHKSP, ZCHKSY, - $ ZCHKSY_ROOK, ZCHKSY_RK, ZCHKTB, ZCHKTP, ZCHKTR, - $ ZCHKTZ, ZDRVGB, ZDRVGE, ZDRVGT, ZDRVHE, ZDRVHE_ROOK, - $ ZDRVHE_RK, ZDRVHE_AA, ZDRVHP, ZDRVLS, ZDRVPB, - $ ZDRVPO, ZDRVPP, ZDRVPT, ZDRVSP, ZDRVSY, ZDRVSY_ROOK, - $ ZDRVSY_RK, ILAVER, ZCHKQRT, ZCHKQRTP, ZCHKLQT, + $ ZCHKSY_ROOK, ZCHKSY_RK, ZCHKSY_AA, ZCHKTB, + $ ZCHKTP, ZCHKTR, ZCHKTZ, ZDRVGB, ZDRVGE, ZDRVGT, + $ ZDRVHE, ZDRVHE_ROOK, ZDRVHE_RK, ZDRVHE_AA, + $ ZDRVHP, ZDRVLS, ZDRVPB, ZDRVPO, ZDRVPP, ZDRVPT, + $ ZDRVSP, ZDRVSY, ZDRVSY_ROOK, ZDRVSY_RK, + $ ZDRVSY_AA, ILAVER, ZCHKQRT, ZCHKQRTP, ZCHKLQT, $ ZCHKLQTP, ZCHKTSQR * .. * .. Scalars in Common .. @@ -828,6 +829,33 @@ WRITE( NOUT, FMT = 9988 )PATH END IF * + ELSE IF( LSAMEN( 2, C2, 'SA' ) ) THEN +* +* SK: symmetric indefinite matrices, +* with bounded Bunch-Kaufman (rook) pivoting algorithm, +* differnet matrix storage format than SR path version. +* + NTYPES = 11 + CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) +* + IF( TSTCHK ) THEN + CALL ZCHKSY_AA( DOTYPE, NN, NVAL, NNB2, NBVAL2, NNS, NSVAL, + $ THRESH, TSTERR, LDA, A( 1, 1 ), A( 1, 2 ), + $ A( 1, 3 ), B( 1, 1 ), B( 1, 2 ), + $ B( 1, 3 ), WORK, RWORK, IWORK, NOUT ) + ELSE + WRITE( NOUT, FMT = 9989 )PATH + END IF +* + IF( TSTDRV ) THEN + CALL ZDRVSY_AA( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, + $ LDA, A( 1, 1 ), A( 1, 2 ), A( 1, 3 ), + $ B( 1, 1 ), B( 1, 2 ), B( 1, 3 ), WORK, + $ RWORK, IWORK, NOUT ) + ELSE + WRITE( NOUT, FMT = 9988 )PATH + END IF +* ELSE IF( LSAMEN( 2, C2, 'SP' ) ) THEN * * SP: symmetric indefinite packed matrices, diff --git a/TESTING/LIN/zchkhe_aa.f b/TESTING/LIN/zchkhe_aa.f index 20d595be..5cc4b9ef 100644 --- a/TESTING/LIN/zchkhe_aa.f +++ b/TESTING/LIN/zchkhe_aa.f @@ -205,13 +205,13 @@ PARAMETER ( NTESTS = 9 ) * .. * .. Local Scalars .. - LOGICAL TRFCON, ZEROT + LOGICAL ZEROT CHARACTER DIST, TYPE, UPLO, XTYPE CHARACTER*3 PATH, MATPATH INTEGER I, I1, I2, IMAT, IN, INB, INFO, IOFF, IRHS, $ IUPLO, IZERO, J, K, KL, KU, LDA, LWORK, MODE, $ N, NB, NERRS, NFAIL, NIMAT, NRHS, NRUN, NT - DOUBLE PRECISION ANORM, CNDNUM, RCOND, RCONDC + DOUBLE PRECISION ANORM, CNDNUM * .. * .. Local Arrays .. CHARACTER UPLOS( 2 ) @@ -224,7 +224,7 @@ * .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, XLAENV, ZERRHE, ZGET04, - $ ZHECON, ZHERFS, ZHET01, ZHETRF_AA, ZHETRI2, + $ ZHECON, ZHERFS, ZHET01_AA, ZHETRF_AA, ZHETRI2, $ ZHETRS_AA, ZLACPY, ZLAIPD, ZLARHS, ZLATB4, $ ZLATMS, ZPOT02, ZPOT03, ZPOT05 * .. @@ -430,10 +430,10 @@ * the block structure of D. AINV is a work array for * block factorization, LWORK is the length of AINV. * - LWORK = ( NB+1 )*LDA + LWORK = MAX( 1, ( NB+1 )*LDA ) SRNAMT = 'ZHETRF_AA' CALL ZHETRF_AA( UPLO, N, AFAC, LDA, IWORK, AINV, - $ LWORK, INFO ) + $ LWORK, INFO ) * * Adjust the expected value of INFO to account for * pivoting. @@ -463,19 +463,11 @@ $ NOUT ) END IF * -* Set the condition estimate flag if the INFO is not 0. -* - IF( INFO.NE.0 ) THEN - TRFCON = .TRUE. - ELSE - TRFCON = .FALSE. - END IF -* *+ TEST 1 * Reconstruct matrix from factors and compute residual. * CALL ZHET01_AA( UPLO, N, A, LDA, AFAC, LDA, IWORK, - $ AINV, LDA, RWORK, RESULT( 1 ) ) + $ AINV, LDA, RWORK, RESULT( 1 ) ) NT = 1 * * @@ -493,10 +485,9 @@ 110 CONTINUE NRUN = NRUN + NT * -* Do only the condition estimate if INFO is not 0. +* Skip solver test if INFO is not 0. * - IF( TRFCON ) THEN - RCONDC = ZERO + IF( INFO.NE.0 ) THEN GO TO 140 END IF * @@ -505,7 +496,7 @@ DO 130 IRHS = 1, NNS NRHS = NSVAL( IRHS ) * -*+ TEST 3 (Using TRS) +*+ TEST 2 (Using TRS) * Solve and compute residual for A * X = B. * * Choose a set of NRHS random solution vectors @@ -518,9 +509,9 @@ CALL ZLACPY( 'Full', N, NRHS, B, LDA, X, LDA ) * SRNAMT = 'ZHETRS_AA' - LWORK = 3*N-2 + LWORK = MAX( 1, 3*N-2 ) CALL ZHETRS_AA( UPLO, N, NRHS, AFAC, LDA, IWORK, - $ X, LDA, WORK, LWORK, INFO ) + $ X, LDA, WORK, LWORK, INFO ) * * Check error code from ZHETRS and handle error. * diff --git a/TESTING/LIN/zchksy_aa.f b/TESTING/LIN/zchksy_aa.f new file mode 100644 index 00000000..4c2cd041 --- /dev/null +++ b/TESTING/LIN/zchksy_aa.f @@ -0,0 +1,572 @@ +*> \brief \b ZCHKSY_AA +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE ZCHKSY_AA( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL, +* THRESH, TSTERR, NMAX, A, AFAC, AINV, B, X, +* XACT, WORK, RWORK, IWORK, NOUT ) +* +* .. Scalar Arguments .. +* LOGICAL TSTERR +* INTEGER NMAX, NN, NNB, NNS, NOUT +* COMPLEX*16 THRESH +* .. +* .. Array Arguments .. +* LOGICAL DOTYPE( * ) +* INTEGER IWORK( * ), NBVAL( * ), NSVAL( * ), NVAL( * ) +* COMPLEX*16 A( * ), AFAC( * ), AINV( * ), B( * ), +* $ RWORK( * ), WORK( * ), X( * ), XACT( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> ZCHKSY_AA tests ZSYTRF_AA, -TRS_AA. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] DOTYPE +*> \verbatim +*> DOTYPE is LOGICAL array, dimension (NTYPES) +*> The matrix types to be used for testing. Matrices of type j +*> (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) = +*> .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used. +*> \endverbatim +*> +*> \param[in] NN +*> \verbatim +*> NN is INTEGER +*> The number of values of N contained in the vector NVAL. +*> \endverbatim +*> +*> \param[in] NVAL +*> \verbatim +*> NVAL is INTEGER array, dimension (NN) +*> The values of the matrix dimension N. +*> \endverbatim +*> +*> \param[in] NNB +*> \verbatim +*> NNB is INTEGER +*> The number of values of NB contained in the vector NBVAL. +*> \endverbatim +*> +*> \param[in] NBVAL +*> \verbatim +*> NBVAL is INTEGER array, dimension (NBVAL) +*> The values of the blocksize NB. +*> \endverbatim +*> +*> \param[in] NNS +*> \verbatim +*> NNS is INTEGER +*> The number of values of NRHS contained in the vector NSVAL. +*> \endverbatim +*> +*> \param[in] NSVAL +*> \verbatim +*> NSVAL is INTEGER array, dimension (NNS) +*> The values of the number of right hand sides NRHS. +*> \endverbatim +*> +*> \param[in] THRESH +*> \verbatim +*> THRESH is COMPLEX*16 +*> The threshold value for the test ratios. A result is +*> included in the output file if RESULT >= THRESH. To have +*> every test ratio printed, use THRESH = 0. +*> \endverbatim +*> +*> \param[in] TSTERR +*> \verbatim +*> TSTERR is LOGICAL +*> Flag that indicates whether error exits are to be tested. +*> \endverbatim +*> +*> \param[in] NMAX +*> \verbatim +*> NMAX is INTEGER +*> The maximum value permitted for N, used in dimensioning the +*> work arrays. +*> \endverbatim +*> +*> \param[out] A +*> \verbatim +*> A is COMPLEX*16 array, dimension (NMAX*NMAX) +*> \endverbatim +*> +*> \param[out] AFAC +*> \verbatim +*> AFAC is COMPLEX*16 array, dimension (NMAX*NMAX) +*> \endverbatim +*> +*> \param[out] AINV +*> \verbatim +*> AINV is COMPLEX*16 array, dimension (NMAX*NMAX) +*> \endverbatim +*> +*> \param[out] B +*> \verbatim +*> B is COMPLEX*16 array, dimension (NMAX*NSMAX) +*> where NSMAX is the largest entry in NSVAL. +*> \endverbatim +*> +*> \param[out] X +*> \verbatim +*> X is COMPLEX*16 array, dimension (NMAX*NSMAX) +*> \endverbatim +*> +*> \param[out] XACT +*> \verbatim +*> XACT is COMPLEX*16 array, dimension (NMAX*NSMAX) +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is COMPLEX*16 array, dimension (NMAX*max(3,NSMAX)) +*> \endverbatim +*> +*> \param[out] RWORK +*> \verbatim +*> RWORK is COMPLEX*16 array, dimension (max(NMAX,2*NSMAX)) +*> \endverbatim +*> +*> \param[out] IWORK +*> \verbatim +*> IWORK is INTEGER array, dimension (2*NMAX) +*> \endverbatim +*> +*> \param[in] NOUT +*> \verbatim +*> NOUT is INTEGER +*> The unit number for output. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2016 +* +* @generated from LIN/dchksy_aa.f, fortran d -> z, Wed Nov 16 21:34:18 2016 +* +*> \ingroup complex16_lin +* +* ===================================================================== + SUBROUTINE ZCHKSY_AA( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL, + $ THRESH, TSTERR, NMAX, A, AFAC, AINV, B, + $ X, XACT, WORK, RWORK, IWORK, NOUT ) +* +* -- LAPACK test routine (version 3.7.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2016 +* + IMPLICIT NONE +* +* .. Scalar Arguments .. + LOGICAL TSTERR + INTEGER NN, NNB, NNS, NMAX, NOUT + DOUBLE PRECISION THRESH +* .. +* .. Array Arguments .. + LOGICAL DOTYPE( * ) + INTEGER IWORK( * ), NBVAL( * ), NSVAL( * ), NVAL( * ) + DOUBLE PRECISION RWORK( * ) + COMPLEX*16 A( * ), AFAC( * ), AINV( * ), B( * ), + $ WORK( * ), X( * ), XACT( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ZERO + PARAMETER ( ZERO = 0.0D+0 ) + COMPLEX*16 CZERO + PARAMETER ( CZERO = 0.0E+0 ) + INTEGER NTYPES + PARAMETER ( NTYPES = 10 ) + INTEGER NTESTS + PARAMETER ( NTESTS = 9 ) +* .. +* .. Local Scalars .. + LOGICAL ZEROT + CHARACTER DIST, TYPE, UPLO, XTYPE + CHARACTER*3 PATH, MATPATH + INTEGER I, I1, I2, IMAT, IN, INB, INFO, IOFF, IRHS, + $ IUPLO, IZERO, J, K, KL, KU, LDA, LWORK, MODE, + $ N, NB, NERRS, NFAIL, NIMAT, NRHS, NRUN, NT + DOUBLE PRECISION ANORM, CNDNUM +* .. +* .. Local Arrays .. + CHARACTER UPLOS( 2 ) + INTEGER ISEED( 4 ), ISEEDY( 4 ) + DOUBLE PRECISION RESULT( NTESTS ) +* .. +* .. External Functions .. + DOUBLE PRECISION DGET06, ZLANSY + EXTERNAL DGET06, ZLANSY +* .. +* .. External Subroutines .. + EXTERNAL ALAERH, ALAHD, ALASUM, ZERRSY, ZGET04, ZLACPY, + $ ZLARHS, ZLATB4, ZLATMS, ZSYT02, DSYT03, DSYT05, + $ DSYCON, ZSYRFS, ZSYT01_AA, ZSYTRF_AA, + $ DSYTRI2, ZSYTRS_AA, XLAENV +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. Scalars in Common .. + LOGICAL LERR, OK + CHARACTER*32 SRNAMT + INTEGER INFOT, NUNIT +* .. +* .. Common blocks .. + COMMON / INFOC / INFOT, NUNIT, OK, LERR + COMMON / SRNAMC / SRNAMT +* .. +* .. Data statements .. + DATA ISEEDY / 1988, 1989, 1990, 1991 / + DATA UPLOS / 'U', 'L' / +* .. +* .. Executable Statements .. +* +* Initialize constants and the random number seed. +* +* Test path +* + PATH( 1: 1 ) = 'Zomplex precision' + PATH( 2: 3 ) = 'SA' +* +* Path to generate matrices +* + MATPATH( 1: 1 ) = 'Zomplex precision' + MATPATH( 2: 3 ) = 'SY' + NRUN = 0 + NFAIL = 0 + NERRS = 0 + DO 10 I = 1, 4 + ISEED( I ) = ISEEDY( I ) + 10 CONTINUE +* +* Test the error exits +* + IF( TSTERR ) + $ CALL ZERRSY( PATH, NOUT ) + INFOT = 0 +* +* Set the minimum block size for which the block routine should +* be used, which will be later returned by ILAENV +* + CALL XLAENV( 2, 2 ) +* +* Do for each value of N in NVAL +* + DO 180 IN = 1, NN + N = NVAL( IN ) + IF( N .GT. NMAX ) THEN + NFAIL = NFAIL + 1 + WRITE(NOUT, 9995) 'M ', N, NMAX + GO TO 180 + END IF + LDA = MAX( N, 1 ) + XTYPE = 'N' + NIMAT = NTYPES + IF( N.LE.0 ) + $ NIMAT = 1 +* + IZERO = 0 +* +* Do for each value of matrix type IMAT +* + DO 170 IMAT = 1, NIMAT +* +* Do the tests only if DOTYPE( IMAT ) is true. +* + IF( .NOT.DOTYPE( IMAT ) ) + $ GO TO 170 +* +* Skip types 3, 4, 5, or 6 if the matrix size is too small. +* + ZEROT = IMAT.GE.3 .AND. IMAT.LE.6 + IF( ZEROT .AND. N.LT.IMAT-2 ) + $ GO TO 170 +* +* Do first for UPLO = 'U', then for UPLO = 'L' +* + DO 160 IUPLO = 1, 2 + UPLO = UPLOS( IUPLO ) +* +* Begin generate the test matrix A. +* +* +* Set up parameters with ZLATB4 for the matrix generator +* based on the type of matrix to be generated. +* + CALL ZLATB4( MATPATH, IMAT, N, N, TYPE, KL, KU, + $ ANORM, MODE, CNDNUM, DIST ) +* +* Generate a matrix with ZLATMS. +* + SRNAMT = 'ZLATMS' + CALL ZLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE, + $ CNDNUM, ANORM, KL, KU, UPLO, A, LDA, WORK, + $ INFO ) +* +* Check error code from ZLATMS and handle error. +* + IF( INFO.NE.0 ) THEN + CALL ALAERH( PATH, 'ZLATMS', INFO, 0, UPLO, N, N, -1, + $ -1, -1, IMAT, NFAIL, NERRS, NOUT ) +* +* Skip all tests for this generated matrix +* + GO TO 160 + END IF +* +* For matrix types 3-6, zero one or more rows and +* columns of the matrix to test that INFO is returned +* correctly. +* + IF( ZEROT ) THEN + IF( IMAT.EQ.3 ) THEN + IZERO = 1 + ELSE IF( IMAT.EQ.4 ) THEN + IZERO = N + ELSE + IZERO = N / 2 + 1 + END IF +* + IF( IMAT.LT.6 ) THEN +* +* Set row and column IZERO to zero. +* + IF( IUPLO.EQ.1 ) THEN + IOFF = ( IZERO-1 )*LDA + DO 20 I = 1, IZERO - 1 + A( IOFF+I ) = CZERO + 20 CONTINUE + IOFF = IOFF + IZERO + DO 30 I = IZERO, N + A( IOFF ) = CZERO + IOFF = IOFF + LDA + 30 CONTINUE + ELSE + IOFF = IZERO + DO 40 I = 1, IZERO - 1 + A( IOFF ) = CZERO + IOFF = IOFF + LDA + 40 CONTINUE + IOFF = IOFF - IZERO + DO 50 I = IZERO, N + A( IOFF+I ) = CZERO + 50 CONTINUE + END IF + ELSE + IF( IUPLO.EQ.1 ) THEN +* +* Set the first IZERO rows and columns to zero. +* + IOFF = 0 + DO 70 J = 1, N + I2 = MIN( J, IZERO ) + DO 60 I = 1, I2 + A( IOFF+I ) = CZERO + 60 CONTINUE + IOFF = IOFF + LDA + 70 CONTINUE + IZERO = 1 + ELSE +* +* Set the last IZERO rows and columns to zero. +* + IOFF = 0 + DO 90 J = 1, N + I1 = MAX( J, IZERO ) + DO 80 I = I1, N + A( IOFF+I ) = CZERO + 80 CONTINUE + IOFF = IOFF + LDA + 90 CONTINUE + END IF + END IF + ELSE + IZERO = 0 + END IF +* +* End generate the test matrix A. +* +* Do for each value of NB in NBVAL +* + DO 150 INB = 1, NNB +* +* Set the optimal blocksize, which will be later +* returned by ILAENV. +* + NB = NBVAL( INB ) + CALL XLAENV( 1, NB ) +* +* Copy the test matrix A into matrix AFAC which +* will be factorized in place. This is needed to +* preserve the test matrix A for subsequent tests. +* + CALL ZLACPY( UPLO, N, N, A, LDA, AFAC, LDA ) +* +* Compute the L*D*L**T or U*D*U**T factorization of the +* matrix. IWORK stores details of the interchanges and +* the block structure of D. AINV is a work array for +* block factorization, LWORK is the length of AINV. +* + SRNAMT = 'ZSYTRF_AA' + LWORK = MAX( 1, N*NB + N ) + CALL ZSYTRF_AA( UPLO, N, AFAC, LDA, IWORK, AINV, + $ LWORK, INFO ) +* +* Adjust the expected value of INFO to account for +* pivoting. +* + IF( IZERO.GT.0 ) THEN + J = 1 + K = IZERO + 100 CONTINUE + IF( J.EQ.K ) THEN + K = IWORK( J ) + ELSE IF( IWORK( J ).EQ.K ) THEN + K = J + END IF + IF( J.LT.K ) THEN + J = J + 1 + GO TO 100 + END IF + ELSE + K = 0 + END IF +* +* Check error code from ZSYTRF and handle error. +* + IF( INFO.NE.K ) THEN + CALL ALAERH( PATH, 'ZSYTRF_AA', INFO, K, UPLO, + $ N, N, -1, -1, NB, IMAT, NFAIL, NERRS, + $ NOUT ) + END IF +* +*+ TEST 1 +* Reconstruct matrix from factors and compute residual. +* + CALL ZSYT01_AA( UPLO, N, A, LDA, AFAC, LDA, IWORK, + $ AINV, LDA, RWORK, RESULT( 1 ) ) + NT = 1 +* +* +* Print information about the tests that did not pass +* the threshold. +* + DO 110 K = 1, NT + IF( RESULT( K ).GE.THRESH ) THEN + IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) + $ CALL ALAHD( NOUT, PATH ) + WRITE( NOUT, FMT = 9999 )UPLO, N, NB, IMAT, K, + $ RESULT( K ) + NFAIL = NFAIL + 1 + END IF + 110 CONTINUE + NRUN = NRUN + NT +* +* Skip solver test if INFO is not 0. +* + IF( INFO.NE.0 ) THEN + GO TO 140 + END IF +* +* Do for each value of NRHS in NSVAL. +* + DO 130 IRHS = 1, NNS + NRHS = NSVAL( IRHS ) +* +*+ TEST 2 (Using TRS) +* Solve and compute residual for A * X = B. +* +* Choose a set of NRHS random solution vectors +* stored in XACT and set up the right hand side B +* + SRNAMT = 'ZLARHS' + CALL ZLARHS( MATPATH, XTYPE, UPLO, ' ', N, N, + $ KL, KU, NRHS, A, LDA, XACT, LDA, + $ B, LDA, ISEED, INFO ) + CALL ZLACPY( 'Full', N, NRHS, B, LDA, X, LDA ) +* + SRNAMT = 'ZSYTRS_AA' + LWORK = MAX( 1, 3*N-2 ) + CALL ZSYTRS_AA( UPLO, N, NRHS, AFAC, LDA, + $ IWORK, X, LDA, WORK, LWORK, + $ INFO ) +* +* Check error code from ZSYTRS and handle error. +* + IF( INFO.NE.0 ) THEN + CALL ALAERH( PATH, 'ZSYTRS_AA', INFO, 0, + $ UPLO, N, N, -1, -1, NRHS, IMAT, + $ NFAIL, NERRS, NOUT ) + END IF +* + CALL ZLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA ) +* +* Compute the residual for the solution +* + CALL ZSYT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK, + $ LDA, RWORK, RESULT( 2 ) ) +* +* +* Print information about the tests that did not pass +* the threshold. +* + DO 120 K = 2, 2 + IF( RESULT( K ).GE.THRESH ) THEN + IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) + $ CALL ALAHD( NOUT, PATH ) + WRITE( NOUT, FMT = 9998 )UPLO, N, NRHS, + $ IMAT, K, RESULT( K ) + NFAIL = NFAIL + 1 + END IF + 120 CONTINUE + NRUN = NRUN + 1 +* +* End do for each value of NRHS in NSVAL. +* + 130 CONTINUE + 140 CONTINUE + 150 CONTINUE + 160 CONTINUE + 170 CONTINUE + 180 CONTINUE +* +* Print a summary of the results. +* + CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS ) +* + 9999 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', NB =', I4, ', type ', + $ I2, ', test ', I2, ', ratio =', G12.5 ) + 9998 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', NRHS=', I3, ', type ', + $ I2, ', test(', I2, ') =', G12.5 ) + 9995 FORMAT( ' Invalid input value: ', A4, '=', I6, '; must be <=', + $ I6 ) + RETURN +* +* End of ZCHKSY_AA +* + END diff --git a/TESTING/LIN/zdrvhe_aa.f b/TESTING/LIN/zdrvhe_aa.f index 3a43cf79..3b59395a 100644 --- a/TESTING/LIN/zdrvhe_aa.f +++ b/TESTING/LIN/zdrvhe_aa.f @@ -9,8 +9,8 @@ * =========== * * SUBROUTINE ZDRVHE_AA( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, NMAX, -* A, AFAC, AINV, B, X, XACT, WORK, RWORK, IWORK, -* NOUT ) +* A, AFAC, AINV, B, X, XACT, WORK, RWORK, IWORK, +* NOUT ) * * .. Scalar Arguments .. * LOGICAL TSTERR @@ -186,9 +186,9 @@ CHARACTER DIST, FACT, TYPE, UPLO, XTYPE CHARACTER*3 MATPATH, PATH INTEGER I, I1, I2, IFACT, IMAT, IN, INFO, IOFF, IUPLO, - $ IZERO, J, K, K1, KL, KU, LDA, LWORK, MODE, N, + $ IZERO, J, K, KL, KU, LDA, LWORK, MODE, N, $ NB, NBMIN, NERRS, NFAIL, NIMAT, NRUN, NT - DOUBLE PRECISION AINVNM, ANORM, CNDNUM, RCOND, RCONDC + DOUBLE PRECISION ANORM, CNDNUM * .. * .. Local Arrays .. CHARACTER FACTS( NFACT ), UPLOS( 2 ) @@ -202,8 +202,8 @@ * .. External Subroutines .. EXTERNAL ALADHD, ALAERH, ALASVM, XLAENV, ZERRVX, ZGET04, $ ZHESV_AA, ZHET01_AA, ZHETRF_AA, - $ ZHETRI2, ZLACPY, ZLAIPD, ZLARHS, ZLATB4, ZLATMS, - $ ZPOT02 + $ ZHETRI2, ZLACPY, ZLAIPD, ZLARHS, ZLATB4, + $ ZLATMS, ZPOT02 * .. * .. Scalars in Common .. LOGICAL LERR, OK @@ -383,43 +383,6 @@ * FACT = FACTS( IFACT ) * -* Compute the condition number for comparison with -* the value returned by ZHESVX. -* - IF( ZEROT ) THEN - IF( IFACT.EQ.1 ) - $ GO TO 150 - RCONDC = ZERO -* - ELSE IF( IFACT.EQ.1 ) THEN -* -* Compute the 1-norm of A. -* - ANORM = ZLANHE( '1', UPLO, N, A, LDA, RWORK ) -* -* Factor the matrix A. -* -c CALL ZLACPY( UPLO, N, N, A, LDA, AFAC, LDA ) -c CALL ZHETRF( UPLO, N, AFAC, LDA, IWORK, WORK, -c $ LWORK, INFO ) -* -* Compute inv(A) and take its norm. -* -c CALL ZLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA ) -c LWORK = (N+NB+1)*(NB+3) -c CALL ZHETRI2( UPLO, N, AINV, LDA, IWORK, WORK, -c $ LWORK, INFO ) -c AINVNM = ZLANHE( '1', UPLO, N, AINV, LDA, RWORK ) -* -* Compute the 1-norm condition number of A. -* -c IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN -c RCONDC = ONE -c ELSE -c RCONDC = ( ONE / ANORM ) / AINVNM -c END IF - END IF -* * Form an exact solution and set the right hand side. * SRNAMT = 'ZLARHS' @@ -483,12 +446,7 @@ c END IF CALL ZLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA ) CALL ZPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK, $ LDA, RWORK, RESULT( 2 ) ) -* -* Check solution from generated exact solution. -* - CALL ZGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, - $ RESULT( 3 ) ) - NT = 3 + NT = 2 * * Print information about the tests that did not pass * the threshold. diff --git a/TESTING/LIN/zdrvsy_aa.f b/TESTING/LIN/zdrvsy_aa.f new file mode 100644 index 00000000..d0a9711b --- /dev/null +++ b/TESTING/LIN/zdrvsy_aa.f @@ -0,0 +1,480 @@ +*> \brief \b ZDRVSY_AA +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE ZDRVSY_AA( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, NMAX, +* A, AFAC, AINV, B, X, XACT, WORK, RWORK, IWORK, +* NOUT ) +* +* .. Scalar Arguments .. +* LOGICAL TSTERR +* INTEGER NMAX, NN, NOUT, NRHS +* DOUBLE PRECISION THRESH +* .. +* .. Array Arguments .. +* LOGICAL DOTYPE( * ) +* INTEGER IWORK( * ), NVAL( * ) +* DOUBLE PRECISION RWORK( * ) +* COMPLEX*16 A( * ), AFAC( * ), AINV( * ), B( * ), +* $ WORK( * ), X( * ), XACT( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> ZDRVSY_AA tests the driver routine ZSYSV_AA. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] DOTYPE +*> \verbatim +*> DOTYPE is LOGICAL array, dimension (NTYPES) +*> The matrix types to be used for testing. Matrices of type j +*> (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) = +*> .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used. +*> \endverbatim +*> +*> \param[in] NN +*> \verbatim +*> NN is INTEGER +*> The number of values of N contained in the vector NVAL. +*> \endverbatim +*> +*> \param[in] NVAL +*> \verbatim +*> NVAL is INTEGER array, dimension (NN) +*> The values of the matrix dimension N. +*> \endverbatim +*> +*> \param[in] NRHS +*> \verbatim +*> NRHS is INTEGER +*> The number of right hand side vectors to be generated for +*> each linear system. +*> \endverbatim +*> +*> \param[in] THRESH +*> \verbatim +*> THRESH is COMPLEX*16 +*> The threshold value for the test ratios. A result is +*> included in the output file if RESULT >= THRESH. To have +*> every test ratio printed, use THRESH = 0. +*> \endverbatim +*> +*> \param[in] TSTERR +*> \verbatim +*> TSTERR is LOGICAL +*> Flag that indicates whether error exits are to be tested. +*> \endverbatim +*> +*> \param[in] NMAX +*> \verbatim +*> NMAX is INTEGER +*> The maximum value permitted for N, used in dimensioning the +*> work arrays. +*> \endverbatim +*> +*> \param[out] A +*> \verbatim +*> A is COMPLEX*16 array, dimension (NMAX*NMAX) +*> \endverbatim +*> +*> \param[out] AFAC +*> \verbatim +*> AFAC is COMPLEX*16 array, dimension (NMAX*NMAX) +*> \endverbatim +*> +*> \param[out] AINV +*> \verbatim +*> AINV is COMPLEX*16 array, dimension (NMAX*NMAX) +*> \endverbatim +*> +*> \param[out] B +*> \verbatim +*> B is COMPLEX*16 array, dimension (NMAX*NRHS) +*> \endverbatim +*> +*> \param[out] X +*> \verbatim +*> X is COMPLEX*16 array, dimension (NMAX*NRHS) +*> \endverbatim +*> +*> \param[out] XACT +*> \verbatim +*> XACT is COMPLEX*16 array, dimension (NMAX*NRHS) +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is COMPLEX*16 array, dimension (NMAX*max(2,NRHS)) +*> \endverbatim +*> +*> \param[out] RWORK +*> \verbatim +*> RWORK is COMPLEX*16 array, dimension (NMAX+2*NRHS) +*> \endverbatim +*> +*> \param[out] IWORK +*> \verbatim +*> IWORK is INTEGER array, dimension (2*NMAX) +*> \endverbatim +*> +*> \param[in] NOUT +*> \verbatim +*> NOUT is INTEGER +*> The unit number for output. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2016 +* +* @generated from LIN/ddrvsy_aa.f, fortran d -> z, Thu Nov 17 12:14:51 2016 +* +*> \ingroup complex16_lin +* +* ===================================================================== + SUBROUTINE ZDRVSY_AA( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, + $ NMAX, A, AFAC, AINV, B, X, XACT, WORK, + $ RWORK, IWORK, NOUT ) +* +* -- LAPACK test routine (version 3.7.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2016 +* +* .. Scalar Arguments .. + LOGICAL TSTERR + INTEGER NMAX, NN, NOUT, NRHS + DOUBLE PRECISION THRESH +* .. +* .. Array Arguments .. + LOGICAL DOTYPE( * ) + INTEGER IWORK( * ), NVAL( * ) + DOUBLE PRECISION RWORK( * ) + COMPLEX*16 A( * ), AFAC( * ), AINV( * ), B( * ), + $ WORK( * ), X( * ), XACT( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ZERO + PARAMETER ( ZERO = 0.0D+0 ) + COMPLEX*16 CZERO + PARAMETER ( CZERO = 0.0E+0 ) + INTEGER NTYPES, NTESTS + PARAMETER ( NTYPES = 10, NTESTS = 3 ) + INTEGER NFACT + PARAMETER ( NFACT = 2 ) +* .. +* .. Local Scalars .. + LOGICAL ZEROT + CHARACTER DIST, FACT, TYPE, UPLO, XTYPE + CHARACTER*3 MATPATH, PATH + INTEGER I, I1, I2, IFACT, IMAT, IN, INFO, IOFF, IUPLO, + $ IZERO, J, K, KL, KU, LDA, LWORK, MODE, N, + $ NB, NBMIN, NERRS, NFAIL, NIMAT, NRUN, NT + DOUBLE PRECISION ANORM, CNDNUM +* .. +* .. Local Arrays .. + CHARACTER FACTS( NFACT ), UPLOS( 2 ) + INTEGER ISEED( 4 ), ISEEDY( 4 ) + DOUBLE PRECISION RESULT( NTESTS ) +* .. +* .. External Functions .. + DOUBLE PRECISION DGET06, ZLANSY + EXTERNAL DGET06, ZLANSY +* .. +* .. External Subroutines .. + EXTERNAL ALADHD, ALAERH, ALASVM, DERRVX, ZGET04, ZLACPY, + $ ZLARHS, ZLASET, ZLATB4, ZLATMS, ZSYT02, DSYT05, + $ ZSYSV_AA, ZSYT01_AA, ZSYTRF_AA, XLAENV +* .. +* .. Scalars in Common .. + LOGICAL LERR, OK + CHARACTER*32 SRNAMT + INTEGER INFOT, NUNIT +* .. +* .. Common blocks .. + COMMON / INFOC / INFOT, NUNIT, OK, LERR + COMMON / SRNAMC / SRNAMT +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. Data statements .. + DATA ISEEDY / 1988, 1989, 1990, 1991 / + DATA UPLOS / 'U', 'L' / , FACTS / 'F', 'N' / +* .. +* .. Executable Statements .. +* +* Initialize constants and the random number seed. +* +* Test path +* + PATH( 1: 1 ) = 'Zomplex precision' + PATH( 2: 3 ) = 'SA' +* +* Path to generate matrices +* + MATPATH( 1: 1 ) = 'Zomplex precision' + MATPATH( 2: 3 ) = 'SY' +* + NRUN = 0 + NFAIL = 0 + NERRS = 0 + DO 10 I = 1, 4 + ISEED( I ) = ISEEDY( I ) + 10 CONTINUE + LWORK = MAX( 2*NMAX, NMAX*NRHS ) +* +* Test the error exits +* + IF( TSTERR ) + $ CALL ZERRVX( PATH, NOUT ) + INFOT = 0 +* +* Set the block size and minimum block size for testing. +* + NB = 1 + NBMIN = 2 + CALL XLAENV( 1, NB ) + CALL XLAENV( 2, NBMIN ) +* +* Do for each value of N in NVAL +* + DO 180 IN = 1, NN + N = NVAL( IN ) + LDA = MAX( N, 1 ) + XTYPE = 'N' + NIMAT = NTYPES + IF( N.LE.0 ) + $ NIMAT = 1 +* + DO 170 IMAT = 1, NIMAT +* +* Do the tests only if DOTYPE( IMAT ) is true. +* + IF( .NOT.DOTYPE( IMAT ) ) + $ GO TO 170 +* +* Skip types 3, 4, 5, or 6 if the matrix size is too small. +* + ZEROT = IMAT.GE.3 .AND. IMAT.LE.6 + IF( ZEROT .AND. N.LT.IMAT-2 ) + $ GO TO 170 +* +* Do first for UPLO = 'U', then for UPLO = 'L' +* + DO 160 IUPLO = 1, 2 + UPLO = UPLOS( IUPLO ) +* +* Set up parameters with ZLATB4 and generate a test matrix +* with ZLATMS. +* + CALL ZLATB4( MATPATH, IMAT, N, N, TYPE, KL, KU, ANORM, + $ MODE, CNDNUM, DIST ) +* + SRNAMT = 'ZLATMS' + CALL ZLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE, + $ CNDNUM, ANORM, KL, KU, UPLO, A, LDA, WORK, + $ INFO ) +* +* Check error code from ZLATMS. +* + IF( INFO.NE.0 ) THEN + CALL ALAERH( PATH, 'ZLATMS', INFO, 0, UPLO, N, N, -1, + $ -1, -1, IMAT, NFAIL, NERRS, NOUT ) + GO TO 160 + END IF +* +* For types 3-6, zero one or more rows and columns of the +* matrix to test that INFO is returned correctly. +* + IF( ZEROT ) THEN + IF( IMAT.EQ.3 ) THEN + IZERO = 1 + ELSE IF( IMAT.EQ.4 ) THEN + IZERO = N + ELSE + IZERO = N / 2 + 1 + END IF +* + IF( IMAT.LT.6 ) THEN +* +* Set row and column IZERO to zero. +* + IF( IUPLO.EQ.1 ) THEN + IOFF = ( IZERO-1 )*LDA + DO 20 I = 1, IZERO - 1 + A( IOFF+I ) = CZERO + 20 CONTINUE + IOFF = IOFF + IZERO + DO 30 I = IZERO, N + A( IOFF ) = CZERO + IOFF = IOFF + LDA + 30 CONTINUE + ELSE + IOFF = IZERO + DO 40 I = 1, IZERO - 1 + A( IOFF ) = CZERO + IOFF = IOFF + LDA + 40 CONTINUE + IOFF = IOFF - IZERO + DO 50 I = IZERO, N + A( IOFF+I ) = CZERO + 50 CONTINUE + END IF + ELSE + IOFF = 0 + IF( IUPLO.EQ.1 ) THEN +* +* Set the first IZERO rows and columns to zero. +* + DO 70 J = 1, N + I2 = MIN( J, IZERO ) + DO 60 I = 1, I2 + A( IOFF+I ) = CZERO + 60 CONTINUE + IOFF = IOFF + LDA + 70 CONTINUE + IZERO = 1 + ELSE +* +* Set the last IZERO rows and columns to zero. +* + DO 90 J = 1, N + I1 = MAX( J, IZERO ) + DO 80 I = I1, N + A( IOFF+I ) = CZERO + 80 CONTINUE + IOFF = IOFF + LDA + 90 CONTINUE + END IF + END IF + ELSE + IZERO = 0 + END IF +* + DO 150 IFACT = 1, NFACT +* +* Do first for FACT = 'F', then for other values. +* + FACT = FACTS( IFACT ) +* +* Form an exact solution and set the right hand side. +* + SRNAMT = 'ZLARHS' + CALL ZLARHS( MATPATH, XTYPE, UPLO, ' ', N, N, KL, KU, + $ NRHS, A, LDA, XACT, LDA, B, LDA, ISEED, + $ INFO ) + XTYPE = 'C' +* +* --- Test ZSYSV_AA --- +* + IF( IFACT.EQ.2 ) THEN + CALL ZLACPY( UPLO, N, N, A, LDA, AFAC, LDA ) + CALL ZLACPY( 'Full', N, NRHS, B, LDA, X, LDA ) +* +* Factor the matrix and solve the system using ZSYSV_AA. +* + SRNAMT = 'ZSYSV_AA' + CALL ZSYSV_AA( UPLO, N, NRHS, AFAC, LDA, IWORK, + $ X, LDA, WORK, LWORK, INFO ) +* +* Adjust the expected value of INFO to account for +* pivoting. +* + IF( IZERO.GT.0 ) THEN + J = 1 + K = IZERO + 100 CONTINUE + IF( J.EQ.K ) THEN + K = IWORK( J ) + ELSE IF( IWORK( J ).EQ.K ) THEN + K = J + END IF + IF( J.LT.K ) THEN + J = J + 1 + GO TO 100 + END IF + ELSE + K = 0 + END IF +* +* Check error code from ZSYSV_AA . +* + IF( INFO.NE.K ) THEN + CALL ALAERH( PATH, 'ZSYSV_AA ', INFO, K, + $ UPLO, N, N, -1, -1, NRHS, + $ IMAT, NFAIL, NERRS, NOUT ) + GO TO 120 + ELSE IF( INFO.NE.0 ) THEN + GO TO 120 + END IF +* +* Reconstruct matrix from factors and compute +* residual. +* + CALL ZSYT01_AA( UPLO, N, A, LDA, AFAC, LDA, + $ IWORK, AINV, LDA, RWORK, + $ RESULT( 1 ) ) +* +* Compute residual of the computed solution. +* + CALL ZLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA ) + CALL ZSYT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK, + $ LDA, RWORK, RESULT( 2 ) ) + NT = 2 +* +* Print information about the tests that did not pass +* the threshold. +* + DO 110 K = 1, NT + IF( RESULT( K ).GE.THRESH ) THEN + IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) + $ CALL ALADHD( NOUT, PATH ) + WRITE( NOUT, FMT = 9999 )'ZSYSV_AA ', + $ UPLO, N, IMAT, K, RESULT( K ) + NFAIL = NFAIL + 1 + END IF + 110 CONTINUE + NRUN = NRUN + NT + 120 CONTINUE + END IF +* + 150 CONTINUE +* + 160 CONTINUE + 170 CONTINUE + 180 CONTINUE +* +* Print a summary of the results. +* + CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS ) +* + 9999 FORMAT( 1X, A, ', UPLO=''', A1, ''', N =', I5, ', type ', I2, + $ ', test ', I2, ', ratio =', G12.5 ) + RETURN +* +* End of ZDRVSY_AA +* + END diff --git a/TESTING/LIN/zerrhe.f b/TESTING/LIN/zerrhe.f index b6304b1c..d15b2cf4 100644 --- a/TESTING/LIN/zerrhe.f +++ b/TESTING/LIN/zerrhe.f @@ -93,7 +93,8 @@ $ ZHETRF_RK, ZHETRF_ROOK, ZHETRF_AA, ZHETRI, $ ZHETRI_3, ZHETRI_3X, ZHETRI_ROOK, ZHETRI2, $ ZHETRI2X, ZHETRS, ZHETRS_3, ZHETRS_ROOK, - $ ZHETRS_AA, ZHPCON, ZHPRFS, ZHPTRF, ZHPTRI, ZHPTRS + $ ZHETRS_AA, ZHPCON, ZHPRFS, ZHPTRF, ZHPTRI, + $ ZHPTRS * .. * .. Scalars in Common .. LOGICAL LERR, OK @@ -489,6 +490,12 @@ INFOT = 4 CALL ZHETRF_AA( 'U', 2, A, 1, IP, W, 4, INFO ) CALL CHKXER( 'ZHETRF_AA', INFOT, NOUT, LERR, OK ) + INFOT = 7 + CALL ZHETRF_AA( 'U', 0, A, 1, IP, W, 0, INFO ) + CALL CHKXER( 'ZHETRF_AA', INFOT, NOUT, LERR, OK ) + INFOT = 7 + CALL ZHETRF_AA( 'U', 0, A, 1, IP, W, -2, INFO ) + CALL CHKXER( 'ZHETRF_AA', INFOT, NOUT, LERR, OK ) * * ZHETRS_AA * @@ -508,6 +515,12 @@ INFOT = 8 CALL ZHETRS_AA( 'U', 2, 1, A, 2, IP, B, 1, W, 1, INFO ) CALL CHKXER( 'ZHETRS_AA', INFOT, NOUT, LERR, OK ) + INFOT = 10 + CALL ZHETRS_AA( 'U', 0, 1, A, 1, IP, B, 1, W, 0, INFO ) + CALL CHKXER( 'ZHETRS_AA', INFOT, NOUT, LERR, OK ) + INFOT = 10 + CALL ZHETRS_AA( 'U', 0, 1, A, 1, IP, B, 1, W, -2, INFO ) + CALL CHKXER( 'ZHETRS_AA', INFOT, NOUT, LERR, OK ) * ELSE IF( LSAMEN( 2, C2, 'HP' ) ) THEN * diff --git a/TESTING/LIN/zerrsy.f b/TESTING/LIN/zerrsy.f index 45e5f0c0..eb8bb628 100644 --- a/TESTING/LIN/zerrsy.f +++ b/TESTING/LIN/zerrsy.f @@ -549,6 +549,50 @@ INFOT = 5 CALL ZSPCON( 'U', 1, A, IP, -ANRM, RCOND, W, INFO ) CALL CHKXER( 'ZSPCON', INFOT, NOUT, LERR, OK ) +* + ELSE IF( LSAMEN( 2, C2, 'SA' ) ) THEN +* +* Test error exits of the routines that use factorization +* of a symmetric indefinite matrix with Aasen's algorithm. +* +* ZSYTRF_AA +* + SRNAMT = 'ZSYTRF_AA' + INFOT = 1 + CALL ZSYTRF_AA( '/', 0, A, 1, IP, W, 1, INFO ) + CALL CHKXER( 'ZSYTRF_AA', INFOT, NOUT, LERR, OK ) + INFOT = 2 + CALL ZSYTRF_AA( 'U', -1, A, 1, IP, W, 1, INFO ) + CALL CHKXER( 'ZSYTRF_AA', INFOT, NOUT, LERR, OK ) + INFOT = 4 + CALL ZSYTRF_AA( 'U', 2, A, 1, IP, W, 4, INFO ) + CALL CHKXER( 'ZSYTRF_AA', INFOT, NOUT, LERR, OK ) + INFOT = 7 + CALL ZSYTRF_AA( 'U', 0, A, 1, IP, W, 0, INFO ) + CALL CHKXER( 'ZSYTRF_AA', INFOT, NOUT, LERR, OK ) + INFOT = 7 + CALL ZSYTRF_AA( 'U', 0, A, 1, IP, W, -2, INFO ) + CALL CHKXER( 'ZSYTRF_AA', INFOT, NOUT, LERR, OK ) +* +* ZSYTRS_AA +* + SRNAMT = 'ZSYTRS_AA' + INFOT = 1 + CALL ZSYTRS_AA( '/', 0, 0, A, 1, IP, B, 1, W, 1, INFO ) + CALL CHKXER( 'ZSYTRS_AA', INFOT, NOUT, LERR, OK ) + INFOT = 2 + CALL ZSYTRS_AA( 'U', -1, 0, A, 1, IP, B, 1, W, 1, INFO ) + CALL CHKXER( 'ZSYTRS_AA', INFOT, NOUT, LERR, OK ) + INFOT = 3 + CALL ZSYTRS_AA( 'U', 0, -1, A, 1, IP, B, 1, W, 1, INFO ) + CALL CHKXER( 'ZSYTRS_AA', INFOT, NOUT, LERR, OK ) + INFOT = 5 + CALL ZSYTRS_AA( 'U', 2, 1, A, 1, IP, B, 2, W, 1, INFO ) + CALL CHKXER( 'ZSYTRS_AA', INFOT, NOUT, LERR, OK ) + INFOT = 8 + CALL ZSYTRS_AA( 'U', 2, 1, A, 2, IP, B, 1, W, 1, INFO ) + CALL CHKXER( 'ZSYTRS_AA', INFOT, NOUT, LERR, OK ) +* END IF * * Print a summary line. diff --git a/TESTING/LIN/zhet01_aa.f b/TESTING/LIN/zhet01_aa.f index d1328c88..c4734fcd 100644 --- a/TESTING/LIN/zhet01_aa.f +++ b/TESTING/LIN/zhet01_aa.f @@ -9,17 +9,17 @@ * =========== * * SUBROUTINE ZHET01_AA( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, -* C, LDC, RWORK, RESID ) +* C, LDC, RWORK, RESID ) * * .. Scalar Arguments .. * CHARACTER UPLO * INTEGER LDA, LDAFAC, LDC, N -* COMPLEX*16 RESID +* DOUBLE PRECISION RESID * .. * .. Array Arguments .. * INTEGER IPIV( * ) -* COMPLEX*16 A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ), -* $ RWORK( * ) +* DOUBLE PRECISION RWORK( * ) +* COMPLEX*16 A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ) * .. * * @@ -123,7 +123,7 @@ * * ===================================================================== SUBROUTINE ZHET01_AA( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, C, - $ LDC, RWORK, RESID ) + $ LDC, RWORK, RESID ) * * -- LAPACK test routine (version 3.7.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -137,8 +137,8 @@ * .. * .. Array Arguments .. INTEGER IPIV( * ) - COMPLEX*16 A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ), - $ RWORK( * ) + DOUBLE PRECISION RWORK( * ) + COMPLEX*16 A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ) * .. * * ===================================================================== @@ -197,40 +197,42 @@ $ LDC+1 ) CALL ZLACGV( N-1, C( 1, 2 ), LDC+1 ) ENDIF - ENDIF * -* Call ZTRMM to form the product U' * D (or L * D ). +* Call ZTRMM to form the product U' * D (or L * D ). * - IF( LSAME( UPLO, 'U' ) ) THEN - CALL ZTRMM( 'Left', UPLO, 'Conjugate transpose', 'Unit', N-1, - $ N, CONE, AFAC( 1, 2 ), LDAFAC, C( 2, 1 ), LDC ) - ELSE - CALL ZTRMM( 'Left', UPLO, 'No transpose', 'Unit', N-1, N, - $ CONE, AFAC( 2, 1 ), LDAFAC, C( 2, 1 ), LDC ) - END IF + IF( LSAME( UPLO, 'U' ) ) THEN + CALL ZTRMM( 'Left', UPLO, 'Conjugate transpose', 'Unit', + $ N-1, N, CONE, AFAC( 1, 2 ), LDAFAC, C( 2, 1 ), + $ LDC ) + ELSE + CALL ZTRMM( 'Left', UPLO, 'No transpose', 'Unit', N-1, N, + $ CONE, AFAC( 2, 1 ), LDAFAC, C( 2, 1 ), LDC ) + END IF * -* Call ZTRMM again to multiply by U (or L ). +* Call ZTRMM again to multiply by U (or L ). * - IF( LSAME( UPLO, 'U' ) ) THEN - CALL ZTRMM( 'Right', UPLO, 'No transpose', 'Unit', N, N-1, - $ CONE, AFAC( 1, 2 ), LDAFAC, C( 1, 2 ), LDC ) - ELSE - CALL ZTRMM( 'Right', UPLO, 'Conjugate transpose', 'Unit', N, - $ N-1, CONE, AFAC( 2, 1 ), LDAFAC, C( 1, 2 ), LDC ) - END IF + IF( LSAME( UPLO, 'U' ) ) THEN + CALL ZTRMM( 'Right', UPLO, 'No transpose', 'Unit', N, N-1, + $ CONE, AFAC( 1, 2 ), LDAFAC, C( 1, 2 ), LDC ) + ELSE + CALL ZTRMM( 'Right', UPLO, 'Conjugate transpose', 'Unit', N, + $ N-1, CONE, AFAC( 2, 1 ), LDAFAC, C( 1, 2 ), + $ LDC ) + END IF +* +* Apply hermitian pivots * -* Apply hermitian pivots -* - DO J = N, 1, -1 - I = IPIV( J ) - IF( I.NE.J ) - $ CALL ZSWAP( N, C( J, 1 ), LDC, C( I, 1 ), LDC ) - END DO - DO J = N, 1, -1 - I = IPIV( J ) - IF( I.NE.J ) - $ CALL ZSWAP( N, C( 1, J ), 1, C( 1, I ), 1 ) - END DO + DO J = N, 1, -1 + I = IPIV( J ) + IF( I.NE.J ) + $ CALL ZSWAP( N, C( J, 1 ), LDC, C( I, 1 ), LDC ) + END DO + DO J = N, 1, -1 + I = IPIV( J ) + IF( I.NE.J ) + $ CALL ZSWAP( N, C( 1, J ), 1, C( 1, I ), 1 ) + END DO + ENDIF * * * Compute the difference C - A . diff --git a/TESTING/LIN/zsyt01_aa.f b/TESTING/LIN/zsyt01_aa.f new file mode 100644 index 00000000..988f4beb --- /dev/null +++ b/TESTING/LIN/zsyt01_aa.f @@ -0,0 +1,265 @@ +*> \brief \b ZSYT01 +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE ZSYT01_AA( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, C, LDC, +* RWORK, RESID ) +* +* .. Scalar Arguments .. +* CHARACTER UPLO +* INTEGER LDA, LDAFAC, LDC, N +* DOUBLE PRECISION RESID +* .. +* .. Array Arguments .. +* INTEGER IPIV( * ) +* DOUBLE PRECISION RWORK( * ) +* COMPLEX*16 A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ), +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> ZSYT01 reconstructs a hermitian indefinite matrix A from its +*> block L*D*L' or U*D*U' factorization and computes the residual +*> norm( C - A ) / ( N * norm(A) * EPS ), +*> where C is the reconstructed matrix and EPS is the machine epsilon. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] UPLO +*> \verbatim +*> UPLO is CHARACTER*1 +*> Specifies whether the upper or lower triangular part of the +*> hermitian matrix A is stored: +*> = 'U': Upper triangular +*> = 'L': Lower triangular +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of rows and columns of the matrix A. N >= 0. +*> \endverbatim +*> +*> \param[in] A +*> \verbatim +*> A is COMPLEX*16 array, dimension (LDA,N) +*> The original hermitian matrix A. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N) +*> \endverbatim +*> +*> \param[in] AFAC +*> \verbatim +*> AFAC is COMPLEX*16 array, dimension (LDAFAC,N) +*> The factored form of the matrix A. AFAC contains the block +*> diagonal matrix D and the multipliers used to obtain the +*> factor L or U from the block L*D*L' or U*D*U' factorization +*> as computed by ZSYTRF. +*> \endverbatim +*> +*> \param[in] LDAFAC +*> \verbatim +*> LDAFAC is INTEGER +*> The leading dimension of the array AFAC. LDAFAC >= max(1,N). +*> \endverbatim +*> +*> \param[in] IPIV +*> \verbatim +*> IPIV is INTEGER array, dimension (N) +*> The pivot indices from ZSYTRF. +*> \endverbatim +*> +*> \param[out] C +*> \verbatim +*> C is COMPLEX*16 array, dimension (LDC,N) +*> \endverbatim +*> +*> \param[in] LDC +*> \verbatim +*> LDC is INTEGER +*> The leading dimension of the array C. LDC >= max(1,N). +*> \endverbatim +*> +*> \param[out] RWORK +*> \verbatim +*> RWORK is COMPLEX*16 array, dimension (N) +*> \endverbatim +*> +*> \param[out] RESID +*> \verbatim +*> RESID is COMPLEX*16 +*> If UPLO = 'L', norm(L*D*L' - A) / ( N * norm(A) * EPS ) +*> If UPLO = 'U', norm(U*D*U' - A) / ( N * norm(A) * EPS ) +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2016 +* +* @generated from LIN/dsyt01_aa.f, fortran d -> z, Thu Nov 17 13:01:50 2016 +* +*> \ingroup complex16_lin +* +* ===================================================================== + SUBROUTINE ZSYT01_AA( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, C, + $ LDC, RWORK, RESID ) +* +* -- LAPACK test routine (version 3.5.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2016 +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER LDA, LDAFAC, LDC, N + DOUBLE PRECISION RESID +* .. +* .. Array Arguments .. + INTEGER IPIV( * ) + COMPLEX*16 A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ) + DOUBLE PRECISION RWORK( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) + COMPLEX*16 CZERO, CONE + PARAMETER ( CZERO = 0.0E+0, CONE = 1.0E+0 ) +* .. +* .. Local Scalars .. + INTEGER I, J + DOUBLE PRECISION ANORM, EPS +* .. +* .. External Functions .. + LOGICAL LSAME + DOUBLE PRECISION DLAMCH, ZLANSY + EXTERNAL LSAME, DLAMCH, ZLANSY +* .. +* .. External Subroutines .. + EXTERNAL ZLASET, ZLAVSY +* .. +* .. Intrinsic Functions .. + INTRINSIC DBLE +* .. +* .. Executable Statements .. +* +* Quick exit if N = 0. +* + IF( N.LE.0 ) THEN + RESID = ZERO + RETURN + END IF +* +* Determine EPS and the norm of A. +* + EPS = DLAMCH( 'Epsilon' ) + ANORM = ZLANSY( '1', UPLO, N, A, LDA, RWORK ) +* +* Initialize C to the tridiagonal matrix T. +* + CALL ZLASET( 'Full', N, N, CZERO, CZERO, C, LDC ) + CALL ZLACPY( 'F', 1, N, AFAC( 1, 1 ), LDAFAC+1, C( 1, 1 ), LDC+1 ) + IF( N.GT.1 ) THEN + IF( LSAME( UPLO, 'U' ) ) THEN + CALL ZLACPY( 'F', 1, N-1, AFAC( 1, 2 ), LDAFAC+1, C( 1, 2 ), + $ LDC+1 ) + CALL ZLACPY( 'F', 1, N-1, AFAC( 1, 2 ), LDAFAC+1, C( 2, 1 ), + $ LDC+1 ) + ELSE + CALL ZLACPY( 'F', 1, N-1, AFAC( 2, 1 ), LDAFAC+1, C( 1, 2 ), + $ LDC+1 ) + CALL ZLACPY( 'F', 1, N-1, AFAC( 2, 1 ), LDAFAC+1, C( 2, 1 ), + $ LDC+1 ) + ENDIF +* +* Call ZTRMM to form the product U' * D (or L * D ). +* + IF( LSAME( UPLO, 'U' ) ) THEN + CALL ZTRMM( 'Left', UPLO, 'Transpose', 'Unit', N-1, N, + $ CONE, AFAC( 1, 2 ), LDAFAC, C( 2, 1 ), LDC ) + ELSE + CALL ZTRMM( 'Left', UPLO, 'No transpose', 'Unit', N-1, N, + $ CONE, AFAC( 2, 1 ), LDAFAC, C( 2, 1 ), LDC ) + END IF +* +* Call ZTRMM again to multiply by U (or L ). +* + IF( LSAME( UPLO, 'U' ) ) THEN + CALL ZTRMM( 'Right', UPLO, 'No transpose', 'Unit', N, N-1, + $ CONE, AFAC( 1, 2 ), LDAFAC, C( 1, 2 ), LDC ) + ELSE + CALL ZTRMM( 'Right', UPLO, 'Transpose', 'Unit', N, N-1, + $ CONE, AFAC( 2, 1 ), LDAFAC, C( 1, 2 ), LDC ) + END IF + ENDIF +* +* Apply symmetric pivots +* + DO J = N, 1, -1 + I = IPIV( J ) + IF( I.NE.J ) + $ CALL ZSWAP( N, C( J, 1 ), LDC, C( I, 1 ), LDC ) + END DO + DO J = N, 1, -1 + I = IPIV( J ) + IF( I.NE.J ) + $ CALL ZSWAP( N, C( 1, J ), 1, C( 1, I ), 1 ) + END DO +* +* +* Compute the difference C - A . +* + IF( LSAME( UPLO, 'U' ) ) THEN + DO J = 1, N + DO I = 1, J + C( I, J ) = C( I, J ) - A( I, J ) + END DO + END DO + ELSE + DO J = 1, N + DO I = J, N + C( I, J ) = C( I, J ) - A( I, J ) + END DO + END DO + END IF +* +* Compute norm( C - A ) / ( N * norm(A) * EPS ) +* + RESID = ZLANSY( '1', UPLO, N, C, LDC, RWORK ) +* + IF( ANORM.LE.ZERO ) THEN + IF( RESID.NE.ZERO ) + $ RESID = ONE / EPS + ELSE + RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS + END IF +* + RETURN +* +* End of ZSYT01 +* + END |