Mercurial > octave-nkf
changeset 5340:15843d76156d
[project @ 2005-05-06 16:26:58 by jwe]
| author | jwe |
|---|---|
| date | Fri, 06 May 2005 16:26:59 +0000 |
| parents | 4266ef7972b2 |
| children | 4bea82210dcd |
| files | libcruft/ChangeLog libcruft/lapack/dlauu2.f libcruft/lapack/dlauum.f libcruft/lapack/dpotri.f libcruft/lapack/zlauu2.f libcruft/lapack/zlauum.f libcruft/lapack/zpotri.f liboctave/Array.cc liboctave/ChangeLog liboctave/CmplxCHOL.cc liboctave/CmplxCHOL.h liboctave/dbleCHOL.cc liboctave/dbleCHOL.h src/ChangeLog src/DLD-FUNCTIONS/chol.cc |
| diffstat | 15 files changed, 1055 insertions(+), 9 deletions(-) [+] |
line wrap: on
line diff
--- a/libcruft/ChangeLog Thu May 05 19:06:26 2005 +0000 +++ b/libcruft/ChangeLog Fri May 06 16:26:59 2005 +0000 @@ -1,3 +1,8 @@ +2005-05-06 John W. Eaton <jwe@octave.org> + + * lapack/dpotri.f, lapack/dlauum.f, lapack/dlauu2.f, + lapack/zpotri.f, lapack/zlauum.f, lapack/zlauu2.f: New files. + 2005-04-08 John W. Eaton <jwe@octave.org> * Makefile.in, Makerules.in (clean, distclean, maintainer-clean):
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/lapack/dlauu2.f Fri May 06 16:26:59 2005 +0000 @@ -0,0 +1,136 @@ + SUBROUTINE DLAUU2( UPLO, N, A, LDA, INFO ) +* +* -- LAPACK auxiliary routine (version 3.0) -- +* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., +* Courant Institute, Argonne National Lab, and Rice University +* February 29, 1992 +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. + DOUBLE PRECISION A( LDA, * ) +* .. +* +* Purpose +* ======= +* +* DLAUU2 computes the product U * U' or L' * L, where the triangular +* factor U or L is stored in the upper or lower triangular part of +* the array A. +* +* If UPLO = 'U' or 'u' then the upper triangle of the result is stored, +* overwriting the factor U in A. +* If UPLO = 'L' or 'l' then the lower triangle of the result is stored, +* overwriting the factor L in A. +* +* This is the unblocked form of the algorithm, calling Level 2 BLAS. +* +* Arguments +* ========= +* +* UPLO (input) CHARACTER*1 +* Specifies whether the triangular factor stored in the array A +* is upper or lower triangular: +* = 'U': Upper triangular +* = 'L': Lower triangular +* +* N (input) INTEGER +* The order of the triangular factor U or L. N >= 0. +* +* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) +* On entry, the triangular factor U or L. +* On exit, if UPLO = 'U', the upper triangle of A is +* overwritten with the upper triangle of the product U * U'; +* if UPLO = 'L', the lower triangle of A is overwritten with +* the lower triangle of the product L' * L. +* +* LDA (input) INTEGER +* The leading dimension of the array A. LDA >= max(1,N). +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -k, the k-th argument had an illegal value +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ONE + PARAMETER ( ONE = 1.0D+0 ) +* .. +* .. Local Scalars .. + LOGICAL UPPER + INTEGER I + DOUBLE PRECISION AII +* .. +* .. External Functions .. + LOGICAL LSAME + DOUBLE PRECISION DDOT + EXTERNAL LSAME, DDOT +* .. +* .. External Subroutines .. + EXTERNAL DGEMV, DSCAL, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX +* .. +* .. 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( 'DLAUU2', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 ) + $ RETURN +* + IF( UPPER ) THEN +* +* Compute the product U * U'. +* + DO 10 I = 1, N + AII = A( I, I ) + IF( I.LT.N ) THEN + A( I, I ) = DDOT( N-I+1, A( I, I ), LDA, A( I, I ), LDA ) + CALL DGEMV( 'No transpose', I-1, N-I, ONE, A( 1, I+1 ), + $ LDA, A( I, I+1 ), LDA, AII, A( 1, I ), 1 ) + ELSE + CALL DSCAL( I, AII, A( 1, I ), 1 ) + END IF + 10 CONTINUE +* + ELSE +* +* Compute the product L' * L. +* + DO 20 I = 1, N + AII = A( I, I ) + IF( I.LT.N ) THEN + A( I, I ) = DDOT( N-I+1, A( I, I ), 1, A( I, I ), 1 ) + CALL DGEMV( 'Transpose', N-I, I-1, ONE, A( I+1, 1 ), LDA, + $ A( I+1, I ), 1, AII, A( I, 1 ), LDA ) + ELSE + CALL DSCAL( I, AII, A( I, 1 ), LDA ) + END IF + 20 CONTINUE + END IF +* + RETURN +* +* End of DLAUU2 +* + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/lapack/dlauum.f Fri May 06 16:26:59 2005 +0000 @@ -0,0 +1,156 @@ + SUBROUTINE DLAUUM( UPLO, N, A, LDA, INFO ) +* +* -- LAPACK auxiliary routine (version 3.0) -- +* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., +* Courant Institute, Argonne National Lab, and Rice University +* February 29, 1992 +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. + DOUBLE PRECISION A( LDA, * ) +* .. +* +* Purpose +* ======= +* +* DLAUUM computes the product U * U' or L' * L, where the triangular +* factor U or L is stored in the upper or lower triangular part of +* the array A. +* +* If UPLO = 'U' or 'u' then the upper triangle of the result is stored, +* overwriting the factor U in A. +* If UPLO = 'L' or 'l' then the lower triangle of the result is stored, +* overwriting the factor L in A. +* +* This is the blocked form of the algorithm, calling Level 3 BLAS. +* +* Arguments +* ========= +* +* UPLO (input) CHARACTER*1 +* Specifies whether the triangular factor stored in the array A +* is upper or lower triangular: +* = 'U': Upper triangular +* = 'L': Lower triangular +* +* N (input) INTEGER +* The order of the triangular factor U or L. N >= 0. +* +* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) +* On entry, the triangular factor U or L. +* On exit, if UPLO = 'U', the upper triangle of A is +* overwritten with the upper triangle of the product U * U'; +* if UPLO = 'L', the lower triangle of A is overwritten with +* the lower triangle of the product L' * L. +* +* LDA (input) INTEGER +* The leading dimension of the array A. LDA >= max(1,N). +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -k, the k-th argument had an illegal value +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ONE + PARAMETER ( ONE = 1.0D+0 ) +* .. +* .. Local Scalars .. + LOGICAL UPPER + INTEGER I, IB, NB +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER ILAENV + EXTERNAL LSAME, ILAENV +* .. +* .. External Subroutines .. + EXTERNAL DGEMM, DLAUU2, DSYRK, DTRMM, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. 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( 'DLAUUM', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 ) + $ RETURN +* +* Determine the block size for this environment. +* + NB = ILAENV( 1, 'DLAUUM', UPLO, N, -1, -1, -1 ) +* + IF( NB.LE.1 .OR. NB.GE.N ) THEN +* +* Use unblocked code +* + CALL DLAUU2( UPLO, N, A, LDA, INFO ) + ELSE +* +* Use blocked code +* + IF( UPPER ) THEN +* +* Compute the product U * U'. +* + DO 10 I = 1, N, NB + IB = MIN( NB, N-I+1 ) + CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Non-unit', + $ I-1, IB, ONE, A( I, I ), LDA, A( 1, I ), + $ LDA ) + CALL DLAUU2( 'Upper', IB, A( I, I ), LDA, INFO ) + IF( I+IB.LE.N ) THEN + CALL DGEMM( 'No transpose', 'Transpose', I-1, IB, + $ N-I-IB+1, ONE, A( 1, I+IB ), LDA, + $ A( I, I+IB ), LDA, ONE, A( 1, I ), LDA ) + CALL DSYRK( 'Upper', 'No transpose', IB, N-I-IB+1, + $ ONE, A( I, I+IB ), LDA, ONE, A( I, I ), + $ LDA ) + END IF + 10 CONTINUE + ELSE +* +* Compute the product L' * L. +* + DO 20 I = 1, N, NB + IB = MIN( NB, N-I+1 ) + CALL DTRMM( 'Left', 'Lower', 'Transpose', 'Non-unit', IB, + $ I-1, ONE, A( I, I ), LDA, A( I, 1 ), LDA ) + CALL DLAUU2( 'Lower', IB, A( I, I ), LDA, INFO ) + IF( I+IB.LE.N ) THEN + CALL DGEMM( 'Transpose', 'No transpose', IB, I-1, + $ N-I-IB+1, ONE, A( I+IB, I ), LDA, + $ A( I+IB, 1 ), LDA, ONE, A( I, 1 ), LDA ) + CALL DSYRK( 'Lower', 'Transpose', IB, N-I-IB+1, ONE, + $ A( I+IB, I ), LDA, ONE, A( I, I ), LDA ) + END IF + 20 CONTINUE + END IF + END IF +* + RETURN +* +* End of DLAUUM +* + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/lapack/dpotri.f Fri May 06 16:26:59 2005 +0000 @@ -0,0 +1,97 @@ + SUBROUTINE DPOTRI( UPLO, N, A, LDA, INFO ) +* +* -- LAPACK routine (version 3.0) -- +* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., +* Courant Institute, Argonne National Lab, and Rice University +* March 31, 1993 +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. + DOUBLE PRECISION A( LDA, * ) +* .. +* +* Purpose +* ======= +* +* DPOTRI computes the inverse of a real symmetric positive definite +* matrix A using the Cholesky factorization A = U**T*U or A = L*L**T +* computed by DPOTRF. +* +* Arguments +* ========= +* +* UPLO (input) CHARACTER*1 +* = 'U': Upper triangle of A is stored; +* = 'L': Lower triangle of A is stored. +* +* N (input) INTEGER +* The order of the matrix A. N >= 0. +* +* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) +* On entry, the triangular factor U or L from the Cholesky +* factorization A = U**T*U or A = L*L**T, as computed by +* DPOTRF. +* On exit, the upper or lower triangle of the (symmetric) +* inverse of A, overwriting the input factor U or L. +* +* LDA (input) INTEGER +* The leading dimension of the array A. LDA >= max(1,N). +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument had an illegal value +* > 0: if INFO = i, the (i,i) element of the factor U or L is +* zero, and the inverse could not be computed. +* +* ===================================================================== +* +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. External Subroutines .. + EXTERNAL DLAUUM, DTRTRI, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX +* .. +* .. Executable Statements .. +* +* Test the input parameters. +* + INFO = 0 + IF( .NOT.LSAME( UPLO, 'U' ) .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( 'DPOTRI', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 ) + $ RETURN +* +* Invert the triangular Cholesky factor U or L. +* + CALL DTRTRI( UPLO, 'Non-unit', N, A, LDA, INFO ) + IF( INFO.GT.0 ) + $ RETURN +* +* Form inv(U)*inv(U)' or inv(L)'*inv(L). +* + CALL DLAUUM( UPLO, N, A, LDA, INFO ) +* + RETURN +* +* End of DPOTRI +* + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/lapack/zlauu2.f Fri May 06 16:26:59 2005 +0000 @@ -0,0 +1,144 @@ + SUBROUTINE ZLAUU2( UPLO, N, A, LDA, INFO ) +* +* -- LAPACK auxiliary routine (version 3.0) -- +* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., +* Courant Institute, Argonne National Lab, and Rice University +* September 30, 1994 +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. + COMPLEX*16 A( LDA, * ) +* .. +* +* Purpose +* ======= +* +* ZLAUU2 computes the product U * U' or L' * L, where the triangular +* factor U or L is stored in the upper or lower triangular part of +* the array A. +* +* If UPLO = 'U' or 'u' then the upper triangle of the result is stored, +* overwriting the factor U in A. +* If UPLO = 'L' or 'l' then the lower triangle of the result is stored, +* overwriting the factor L in A. +* +* This is the unblocked form of the algorithm, calling Level 2 BLAS. +* +* Arguments +* ========= +* +* UPLO (input) CHARACTER*1 +* Specifies whether the triangular factor stored in the array A +* is upper or lower triangular: +* = 'U': Upper triangular +* = 'L': Lower triangular +* +* N (input) INTEGER +* The order of the triangular factor U or L. N >= 0. +* +* A (input/output) COMPLEX*16 array, dimension (LDA,N) +* On entry, the triangular factor U or L. +* On exit, if UPLO = 'U', the upper triangle of A is +* overwritten with the upper triangle of the product U * U'; +* if UPLO = 'L', the lower triangle of A is overwritten with +* the lower triangle of the product L' * L. +* +* LDA (input) INTEGER +* The leading dimension of the array A. LDA >= max(1,N). +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -k, the k-th argument had an illegal value +* +* ===================================================================== +* +* .. Parameters .. + COMPLEX*16 ONE + PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) ) +* .. +* .. Local Scalars .. + LOGICAL UPPER + INTEGER I + DOUBLE PRECISION AII +* .. +* .. External Functions .. + LOGICAL LSAME + COMPLEX*16 ZDOTC + EXTERNAL LSAME, ZDOTC +* .. +* .. External Subroutines .. + EXTERNAL XERBLA, ZDSCAL, ZGEMV, ZLACGV +* .. +* .. Intrinsic Functions .. + INTRINSIC DBLE, DCMPLX, MAX +* .. +* .. 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( 'ZLAUU2', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 ) + $ RETURN +* + IF( UPPER ) THEN +* +* Compute the product U * U'. +* + DO 10 I = 1, N + AII = A( I, I ) + IF( I.LT.N ) THEN + A( I, I ) = AII*AII + DBLE( ZDOTC( N-I, A( I, I+1 ), LDA, + $ A( I, I+1 ), LDA ) ) + CALL ZLACGV( N-I, A( I, I+1 ), LDA ) + CALL ZGEMV( 'No transpose', I-1, N-I, ONE, A( 1, I+1 ), + $ LDA, A( I, I+1 ), LDA, DCMPLX( AII ), + $ A( 1, I ), 1 ) + CALL ZLACGV( N-I, A( I, I+1 ), LDA ) + ELSE + CALL ZDSCAL( I, AII, A( 1, I ), 1 ) + END IF + 10 CONTINUE +* + ELSE +* +* Compute the product L' * L. +* + DO 20 I = 1, N + AII = A( I, I ) + IF( I.LT.N ) THEN + A( I, I ) = AII*AII + DBLE( ZDOTC( N-I, A( I+1, I ), 1, + $ A( I+1, I ), 1 ) ) + CALL ZLACGV( I-1, A( I, 1 ), LDA ) + CALL ZGEMV( 'Conjugate transpose', N-I, I-1, ONE, + $ A( I+1, 1 ), LDA, A( I+1, I ), 1, + $ DCMPLX( AII ), A( I, 1 ), LDA ) + CALL ZLACGV( I-1, A( I, 1 ), LDA ) + ELSE + CALL ZDSCAL( I, AII, A( I, 1 ), LDA ) + END IF + 20 CONTINUE + END IF +* + RETURN +* +* End of ZLAUU2 +* + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/lapack/zlauum.f Fri May 06 16:26:59 2005 +0000 @@ -0,0 +1,161 @@ + SUBROUTINE ZLAUUM( UPLO, N, A, LDA, INFO ) +* +* -- LAPACK auxiliary routine (version 3.0) -- +* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., +* Courant Institute, Argonne National Lab, and Rice University +* September 30, 1994 +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. + COMPLEX*16 A( LDA, * ) +* .. +* +* Purpose +* ======= +* +* ZLAUUM computes the product U * U' or L' * L, where the triangular +* factor U or L is stored in the upper or lower triangular part of +* the array A. +* +* If UPLO = 'U' or 'u' then the upper triangle of the result is stored, +* overwriting the factor U in A. +* If UPLO = 'L' or 'l' then the lower triangle of the result is stored, +* overwriting the factor L in A. +* +* This is the blocked form of the algorithm, calling Level 3 BLAS. +* +* Arguments +* ========= +* +* UPLO (input) CHARACTER*1 +* Specifies whether the triangular factor stored in the array A +* is upper or lower triangular: +* = 'U': Upper triangular +* = 'L': Lower triangular +* +* N (input) INTEGER +* The order of the triangular factor U or L. N >= 0. +* +* A (input/output) COMPLEX*16 array, dimension (LDA,N) +* On entry, the triangular factor U or L. +* On exit, if UPLO = 'U', the upper triangle of A is +* overwritten with the upper triangle of the product U * U'; +* if UPLO = 'L', the lower triangle of A is overwritten with +* the lower triangle of the product L' * L. +* +* LDA (input) INTEGER +* The leading dimension of the array A. LDA >= max(1,N). +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -k, the k-th argument had an illegal value +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ONE + PARAMETER ( ONE = 1.0D+0 ) + COMPLEX*16 CONE + PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) ) +* .. +* .. Local Scalars .. + LOGICAL UPPER + INTEGER I, IB, NB +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER ILAENV + EXTERNAL LSAME, ILAENV +* .. +* .. External Subroutines .. + EXTERNAL XERBLA, ZGEMM, ZHERK, ZLAUU2, ZTRMM +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. 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( 'ZLAUUM', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 ) + $ RETURN +* +* Determine the block size for this environment. +* + NB = ILAENV( 1, 'ZLAUUM', UPLO, N, -1, -1, -1 ) +* + IF( NB.LE.1 .OR. NB.GE.N ) THEN +* +* Use unblocked code +* + CALL ZLAUU2( UPLO, N, A, LDA, INFO ) + ELSE +* +* Use blocked code +* + IF( UPPER ) THEN +* +* Compute the product U * U'. +* + DO 10 I = 1, N, NB + IB = MIN( NB, N-I+1 ) + CALL ZTRMM( 'Right', 'Upper', 'Conjugate transpose', + $ 'Non-unit', I-1, IB, CONE, A( I, I ), LDA, + $ A( 1, I ), LDA ) + CALL ZLAUU2( 'Upper', IB, A( I, I ), LDA, INFO ) + IF( I+IB.LE.N ) THEN + CALL ZGEMM( 'No transpose', 'Conjugate transpose', + $ I-1, IB, N-I-IB+1, CONE, A( 1, I+IB ), + $ LDA, A( I, I+IB ), LDA, CONE, A( 1, I ), + $ LDA ) + CALL ZHERK( 'Upper', 'No transpose', IB, N-I-IB+1, + $ ONE, A( I, I+IB ), LDA, ONE, A( I, I ), + $ LDA ) + END IF + 10 CONTINUE + ELSE +* +* Compute the product L' * L. +* + DO 20 I = 1, N, NB + IB = MIN( NB, N-I+1 ) + CALL ZTRMM( 'Left', 'Lower', 'Conjugate transpose', + $ 'Non-unit', IB, I-1, CONE, A( I, I ), LDA, + $ A( I, 1 ), LDA ) + CALL ZLAUU2( 'Lower', IB, A( I, I ), LDA, INFO ) + IF( I+IB.LE.N ) THEN + CALL ZGEMM( 'Conjugate transpose', 'No transpose', IB, + $ I-1, N-I-IB+1, CONE, A( I+IB, I ), LDA, + $ A( I+IB, 1 ), LDA, CONE, A( I, 1 ), LDA ) + CALL ZHERK( 'Lower', 'Conjugate transpose', IB, + $ N-I-IB+1, ONE, A( I+IB, I ), LDA, ONE, + $ A( I, I ), LDA ) + END IF + 20 CONTINUE + END IF + END IF +* + RETURN +* +* End of ZLAUUM +* + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/lapack/zpotri.f Fri May 06 16:26:59 2005 +0000 @@ -0,0 +1,97 @@ + SUBROUTINE ZPOTRI( UPLO, N, A, LDA, INFO ) +* +* -- LAPACK routine (version 3.0) -- +* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., +* Courant Institute, Argonne National Lab, and Rice University +* March 31, 1993 +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. + COMPLEX*16 A( LDA, * ) +* .. +* +* Purpose +* ======= +* +* ZPOTRI computes the inverse of a complex Hermitian positive definite +* matrix A using the Cholesky factorization A = U**H*U or A = L*L**H +* computed by ZPOTRF. +* +* Arguments +* ========= +* +* UPLO (input) CHARACTER*1 +* = 'U': Upper triangle of A is stored; +* = 'L': Lower triangle of A is stored. +* +* N (input) INTEGER +* The order of the matrix A. N >= 0. +* +* A (input/output) COMPLEX*16 array, dimension (LDA,N) +* On entry, the triangular factor U or L from the Cholesky +* factorization A = U**H*U or A = L*L**H, as computed by +* ZPOTRF. +* On exit, the upper or lower triangle of the (Hermitian) +* inverse of A, overwriting the input factor U or L. +* +* LDA (input) INTEGER +* The leading dimension of the array A. LDA >= max(1,N). +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument had an illegal value +* > 0: if INFO = i, the (i,i) element of the factor U or L is +* zero, and the inverse could not be computed. +* +* ===================================================================== +* +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. External Subroutines .. + EXTERNAL XERBLA, ZLAUUM, ZTRTRI +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX +* .. +* .. Executable Statements .. +* +* Test the input parameters. +* + INFO = 0 + IF( .NOT.LSAME( UPLO, 'U' ) .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( 'ZPOTRI', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 ) + $ RETURN +* +* Invert the triangular Cholesky factor U or L. +* + CALL ZTRTRI( UPLO, 'Non-unit', N, A, LDA, INFO ) + IF( INFO.GT.0 ) + $ RETURN +* +* Form inv(U)*inv(U)' or inv(L)'*inv(L). +* + CALL ZLAUUM( UPLO, N, A, LDA, INFO ) +* + RETURN +* +* End of ZPOTRI +* + END
--- a/liboctave/Array.cc Thu May 05 19:06:26 2005 +0000 +++ b/liboctave/Array.cc Fri May 06 16:26:59 2005 +0000 @@ -496,7 +496,7 @@ increment_index (old_idx, dv); } - chop_trailing_singletons (); + retval.chop_trailing_singletons (); return retval; }
--- a/liboctave/ChangeLog Thu May 05 19:06:26 2005 +0000 +++ b/liboctave/ChangeLog Fri May 06 16:26:59 2005 +0000 @@ -1,3 +1,17 @@ +2005-05-06 John W. Eaton <jwe@octave.org> + + * dbleCHOL.cc (CHOL::init): Use xelem instead of elem for indexing + chol_mat. + (chol2mat_internal, chol2mat, CHOL::inverse): New functions. + * dbleCHOL.h (chol2mat_internal, chol2mat, CHOL::inverse): + Provide decls. + + * CmplxChol.cc (ComplexCHOL::init): Use xelem instead of elem for + indexing chol_mat. + (chol2mat_internal, chol2mat, ComplexCHOL::inverse): New functions. + * CmplxCHOL.h (chol2mat_internal, chol2mat, CmplxCHOL::inverse): + Provide decls. + 2005-05-05 John W. Eaton <jwe@octave.org> * Array.cc (Array<T>::permute): Call chop_trailing_singletons on
--- a/liboctave/CmplxCHOL.cc Thu May 05 19:06:26 2005 +0000 +++ b/liboctave/CmplxCHOL.cc Fri May 06 16:26:59 2005 +0000 @@ -32,8 +32,12 @@ extern "C" { F77_RET_T - F77_FUNC (zpotrf, ZPOTRF) (F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, Complex*, const octave_idx_type&, octave_idx_type& + F77_FUNC (zpotrf, ZPOTRF) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, + Complex*, const octave_idx_type&, octave_idx_type& + F77_CHAR_ARG_LEN_DECL); + F77_RET_T + F77_FUNC (zpotri, ZPOTRI) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, + Complex*, const octave_idx_type&, octave_idx_type& F77_CHAR_ARG_LEN_DECL); } @@ -69,12 +73,65 @@ if (n > 1) for (octave_idx_type j = 0; j < a_nc; j++) for (octave_idx_type i = j+1; i < a_nr; i++) - chol_mat.elem (i, j) = 0.0; + chol_mat.xelem (i, j) = 0.0; } return info; } +static ComplexMatrix +chol2inv_internal (const ComplexMatrix& r) +{ + ComplexMatrix retval; + + octave_idx_type r_nr = r.rows (); + octave_idx_type r_nc = r.cols (); + + if (r_nr == r_nc) + { + octave_idx_type n = r_nc; + octave_idx_type info; + + ComplexMatrix tmp = r; + + F77_XFCN (zpotri, ZPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n, + tmp.fortran_vec (), n, info + F77_CHAR_ARG_LEN (1))); + + if (f77_exception_encountered) + (*current_liboctave_error_handler) ("unrecoverable error in zpotri"); + else + { + // If someone thinks of a more graceful way of doing this (or + // faster for that matter :-)), please let me know! + + if (n > 1) + for (octave_idx_type j = 0; j < r_nc; j++) + for (octave_idx_type i = j+1; i < r_nr; i++) + tmp.xelem (i, j) = std::conj (tmp.xelem (j, i)); + + retval = tmp; + } + } + else + (*current_liboctave_error_handler) ("chol2inv requires square matrix"); + + return retval; +} + +// Compute the inverse of a matrix using the Cholesky factorization. +ComplexMatrix +ComplexCHOL::inverse (void) const +{ + return chol2inv_internal (chol_mat); +} + +ComplexMatrix +chol2inv (const ComplexMatrix& r) +{ + return chol2inv_internal (r); +} + /* ;;; Local Variables: *** ;;; mode: C++ ***
--- a/liboctave/CmplxCHOL.h Thu May 05 19:06:26 2005 +0000 +++ b/liboctave/CmplxCHOL.h Fri May 06 16:26:59 2005 +0000 @@ -53,10 +53,9 @@ return *this; } - ComplexMatrix chol_matrix (void) const - { - return chol_mat; - } + ComplexMatrix chol_matrix (void) const { return chol_mat; } + + ComplexMatrix inverse (void) const; friend std::ostream& operator << (std::ostream& os, const ComplexCHOL& a); @@ -67,6 +66,8 @@ octave_idx_type init (const ComplexMatrix& a); }; +ComplexMatrix chol2inv (const ComplexMatrix& r); + #endif /*
--- a/liboctave/dbleCHOL.cc Thu May 05 19:06:26 2005 +0000 +++ b/liboctave/dbleCHOL.cc Fri May 06 16:26:59 2005 +0000 @@ -35,6 +35,11 @@ F77_FUNC (dpotrf, DPOTRF) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, double*, const octave_idx_type&, octave_idx_type& F77_CHAR_ARG_LEN_DECL); + + F77_RET_T + F77_FUNC (dpotri, DPOTRI) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, + double*, const octave_idx_type&, octave_idx_type& + F77_CHAR_ARG_LEN_DECL); } octave_idx_type @@ -69,12 +74,65 @@ if (n > 1) for (octave_idx_type j = 0; j < a_nc; j++) for (octave_idx_type i = j+1; i < a_nr; i++) - chol_mat.elem (i, j) = 0.0; + chol_mat.xelem (i, j) = 0.0; } return info; } +static Matrix +chol2inv_internal (const Matrix& r) +{ + Matrix retval; + + octave_idx_type r_nr = r.rows (); + octave_idx_type r_nc = r.cols (); + + if (r_nr == r_nc) + { + octave_idx_type n = r_nc; + octave_idx_type info; + + Matrix tmp = r; + + F77_XFCN (dpotri, DPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n, + tmp.fortran_vec (), n, info + F77_CHAR_ARG_LEN (1))); + + if (f77_exception_encountered) + (*current_liboctave_error_handler) ("unrecoverable error in dpotri"); + else + { + // If someone thinks of a more graceful way of doing this (or + // faster for that matter :-)), please let me know! + + if (n > 1) + for (octave_idx_type j = 0; j < r_nc; j++) + for (octave_idx_type i = j+1; i < r_nr; i++) + tmp.xelem (i, j) = tmp.xelem (j, i); + + retval = tmp; + } + } + else + (*current_liboctave_error_handler) ("chol2inv requires square matrix"); + + return retval; +} + +// Compute the inverse of a matrix using the Cholesky factorization. +Matrix +CHOL::inverse (void) const +{ + return chol2inv_internal (chol_mat); +} + +Matrix +chol2inv (const Matrix& r) +{ + return chol2inv_internal (r); +} + /* ;;; Local Variables: *** ;;; mode: C++ ***
--- a/liboctave/dbleCHOL.h Thu May 05 19:06:26 2005 +0000 +++ b/liboctave/dbleCHOL.h Fri May 06 16:26:59 2005 +0000 @@ -51,6 +51,9 @@ Matrix chol_matrix (void) const { return chol_mat; } + // Compute the inverse of a matrix using the Cholesky factorization. + Matrix inverse (void) const; + friend std::ostream& operator << (std::ostream& os, const CHOL& a); private: @@ -60,6 +63,8 @@ octave_idx_type init (const Matrix& a); }; +Matrix chol2inv (const Matrix& r); + #endif /*
--- a/src/ChangeLog Thu May 05 19:06:26 2005 +0000 +++ b/src/ChangeLog Fri May 06 16:26:59 2005 +0000 @@ -1,3 +1,8 @@ +2005-05-06 John W. Eaton <jwe@octave.org> + + * DLD-FUNCTIONS/chol.cc (Fcholinv): New function. + (Fchol2inv): New function. + 2005-05-05 Keith Goodman <kwgoodman@gmail.com> * ov-usr-fcn.cc (Fnargout, Fnargin): Update doc strings.
--- a/src/DLD-FUNCTIONS/chol.cc Thu May 05 19:06:26 2005 +0000 +++ b/src/DLD-FUNCTIONS/chol.cc Fri May 06 16:26:59 2005 +0000 @@ -51,6 +51,7 @@ r' * r = a.\n\ @end example\n\ @end ifinfo\n\ +@seealso{cholinv, chol2inv}\n\ @end deftypefn") { octave_value_list retval; @@ -117,6 +118,115 @@ return retval; } +DEFUN_DLD (cholinv, args, nargout, + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {} cholinv (@var{a})\n\ +Use the Cholesky factorization to compute the inverse of of the\n\ +symmetric positive definite matrix @var{a}.\n\ +@seealso{chol, chol2inv}\n\ +@end deftypefn") +{ + octave_value retval; + + int nargin = args.length (); + + if (nargin == 1) + { + octave_value arg = args(0); + + octave_idx_type nr = arg.rows (); + octave_idx_type nc = arg.columns (); + + if (nr == 0 || nc == 0) + retval = Matrix (); + else + { + if (arg.is_real_type ()) + { + Matrix m = arg.matrix_value (); + + if (! error_state) + { + octave_idx_type info; + CHOL chol (m, info); + if (info == 0) + retval = chol.inverse (); + else + error ("cholinv: matrix not positive definite"); + } + } + else if (arg.is_complex_type ()) + { + ComplexMatrix m = arg.complex_matrix_value (); + + if (! error_state) + { + octave_idx_type info; + ComplexCHOL chol (m, info); + if (info == 0) + retval = chol.inverse (); + else + error ("cholinv: matrix not positive definite"); + } + } + else + gripe_wrong_type_arg ("chol", arg); + } + } + else + print_usage ("chol"); + + return retval; +} + +DEFUN_DLD (chol2inv, args, nargout, + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {} chol2inv (@var{r})\n\ +Invert a symmetric, positive definite square matrix from its Cholesky\n\ +decomposition, @var{r}. Note that no check is performed to ensure\n\ +that @var{r} is actually a Cholesky factor.\n\ +@seealso{chol, cholinv}\n\ +@end deftypefn") +{ + octave_value retval; + + int nargin = args.length (); + + if (nargin == 1) + { + octave_value arg = args(0); + + octave_idx_type nr = arg.rows (); + octave_idx_type nc = arg.columns (); + + if (nr == 0 || nc == 0) + retval = Matrix (); + else + { + if (arg.is_real_type ()) + { + Matrix r = arg.matrix_value (); + + if (! error_state) + retval = chol2inv (r); + } + else if (arg.is_complex_type ()) + { + ComplexMatrix r = arg.complex_matrix_value (); + + if (! error_state) + retval = chol2inv (r); + } + else + gripe_wrong_type_arg ("chol2inv", arg); + } + } + else + print_usage ("chol2inv"); + + return retval; +} + /* ;;; Local Variables: *** ;;; mode: C++ ***