From e5c236b03ae6f2c2af73a3331531825ca6b16e70 Mon Sep 17 00:00:00 2001 From: langou Date: Sun, 15 Nov 2015 16:49:46 +0000 Subject: [PATCH] integration of xPOTRF2 --- SRC/Makefile | 5 +- SRC/cpotrf.f | 8 +- SRC/cpotrf2.f | 246 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SRC/dpotrf.f | 8 +- SRC/dpotrf2.f | 237 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ SRC/spotrf.f | 8 +- SRC/spotrf2.f | 237 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ SRC/zpotrf.f | 8 +- SRC/zpotrf2.f | 241 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 9 files changed, 981 insertions(+), 17 deletions(-) create mode 100644 SRC/cpotrf2.f create mode 100644 SRC/dpotrf2.f create mode 100644 SRC/spotrf2.f create mode 100644 SRC/zpotrf2.f diff --git a/SRC/Makefile b/SRC/Makefile index 671b457..d78c6cf 100644 --- a/SRC/Makefile +++ b/SRC/Makefile @@ -97,7 +97,7 @@ DZLAUX = \ ../INSTALL/dlamch.o ../INSTALL/dsecnd_$(TIMER).o SLASRC = \ - sbdsvdx.o \ + sbdsvdx.o spotrf2.o \ sgbbrd.o sgbcon.o sgbequ.o sgbrfs.o sgbsv.o \ sgbsvx.o sgbtf2.o sgbtrf.o sgbtrs.o sgebak.o sgebal.o sgebd2.o \ sgebrd.o sgecon.o sgeequ.o sgees.o sgeesx.o sgeev.o sgeevx.o \ @@ -174,6 +174,7 @@ SXLASRC = sgesvxx.o sgerfsx.o sla_gerfsx_extended.o sla_geamv.o \ endif CLASRC = \ + cpotrf2.o \ cbdsqr.o cgbbrd.o cgbcon.o cgbequ.o cgbrfs.o cgbsv.o cgbsvx.o \ cgbtf2.o cgbtrf.o cgbtrs.o cgebak.o cgebal.o cgebd2.o cgebrd.o \ cgecon.o cgeequ.o cgees.o cgeesx.o cgeev.o cgeevx.o \ @@ -261,6 +262,7 @@ endif ZCLASRC = cpotrs.o cgetrs.o cpotrf.o cgetrf.o DLASRC = \ + dpotrf2.o \ dbdsvdx.o \ dgbbrd.o dgbcon.o dgbequ.o dgbrfs.o dgbsv.o \ dgbsvx.o dgbtf2.o dgbtrf.o dgbtrs.o dgebak.o dgebal.o dgebd2.o \ @@ -337,6 +339,7 @@ DXLASRC = dgesvxx.o dgerfsx.o dla_gerfsx_extended.o dla_geamv.o \ endif ZLASRC = \ + zpotrf2.o \ zbdsqr.o zgbbrd.o zgbcon.o zgbequ.o zgbrfs.o zgbsv.o zgbsvx.o \ zgbtf2.o zgbtrf.o zgbtrs.o zgebak.o zgebal.o zgebd2.o zgebrd.o \ zgecon.o zgeequ.o zgees.o zgeesx.o zgeev.o zgeevx.o \ diff --git a/SRC/cpotrf.f b/SRC/cpotrf.f index 106652b..ec91495 100644 --- a/SRC/cpotrf.f +++ b/SRC/cpotrf.f @@ -137,7 +137,7 @@ EXTERNAL LSAME, ILAENV * .. * .. External Subroutines .. - EXTERNAL CGEMM, CHERK, CPOTF2, CTRSM, XERBLA + EXTERNAL CGEMM, CHERK, CPOTRF2, CTRSM, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN @@ -172,7 +172,7 @@ * * Use unblocked code. * - CALL CPOTF2( UPLO, N, A, LDA, INFO ) + CALL CPOTRF2( UPLO, N, A, LDA, INFO ) ELSE * * Use blocked code. @@ -189,7 +189,7 @@ JB = MIN( NB, N-J+1 ) CALL CHERK( 'Upper', 'Conjugate transpose', JB, J-1, $ -ONE, A( 1, J ), LDA, ONE, A( J, J ), LDA ) - CALL CPOTF2( 'Upper', JB, A( J, J ), LDA, INFO ) + CALL CPOTRF2( 'Upper', JB, A( J, J ), LDA, INFO ) IF( INFO.NE.0 ) $ GO TO 30 IF( J+JB.LE.N ) THEN @@ -218,7 +218,7 @@ JB = MIN( NB, N-J+1 ) CALL CHERK( 'Lower', 'No transpose', JB, J-1, -ONE, $ A( J, 1 ), LDA, ONE, A( J, J ), LDA ) - CALL CPOTF2( 'Lower', JB, A( J, J ), LDA, INFO ) + CALL CPOTRF2( 'Lower', JB, A( J, J ), LDA, INFO ) IF( INFO.NE.0 ) $ GO TO 30 IF( J+JB.LE.N ) THEN diff --git a/SRC/cpotrf2.f b/SRC/cpotrf2.f new file mode 100644 index 0000000..6ab06a6 --- /dev/null +++ b/SRC/cpotrf2.f @@ -0,0 +1,246 @@ +*> \brief \b CPOTRF2 +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* RECURSIVE SUBROUTINE CPOTRF2( UPLO, N, A, LDA, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER UPLO +* INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. +* COMPLEX A( LDA, * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> CPOTRF2 computes the Cholesky factorization of a real symmetric +*> positive definite matrix A using the recursive algorithm. +*> +*> The factorization has the form +*> A = U**H * U, if UPLO = 'U', or +*> A = L * L**H, if UPLO = 'L', +*> where U is an upper triangular matrix and L is lower triangular. +*> +*> This is the recursive version of the algorithm. It divides +*> the matrix into four submatrices: +*> +*> [ A11 | A12 ] where A11 is n1 by n1 and A22 is n2 by n2 +*> A = [ -----|----- ] with n1 = n/2 +*> [ A21 | A22 ] n2 = n-n1 +*> +*> The subroutine calls itself to factor A11. Update and scale A21 +*> or A12, update A22 then calls itself to factor A22. +*> +*> \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 DOUBLE PRECISION 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 factor U or L from the Cholesky +*> factorization A = U**H*U or A = L*L**H. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N). +*> \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, the leading minor of order i is not +*> positive definite, and the factorization could not be +*> completed. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2015 +* +*> \ingroup complexPOcomputational +* +* ===================================================================== + RECURSIVE SUBROUTINE CPOTRF2( UPLO, N, A, LDA, INFO ) +* +* -- LAPACK computational routine (version 3.6.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2015 +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. + COMPLEX A( LDA, * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + REAL ONE, ZERO + PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) + COMPLEX CONE + PARAMETER ( CONE = (1.0E+0, 0.0E+0) ) +* .. +* .. Local Scalars .. + LOGICAL UPPER + INTEGER N1, N2, IINFO + REAL AJJ +* .. +* .. External Functions .. + LOGICAL LSAME, SISNAN + EXTERNAL LSAME, SISNAN +* .. +* .. External Subroutines .. + EXTERNAL CHERK, CTRSM, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, REAL, SQRT +* .. +* .. Executable Statements .. +* +* Test the input parameters +* + INFO = 0 + UPPER = LSAME( UPLO, 'U' ) + 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 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'CPOTRF2', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 ) + $ RETURN +* +* N=1 case +* + IF( N.EQ.1 ) THEN +* +* Test for non-positive-definiteness +* + AJJ = REAL( A( 1, 1 ) ) + IF( AJJ.LE.ZERO.OR.SISNAN( AJJ ) ) THEN + INFO = 1 + RETURN + END IF +* +* Factor +* + A( 1, 1 ) = SQRT( AJJ ) +* +* Use recursive code +* + ELSE + N1 = N/2 + N2 = N-N1 +* +* Factor A11 +* + CALL CPOTRF2( UPLO, N1, A( 1, 1 ), LDA, IINFO ) + IF ( IINFO.NE.0 ) THEN + INFO = IINFO + RETURN + END IF +* +* Compute the Cholesky factorization A = U**H*U +* + IF( UPPER ) THEN +* +* Update and scale A12 +* + CALL CTRSM( 'L', 'U', 'C', 'N', N1, N2, CONE, + $ A( 1, 1 ), LDA, A( 1, N1+1 ), LDA ) +* +* Update and factor A22 +* + CALL CHERK( UPLO, 'C', N2, N1, -ONE, A( 1, N1+1 ), LDA, + $ ONE, A( N1+1, N1+1 ), LDA ) +* + CALL CPOTRF2( UPLO, N2, A( N1+1, N1+1 ), LDA, IINFO ) +* + IF ( IINFO.NE.0 ) THEN + INFO = IINFO + N1 + RETURN + END IF +* +* Compute the Cholesky factorization A = L*L**H +* + ELSE +* +* Update and scale A21 +* + CALL CTRSM( 'R', 'L', 'C', 'N', N2, N1, CONE, + $ A( 1, 1 ), LDA, A( N1+1, 1 ), LDA ) +* +* Update and factor A22 +* + CALL CHERK( UPLO, 'N', N2, N1, -ONE, A( N1+1, 1 ), LDA, + $ ONE, A( N1+1, N1+1 ), LDA ) +* + CALL CPOTRF2( UPLO, N2, A( N1+1, N1+1 ), LDA, IINFO ) +* + IF ( IINFO.NE.0 ) THEN + INFO = IINFO + N1 + RETURN + END IF +* + END IF + END IF + RETURN +* +* End of CPOTRF2 +* + END diff --git a/SRC/dpotrf.f b/SRC/dpotrf.f index 3457230..6e4865c 100644 --- a/SRC/dpotrf.f +++ b/SRC/dpotrf.f @@ -136,7 +136,7 @@ EXTERNAL LSAME, ILAENV * .. * .. External Subroutines .. - EXTERNAL DGEMM, DPOTF2, DSYRK, DTRSM, XERBLA + EXTERNAL DGEMM, DPOTRF2, DSYRK, DTRSM, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN @@ -171,7 +171,7 @@ * * Use unblocked code. * - CALL DPOTF2( UPLO, N, A, LDA, INFO ) + CALL DPOTRF2( UPLO, N, A, LDA, INFO ) ELSE * * Use blocked code. @@ -188,7 +188,7 @@ JB = MIN( NB, N-J+1 ) CALL DSYRK( 'Upper', 'Transpose', JB, J-1, -ONE, $ A( 1, J ), LDA, ONE, A( J, J ), LDA ) - CALL DPOTF2( 'Upper', JB, A( J, J ), LDA, INFO ) + CALL DPOTRF2( 'Upper', JB, A( J, J ), LDA, INFO ) IF( INFO.NE.0 ) $ GO TO 30 IF( J+JB.LE.N ) THEN @@ -216,7 +216,7 @@ JB = MIN( NB, N-J+1 ) CALL DSYRK( 'Lower', 'No transpose', JB, J-1, -ONE, $ A( J, 1 ), LDA, ONE, A( J, J ), LDA ) - CALL DPOTF2( 'Lower', JB, A( J, J ), LDA, INFO ) + CALL DPOTRF2( 'Lower', JB, A( J, J ), LDA, INFO ) IF( INFO.NE.0 ) $ GO TO 30 IF( J+JB.LE.N ) THEN diff --git a/SRC/dpotrf2.f b/SRC/dpotrf2.f new file mode 100644 index 0000000..751ff76 --- /dev/null +++ b/SRC/dpotrf2.f @@ -0,0 +1,237 @@ +*> \brief \b DPOTRF2 +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* RECURSIVE SUBROUTINE DPOTRF2( UPLO, N, A, LDA, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER UPLO +* INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. +* REAL A( LDA, * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DPOTRF2 computes the Cholesky factorization of a real symmetric +*> positive definite matrix A using the recursive algorithm. +*> +*> The factorization has the form +*> A = U**T * U, if UPLO = 'U', or +*> A = L * L**T, if UPLO = 'L', +*> where U is an upper triangular matrix and L is lower triangular. +*> +*> This is the recursive version of the algorithm. It divides +*> the matrix into four submatrices: +*> +*> [ A11 | A12 ] where A11 is n1 by n1 and A22 is n2 by n2 +*> A = [ -----|----- ] with n1 = n/2 +*> [ A21 | A22 ] n2 = n-n1 +*> +*> The subroutine calls itself to factor A11. Update and scale A21 +*> or A12, update A22 then calls itself to factor A22. +*> +*> \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 DOUBLE PRECISION 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 factor U or L from the Cholesky +*> factorization A = U**T*U or A = L*L**T. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N). +*> \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, the leading minor of order i is not +*> positive definite, and the factorization could not be +*> completed. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2015 +* +*> \ingroup doublePOcomputational +* +* ===================================================================== + RECURSIVE SUBROUTINE DPOTRF2( UPLO, N, A, LDA, INFO ) +* +* -- LAPACK computational routine (version 3.6.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2015 +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. + DOUBLE PRECISION A( LDA, * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ONE, ZERO + PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) +* .. +* .. Local Scalars .. + LOGICAL UPPER + INTEGER N1, N2, IINFO +* .. +* .. External Functions .. + LOGICAL LSAME, DISNAN + EXTERNAL LSAME, DISNAN +* .. +* .. External Subroutines .. + EXTERNAL DSYRK, DTRSM, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, SQRT +* .. +* .. Executable Statements .. +* +* Test the input parameters +* + INFO = 0 + UPPER = LSAME( UPLO, 'U' ) + 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 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DPOTRF2', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 ) + $ RETURN +* +* N=1 case +* + IF( N.EQ.1 ) THEN +* +* Test for non-positive-definiteness +* + IF( A( 1, 1 ).LE.ZERO.OR.DISNAN( A( 1, 1 ) ) ) THEN + INFO = 1 + RETURN + END IF +* +* Factor +* + A( 1, 1 ) = SQRT( A( 1, 1 ) ) +* +* Use recursive code +* + ELSE + N1 = N/2 + N2 = N-N1 +* +* Factor A11 +* + CALL DPOTRF2( UPLO, N1, A( 1, 1 ), LDA, IINFO ) + IF ( IINFO.NE.0 ) THEN + INFO = IINFO + RETURN + END IF +* +* Compute the Cholesky factorization A = U**T*U +* + IF( UPPER ) THEN +* +* Update and scale A12 +* + CALL DTRSM( 'L', 'U', 'T', 'N', N1, N2, ONE, + $ A( 1, 1 ), LDA, A( 1, N1+1 ), LDA ) +* +* Update and factor A22 +* + CALL DSYRK( UPLO, 'T', N2, N1, -ONE, A( 1, N1+1 ), LDA, + $ ONE, A( N1+1, N1+1 ), LDA ) + CALL DPOTRF2( UPLO, N2, A( N1+1, N1+1 ), LDA, IINFO ) + IF ( IINFO.NE.0 ) THEN + INFO = IINFO + N1 + RETURN + END IF +* +* Compute the Cholesky factorization A = L*L**T +* + ELSE +* +* Update and scale A21 +* + CALL DTRSM( 'R', 'L', 'T', 'N', N2, N1, ONE, + $ A( 1, 1 ), LDA, A( N1+1, 1 ), LDA ) +* +* Update and factor A22 +* + CALL DSYRK( UPLO, 'N', N2, N1, -ONE, A( N1+1, 1 ), LDA, + $ ONE, A( N1+1, N1+1 ), LDA ) + CALL DPOTRF2( UPLO, N2, A( N1+1, N1+1 ), LDA, IINFO ) + IF ( IINFO.NE.0 ) THEN + INFO = IINFO + N1 + RETURN + END IF + END IF + END IF + RETURN +* +* End of DPOTRF2 +* + END diff --git a/SRC/spotrf.f b/SRC/spotrf.f index f9ea451..6a1f0c1 100644 --- a/SRC/spotrf.f +++ b/SRC/spotrf.f @@ -136,7 +136,7 @@ EXTERNAL LSAME, ILAENV * .. * .. External Subroutines .. - EXTERNAL SGEMM, SPOTF2, SSYRK, STRSM, XERBLA + EXTERNAL SGEMM, SPOTRF2, SSYRK, STRSM, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN @@ -171,7 +171,7 @@ * * Use unblocked code. * - CALL SPOTF2( UPLO, N, A, LDA, INFO ) + CALL SPOTRF2( UPLO, N, A, LDA, INFO ) ELSE * * Use blocked code. @@ -188,7 +188,7 @@ JB = MIN( NB, N-J+1 ) CALL SSYRK( 'Upper', 'Transpose', JB, J-1, -ONE, $ A( 1, J ), LDA, ONE, A( J, J ), LDA ) - CALL SPOTF2( 'Upper', JB, A( J, J ), LDA, INFO ) + CALL SPOTRF2( 'Upper', JB, A( J, J ), LDA, INFO ) IF( INFO.NE.0 ) $ GO TO 30 IF( J+JB.LE.N ) THEN @@ -216,7 +216,7 @@ JB = MIN( NB, N-J+1 ) CALL SSYRK( 'Lower', 'No transpose', JB, J-1, -ONE, $ A( J, 1 ), LDA, ONE, A( J, J ), LDA ) - CALL SPOTF2( 'Lower', JB, A( J, J ), LDA, INFO ) + CALL SPOTRF2( 'Lower', JB, A( J, J ), LDA, INFO ) IF( INFO.NE.0 ) $ GO TO 30 IF( J+JB.LE.N ) THEN diff --git a/SRC/spotrf2.f b/SRC/spotrf2.f new file mode 100644 index 0000000..dfdf16e --- /dev/null +++ b/SRC/spotrf2.f @@ -0,0 +1,237 @@ +*> \brief \b SPOTRF2 +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* RECURSIVE SUBROUTINE SPOTRF2( UPLO, N, A, LDA, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER UPLO +* INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. +* REAL A( LDA, * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> SPOTRF2 computes the Cholesky factorization of a real symmetric +*> positive definite matrix A using the recursive algorithm. +*> +*> The factorization has the form +*> A = U**T * U, if UPLO = 'U', or +*> A = L * L**T, if UPLO = 'L', +*> where U is an upper triangular matrix and L is lower triangular. +*> +*> This is the recursive version of the algorithm. It divides +*> the matrix into four submatrices: +*> +*> [ A11 | A12 ] where A11 is n1 by n1 and A22 is n2 by n2 +*> A = [ -----|----- ] with n1 = n/2 +*> [ A21 | A22 ] n2 = n-n1 +*> +*> The subroutine calls itself to factor A11. Update and scale A21 +*> or A12, update A22 then call itself to factor A22. +*> +*> \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, if INFO = 0, the factor U or L from the Cholesky +*> factorization A = U**T*U or A = L*L**T. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N). +*> \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, the leading minor of order i is not +*> positive definite, and the factorization could not be +*> completed. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2015 +* +*> \ingroup realPOcomputational +* +* ===================================================================== + RECURSIVE SUBROUTINE SPOTRF2( UPLO, N, A, LDA, INFO ) +* +* -- LAPACK computational routine (version 3.6.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2015 +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. + REAL A( LDA, * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + REAL ONE, ZERO + PARAMETER ( ONE = 1.0E+0, ZERO=0.0E+0 ) +* .. +* .. Local Scalars .. + LOGICAL UPPER + INTEGER N1, N2, IINFO +* .. +* .. External Functions .. + LOGICAL LSAME, SISNAN + EXTERNAL LSAME, SISNAN +* .. +* .. External Subroutines .. + EXTERNAL SGEMM, SSYRK, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, SQRT +* .. +* .. Executable Statements .. +* +* Test the input parameters +* + INFO = 0 + UPPER = LSAME( UPLO, 'U' ) + 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 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'SPOTRF2', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 ) + $ RETURN +* +* N=1 case +* + IF( N.EQ.1 ) THEN +* +* Test for non-positive-definiteness +* + IF( A( 1, 1 ).LE.ZERO.OR.SISNAN( A( 1, 1 ) ) ) THEN + INFO = 1 + RETURN + END IF +* +* Factor +* + A( 1, 1 ) = SQRT( A( 1, 1 ) ) +* +* Use recursive code +* + ELSE + N1 = N/2 + N2 = N-N1 +* +* Factor A11 +* + CALL SPOTRF2( UPLO, N1, A( 1, 1 ), LDA, IINFO ) + IF ( IINFO.NE.0 ) THEN + INFO = IINFO + RETURN + END IF +* +* Compute the Cholesky factorization A = U**T*U +* + IF( UPPER ) THEN +* +* Update and scale A12 +* + CALL STRSM( 'L', 'U', 'T', 'N', N1, N2, ONE, + $ A( 1, 1 ), LDA, A( 1, N1+1 ), LDA ) +* +* Update and factor A22 +* + CALL SSYRK( UPLO, 'T', N2, N1, -ONE, A( 1, N1+1 ), LDA, + $ ONE, A( N1+1, N1+1 ), LDA ) + CALL SPOTRF2( UPLO, N2, A( N1+1, N1+1 ), LDA, IINFO ) + IF ( IINFO.NE.0 ) THEN + INFO = IINFO + N1 + RETURN + END IF +* +* Compute the Cholesky factorization A = L*L**T +* + ELSE +* +* Update and scale A21 +* + CALL STRSM( 'R', 'L', 'T', 'N', N2, N1, ONE, + $ A( 1, 1 ), LDA, A( N1+1, 1 ), LDA ) +* +* Update and factor A22 +* + CALL SSYRK( UPLO, 'N', N2, N1, -ONE, A( N1+1, 1 ), LDA, + $ ONE, A( N1+1, N1+1 ), LDA ) + CALL SPOTRF2( UPLO, N2, A( N1+1, N1+1 ), LDA, IINFO ) + IF ( IINFO.NE.0 ) THEN + INFO = IINFO + N1 + RETURN + END IF + END IF + END IF + RETURN +* +* End of SPOTRF2 +* + END diff --git a/SRC/zpotrf.f b/SRC/zpotrf.f index 0d8055f..216d7d3 100644 --- a/SRC/zpotrf.f +++ b/SRC/zpotrf.f @@ -137,7 +137,7 @@ EXTERNAL LSAME, ILAENV * .. * .. External Subroutines .. - EXTERNAL XERBLA, ZGEMM, ZHERK, ZPOTF2, ZTRSM + EXTERNAL XERBLA, ZGEMM, ZHERK, ZPOTRF2, ZTRSM * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN @@ -172,7 +172,7 @@ * * Use unblocked code. * - CALL ZPOTF2( UPLO, N, A, LDA, INFO ) + CALL ZPOTRF2( UPLO, N, A, LDA, INFO ) ELSE * * Use blocked code. @@ -189,7 +189,7 @@ JB = MIN( NB, N-J+1 ) CALL ZHERK( 'Upper', 'Conjugate transpose', JB, J-1, $ -ONE, A( 1, J ), LDA, ONE, A( J, J ), LDA ) - CALL ZPOTF2( 'Upper', JB, A( J, J ), LDA, INFO ) + CALL ZPOTRF2( 'Upper', JB, A( J, J ), LDA, INFO ) IF( INFO.NE.0 ) $ GO TO 30 IF( J+JB.LE.N ) THEN @@ -218,7 +218,7 @@ JB = MIN( NB, N-J+1 ) CALL ZHERK( 'Lower', 'No transpose', JB, J-1, -ONE, $ A( J, 1 ), LDA, ONE, A( J, J ), LDA ) - CALL ZPOTF2( 'Lower', JB, A( J, J ), LDA, INFO ) + CALL ZPOTRF2( 'Lower', JB, A( J, J ), LDA, INFO ) IF( INFO.NE.0 ) $ GO TO 30 IF( J+JB.LE.N ) THEN diff --git a/SRC/zpotrf2.f b/SRC/zpotrf2.f new file mode 100644 index 0000000..c2a0829 --- /dev/null +++ b/SRC/zpotrf2.f @@ -0,0 +1,241 @@ +*> \brief \b ZPOTRF2 +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* RECURSIVE SUBROUTINE ZPOTRF2( UPLO, N, A, LDA, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER UPLO +* INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. +* COMPLEX*16 A( LDA, * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> ZPOTRF2 computes the Cholesky factorization of a real symmetric +*> positive definite matrix A using the recursive algorithm. +*> +*> The factorization has the form +*> A = U**H * U, if UPLO = 'U', or +*> A = L * L**H, if UPLO = 'L', +*> where U is an upper triangular matrix and L is lower triangular. +*> +*> This is the recursive version of the algorithm. It divides +*> the matrix into four submatrices: +*> +*> [ A11 | A12 ] where A11 is n1 by n1 and A22 is n2 by n2 +*> A = [ -----|----- ] with n1 = n/2 +*> [ A21 | A22 ] n2 = n-n1 +*> +*> The subroutine calls itself to factor A11. Update and scale A21 +*> or A12, update A22 then call itself to factor A22. +*> +*> \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, if INFO = 0, the factor U or L from the Cholesky +*> factorization A = U**H*U or A = L*L**H. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N). +*> \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, the leading minor of order i is not +*> positive definite, and the factorization could not be +*> completed. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2015 +* +*> \ingroup complex16POcomputational +* +* ===================================================================== + RECURSIVE SUBROUTINE ZPOTRF2( UPLO, N, A, LDA, INFO ) +* +* -- LAPACK computational routine (version 3.6.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2015 +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. + COMPLEX*16 A( LDA, * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ONE, ZERO + PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) + COMPLEX*16 CONE + PARAMETER ( CONE = (1.0D+0, 0.0D+0) ) +* .. +* .. Local Scalars .. + LOGICAL UPPER + INTEGER N1, N2, IINFO + DOUBLE PRECISION AJJ +* .. +* .. External Functions .. + LOGICAL LSAME, DISNAN + EXTERNAL LSAME, DISNAN +* .. +* .. External Subroutines .. + EXTERNAL ZHERK, ZTRSM, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, DBLE, SQRT +* .. +* .. Executable Statements .. +* +* Test the input parameters +* + INFO = 0 + UPPER = LSAME( UPLO, 'U' ) + 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 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'ZPOTRF2', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 ) + $ RETURN +* +* N=1 case +* + IF( N.EQ.1 ) THEN +* +* Test for non-positive-definiteness +* + AJJ = DBLE( A( 1, 1 ) ) + IF( AJJ.LE.ZERO.OR.DISNAN( AJJ ) ) THEN + INFO = 1 + RETURN + END IF +* +* Factor +* + A( 1, 1 ) = SQRT( AJJ ) +* +* Use recursive code +* + ELSE + N1 = N/2 + N2 = N-N1 +* +* Factor A11 +* + CALL ZPOTRF2( UPLO, N1, A( 1, 1 ), LDA, IINFO ) + IF ( IINFO.NE.0 ) THEN + INFO = IINFO + RETURN + END IF +* +* Compute the Cholesky factorization A = U**H*U +* + IF( UPPER ) THEN +* +* Update and scale A12 +* + CALL ZTRSM( 'L', 'U', 'C', 'N', N1, N2, CONE, + $ A( 1, 1 ), LDA, A( 1, N1+1 ), LDA ) +* +* Update and factor A22 +* + CALL ZHERK( UPLO, 'C', N2, N1, -ONE, A( 1, N1+1 ), LDA, + $ ONE, A( N1+1, N1+1 ), LDA ) + CALL ZPOTRF2( UPLO, N2, A( N1+1, N1+1 ), LDA, IINFO ) + IF ( IINFO.NE.0 ) THEN + INFO = IINFO + N1 + RETURN + END IF +* +* Compute the Cholesky factorization A = L*L**H +* + ELSE +* +* Update and scale A21 +* + CALL ZTRSM( 'R', 'L', 'C', 'N', N2, N1, CONE, + $ A( 1, 1 ), LDA, A( N1+1, 1 ), LDA ) +* +* Update and factor A22 +* + CALL ZHERK( UPLO, 'N', N2, N1, -ONE, A( N1+1, 1 ), LDA, + $ ONE, A( N1+1, N1+1 ), LDA ) + CALL ZPOTRF2( UPLO, N2, A( N1+1, N1+1 ), LDA, IINFO ) + IF ( IINFO.NE.0 ) THEN + INFO = IINFO + N1 + RETURN + END IF + END IF + END IF + RETURN +* +* End of ZPOTRF2 +* + END -- 2.7.4