Mercurial > octave-nkf
changeset 4329:d53c33d93440
[project @ 2003-02-18 20:00:48 by jwe]
line wrap: on
line diff
--- a/ChangeLog Sun Feb 16 04:29:00 2003 +0000 +++ b/ChangeLog Tue Feb 18 20:08:20 2003 +0000 @@ -1,3 +1,7 @@ +2003-02-18 David Bateman <dbateman@free.fr> + + * configure.in: Eliminate linpack + 2003-02-15 John W. Eaton <jwe@bevo.che.wisc.edu> * configure.in: Check for mkstemp too.
--- a/configure.in Sun Feb 16 04:29:00 2003 +0000 +++ b/configure.in Tue Feb 18 20:08:20 2003 +0000 @@ -22,7 +22,7 @@ ### 02111-1307, USA. AC_INIT -AC_REVISION($Revision: 1.411 $) +AC_REVISION($Revision: 1.412 $) AC_PREREQ(2.52) AC_CONFIG_SRCDIR([src/octave.cc]) AC_CONFIG_HEADER(config.h) @@ -1374,17 +1374,15 @@ doc/Makefile doc/faq/Makefile doc/interpreter/Makefile \ doc/liboctave/Makefile doc/refcard/Makefile emacs/Makefile \ examples/Makefile liboctave/Makefile src/Makefile \ - libcruft/Makefile libcruft/Makerules \ - libcruft/amos/Makefile libcruft/blas/Makefile \ - libcruft/daspk/Makefile libcruft/dasrt/Makefile \ - libcruft/dassl/Makefile libcruft/fftpack/Makefile \ - libcruft/lapack/Makefile libcruft/linpack/Makefile \ + libcruft/Makefile libcruft/Makerules libcruft/amos/Makefile \ + libcruft/blas/Makefile libcruft/daspk/Makefile \ + libcruft/dasrt/Makefile libcruft/dassl/Makefile \ + libcruft/fftpack/Makefile libcruft/lapack/Makefile \ libcruft/minpack/Makefile libcruft/misc/Makefile \ libcruft/odepack/Makefile libcruft/odessa/Makefile \ - libcruft/ordered-qz/Makefile \ - libcruft/quadpack/Makefile libcruft/ranlib/Makefile \ - libcruft/slatec-fn/Makefile libcruft/slatec-err/Makefile \ - libcruft/villad/Makefile \ + libcruft/ordered-qz/Makefile libcruft/quadpack/Makefile \ + libcruft/ranlib/Makefile libcruft/slatec-fn/Makefile \ + libcruft/slatec-err/Makefile libcruft/villad/Makefile \ libcruft/blas-xtra/Makefile libcruft/lapack-xtra/Makefile]) AC_OUTPUT
--- a/libcruft/ChangeLog Sun Feb 16 04:29:00 2003 +0000 +++ b/libcruft/ChangeLog Tue Feb 18 20:08:20 2003 +0000 @@ -1,3 +1,37 @@ +2003-02-18 John W. Eaton <jwe@bevo.che.wisc.edu> + + * libcruft/blas/dtbsv.f: New file. + * libcruft/lapack/dlatrs.f, libcruft/lapack/dtrti2.f, + libcruft/lapack/dtrtri.f, libcruft/lapack/ztrti2.f, + libcruft/lapack/ztrtri.f: New files. + +2003-02-04 David Bateman <dbateman@free.fr> + + * libcruft/Makefile.in (CRUFT_DIRS): Remove linpack from list. + + * libcruft/linpackdgbfa.f, libcruft/linpackdgbsl.f + libcruft/linpackdgeco.f, libcruft/linpackdgedi.f, + libcruft/linpackdgefa.f, libcruft/linpackdgesl.f, + libcruft/linpackspofa.f, libcruft/linpackzgeco.f, + libcruft/linpackzgedi.f, libcruft/linpackzgefa.f, + libcruft/linpackzgesl.f: Delete. + + * libcruft/dassl/ddajac.f, libcruft/dassl/ddaslv.f: Use DGxTRF and + DGxTRS instead of DGxFA and DGxSL. + * libcruft/daspk/ddaspk.f, libcruft/daspk/dmatd.f, + libcruft/daspk/dslvd.f: Likewise. + * libcruft/odepack/lsode.f, libcruft/odepack/prepj.f, + libcruft/odepack/solsy.f: Likewise. + * libcruft/odessa/odessa.f, libcruft/odessa/odessa_prepj.f, + libcruft/odessa/odessa_solsy.f: Likewise. + * libcrudt/ranlib/setgmn.f: Use SPOTRF instead of SPOFA. + + * libcruft/lapack/dgbtf2.f, libcruft/lapack/dgbtrf.f, + libcruft/lapack/dgbtrs.f, libcruft/lapack/dgecon.f, + libcruft/lapack/dgetri.f, libcruft/lapack/spotf2.f, + libcruft/lapack/spotrf.f, libcruft/lapack/zgecon.f, + libcruft/lapack/zgetri.f: New files. + 2003-01-22 John W. Eaton <jwe@bevo.che.wisc.edu> * misc/quit.h (BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE_1,
--- a/libcruft/Makefile.in Sun Feb 16 04:29:00 2003 +0000 +++ b/libcruft/Makefile.in Tue Feb 18 20:08:20 2003 +0000 @@ -29,7 +29,7 @@ # dirname is substituted by configure and may be the empty string. CRUFT_DIRS = amos @BLAS_DIR@ blas-xtra daspk dasrt dassl \ - @FFT_DIR@ @LAPACK_DIR@ lapack-xtra linpack minpack \ + @FFT_DIR@ @LAPACK_DIR@ lapack-xtra minpack \ misc odepack odessa ordered-qz quadpack ranlib \ slatec-err slatec-fn villad
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/blas/dtbsv.f Tue Feb 18 20:08:20 2003 +0000 @@ -0,0 +1,346 @@ + SUBROUTINE DTBSV ( UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX ) +* .. Scalar Arguments .. + INTEGER INCX, K, LDA, N + CHARACTER*1 DIAG, TRANS, UPLO +* .. Array Arguments .. + DOUBLE PRECISION A( LDA, * ), X( * ) +* .. +* +* Purpose +* ======= +* +* DTBSV solves one of the systems of equations +* +* A*x = b, or A'*x = b, +* +* where b and x are n element vectors and A is an n by n unit, or +* non-unit, upper or lower triangular band matrix, with ( k + 1 ) +* diagonals. +* +* No test for singularity or near-singularity is included in this +* routine. Such tests must be performed before calling this routine. +* +* Parameters +* ========== +* +* UPLO - CHARACTER*1. +* On entry, UPLO specifies whether the matrix is an upper or +* lower triangular matrix as follows: +* +* UPLO = 'U' or 'u' A is an upper triangular matrix. +* +* UPLO = 'L' or 'l' A is a lower triangular matrix. +* +* Unchanged on exit. +* +* TRANS - CHARACTER*1. +* On entry, TRANS specifies the equations to be solved as +* follows: +* +* TRANS = 'N' or 'n' A*x = b. +* +* TRANS = 'T' or 't' A'*x = b. +* +* TRANS = 'C' or 'c' A'*x = b. +* +* Unchanged on exit. +* +* DIAG - CHARACTER*1. +* On entry, DIAG specifies whether or not A is unit +* triangular as follows: +* +* DIAG = 'U' or 'u' A is assumed to be unit triangular. +* +* DIAG = 'N' or 'n' A is not assumed to be unit +* triangular. +* +* Unchanged on exit. +* +* N - INTEGER. +* On entry, N specifies the order of the matrix A. +* N must be at least zero. +* Unchanged on exit. +* +* K - INTEGER. +* On entry with UPLO = 'U' or 'u', K specifies the number of +* super-diagonals of the matrix A. +* On entry with UPLO = 'L' or 'l', K specifies the number of +* sub-diagonals of the matrix A. +* K must satisfy 0 .le. K. +* Unchanged on exit. +* +* A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). +* Before entry with UPLO = 'U' or 'u', the leading ( k + 1 ) +* by n part of the array A must contain the upper triangular +* band part of the matrix of coefficients, supplied column by +* column, with the leading diagonal of the matrix in row +* ( k + 1 ) of the array, the first super-diagonal starting at +* position 2 in row k, and so on. The top left k by k triangle +* of the array A is not referenced. +* The following program segment will transfer an upper +* triangular band matrix from conventional full matrix storage +* to band storage: +* +* DO 20, J = 1, N +* M = K + 1 - J +* DO 10, I = MAX( 1, J - K ), J +* A( M + I, J ) = matrix( I, J ) +* 10 CONTINUE +* 20 CONTINUE +* +* Before entry with UPLO = 'L' or 'l', the leading ( k + 1 ) +* by n part of the array A must contain the lower triangular +* band part of the matrix of coefficients, supplied column by +* column, with the leading diagonal of the matrix in row 1 of +* the array, the first sub-diagonal starting at position 1 in +* row 2, and so on. The bottom right k by k triangle of the +* array A is not referenced. +* The following program segment will transfer a lower +* triangular band matrix from conventional full matrix storage +* to band storage: +* +* DO 20, J = 1, N +* M = 1 - J +* DO 10, I = J, MIN( N, J + K ) +* A( M + I, J ) = matrix( I, J ) +* 10 CONTINUE +* 20 CONTINUE +* +* Note that when DIAG = 'U' or 'u' the elements of the array A +* corresponding to the diagonal elements of the matrix are not +* referenced, but are assumed to be unity. +* Unchanged on exit. +* +* LDA - INTEGER. +* On entry, LDA specifies the first dimension of A as declared +* in the calling (sub) program. LDA must be at least +* ( k + 1 ). +* Unchanged on exit. +* +* X - DOUBLE PRECISION array of dimension at least +* ( 1 + ( n - 1 )*abs( INCX ) ). +* Before entry, the incremented array X must contain the n +* element right-hand side vector b. On exit, X is overwritten +* with the solution vector x. +* +* INCX - INTEGER. +* On entry, INCX specifies the increment for the elements of +* X. INCX must not be zero. +* Unchanged on exit. +* +* +* Level 2 Blas routine. +* +* -- Written on 22-October-1986. +* Jack Dongarra, Argonne National Lab. +* Jeremy Du Croz, Nag Central Office. +* Sven Hammarling, Nag Central Office. +* Richard Hanson, Sandia National Labs. +* +* +* .. Parameters .. + DOUBLE PRECISION ZERO + PARAMETER ( ZERO = 0.0D+0 ) +* .. Local Scalars .. + DOUBLE PRECISION TEMP + INTEGER I, INFO, IX, J, JX, KPLUS1, KX, L + LOGICAL NOUNIT +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. External Subroutines .. + EXTERNAL XERBLA +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. Executable Statements .. +* +* Test the input parameters. +* + INFO = 0 + IF ( .NOT.LSAME( UPLO , 'U' ).AND. + $ .NOT.LSAME( UPLO , 'L' ) )THEN + INFO = 1 + ELSE IF( .NOT.LSAME( TRANS, 'N' ).AND. + $ .NOT.LSAME( TRANS, 'T' ).AND. + $ .NOT.LSAME( TRANS, 'C' ) )THEN + INFO = 2 + ELSE IF( .NOT.LSAME( DIAG , 'U' ).AND. + $ .NOT.LSAME( DIAG , 'N' ) )THEN + INFO = 3 + ELSE IF( N.LT.0 )THEN + INFO = 4 + ELSE IF( K.LT.0 )THEN + INFO = 5 + ELSE IF( LDA.LT.( K + 1 ) )THEN + INFO = 7 + ELSE IF( INCX.EQ.0 )THEN + INFO = 9 + END IF + IF( INFO.NE.0 )THEN + CALL XERBLA( 'DTBSV ', INFO ) + RETURN + END IF +* +* Quick return if possible. +* + IF( N.EQ.0 ) + $ RETURN +* + NOUNIT = LSAME( DIAG, 'N' ) +* +* Set up the start point in X if the increment is not unity. This +* will be ( N - 1 )*INCX too small for descending loops. +* + IF( INCX.LE.0 )THEN + KX = 1 - ( N - 1 )*INCX + ELSE IF( INCX.NE.1 )THEN + KX = 1 + END IF +* +* Start the operations. In this version the elements of A are +* accessed by sequentially with one pass through A. +* + IF( LSAME( TRANS, 'N' ) )THEN +* +* Form x := inv( A )*x. +* + IF( LSAME( UPLO, 'U' ) )THEN + KPLUS1 = K + 1 + IF( INCX.EQ.1 )THEN + DO 20, J = N, 1, -1 + IF( X( J ).NE.ZERO )THEN + L = KPLUS1 - J + IF( NOUNIT ) + $ X( J ) = X( J )/A( KPLUS1, J ) + TEMP = X( J ) + DO 10, I = J - 1, MAX( 1, J - K ), -1 + X( I ) = X( I ) - TEMP*A( L + I, J ) + 10 CONTINUE + END IF + 20 CONTINUE + ELSE + KX = KX + ( N - 1 )*INCX + JX = KX + DO 40, J = N, 1, -1 + KX = KX - INCX + IF( X( JX ).NE.ZERO )THEN + IX = KX + L = KPLUS1 - J + IF( NOUNIT ) + $ X( JX ) = X( JX )/A( KPLUS1, J ) + TEMP = X( JX ) + DO 30, I = J - 1, MAX( 1, J - K ), -1 + X( IX ) = X( IX ) - TEMP*A( L + I, J ) + IX = IX - INCX + 30 CONTINUE + END IF + JX = JX - INCX + 40 CONTINUE + END IF + ELSE + IF( INCX.EQ.1 )THEN + DO 60, J = 1, N + IF( X( J ).NE.ZERO )THEN + L = 1 - J + IF( NOUNIT ) + $ X( J ) = X( J )/A( 1, J ) + TEMP = X( J ) + DO 50, I = J + 1, MIN( N, J + K ) + X( I ) = X( I ) - TEMP*A( L + I, J ) + 50 CONTINUE + END IF + 60 CONTINUE + ELSE + JX = KX + DO 80, J = 1, N + KX = KX + INCX + IF( X( JX ).NE.ZERO )THEN + IX = KX + L = 1 - J + IF( NOUNIT ) + $ X( JX ) = X( JX )/A( 1, J ) + TEMP = X( JX ) + DO 70, I = J + 1, MIN( N, J + K ) + X( IX ) = X( IX ) - TEMP*A( L + I, J ) + IX = IX + INCX + 70 CONTINUE + END IF + JX = JX + INCX + 80 CONTINUE + END IF + END IF + ELSE +* +* Form x := inv( A')*x. +* + IF( LSAME( UPLO, 'U' ) )THEN + KPLUS1 = K + 1 + IF( INCX.EQ.1 )THEN + DO 100, J = 1, N + TEMP = X( J ) + L = KPLUS1 - J + DO 90, I = MAX( 1, J - K ), J - 1 + TEMP = TEMP - A( L + I, J )*X( I ) + 90 CONTINUE + IF( NOUNIT ) + $ TEMP = TEMP/A( KPLUS1, J ) + X( J ) = TEMP + 100 CONTINUE + ELSE + JX = KX + DO 120, J = 1, N + TEMP = X( JX ) + IX = KX + L = KPLUS1 - J + DO 110, I = MAX( 1, J - K ), J - 1 + TEMP = TEMP - A( L + I, J )*X( IX ) + IX = IX + INCX + 110 CONTINUE + IF( NOUNIT ) + $ TEMP = TEMP/A( KPLUS1, J ) + X( JX ) = TEMP + JX = JX + INCX + IF( J.GT.K ) + $ KX = KX + INCX + 120 CONTINUE + END IF + ELSE + IF( INCX.EQ.1 )THEN + DO 140, J = N, 1, -1 + TEMP = X( J ) + L = 1 - J + DO 130, I = MIN( N, J + K ), J + 1, -1 + TEMP = TEMP - A( L + I, J )*X( I ) + 130 CONTINUE + IF( NOUNIT ) + $ TEMP = TEMP/A( 1, J ) + X( J ) = TEMP + 140 CONTINUE + ELSE + KX = KX + ( N - 1 )*INCX + JX = KX + DO 160, J = N, 1, -1 + TEMP = X( JX ) + IX = KX + L = 1 - J + DO 150, I = MIN( N, J + K ), J + 1, -1 + TEMP = TEMP - A( L + I, J )*X( IX ) + IX = IX - INCX + 150 CONTINUE + IF( NOUNIT ) + $ TEMP = TEMP/A( 1, J ) + X( JX ) = TEMP + JX = JX - INCX + IF( ( N - J ).GE.K ) + $ KX = KX - INCX + 160 CONTINUE + END IF + END IF + END IF +* + RETURN +* +* End of DTBSV . +* + END
--- a/libcruft/daspk/ddaspk.f Sun Feb 16 04:29:00 2003 +0000 +++ b/libcruft/daspk/ddaspk.f Tue Feb 18 20:08:20 2003 +0000 @@ -1326,7 +1326,7 @@ C DORTH performs orthogonalization of Krylov basis vectors. C DHEQR performs QR factorization of Hessenberg matrix. C DHELS finds least-squares solution of Hessenberg linear system. -C DGEFA, DGESL, DGBFA, DGBSL are LINPACK routines for solving +C DGETRF, DGETRS, DGBTRF, DGBTRS are LAPACK routines for solving C linear systems (dense or band direct methods). C DAXPY, DCOPY, DDOT, DNRM2, DSCAL are Basic Linear Algebra (BLAS) C routines.
--- a/libcruft/daspk/dmatd.f Sun Feb 16 04:29:00 2003 +0000 +++ b/libcruft/daspk/dmatd.f Tue Feb 18 20:08:20 2003 +0000 @@ -56,7 +56,7 @@ C They are not altered by DMATD. C----------------------------------------------------------------------- C***ROUTINES CALLED -C JACD, RES, DGEFA, DGBFA +C JACD, RES, DGETRF, DGBTRF C C***END PROLOGUE DMATD C @@ -111,7 +111,7 @@ C C Do dense-matrix LU decomposition on J. C -230 CALL DGEFA(WM,NEQ,NEQ,IWM(LIPVT),IER) +230 CALL DGETRF( NEQ, NEQ, WM, NEQ, IWM(LIPVT), IER) RETURN C C @@ -175,7 +175,8 @@ C C Do LU decomposition of banded J. C -550 CALL DGBFA (WM,MEBAND,NEQ,IWM(LML),IWM(LMU),IWM(LIPVT),IER) +550 CALL DGBTRF(NEQ, NEQ, IWM(LML), IWM(LMU), WM, MEBAND, + * IWM(LIPVT), IER) RETURN C C------END OF SUBROUTINE DMATD------------------------------------------
--- a/libcruft/daspk/dslvd.f Sun Feb 16 04:29:00 2003 +0000 +++ b/libcruft/daspk/dslvd.f Tue Feb 18 20:08:20 2003 +0000 @@ -18,11 +18,11 @@ C Real matrix information and real temporary storage C is stored in the array WM. C Integer matrix information is stored in the array IWM. -C For a dense matrix, the LINPACK routine DGESL is called. -C For a banded matrix, the LINPACK routine DGBSL is called. +C For a dense matrix, the LAPACK routine DGETRS is called. +C For a banded matrix, the LAPACK routine DGBTRS is called. C----------------------------------------------------------------------- C***ROUTINES CALLED -C DGESL, DGBSL +C DGETRS, DGBTRS C C***END PROLOGUE DSLVD C @@ -38,7 +38,7 @@ C C Dense matrix. C -100 CALL DGESL(WM,NEQ,NEQ,IWM(LIPVT),DELTA,0) +100 CALL DGETRS('N', NEQ, 1, WM, NEQ, IWM(LIPVT), DELTA, NEQ, INLPCK) RETURN C C Dummy section for MTYPE=3. @@ -49,8 +49,8 @@ C Banded matrix. C 400 MEBAND=2*IWM(LML)+IWM(LMU)+1 - CALL DGBSL(WM,MEBAND,NEQ,IWM(LML), - * IWM(LMU),IWM(LIPVT),DELTA,0) + CALL DGBTRS('N', NEQ, IWM(LML), IWM(LMU), 1, WM, MEBAND, + * IWM(LIPVT), DELTA, NEQ, INLPCK) RETURN C C------END OF SUBROUTINE DSLVD------------------------------------------
--- a/libcruft/dassl/ddajac.f Sun Feb 16 04:29:00 2003 +0000 +++ b/libcruft/dassl/ddajac.f Tue Feb 18 20:08:20 2003 +0000 @@ -44,7 +44,7 @@ C TO EVALUATE THE ITERATION MATRIX (THIS ROUTINE C IS ONLY USED IF IWM(MTYPE) IS 1 OR 4) C----------------------------------------------------------------------- -C***ROUTINES CALLED DGBFA, DGEFA +C***ROUTINES CALLED DGBTRF, DGETRF C***REVISION HISTORY (YYMMDD) C 830315 DATE WRITTEN C 901009 Finished conversion to SLATEC 4.0 format (F.N.Fritsch) @@ -53,6 +53,7 @@ C 901026 Added explicit declarations for all variables and minor C cosmetic changes to prologue. (FNF) C 901101 Corrected PURPOSE. (FNF) +C 020204 Convert to use LAPACK C***END PROLOGUE DDAJAC C INTEGER NEQ, IER, IWM(*), IRES, IPAR(*), NTEMP @@ -61,7 +62,7 @@ * UROUND, RPAR(*) EXTERNAL RES, JAC C - EXTERNAL DGBFA, DGEFA + EXTERNAL DGBTRF, DGETRF C INTEGER I, I1, I2, II, IPSAVE, ISAVE, J, K, L, LENPD, LIPVT, * LML, LMTYPE, LMU, MBA, MBAND, MEB1, MEBAND, MSAVE, MTYPE, N, @@ -113,7 +114,7 @@ C C C DO DENSE-MATRIX LU DECOMPOSITION ON PD -230 CALL DGEFA(WM(NPD),NEQ,NEQ,IWM(LIPVT),IER) +230 CALL DGETRF( NEQ, NEQ, WM(NPD), NEQ, IWM(LIPVT), IER) RETURN C C @@ -170,8 +171,8 @@ C C C DO LU DECOMPOSITION OF BANDED PD -550 CALL DGBFA(WM(NPD),MEBAND,NEQ, - * IWM(LML),IWM(LMU),IWM(LIPVT),IER) +550 CALL DGBTRF(NEQ, NEQ, IWM(LML), IWM(LMU), WM(NPD), MEBAND, + * IWM(LIPVT), IER) RETURN C------END OF SUBROUTINE DDAJAC------ END
--- a/libcruft/dassl/ddaslv.f Sun Feb 16 04:29:00 2003 +0000 +++ b/libcruft/dassl/ddaslv.f Tue Feb 18 20:08:20 2003 +0000 @@ -13,26 +13,27 @@ C REAL INFORMATION ARE STORED IN THE ARRAY WM. C INTEGER MATRIX INFORMATION IS STORED IN C THE ARRAY IWM. -C FOR A DENSE MATRIX, THE LINPACK ROUTINE -C DGESL IS CALLED. -C FOR A BANDED MATRIX,THE LINPACK ROUTINE -C DGBSL IS CALLED. +C FOR A DENSE MATRIX, THE LAPACK ROUTINE +C DGETRS IS CALLED. +C FOR A BANDED MATRIX,THE LAPACK ROUTINE +C DGBTRS IS CALLED. C----------------------------------------------------------------------- -C***ROUTINES CALLED DGBSL, DGESL +C***ROUTINES CALLED DGBTRS, DGETRF C***REVISION HISTORY (YYMMDD) C 830315 DATE WRITTEN C 901009 Finished conversion to SLATEC 4.0 format (F.N.Fritsch) C 901019 Merged changes made by C. Ulrich with SLATEC 4.0 format. C 901026 Added explicit declarations for all variables and minor C cosmetic changes to prologue. (FNF) +C 020204 Convert to use LAPACK C***END PROLOGUE DDASLV C INTEGER NEQ, IWM(*) DOUBLE PRECISION DELTA(*), WM(*) C - EXTERNAL DGBSL, DGESL + EXTERNAL DGBTRS, DGETRS C - INTEGER LIPVT, LML, LMU, LMTYPE, MEBAND, MTYPE, NPD + INTEGER LIPVT, LML, LMU, LMTYPE, MEBAND, MTYPE, NPD, INLAPACK PARAMETER (NPD=1) PARAMETER (LML=1) PARAMETER (LMU=2) @@ -44,7 +45,8 @@ GO TO(100,100,300,400,400),MTYPE C C DENSE MATRIX -100 CALL DGESL(WM(NPD),NEQ,NEQ,IWM(LIPVT),DELTA,0) +100 CALL DGETRS('N', NEQ, 1, WM(NPD), NEQ, IWM(LIPVT), DELTA, NEQ, + * INLPCK) RETURN C C DUMMY SECTION FOR MTYPE=3 @@ -53,8 +55,8 @@ C C BANDED MATRIX 400 MEBAND=2*IWM(LML)+IWM(LMU)+1 - CALL DGBSL(WM(NPD),MEBAND,NEQ,IWM(LML), - * IWM(LMU),IWM(LIPVT),DELTA,0) + CALL DGBTRS ('N', NEQ, IWM(LML), IWM(LMU), 1, WM(NPD), MEBAND, + * IWM(LIPVT), DELTA, NEQ, INLPCK); RETURN C------END OF SUBROUTINE DDASLV------ END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/lapack/dgbtf2.f Tue Feb 18 20:08:20 2003 +0000 @@ -0,0 +1,203 @@ + SUBROUTINE DGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO ) +* +* -- LAPACK 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 .. + INTEGER INFO, KL, KU, LDAB, M, N +* .. +* .. Array Arguments .. + INTEGER IPIV( * ) + DOUBLE PRECISION AB( LDAB, * ) +* .. +* +* Purpose +* ======= +* +* DGBTF2 computes an LU factorization of a real m-by-n band matrix A +* using partial pivoting with row interchanges. +* +* This is the unblocked version of the algorithm, calling Level 2 BLAS. +* +* Arguments +* ========= +* +* M (input) INTEGER +* The number of rows of the matrix A. M >= 0. +* +* N (input) INTEGER +* The number of columns of the matrix A. N >= 0. +* +* KL (input) INTEGER +* The number of subdiagonals within the band of A. KL >= 0. +* +* KU (input) INTEGER +* The number of superdiagonals within the band of A. KU >= 0. +* +* AB (input/output) DOUBLE PRECISION array, dimension (LDAB,N) +* On entry, the matrix A in band storage, in rows KL+1 to +* 2*KL+KU+1; rows 1 to KL of the array need not be set. +* The j-th column of A is stored in the j-th column of the +* array AB as follows: +* AB(kl+ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl) +* +* On exit, details of the factorization: U is stored as an +* upper triangular band matrix with KL+KU superdiagonals in +* rows 1 to KL+KU+1, and the multipliers used during the +* factorization are stored in rows KL+KU+2 to 2*KL+KU+1. +* See below for further details. +* +* LDAB (input) INTEGER +* The leading dimension of the array AB. LDAB >= 2*KL+KU+1. +* +* IPIV (output) INTEGER array, dimension (min(M,N)) +* The pivot indices; for 1 <= i <= min(M,N), row i of the +* matrix was interchanged with row IPIV(i). +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument had an illegal value +* > 0: if INFO = +i, U(i,i) is exactly zero. The factorization +* has been completed, but the factor U is exactly +* singular, and division by zero will occur if it is used +* to solve a system of equations. +* +* Further Details +* =============== +* +* The band storage scheme is illustrated by the following example, when +* M = N = 6, KL = 2, KU = 1: +* +* On entry: On exit: +* +* * * * + + + * * * u14 u25 u36 +* * * + + + + * * u13 u24 u35 u46 +* * a12 a23 a34 a45 a56 * u12 u23 u34 u45 u56 +* a11 a22 a33 a44 a55 a66 u11 u22 u33 u44 u55 u66 +* a21 a32 a43 a54 a65 * m21 m32 m43 m54 m65 * +* a31 a42 a53 a64 * * m31 m42 m53 m64 * * +* +* Array elements marked * are not used by the routine; elements marked +* + need not be set on entry, but are required by the routine to store +* elements of U, because of fill-in resulting from the row +* interchanges. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ONE, ZERO + PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) +* .. +* .. Local Scalars .. + INTEGER I, J, JP, JU, KM, KV +* .. +* .. External Functions .. + INTEGER IDAMAX + EXTERNAL IDAMAX +* .. +* .. External Subroutines .. + EXTERNAL DGER, DSCAL, DSWAP, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. Executable Statements .. +* +* KV is the number of superdiagonals in the factor U, allowing for +* fill-in. +* + KV = KU + KL +* +* Test the input parameters. +* + INFO = 0 + IF( M.LT.0 ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( KL.LT.0 ) THEN + INFO = -3 + ELSE IF( KU.LT.0 ) THEN + INFO = -4 + ELSE IF( LDAB.LT.KL+KV+1 ) THEN + INFO = -6 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DGBTF2', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( M.EQ.0 .OR. N.EQ.0 ) + $ RETURN +* +* Gaussian elimination with partial pivoting +* +* Set fill-in elements in columns KU+2 to KV to zero. +* + DO 20 J = KU + 2, MIN( KV, N ) + DO 10 I = KV - J + 2, KL + AB( I, J ) = ZERO + 10 CONTINUE + 20 CONTINUE +* +* JU is the index of the last column affected by the current stage +* of the factorization. +* + JU = 1 +* + DO 40 J = 1, MIN( M, N ) +* +* Set fill-in elements in column J+KV to zero. +* + IF( J+KV.LE.N ) THEN + DO 30 I = 1, KL + AB( I, J+KV ) = ZERO + 30 CONTINUE + END IF +* +* Find pivot and test for singularity. KM is the number of +* subdiagonal elements in the current column. +* + KM = MIN( KL, M-J ) + JP = IDAMAX( KM+1, AB( KV+1, J ), 1 ) + IPIV( J ) = JP + J - 1 + IF( AB( KV+JP, J ).NE.ZERO ) THEN + JU = MAX( JU, MIN( J+KU+JP-1, N ) ) +* +* Apply interchange to columns J to JU. +* + IF( JP.NE.1 ) + $ CALL DSWAP( JU-J+1, AB( KV+JP, J ), LDAB-1, + $ AB( KV+1, J ), LDAB-1 ) +* + IF( KM.GT.0 ) THEN +* +* Compute multipliers. +* + CALL DSCAL( KM, ONE / AB( KV+1, J ), AB( KV+2, J ), 1 ) +* +* Update trailing submatrix within the band. +* + IF( JU.GT.J ) + $ CALL DGER( KM, JU-J, -ONE, AB( KV+2, J ), 1, + $ AB( KV, J+1 ), LDAB-1, AB( KV+1, J+1 ), + $ LDAB-1 ) + END IF + ELSE +* +* If pivot is zero, set INFO to the index of the pivot +* unless a zero pivot has already been found. +* + IF( INFO.EQ.0 ) + $ INFO = J + END IF + 40 CONTINUE + RETURN +* +* End of DGBTF2 +* + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/lapack/dgbtrf.f Tue Feb 18 20:08:20 2003 +0000 @@ -0,0 +1,442 @@ + SUBROUTINE DGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) +* +* -- LAPACK 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 .. + INTEGER INFO, KL, KU, LDAB, M, N +* .. +* .. Array Arguments .. + INTEGER IPIV( * ) + DOUBLE PRECISION AB( LDAB, * ) +* .. +* +* Purpose +* ======= +* +* DGBTRF computes an LU factorization of a real m-by-n band matrix A +* using partial pivoting with row interchanges. +* +* This is the blocked version of the algorithm, calling Level 3 BLAS. +* +* Arguments +* ========= +* +* M (input) INTEGER +* The number of rows of the matrix A. M >= 0. +* +* N (input) INTEGER +* The number of columns of the matrix A. N >= 0. +* +* KL (input) INTEGER +* The number of subdiagonals within the band of A. KL >= 0. +* +* KU (input) INTEGER +* The number of superdiagonals within the band of A. KU >= 0. +* +* AB (input/output) DOUBLE PRECISION array, dimension (LDAB,N) +* On entry, the matrix A in band storage, in rows KL+1 to +* 2*KL+KU+1; rows 1 to KL of the array need not be set. +* The j-th column of A is stored in the j-th column of the +* array AB as follows: +* AB(kl+ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl) +* +* On exit, details of the factorization: U is stored as an +* upper triangular band matrix with KL+KU superdiagonals in +* rows 1 to KL+KU+1, and the multipliers used during the +* factorization are stored in rows KL+KU+2 to 2*KL+KU+1. +* See below for further details. +* +* LDAB (input) INTEGER +* The leading dimension of the array AB. LDAB >= 2*KL+KU+1. +* +* IPIV (output) INTEGER array, dimension (min(M,N)) +* The pivot indices; for 1 <= i <= min(M,N), row i of the +* matrix was interchanged with row IPIV(i). +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument had an illegal value +* > 0: if INFO = +i, U(i,i) is exactly zero. The factorization +* has been completed, but the factor U is exactly +* singular, and division by zero will occur if it is used +* to solve a system of equations. +* +* Further Details +* =============== +* +* The band storage scheme is illustrated by the following example, when +* M = N = 6, KL = 2, KU = 1: +* +* On entry: On exit: +* +* * * * + + + * * * u14 u25 u36 +* * * + + + + * * u13 u24 u35 u46 +* * a12 a23 a34 a45 a56 * u12 u23 u34 u45 u56 +* a11 a22 a33 a44 a55 a66 u11 u22 u33 u44 u55 u66 +* a21 a32 a43 a54 a65 * m21 m32 m43 m54 m65 * +* a31 a42 a53 a64 * * m31 m42 m53 m64 * * +* +* Array elements marked * are not used by the routine; elements marked +* + need not be set on entry, but are required by the routine to store +* elements of U because of fill-in resulting from the row interchanges. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ONE, ZERO + PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) + INTEGER NBMAX, LDWORK + PARAMETER ( NBMAX = 64, LDWORK = NBMAX+1 ) +* .. +* .. Local Scalars .. + INTEGER I, I2, I3, II, IP, J, J2, J3, JB, JJ, JM, JP, + $ JU, K2, KM, KV, NB, NW + DOUBLE PRECISION TEMP +* .. +* .. Local Arrays .. + DOUBLE PRECISION WORK13( LDWORK, NBMAX ), + $ WORK31( LDWORK, NBMAX ) +* .. +* .. External Functions .. + INTEGER IDAMAX, ILAENV + EXTERNAL IDAMAX, ILAENV +* .. +* .. External Subroutines .. + EXTERNAL DCOPY, DGBTF2, DGEMM, DGER, DLASWP, DSCAL, + $ DSWAP, DTRSM, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. Executable Statements .. +* +* KV is the number of superdiagonals in the factor U, allowing for +* fill-in +* + KV = KU + KL +* +* Test the input parameters. +* + INFO = 0 + IF( M.LT.0 ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( KL.LT.0 ) THEN + INFO = -3 + ELSE IF( KU.LT.0 ) THEN + INFO = -4 + ELSE IF( LDAB.LT.KL+KV+1 ) THEN + INFO = -6 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DGBTRF', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( M.EQ.0 .OR. N.EQ.0 ) + $ RETURN +* +* Determine the block size for this environment +* + NB = ILAENV( 1, 'DGBTRF', ' ', M, N, KL, KU ) +* +* The block size must not exceed the limit set by the size of the +* local arrays WORK13 and WORK31. +* + NB = MIN( NB, NBMAX ) +* + IF( NB.LE.1 .OR. NB.GT.KL ) THEN +* +* Use unblocked code +* + CALL DGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO ) + ELSE +* +* Use blocked code +* +* Zero the superdiagonal elements of the work array WORK13 +* + DO 20 J = 1, NB + DO 10 I = 1, J - 1 + WORK13( I, J ) = ZERO + 10 CONTINUE + 20 CONTINUE +* +* Zero the subdiagonal elements of the work array WORK31 +* + DO 40 J = 1, NB + DO 30 I = J + 1, NB + WORK31( I, J ) = ZERO + 30 CONTINUE + 40 CONTINUE +* +* Gaussian elimination with partial pivoting +* +* Set fill-in elements in columns KU+2 to KV to zero +* + DO 60 J = KU + 2, MIN( KV, N ) + DO 50 I = KV - J + 2, KL + AB( I, J ) = ZERO + 50 CONTINUE + 60 CONTINUE +* +* JU is the index of the last column affected by the current +* stage of the factorization +* + JU = 1 +* + DO 180 J = 1, MIN( M, N ), NB + JB = MIN( NB, MIN( M, N )-J+1 ) +* +* The active part of the matrix is partitioned +* +* A11 A12 A13 +* A21 A22 A23 +* A31 A32 A33 +* +* Here A11, A21 and A31 denote the current block of JB columns +* which is about to be factorized. The number of rows in the +* partitioning are JB, I2, I3 respectively, and the numbers +* of columns are JB, J2, J3. The superdiagonal elements of A13 +* and the subdiagonal elements of A31 lie outside the band. +* + I2 = MIN( KL-JB, M-J-JB+1 ) + I3 = MIN( JB, M-J-KL+1 ) +* +* J2 and J3 are computed after JU has been updated. +* +* Factorize the current block of JB columns +* + DO 80 JJ = J, J + JB - 1 +* +* Set fill-in elements in column JJ+KV to zero +* + IF( JJ+KV.LE.N ) THEN + DO 70 I = 1, KL + AB( I, JJ+KV ) = ZERO + 70 CONTINUE + END IF +* +* Find pivot and test for singularity. KM is the number of +* subdiagonal elements in the current column. +* + KM = MIN( KL, M-JJ ) + JP = IDAMAX( KM+1, AB( KV+1, JJ ), 1 ) + IPIV( JJ ) = JP + JJ - J + IF( AB( KV+JP, JJ ).NE.ZERO ) THEN + JU = MAX( JU, MIN( JJ+KU+JP-1, N ) ) + IF( JP.NE.1 ) THEN +* +* Apply interchange to columns J to J+JB-1 +* + IF( JP+JJ-1.LT.J+KL ) THEN +* + CALL DSWAP( JB, AB( KV+1+JJ-J, J ), LDAB-1, + $ AB( KV+JP+JJ-J, J ), LDAB-1 ) + ELSE +* +* The interchange affects columns J to JJ-1 of A31 +* which are stored in the work array WORK31 +* + CALL DSWAP( JJ-J, AB( KV+1+JJ-J, J ), LDAB-1, + $ WORK31( JP+JJ-J-KL, 1 ), LDWORK ) + CALL DSWAP( J+JB-JJ, AB( KV+1, JJ ), LDAB-1, + $ AB( KV+JP, JJ ), LDAB-1 ) + END IF + END IF +* +* Compute multipliers +* + CALL DSCAL( KM, ONE / AB( KV+1, JJ ), AB( KV+2, JJ ), + $ 1 ) +* +* Update trailing submatrix within the band and within +* the current block. JM is the index of the last column +* which needs to be updated. +* + JM = MIN( JU, J+JB-1 ) + IF( JM.GT.JJ ) + $ CALL DGER( KM, JM-JJ, -ONE, AB( KV+2, JJ ), 1, + $ AB( KV, JJ+1 ), LDAB-1, + $ AB( KV+1, JJ+1 ), LDAB-1 ) + ELSE +* +* If pivot is zero, set INFO to the index of the pivot +* unless a zero pivot has already been found. +* + IF( INFO.EQ.0 ) + $ INFO = JJ + END IF +* +* Copy current column of A31 into the work array WORK31 +* + NW = MIN( JJ-J+1, I3 ) + IF( NW.GT.0 ) + $ CALL DCOPY( NW, AB( KV+KL+1-JJ+J, JJ ), 1, + $ WORK31( 1, JJ-J+1 ), 1 ) + 80 CONTINUE + IF( J+JB.LE.N ) THEN +* +* Apply the row interchanges to the other blocks. +* + J2 = MIN( JU-J+1, KV ) - JB + J3 = MAX( 0, JU-J-KV+1 ) +* +* Use DLASWP to apply the row interchanges to A12, A22, and +* A32. +* + CALL DLASWP( J2, AB( KV+1-JB, J+JB ), LDAB-1, 1, JB, + $ IPIV( J ), 1 ) +* +* Adjust the pivot indices. +* + DO 90 I = J, J + JB - 1 + IPIV( I ) = IPIV( I ) + J - 1 + 90 CONTINUE +* +* Apply the row interchanges to A13, A23, and A33 +* columnwise. +* + K2 = J - 1 + JB + J2 + DO 110 I = 1, J3 + JJ = K2 + I + DO 100 II = J + I - 1, J + JB - 1 + IP = IPIV( II ) + IF( IP.NE.II ) THEN + TEMP = AB( KV+1+II-JJ, JJ ) + AB( KV+1+II-JJ, JJ ) = AB( KV+1+IP-JJ, JJ ) + AB( KV+1+IP-JJ, JJ ) = TEMP + END IF + 100 CONTINUE + 110 CONTINUE +* +* Update the relevant part of the trailing submatrix +* + IF( J2.GT.0 ) THEN +* +* Update A12 +* + CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Unit', + $ JB, J2, ONE, AB( KV+1, J ), LDAB-1, + $ AB( KV+1-JB, J+JB ), LDAB-1 ) +* + IF( I2.GT.0 ) THEN +* +* Update A22 +* + CALL DGEMM( 'No transpose', 'No transpose', I2, J2, + $ JB, -ONE, AB( KV+1+JB, J ), LDAB-1, + $ AB( KV+1-JB, J+JB ), LDAB-1, ONE, + $ AB( KV+1, J+JB ), LDAB-1 ) + END IF +* + IF( I3.GT.0 ) THEN +* +* Update A32 +* + CALL DGEMM( 'No transpose', 'No transpose', I3, J2, + $ JB, -ONE, WORK31, LDWORK, + $ AB( KV+1-JB, J+JB ), LDAB-1, ONE, + $ AB( KV+KL+1-JB, J+JB ), LDAB-1 ) + END IF + END IF +* + IF( J3.GT.0 ) THEN +* +* Copy the lower triangle of A13 into the work array +* WORK13 +* + DO 130 JJ = 1, J3 + DO 120 II = JJ, JB + WORK13( II, JJ ) = AB( II-JJ+1, JJ+J+KV-1 ) + 120 CONTINUE + 130 CONTINUE +* +* Update A13 in the work array +* + CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Unit', + $ JB, J3, ONE, AB( KV+1, J ), LDAB-1, + $ WORK13, LDWORK ) +* + IF( I2.GT.0 ) THEN +* +* Update A23 +* + CALL DGEMM( 'No transpose', 'No transpose', I2, J3, + $ JB, -ONE, AB( KV+1+JB, J ), LDAB-1, + $ WORK13, LDWORK, ONE, AB( 1+JB, J+KV ), + $ LDAB-1 ) + END IF +* + IF( I3.GT.0 ) THEN +* +* Update A33 +* + CALL DGEMM( 'No transpose', 'No transpose', I3, J3, + $ JB, -ONE, WORK31, LDWORK, WORK13, + $ LDWORK, ONE, AB( 1+KL, J+KV ), LDAB-1 ) + END IF +* +* Copy the lower triangle of A13 back into place +* + DO 150 JJ = 1, J3 + DO 140 II = JJ, JB + AB( II-JJ+1, JJ+J+KV-1 ) = WORK13( II, JJ ) + 140 CONTINUE + 150 CONTINUE + END IF + ELSE +* +* Adjust the pivot indices. +* + DO 160 I = J, J + JB - 1 + IPIV( I ) = IPIV( I ) + J - 1 + 160 CONTINUE + END IF +* +* Partially undo the interchanges in the current block to +* restore the upper triangular form of A31 and copy the upper +* triangle of A31 back into place +* + DO 170 JJ = J + JB - 1, J, -1 + JP = IPIV( JJ ) - JJ + 1 + IF( JP.NE.1 ) THEN +* +* Apply interchange to columns J to JJ-1 +* + IF( JP+JJ-1.LT.J+KL ) THEN +* +* The interchange does not affect A31 +* + CALL DSWAP( JJ-J, AB( KV+1+JJ-J, J ), LDAB-1, + $ AB( KV+JP+JJ-J, J ), LDAB-1 ) + ELSE +* +* The interchange does affect A31 +* + CALL DSWAP( JJ-J, AB( KV+1+JJ-J, J ), LDAB-1, + $ WORK31( JP+JJ-J-KL, 1 ), LDWORK ) + END IF + END IF +* +* Copy the current column of A31 back into place +* + NW = MIN( I3, JJ-J+1 ) + IF( NW.GT.0 ) + $ CALL DCOPY( NW, WORK31( 1, JJ-J+1 ), 1, + $ AB( KV+KL+1-JJ+J, JJ ), 1 ) + 170 CONTINUE + 180 CONTINUE + END IF +* + RETURN +* +* End of DGBTRF +* + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/lapack/dgbtrs.f Tue Feb 18 20:08:20 2003 +0000 @@ -0,0 +1,187 @@ + SUBROUTINE DGBTRS( TRANS, N, KL, KU, NRHS, AB, LDAB, IPIV, B, LDB, + $ 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 TRANS + INTEGER INFO, KL, KU, LDAB, LDB, N, NRHS +* .. +* .. Array Arguments .. + INTEGER IPIV( * ) + DOUBLE PRECISION AB( LDAB, * ), B( LDB, * ) +* .. +* +* Purpose +* ======= +* +* DGBTRS solves a system of linear equations +* A * X = B or A' * X = B +* with a general band matrix A using the LU factorization computed +* by DGBTRF. +* +* Arguments +* ========= +* +* TRANS (input) CHARACTER*1 +* Specifies the form of the system of equations. +* = 'N': A * X = B (No transpose) +* = 'T': A'* X = B (Transpose) +* = 'C': A'* X = B (Conjugate transpose = Transpose) +* +* N (input) INTEGER +* The order of the matrix A. N >= 0. +* +* KL (input) INTEGER +* The number of subdiagonals within the band of A. KL >= 0. +* +* KU (input) INTEGER +* The number of superdiagonals within the band of A. KU >= 0. +* +* NRHS (input) INTEGER +* The number of right hand sides, i.e., the number of columns +* of the matrix B. NRHS >= 0. +* +* AB (input) DOUBLE PRECISION array, dimension (LDAB,N) +* Details of the LU factorization of the band matrix A, as +* computed by DGBTRF. U is stored as an upper triangular band +* matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and +* the multipliers used during the factorization are stored in +* rows KL+KU+2 to 2*KL+KU+1. +* +* LDAB (input) INTEGER +* The leading dimension of the array AB. LDAB >= 2*KL+KU+1. +* +* IPIV (input) INTEGER array, dimension (N) +* The pivot indices; for 1 <= i <= N, row i of the matrix was +* interchanged with row IPIV(i). +* +* B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) +* On entry, the right hand side matrix B. +* On exit, the solution matrix X. +* +* LDB (input) INTEGER +* The leading dimension of the array B. LDB >= max(1,N). +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument had an illegal value +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ONE + PARAMETER ( ONE = 1.0D+0 ) +* .. +* .. Local Scalars .. + LOGICAL LNOTI, NOTRAN + INTEGER I, J, KD, L, LM +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. External Subroutines .. + EXTERNAL DGEMV, DGER, DSWAP, DTBSV, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. Executable Statements .. +* +* Test the input parameters. +* + INFO = 0 + NOTRAN = LSAME( TRANS, 'N' ) + IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT. + $ LSAME( TRANS, 'C' ) ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( KL.LT.0 ) THEN + INFO = -3 + ELSE IF( KU.LT.0 ) THEN + INFO = -4 + ELSE IF( NRHS.LT.0 ) THEN + INFO = -5 + ELSE IF( LDAB.LT.( 2*KL+KU+1 ) ) THEN + INFO = -7 + ELSE IF( LDB.LT.MAX( 1, N ) ) THEN + INFO = -10 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DGBTRS', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 .OR. NRHS.EQ.0 ) + $ RETURN +* + KD = KU + KL + 1 + LNOTI = KL.GT.0 +* + IF( NOTRAN ) THEN +* +* Solve A*X = B. +* +* Solve L*X = B, overwriting B with X. +* +* L is represented as a product of permutations and unit lower +* triangular matrices L = P(1) * L(1) * ... * P(n-1) * L(n-1), +* where each transformation L(i) is a rank-one modification of +* the identity matrix. +* + IF( LNOTI ) THEN + DO 10 J = 1, N - 1 + LM = MIN( KL, N-J ) + L = IPIV( J ) + IF( L.NE.J ) + $ CALL DSWAP( NRHS, B( L, 1 ), LDB, B( J, 1 ), LDB ) + CALL DGER( LM, NRHS, -ONE, AB( KD+1, J ), 1, B( J, 1 ), + $ LDB, B( J+1, 1 ), LDB ) + 10 CONTINUE + END IF +* + DO 20 I = 1, NRHS +* +* Solve U*X = B, overwriting B with X. +* + CALL DTBSV( 'Upper', 'No transpose', 'Non-unit', N, KL+KU, + $ AB, LDAB, B( 1, I ), 1 ) + 20 CONTINUE +* + ELSE +* +* Solve A'*X = B. +* + DO 30 I = 1, NRHS +* +* Solve U'*X = B, overwriting B with X. +* + CALL DTBSV( 'Upper', 'Transpose', 'Non-unit', N, KL+KU, AB, + $ LDAB, B( 1, I ), 1 ) + 30 CONTINUE +* +* Solve L'*X = B, overwriting B with X. +* + IF( LNOTI ) THEN + DO 40 J = N - 1, 1, -1 + LM = MIN( KL, N-J ) + CALL DGEMV( 'Transpose', LM, NRHS, -ONE, B( J+1, 1 ), + $ LDB, AB( KD+1, J ), 1, ONE, B( J, 1 ), LDB ) + L = IPIV( J ) + IF( L.NE.J ) + $ CALL DSWAP( NRHS, B( L, 1 ), LDB, B( J, 1 ), LDB ) + 40 CONTINUE + END IF + END IF + RETURN +* +* End of DGBTRS +* + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/lapack/dgecon.f Tue Feb 18 20:08:20 2003 +0000 @@ -0,0 +1,181 @@ + SUBROUTINE DGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, IWORK, + $ INFO ) +* +* -- LAPACK 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 NORM + INTEGER INFO, LDA, N + DOUBLE PRECISION ANORM, RCOND +* .. +* .. Array Arguments .. + INTEGER IWORK( * ) + DOUBLE PRECISION A( LDA, * ), WORK( * ) +* .. +* +* Purpose +* ======= +* +* DGECON estimates the reciprocal of the condition number of a general +* real matrix A, in either the 1-norm or the infinity-norm, using +* the LU factorization computed by DGETRF. +* +* An estimate is obtained for norm(inv(A)), and the reciprocal of the +* condition number is computed as +* RCOND = 1 / ( norm(A) * norm(inv(A)) ). +* +* Arguments +* ========= +* +* NORM (input) CHARACTER*1 +* Specifies whether the 1-norm condition number or the +* infinity-norm condition number is required: +* = '1' or 'O': 1-norm; +* = 'I': Infinity-norm. +* +* N (input) INTEGER +* The order of the matrix A. N >= 0. +* +* A (input) DOUBLE PRECISION array, dimension (LDA,N) +* The factors L and U from the factorization A = P*L*U +* as computed by DGETRF. +* +* LDA (input) INTEGER +* The leading dimension of the array A. LDA >= max(1,N). +* +* ANORM (input) DOUBLE PRECISION +* If NORM = '1' or 'O', the 1-norm of the original matrix A. +* If NORM = 'I', the infinity-norm of the original matrix A. +* +* RCOND (output) DOUBLE PRECISION +* The reciprocal of the condition number of the matrix A, +* computed as RCOND = 1/(norm(A) * norm(inv(A))). +* +* WORK (workspace) DOUBLE PRECISION array, dimension (4*N) +* +* IWORK (workspace) INTEGER array, dimension (N) +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument had an illegal value +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ONE, ZERO + PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) +* .. +* .. Local Scalars .. + LOGICAL ONENRM + CHARACTER NORMIN + INTEGER IX, KASE, KASE1 + DOUBLE PRECISION AINVNM, SCALE, SL, SMLNUM, SU +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER IDAMAX + DOUBLE PRECISION DLAMCH + EXTERNAL LSAME, IDAMAX, DLAMCH +* .. +* .. External Subroutines .. + EXTERNAL DLACON, DLATRS, DRSCL, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, MAX +* .. +* .. Executable Statements .. +* +* Test the input parameters. +* + INFO = 0 + ONENRM = NORM.EQ.'1' .OR. LSAME( NORM, 'O' ) + IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM, 'I' ) ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -4 + ELSE IF( ANORM.LT.ZERO ) THEN + INFO = -5 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DGECON', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + RCOND = ZERO + IF( N.EQ.0 ) THEN + RCOND = ONE + RETURN + ELSE IF( ANORM.EQ.ZERO ) THEN + RETURN + END IF +* + SMLNUM = DLAMCH( 'Safe minimum' ) +* +* Estimate the norm of inv(A). +* + AINVNM = ZERO + NORMIN = 'N' + IF( ONENRM ) THEN + KASE1 = 1 + ELSE + KASE1 = 2 + END IF + KASE = 0 + 10 CONTINUE + CALL DLACON( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE ) + IF( KASE.NE.0 ) THEN + IF( KASE.EQ.KASE1 ) THEN +* +* Multiply by inv(L). +* + CALL DLATRS( 'Lower', 'No transpose', 'Unit', NORMIN, N, A, + $ LDA, WORK, SL, WORK( 2*N+1 ), INFO ) +* +* Multiply by inv(U). +* + CALL DLATRS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N, + $ A, LDA, WORK, SU, WORK( 3*N+1 ), INFO ) + ELSE +* +* Multiply by inv(U'). +* + CALL DLATRS( 'Upper', 'Transpose', 'Non-unit', NORMIN, N, A, + $ LDA, WORK, SU, WORK( 3*N+1 ), INFO ) +* +* Multiply by inv(L'). +* + CALL DLATRS( 'Lower', 'Transpose', 'Unit', NORMIN, N, A, + $ LDA, WORK, SL, WORK( 2*N+1 ), INFO ) + END IF +* +* Divide X by 1/(SL*SU) if doing so will not cause overflow. +* + SCALE = SL*SU + NORMIN = 'Y' + IF( SCALE.NE.ONE ) THEN + IX = IDAMAX( N, WORK, 1 ) + IF( SCALE.LT.ABS( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO ) + $ GO TO 20 + CALL DRSCL( N, SCALE, WORK, 1 ) + END IF + GO TO 10 + END IF +* +* Compute the estimate of the reciprocal condition number. +* + IF( AINVNM.NE.ZERO ) + $ RCOND = ( ONE / AINVNM ) / ANORM +* + 20 CONTINUE + RETURN +* +* End of DGECON +* + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/lapack/dgetri.f Tue Feb 18 20:08:20 2003 +0000 @@ -0,0 +1,193 @@ + SUBROUTINE DGETRI( N, A, LDA, IPIV, WORK, LWORK, INFO ) +* +* -- LAPACK routine (version 3.0) -- +* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., +* Courant Institute, Argonne National Lab, and Rice University +* June 30, 1999 +* +* .. Scalar Arguments .. + INTEGER INFO, LDA, LWORK, N +* .. +* .. Array Arguments .. + INTEGER IPIV( * ) + DOUBLE PRECISION A( LDA, * ), WORK( * ) +* .. +* +* Purpose +* ======= +* +* DGETRI computes the inverse of a matrix using the LU factorization +* computed by DGETRF. +* +* This method inverts U and then computes inv(A) by solving the system +* inv(A)*L = inv(U) for inv(A). +* +* Arguments +* ========= +* +* N (input) INTEGER +* The order of the matrix A. N >= 0. +* +* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) +* On entry, the factors L and U from the factorization +* A = P*L*U as computed by DGETRF. +* On exit, if INFO = 0, the inverse of the original matrix A. +* +* LDA (input) INTEGER +* The leading dimension of the array A. LDA >= max(1,N). +* +* IPIV (input) INTEGER array, dimension (N) +* The pivot indices from DGETRF; for 1<=i<=N, row i of the +* matrix was interchanged with row IPIV(i). +* +* WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) +* On exit, if INFO=0, then WORK(1) returns the optimal LWORK. +* +* LWORK (input) INTEGER +* The dimension of the array WORK. LWORK >= max(1,N). +* For optimal performance LWORK >= N*NB, where NB is +* the optimal blocksize returned by ILAENV. +* +* 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. +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument had an illegal value +* > 0: if INFO = i, U(i,i) is exactly zero; the matrix is +* singular and its inverse could not be computed. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) +* .. +* .. Local Scalars .. + LOGICAL LQUERY + INTEGER I, IWS, J, JB, JJ, JP, LDWORK, LWKOPT, NB, + $ NBMIN, NN +* .. +* .. External Functions .. + INTEGER ILAENV + EXTERNAL ILAENV +* .. +* .. External Subroutines .. + EXTERNAL DGEMM, DGEMV, DSWAP, DTRSM, DTRTRI, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. Executable Statements .. +* +* Test the input parameters. +* + INFO = 0 + NB = ILAENV( 1, 'DGETRI', ' ', N, -1, -1, -1 ) + LWKOPT = N*NB + WORK( 1 ) = LWKOPT + LQUERY = ( LWORK.EQ.-1 ) + IF( N.LT.0 ) THEN + INFO = -1 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -3 + ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN + INFO = -6 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DGETRI', -INFO ) + RETURN + ELSE IF( LQUERY ) THEN + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 ) + $ RETURN +* +* Form inv(U). If INFO > 0 from DTRTRI, then U is singular, +* and the inverse is not computed. +* + CALL DTRTRI( 'Upper', 'Non-unit', N, A, LDA, INFO ) + IF( INFO.GT.0 ) + $ RETURN +* + NBMIN = 2 + LDWORK = N + IF( NB.GT.1 .AND. NB.LT.N ) THEN + IWS = MAX( LDWORK*NB, 1 ) + IF( LWORK.LT.IWS ) THEN + NB = LWORK / LDWORK + NBMIN = MAX( 2, ILAENV( 2, 'DGETRI', ' ', N, -1, -1, -1 ) ) + END IF + ELSE + IWS = N + END IF +* +* Solve the equation inv(A)*L = inv(U) for inv(A). +* + IF( NB.LT.NBMIN .OR. NB.GE.N ) THEN +* +* Use unblocked code. +* + DO 20 J = N, 1, -1 +* +* Copy current column of L to WORK and replace with zeros. +* + DO 10 I = J + 1, N + WORK( I ) = A( I, J ) + A( I, J ) = ZERO + 10 CONTINUE +* +* Compute current column of inv(A). +* + IF( J.LT.N ) + $ CALL DGEMV( 'No transpose', N, N-J, -ONE, A( 1, J+1 ), + $ LDA, WORK( J+1 ), 1, ONE, A( 1, J ), 1 ) + 20 CONTINUE + ELSE +* +* Use blocked code. +* + NN = ( ( N-1 ) / NB )*NB + 1 + DO 50 J = NN, 1, -NB + JB = MIN( NB, N-J+1 ) +* +* Copy current block column of L to WORK and replace with +* zeros. +* + DO 40 JJ = J, J + JB - 1 + DO 30 I = JJ + 1, N + WORK( I+( JJ-J )*LDWORK ) = A( I, JJ ) + A( I, JJ ) = ZERO + 30 CONTINUE + 40 CONTINUE +* +* Compute current block column of inv(A). +* + IF( J+JB.LE.N ) + $ CALL DGEMM( 'No transpose', 'No transpose', N, JB, + $ N-J-JB+1, -ONE, A( 1, J+JB ), LDA, + $ WORK( J+JB ), LDWORK, ONE, A( 1, J ), LDA ) + CALL DTRSM( 'Right', 'Lower', 'No transpose', 'Unit', N, JB, + $ ONE, WORK( J ), LDWORK, A( 1, J ), LDA ) + 50 CONTINUE + END IF +* +* Apply column interchanges. +* + DO 60 J = N - 1, 1, -1 + JP = IPIV( J ) + IF( JP.NE.J ) + $ CALL DSWAP( N, A( 1, J ), 1, A( 1, JP ), 1 ) + 60 CONTINUE +* + WORK( 1 ) = IWS + RETURN +* +* End of DGETRI +* + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/lapack/dlatrs.f Tue Feb 18 20:08:20 2003 +0000 @@ -0,0 +1,702 @@ + SUBROUTINE DLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, + $ CNORM, INFO ) +* +* -- LAPACK auxiliary routine (version 3.0) -- +* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., +* Courant Institute, Argonne National Lab, and Rice University +* June 30, 1992 +* +* .. Scalar Arguments .. + CHARACTER DIAG, NORMIN, TRANS, UPLO + INTEGER INFO, LDA, N + DOUBLE PRECISION SCALE +* .. +* .. Array Arguments .. + DOUBLE PRECISION A( LDA, * ), CNORM( * ), X( * ) +* .. +* +* Purpose +* ======= +* +* DLATRS solves one of the triangular systems +* +* A *x = s*b or A'*x = s*b +* +* with scaling to prevent overflow. Here A is an upper or lower +* triangular matrix, A' denotes the transpose of A, x and b are +* n-element vectors, and s is a scaling factor, usually less than +* or equal to 1, chosen so that the components of x will be less than +* the overflow threshold. If the unscaled problem will not cause +* overflow, the Level 2 BLAS routine DTRSV is called. If the matrix A +* is singular (A(j,j) = 0 for some j), then s is set to 0 and a +* non-trivial solution to A*x = 0 is returned. +* +* Arguments +* ========= +* +* UPLO (input) CHARACTER*1 +* Specifies whether the matrix A is upper or lower triangular. +* = 'U': Upper triangular +* = 'L': Lower triangular +* +* TRANS (input) CHARACTER*1 +* Specifies the operation applied to A. +* = 'N': Solve A * x = s*b (No transpose) +* = 'T': Solve A'* x = s*b (Transpose) +* = 'C': Solve A'* x = s*b (Conjugate transpose = Transpose) +* +* DIAG (input) CHARACTER*1 +* Specifies whether or not the matrix A is unit triangular. +* = 'N': Non-unit triangular +* = 'U': Unit triangular +* +* NORMIN (input) CHARACTER*1 +* Specifies whether CNORM has been set or not. +* = 'Y': CNORM contains the column norms on entry +* = 'N': CNORM is not set on entry. On exit, the norms will +* be computed and stored in CNORM. +* +* N (input) INTEGER +* The order of the matrix A. N >= 0. +* +* A (input) DOUBLE PRECISION array, dimension (LDA,N) +* The triangular matrix A. If UPLO = 'U', the leading n by n +* upper triangular part of the array A contains the upper +* triangular matrix, and the strictly lower triangular part of +* A is not referenced. If UPLO = 'L', the leading n by n lower +* triangular part of the array A contains the lower triangular +* matrix, and the strictly upper triangular part of A is not +* referenced. If DIAG = 'U', the diagonal elements of A are +* also not referenced and are assumed to be 1. +* +* LDA (input) INTEGER +* The leading dimension of the array A. LDA >= max (1,N). +* +* X (input/output) DOUBLE PRECISION array, dimension (N) +* On entry, the right hand side b of the triangular system. +* On exit, X is overwritten by the solution vector x. +* +* SCALE (output) DOUBLE PRECISION +* The scaling factor s for the triangular system +* A * x = s*b or A'* x = s*b. +* If SCALE = 0, the matrix A is singular or badly scaled, and +* the vector x is an exact or approximate solution to A*x = 0. +* +* CNORM (input or output) DOUBLE PRECISION array, dimension (N) +* +* If NORMIN = 'Y', CNORM is an input argument and CNORM(j) +* contains the norm of the off-diagonal part of the j-th column +* of A. If TRANS = 'N', CNORM(j) must be greater than or equal +* to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j) +* must be greater than or equal to the 1-norm. +* +* If NORMIN = 'N', CNORM is an output argument and CNORM(j) +* returns the 1-norm of the offdiagonal part of the j-th column +* of A. +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -k, the k-th argument had an illegal value +* +* Further Details +* ======= ======= +* +* A rough bound on x is computed; if that is less than overflow, DTRSV +* is called, otherwise, specific code is used which checks for possible +* overflow or divide-by-zero at every operation. +* +* A columnwise scheme is used for solving A*x = b. The basic algorithm +* if A is lower triangular is +* +* x[1:n] := b[1:n] +* for j = 1, ..., n +* x(j) := x(j) / A(j,j) +* x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j] +* end +* +* Define bounds on the components of x after j iterations of the loop: +* M(j) = bound on x[1:j] +* G(j) = bound on x[j+1:n] +* Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}. +* +* Then for iteration j+1 we have +* M(j+1) <= G(j) / | A(j+1,j+1) | +* G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] | +* <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | ) +* +* where CNORM(j+1) is greater than or equal to the infinity-norm of +* column j+1 of A, not counting the diagonal. Hence +* +* G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | ) +* 1<=i<=j +* and +* +* |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| ) +* 1<=i< j +* +* Since |x(j)| <= M(j), we use the Level 2 BLAS routine DTRSV if the +* reciprocal of the largest M(j), j=1,..,n, is larger than +* max(underflow, 1/overflow). +* +* The bound on x(j) is also used to determine when a step in the +* columnwise method can be performed without fear of overflow. If +* the computed bound is greater than a large constant, x is scaled to +* prevent overflow, but if the bound overflows, x is set to 0, x(j) to +* 1, and scale to 0, and a non-trivial solution to A*x = 0 is found. +* +* Similarly, a row-wise scheme is used to solve A'*x = b. The basic +* algorithm for A upper triangular is +* +* for j = 1, ..., n +* x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j) +* end +* +* We simultaneously compute two bounds +* G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j +* M(j) = bound on x(i), 1<=i<=j +* +* The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we +* add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1. +* Then the bound on x(j) is +* +* M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) | +* +* <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| ) +* 1<=i<=j +* +* and we can safely call DTRSV if 1/M(n) and 1/G(n) are both greater +* than max(underflow, 1/overflow). +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ZERO, HALF, ONE + PARAMETER ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0 ) +* .. +* .. Local Scalars .. + LOGICAL NOTRAN, NOUNIT, UPPER + INTEGER I, IMAX, J, JFIRST, JINC, JLAST + DOUBLE PRECISION BIGNUM, GROW, REC, SMLNUM, SUMJ, TJJ, TJJS, + $ TMAX, TSCAL, USCAL, XBND, XJ, XMAX +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER IDAMAX + DOUBLE PRECISION DASUM, DDOT, DLAMCH + EXTERNAL LSAME, IDAMAX, DASUM, DDOT, DLAMCH +* .. +* .. External Subroutines .. + EXTERNAL DAXPY, DSCAL, DTRSV, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, MAX, MIN +* .. +* .. Executable Statements .. +* + INFO = 0 + UPPER = LSAME( UPLO, 'U' ) + NOTRAN = LSAME( TRANS, 'N' ) + NOUNIT = LSAME( DIAG, 'N' ) +* +* Test the input parameters. +* + IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN + INFO = -1 + ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT. + $ LSAME( TRANS, 'C' ) ) THEN + INFO = -2 + ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN + INFO = -3 + ELSE IF( .NOT.LSAME( NORMIN, 'Y' ) .AND. .NOT. + $ LSAME( NORMIN, 'N' ) ) THEN + INFO = -4 + ELSE IF( N.LT.0 ) THEN + INFO = -5 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -7 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DLATRS', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 ) + $ RETURN +* +* Determine machine dependent parameters to control overflow. +* + SMLNUM = DLAMCH( 'Safe minimum' ) / DLAMCH( 'Precision' ) + BIGNUM = ONE / SMLNUM + SCALE = ONE +* + IF( LSAME( NORMIN, 'N' ) ) THEN +* +* Compute the 1-norm of each column, not including the diagonal. +* + IF( UPPER ) THEN +* +* A is upper triangular. +* + DO 10 J = 1, N + CNORM( J ) = DASUM( J-1, A( 1, J ), 1 ) + 10 CONTINUE + ELSE +* +* A is lower triangular. +* + DO 20 J = 1, N - 1 + CNORM( J ) = DASUM( N-J, A( J+1, J ), 1 ) + 20 CONTINUE + CNORM( N ) = ZERO + END IF + END IF +* +* Scale the column norms by TSCAL if the maximum element in CNORM is +* greater than BIGNUM. +* + IMAX = IDAMAX( N, CNORM, 1 ) + TMAX = CNORM( IMAX ) + IF( TMAX.LE.BIGNUM ) THEN + TSCAL = ONE + ELSE + TSCAL = ONE / ( SMLNUM*TMAX ) + CALL DSCAL( N, TSCAL, CNORM, 1 ) + END IF +* +* Compute a bound on the computed solution vector to see if the +* Level 2 BLAS routine DTRSV can be used. +* + J = IDAMAX( N, X, 1 ) + XMAX = ABS( X( J ) ) + XBND = XMAX + IF( NOTRAN ) THEN +* +* Compute the growth in A * x = b. +* + IF( UPPER ) THEN + JFIRST = N + JLAST = 1 + JINC = -1 + ELSE + JFIRST = 1 + JLAST = N + JINC = 1 + END IF +* + IF( TSCAL.NE.ONE ) THEN + GROW = ZERO + GO TO 50 + END IF +* + IF( NOUNIT ) THEN +* +* A is non-unit triangular. +* +* Compute GROW = 1/G(j) and XBND = 1/M(j). +* Initially, G(0) = max{x(i), i=1,...,n}. +* + GROW = ONE / MAX( XBND, SMLNUM ) + XBND = GROW + DO 30 J = JFIRST, JLAST, JINC +* +* Exit the loop if the growth factor is too small. +* + IF( GROW.LE.SMLNUM ) + $ GO TO 50 +* +* M(j) = G(j-1) / abs(A(j,j)) +* + TJJ = ABS( A( J, J ) ) + XBND = MIN( XBND, MIN( ONE, TJJ )*GROW ) + IF( TJJ+CNORM( J ).GE.SMLNUM ) THEN +* +* G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) ) +* + GROW = GROW*( TJJ / ( TJJ+CNORM( J ) ) ) + ELSE +* +* G(j) could overflow, set GROW to 0. +* + GROW = ZERO + END IF + 30 CONTINUE + GROW = XBND + ELSE +* +* A is unit triangular. +* +* Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}. +* + GROW = MIN( ONE, ONE / MAX( XBND, SMLNUM ) ) + DO 40 J = JFIRST, JLAST, JINC +* +* Exit the loop if the growth factor is too small. +* + IF( GROW.LE.SMLNUM ) + $ GO TO 50 +* +* G(j) = G(j-1)*( 1 + CNORM(j) ) +* + GROW = GROW*( ONE / ( ONE+CNORM( J ) ) ) + 40 CONTINUE + END IF + 50 CONTINUE +* + ELSE +* +* Compute the growth in A' * x = b. +* + IF( UPPER ) THEN + JFIRST = 1 + JLAST = N + JINC = 1 + ELSE + JFIRST = N + JLAST = 1 + JINC = -1 + END IF +* + IF( TSCAL.NE.ONE ) THEN + GROW = ZERO + GO TO 80 + END IF +* + IF( NOUNIT ) THEN +* +* A is non-unit triangular. +* +* Compute GROW = 1/G(j) and XBND = 1/M(j). +* Initially, M(0) = max{x(i), i=1,...,n}. +* + GROW = ONE / MAX( XBND, SMLNUM ) + XBND = GROW + DO 60 J = JFIRST, JLAST, JINC +* +* Exit the loop if the growth factor is too small. +* + IF( GROW.LE.SMLNUM ) + $ GO TO 80 +* +* G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) ) +* + XJ = ONE + CNORM( J ) + GROW = MIN( GROW, XBND / XJ ) +* +* M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j)) +* + TJJ = ABS( A( J, J ) ) + IF( XJ.GT.TJJ ) + $ XBND = XBND*( TJJ / XJ ) + 60 CONTINUE + GROW = MIN( GROW, XBND ) + ELSE +* +* A is unit triangular. +* +* Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}. +* + GROW = MIN( ONE, ONE / MAX( XBND, SMLNUM ) ) + DO 70 J = JFIRST, JLAST, JINC +* +* Exit the loop if the growth factor is too small. +* + IF( GROW.LE.SMLNUM ) + $ GO TO 80 +* +* G(j) = ( 1 + CNORM(j) )*G(j-1) +* + XJ = ONE + CNORM( J ) + GROW = GROW / XJ + 70 CONTINUE + END IF + 80 CONTINUE + END IF +* + IF( ( GROW*TSCAL ).GT.SMLNUM ) THEN +* +* Use the Level 2 BLAS solve if the reciprocal of the bound on +* elements of X is not too small. +* + CALL DTRSV( UPLO, TRANS, DIAG, N, A, LDA, X, 1 ) + ELSE +* +* Use a Level 1 BLAS solve, scaling intermediate results. +* + IF( XMAX.GT.BIGNUM ) THEN +* +* Scale X so that its components are less than or equal to +* BIGNUM in absolute value. +* + SCALE = BIGNUM / XMAX + CALL DSCAL( N, SCALE, X, 1 ) + XMAX = BIGNUM + END IF +* + IF( NOTRAN ) THEN +* +* Solve A * x = b +* + DO 110 J = JFIRST, JLAST, JINC +* +* Compute x(j) = b(j) / A(j,j), scaling x if necessary. +* + XJ = ABS( X( J ) ) + IF( NOUNIT ) THEN + TJJS = A( J, J )*TSCAL + ELSE + TJJS = TSCAL + IF( TSCAL.EQ.ONE ) + $ GO TO 100 + END IF + TJJ = ABS( TJJS ) + IF( TJJ.GT.SMLNUM ) THEN +* +* abs(A(j,j)) > SMLNUM: +* + IF( TJJ.LT.ONE ) THEN + IF( XJ.GT.TJJ*BIGNUM ) THEN +* +* Scale x by 1/b(j). +* + REC = ONE / XJ + CALL DSCAL( N, REC, X, 1 ) + SCALE = SCALE*REC + XMAX = XMAX*REC + END IF + END IF + X( J ) = X( J ) / TJJS + XJ = ABS( X( J ) ) + ELSE IF( TJJ.GT.ZERO ) THEN +* +* 0 < abs(A(j,j)) <= SMLNUM: +* + IF( XJ.GT.TJJ*BIGNUM ) THEN +* +* Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM +* to avoid overflow when dividing by A(j,j). +* + REC = ( TJJ*BIGNUM ) / XJ + IF( CNORM( J ).GT.ONE ) THEN +* +* Scale by 1/CNORM(j) to avoid overflow when +* multiplying x(j) times column j. +* + REC = REC / CNORM( J ) + END IF + CALL DSCAL( N, REC, X, 1 ) + SCALE = SCALE*REC + XMAX = XMAX*REC + END IF + X( J ) = X( J ) / TJJS + XJ = ABS( X( J ) ) + ELSE +* +* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and +* scale = 0, and compute a solution to A*x = 0. +* + DO 90 I = 1, N + X( I ) = ZERO + 90 CONTINUE + X( J ) = ONE + XJ = ONE + SCALE = ZERO + XMAX = ZERO + END IF + 100 CONTINUE +* +* Scale x if necessary to avoid overflow when adding a +* multiple of column j of A. +* + IF( XJ.GT.ONE ) THEN + REC = ONE / XJ + IF( CNORM( J ).GT.( BIGNUM-XMAX )*REC ) THEN +* +* Scale x by 1/(2*abs(x(j))). +* + REC = REC*HALF + CALL DSCAL( N, REC, X, 1 ) + SCALE = SCALE*REC + END IF + ELSE IF( XJ*CNORM( J ).GT.( BIGNUM-XMAX ) ) THEN +* +* Scale x by 1/2. +* + CALL DSCAL( N, HALF, X, 1 ) + SCALE = SCALE*HALF + END IF +* + IF( UPPER ) THEN + IF( J.GT.1 ) THEN +* +* Compute the update +* x(1:j-1) := x(1:j-1) - x(j) * A(1:j-1,j) +* + CALL DAXPY( J-1, -X( J )*TSCAL, A( 1, J ), 1, X, + $ 1 ) + I = IDAMAX( J-1, X, 1 ) + XMAX = ABS( X( I ) ) + END IF + ELSE + IF( J.LT.N ) THEN +* +* Compute the update +* x(j+1:n) := x(j+1:n) - x(j) * A(j+1:n,j) +* + CALL DAXPY( N-J, -X( J )*TSCAL, A( J+1, J ), 1, + $ X( J+1 ), 1 ) + I = J + IDAMAX( N-J, X( J+1 ), 1 ) + XMAX = ABS( X( I ) ) + END IF + END IF + 110 CONTINUE +* + ELSE +* +* Solve A' * x = b +* + DO 160 J = JFIRST, JLAST, JINC +* +* Compute x(j) = b(j) - sum A(k,j)*x(k). +* k<>j +* + XJ = ABS( X( J ) ) + USCAL = TSCAL + REC = ONE / MAX( XMAX, ONE ) + IF( CNORM( J ).GT.( BIGNUM-XJ )*REC ) THEN +* +* If x(j) could overflow, scale x by 1/(2*XMAX). +* + REC = REC*HALF + IF( NOUNIT ) THEN + TJJS = A( J, J )*TSCAL + ELSE + TJJS = TSCAL + END IF + TJJ = ABS( TJJS ) + IF( TJJ.GT.ONE ) THEN +* +* Divide by A(j,j) when scaling x if A(j,j) > 1. +* + REC = MIN( ONE, REC*TJJ ) + USCAL = USCAL / TJJS + END IF + IF( REC.LT.ONE ) THEN + CALL DSCAL( N, REC, X, 1 ) + SCALE = SCALE*REC + XMAX = XMAX*REC + END IF + END IF +* + SUMJ = ZERO + IF( USCAL.EQ.ONE ) THEN +* +* If the scaling needed for A in the dot product is 1, +* call DDOT to perform the dot product. +* + IF( UPPER ) THEN + SUMJ = DDOT( J-1, A( 1, J ), 1, X, 1 ) + ELSE IF( J.LT.N ) THEN + SUMJ = DDOT( N-J, A( J+1, J ), 1, X( J+1 ), 1 ) + END IF + ELSE +* +* Otherwise, use in-line code for the dot product. +* + IF( UPPER ) THEN + DO 120 I = 1, J - 1 + SUMJ = SUMJ + ( A( I, J )*USCAL )*X( I ) + 120 CONTINUE + ELSE IF( J.LT.N ) THEN + DO 130 I = J + 1, N + SUMJ = SUMJ + ( A( I, J )*USCAL )*X( I ) + 130 CONTINUE + END IF + END IF +* + IF( USCAL.EQ.TSCAL ) THEN +* +* Compute x(j) := ( x(j) - sumj ) / A(j,j) if 1/A(j,j) +* was not used to scale the dotproduct. +* + X( J ) = X( J ) - SUMJ + XJ = ABS( X( J ) ) + IF( NOUNIT ) THEN + TJJS = A( J, J )*TSCAL + ELSE + TJJS = TSCAL + IF( TSCAL.EQ.ONE ) + $ GO TO 150 + END IF +* +* Compute x(j) = x(j) / A(j,j), scaling if necessary. +* + TJJ = ABS( TJJS ) + IF( TJJ.GT.SMLNUM ) THEN +* +* abs(A(j,j)) > SMLNUM: +* + IF( TJJ.LT.ONE ) THEN + IF( XJ.GT.TJJ*BIGNUM ) THEN +* +* Scale X by 1/abs(x(j)). +* + REC = ONE / XJ + CALL DSCAL( N, REC, X, 1 ) + SCALE = SCALE*REC + XMAX = XMAX*REC + END IF + END IF + X( J ) = X( J ) / TJJS + ELSE IF( TJJ.GT.ZERO ) THEN +* +* 0 < abs(A(j,j)) <= SMLNUM: +* + IF( XJ.GT.TJJ*BIGNUM ) THEN +* +* Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM. +* + REC = ( TJJ*BIGNUM ) / XJ + CALL DSCAL( N, REC, X, 1 ) + SCALE = SCALE*REC + XMAX = XMAX*REC + END IF + X( J ) = X( J ) / TJJS + ELSE +* +* A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and +* scale = 0, and compute a solution to A'*x = 0. +* + DO 140 I = 1, N + X( I ) = ZERO + 140 CONTINUE + X( J ) = ONE + SCALE = ZERO + XMAX = ZERO + END IF + 150 CONTINUE + ELSE +* +* Compute x(j) := x(j) / A(j,j) - sumj if the dot +* product has already been divided by 1/A(j,j). +* + X( J ) = X( J ) / TJJS - SUMJ + END IF + XMAX = MAX( XMAX, ABS( X( J ) ) ) + 160 CONTINUE + END IF + SCALE = SCALE / TSCAL + END IF +* +* Scale the column norms by 1/TSCAL for return. +* + IF( TSCAL.NE.ONE ) THEN + CALL DSCAL( N, ONE / TSCAL, CNORM, 1 ) + END IF +* + RETURN +* +* End of DLATRS +* + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/lapack/dtrti2.f Tue Feb 18 20:08:20 2003 +0000 @@ -0,0 +1,147 @@ + SUBROUTINE DTRTI2( UPLO, DIAG, 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 +* February 29, 1992 +* +* .. Scalar Arguments .. + CHARACTER DIAG, UPLO + INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. + DOUBLE PRECISION A( LDA, * ) +* .. +* +* Purpose +* ======= +* +* DTRTI2 computes the inverse of a real upper or lower triangular +* matrix. +* +* This is the Level 2 BLAS version of the algorithm. +* +* Arguments +* ========= +* +* UPLO (input) CHARACTER*1 +* Specifies whether the matrix A is upper or lower triangular. +* = 'U': Upper triangular +* = 'L': Lower triangular +* +* DIAG (input) CHARACTER*1 +* Specifies whether or not the matrix A is unit triangular. +* = 'N': Non-unit triangular +* = 'U': Unit triangular +* +* N (input) INTEGER +* The order of the matrix A. N >= 0. +* +* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) +* On entry, the triangular matrix A. If UPLO = 'U', the +* leading n by n upper triangular part of the array A contains +* the upper triangular matrix, and the strictly lower +* triangular part of A is not referenced. If UPLO = 'L', the +* leading n by n lower triangular part of the array A contains +* the lower triangular matrix, and the strictly upper +* triangular part of A is not referenced. If DIAG = 'U', the +* diagonal elements of A are also not referenced and are +* assumed to be 1. +* +* On exit, the (triangular) inverse of the original matrix, in +* the same storage format. +* +* 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 NOUNIT, UPPER + INTEGER J + DOUBLE PRECISION AJJ +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. External Subroutines .. + EXTERNAL DSCAL, DTRMV, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX +* .. +* .. Executable Statements .. +* +* Test the input parameters. +* + INFO = 0 + UPPER = LSAME( UPLO, 'U' ) + NOUNIT = LSAME( DIAG, 'N' ) + IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN + INFO = -1 + ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN + INFO = -2 + ELSE IF( N.LT.0 ) THEN + INFO = -3 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -5 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DTRTI2', -INFO ) + RETURN + END IF +* + IF( UPPER ) THEN +* +* Compute inverse of upper triangular matrix. +* + DO 10 J = 1, N + IF( NOUNIT ) THEN + A( J, J ) = ONE / A( J, J ) + AJJ = -A( J, J ) + ELSE + AJJ = -ONE + END IF +* +* Compute elements 1:j-1 of j-th column. +* + CALL DTRMV( 'Upper', 'No transpose', DIAG, J-1, A, LDA, + $ A( 1, J ), 1 ) + CALL DSCAL( J-1, AJJ, A( 1, J ), 1 ) + 10 CONTINUE + ELSE +* +* Compute inverse of lower triangular matrix. +* + DO 20 J = N, 1, -1 + IF( NOUNIT ) THEN + A( J, J ) = ONE / A( J, J ) + AJJ = -A( J, J ) + ELSE + AJJ = -ONE + END IF + IF( J.LT.N ) THEN +* +* Compute elements j+1:n of j-th column. +* + CALL DTRMV( 'Lower', 'No transpose', DIAG, N-J, + $ A( J+1, J+1 ), LDA, A( J+1, J ), 1 ) + CALL DSCAL( N-J, AJJ, A( J+1, J ), 1 ) + END IF + 20 CONTINUE + END IF +* + RETURN +* +* End of DTRTI2 +* + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/lapack/dtrtri.f Tue Feb 18 20:08:20 2003 +0000 @@ -0,0 +1,177 @@ + SUBROUTINE DTRTRI( UPLO, DIAG, 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 DIAG, UPLO + INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. + DOUBLE PRECISION A( LDA, * ) +* .. +* +* Purpose +* ======= +* +* DTRTRI computes the inverse of a real upper or lower triangular +* matrix A. +* +* This is the Level 3 BLAS version of the algorithm. +* +* Arguments +* ========= +* +* UPLO (input) CHARACTER*1 +* = 'U': A is upper triangular; +* = 'L': A is lower triangular. +* +* DIAG (input) CHARACTER*1 +* = 'N': A is non-unit triangular; +* = 'U': A is unit triangular. +* +* N (input) INTEGER +* The order of the matrix A. N >= 0. +* +* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) +* On entry, the triangular matrix A. If UPLO = 'U', the +* leading N-by-N upper triangular part of the array A contains +* the upper triangular matrix, and the strictly lower +* triangular part of A is not referenced. If UPLO = 'L', the +* leading N-by-N lower triangular part of the array A contains +* the lower triangular matrix, and the strictly upper +* triangular part of A is not referenced. If DIAG = 'U', the +* diagonal elements of A are also not referenced and are +* assumed to be 1. +* On exit, the (triangular) inverse of the original matrix, in +* the same storage format. +* +* 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, A(i,i) is exactly zero. The triangular +* matrix is singular and its inverse can not be computed. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ONE, ZERO + PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) +* .. +* .. Local Scalars .. + LOGICAL NOUNIT, UPPER + INTEGER J, JB, NB, NN +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER ILAENV + EXTERNAL LSAME, ILAENV +* .. +* .. External Subroutines .. + EXTERNAL DTRMM, DTRSM, DTRTI2, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. Executable Statements .. +* +* Test the input parameters. +* + INFO = 0 + UPPER = LSAME( UPLO, 'U' ) + NOUNIT = LSAME( DIAG, 'N' ) + IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN + INFO = -1 + ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN + INFO = -2 + ELSE IF( N.LT.0 ) THEN + INFO = -3 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -5 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DTRTRI', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 ) + $ RETURN +* +* Check for singularity if non-unit. +* + IF( NOUNIT ) THEN + DO 10 INFO = 1, N + IF( A( INFO, INFO ).EQ.ZERO ) + $ RETURN + 10 CONTINUE + INFO = 0 + END IF +* +* Determine the block size for this environment. +* + NB = ILAENV( 1, 'DTRTRI', UPLO // DIAG, N, -1, -1, -1 ) + IF( NB.LE.1 .OR. NB.GE.N ) THEN +* +* Use unblocked code +* + CALL DTRTI2( UPLO, DIAG, N, A, LDA, INFO ) + ELSE +* +* Use blocked code +* + IF( UPPER ) THEN +* +* Compute inverse of upper triangular matrix +* + DO 20 J = 1, N, NB + JB = MIN( NB, N-J+1 ) +* +* Compute rows 1:j-1 of current block column +* + CALL DTRMM( 'Left', 'Upper', 'No transpose', DIAG, J-1, + $ JB, ONE, A, LDA, A( 1, J ), LDA ) + CALL DTRSM( 'Right', 'Upper', 'No transpose', DIAG, J-1, + $ JB, -ONE, A( J, J ), LDA, A( 1, J ), LDA ) +* +* Compute inverse of current diagonal block +* + CALL DTRTI2( 'Upper', DIAG, JB, A( J, J ), LDA, INFO ) + 20 CONTINUE + ELSE +* +* Compute inverse of lower triangular matrix +* + NN = ( ( N-1 ) / NB )*NB + 1 + DO 30 J = NN, 1, -NB + JB = MIN( NB, N-J+1 ) + IF( J+JB.LE.N ) THEN +* +* Compute rows j+jb:n of current block column +* + CALL DTRMM( 'Left', 'Lower', 'No transpose', DIAG, + $ N-J-JB+1, JB, ONE, A( J+JB, J+JB ), LDA, + $ A( J+JB, J ), LDA ) + CALL DTRSM( 'Right', 'Lower', 'No transpose', DIAG, + $ N-J-JB+1, JB, -ONE, A( J, J ), LDA, + $ A( J+JB, J ), LDA ) + END IF +* +* Compute inverse of current diagonal block +* + CALL DTRTI2( 'Lower', DIAG, JB, A( J, J ), LDA, INFO ) + 30 CONTINUE + END IF + END IF +* + RETURN +* +* End of DTRTRI +* + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/lapack/spotf2.f Tue Feb 18 20:08:20 2003 +0000 @@ -0,0 +1,168 @@ + SUBROUTINE SPOTF2( 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 +* February 29, 1992 +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. + REAL A( LDA, * ) +* .. +* +* Purpose +* ======= +* +* SPOTF2 computes the Cholesky factorization of a real symmetric +* positive definite matrix A. +* +* The factorization has the form +* A = U' * U , if UPLO = 'U', or +* A = L * L', if UPLO = 'L', +* where U is an upper triangular matrix and L is lower triangular. +* +* This is the unblocked version of the algorithm, calling Level 2 BLAS. +* +* Arguments +* ========= +* +* UPLO (input) CHARACTER*1 +* Specifies whether the upper or lower triangular part of the +* symmetric matrix A is stored. +* = 'U': Upper triangular +* = 'L': Lower triangular +* +* N (input) INTEGER +* The order of the matrix A. N >= 0. +* +* A (input/output) 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'*U or A = 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 +* > 0: if INFO = k, the leading minor of order k is not +* positive definite, and the factorization could not be +* completed. +* +* ===================================================================== +* +* .. Parameters .. + REAL ONE, ZERO + PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) +* .. +* .. Local Scalars .. + LOGICAL UPPER + INTEGER J + REAL AJJ +* .. +* .. External Functions .. + LOGICAL LSAME + REAL SDOT + EXTERNAL LSAME, SDOT +* .. +* .. External Subroutines .. + EXTERNAL SGEMV, SSCAL, 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( 'SPOTF2', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 ) + $ RETURN +* + IF( UPPER ) THEN +* +* Compute the Cholesky factorization A = U'*U. +* + DO 10 J = 1, N +* +* Compute U(J,J) and test for non-positive-definiteness. +* + AJJ = A( J, J ) - SDOT( J-1, A( 1, J ), 1, A( 1, J ), 1 ) + IF( AJJ.LE.ZERO ) THEN + A( J, J ) = AJJ + GO TO 30 + END IF + AJJ = SQRT( AJJ ) + A( J, J ) = AJJ +* +* Compute elements J+1:N of row J. +* + IF( J.LT.N ) THEN + CALL SGEMV( 'Transpose', J-1, N-J, -ONE, A( 1, J+1 ), + $ LDA, A( 1, J ), 1, ONE, A( J, J+1 ), LDA ) + CALL SSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA ) + END IF + 10 CONTINUE + ELSE +* +* Compute the Cholesky factorization A = L*L'. +* + DO 20 J = 1, N +* +* Compute L(J,J) and test for non-positive-definiteness. +* + AJJ = A( J, J ) - SDOT( J-1, A( J, 1 ), LDA, A( J, 1 ), + $ LDA ) + IF( AJJ.LE.ZERO ) THEN + A( J, J ) = AJJ + GO TO 30 + END IF + AJJ = SQRT( AJJ ) + A( J, J ) = AJJ +* +* Compute elements J+1:N of column J. +* + IF( J.LT.N ) THEN + CALL SGEMV( 'No transpose', N-J, J-1, -ONE, A( J+1, 1 ), + $ LDA, A( J, 1 ), LDA, ONE, A( J+1, J ), 1 ) + CALL SSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 ) + END IF + 20 CONTINUE + END IF + GO TO 40 +* + 30 CONTINUE + INFO = J +* + 40 CONTINUE + RETURN +* +* End of SPOTF2 +* + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/lapack/spotrf.f Tue Feb 18 20:08:20 2003 +0000 @@ -0,0 +1,184 @@ + SUBROUTINE SPOTRF( 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 .. + REAL A( LDA, * ) +* .. +* +* Purpose +* ======= +* +* SPOTRF computes the Cholesky factorization of a real symmetric +* positive definite matrix A. +* +* 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 block version of the algorithm, calling Level 3 BLAS. +* +* 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) 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. +* +* 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 leading minor of order i is not +* positive definite, and the factorization could not be +* completed. +* +* ===================================================================== +* +* .. Parameters .. + REAL ONE + PARAMETER ( ONE = 1.0E+0 ) +* .. +* .. Local Scalars .. + LOGICAL UPPER + INTEGER J, JB, NB +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER ILAENV + EXTERNAL LSAME, ILAENV +* .. +* .. External Subroutines .. + EXTERNAL SGEMM, SPOTF2, SSYRK, STRSM, 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( 'SPOTRF', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 ) + $ RETURN +* +* Determine the block size for this environment. +* + NB = ILAENV( 1, 'SPOTRF', UPLO, N, -1, -1, -1 ) + IF( NB.LE.1 .OR. NB.GE.N ) THEN +* +* Use unblocked code. +* + CALL SPOTF2( UPLO, N, A, LDA, INFO ) + ELSE +* +* Use blocked code. +* + IF( UPPER ) THEN +* +* Compute the Cholesky factorization A = U'*U. +* + DO 10 J = 1, N, NB +* +* Update and factorize the current diagonal block and test +* for non-positive-definiteness. +* + 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 ) + IF( INFO.NE.0 ) + $ GO TO 30 + IF( J+JB.LE.N ) THEN +* +* Compute the current block row. +* + CALL SGEMM( 'Transpose', 'No transpose', JB, N-J-JB+1, + $ J-1, -ONE, A( 1, J ), LDA, A( 1, J+JB ), + $ LDA, ONE, A( J, J+JB ), LDA ) + CALL STRSM( 'Left', 'Upper', 'Transpose', 'Non-unit', + $ JB, N-J-JB+1, ONE, A( J, J ), LDA, + $ A( J, J+JB ), LDA ) + END IF + 10 CONTINUE +* + ELSE +* +* Compute the Cholesky factorization A = L*L'. +* + DO 20 J = 1, N, NB +* +* Update and factorize the current diagonal block and test +* for non-positive-definiteness. +* + 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 ) + IF( INFO.NE.0 ) + $ GO TO 30 + IF( J+JB.LE.N ) THEN +* +* Compute the current block column. +* + CALL SGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB, + $ J-1, -ONE, A( J+JB, 1 ), LDA, A( J, 1 ), + $ LDA, ONE, A( J+JB, J ), LDA ) + CALL STRSM( 'Right', 'Lower', 'Transpose', 'Non-unit', + $ N-J-JB+1, JB, ONE, A( J, J ), LDA, + $ A( J+JB, J ), LDA ) + END IF + 20 CONTINUE + END IF + END IF + GO TO 40 +* + 30 CONTINUE + INFO = INFO + J - 1 +* + 40 CONTINUE + RETURN +* +* End of SPOTRF +* + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/lapack/zgecon.f Tue Feb 18 20:08:20 2003 +0000 @@ -0,0 +1,189 @@ + SUBROUTINE ZGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, RWORK, + $ 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 NORM + INTEGER INFO, LDA, N + DOUBLE PRECISION ANORM, RCOND +* .. +* .. Array Arguments .. + DOUBLE PRECISION RWORK( * ) + COMPLEX*16 A( LDA, * ), WORK( * ) +* .. +* +* Purpose +* ======= +* +* ZGECON estimates the reciprocal of the condition number of a general +* complex matrix A, in either the 1-norm or the infinity-norm, using +* the LU factorization computed by ZGETRF. +* +* An estimate is obtained for norm(inv(A)), and the reciprocal of the +* condition number is computed as +* RCOND = 1 / ( norm(A) * norm(inv(A)) ). +* +* Arguments +* ========= +* +* NORM (input) CHARACTER*1 +* Specifies whether the 1-norm condition number or the +* infinity-norm condition number is required: +* = '1' or 'O': 1-norm; +* = 'I': Infinity-norm. +* +* N (input) INTEGER +* The order of the matrix A. N >= 0. +* +* A (input) COMPLEX*16 array, dimension (LDA,N) +* The factors L and U from the factorization A = P*L*U +* as computed by ZGETRF. +* +* LDA (input) INTEGER +* The leading dimension of the array A. LDA >= max(1,N). +* +* ANORM (input) DOUBLE PRECISION +* If NORM = '1' or 'O', the 1-norm of the original matrix A. +* If NORM = 'I', the infinity-norm of the original matrix A. +* +* RCOND (output) DOUBLE PRECISION +* The reciprocal of the condition number of the matrix A, +* computed as RCOND = 1/(norm(A) * norm(inv(A))). +* +* WORK (workspace) COMPLEX*16 array, dimension (2*N) +* +* RWORK (workspace) DOUBLE PRECISION array, dimension (2*N) +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument had an illegal value +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ONE, ZERO + PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) +* .. +* .. Local Scalars .. + LOGICAL ONENRM + CHARACTER NORMIN + INTEGER IX, KASE, KASE1 + DOUBLE PRECISION AINVNM, SCALE, SL, SMLNUM, SU + COMPLEX*16 ZDUM +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER IZAMAX + DOUBLE PRECISION DLAMCH + EXTERNAL LSAME, IZAMAX, DLAMCH +* .. +* .. External Subroutines .. + EXTERNAL XERBLA, ZDRSCL, ZLACON, ZLATRS +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, DBLE, DIMAG, MAX +* .. +* .. Statement Functions .. + DOUBLE PRECISION CABS1 +* .. +* .. Statement Function definitions .. + CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) ) +* .. +* .. Executable Statements .. +* +* Test the input parameters. +* + INFO = 0 + ONENRM = NORM.EQ.'1' .OR. LSAME( NORM, 'O' ) + IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM, 'I' ) ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -4 + ELSE IF( ANORM.LT.ZERO ) THEN + INFO = -5 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'ZGECON', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + RCOND = ZERO + IF( N.EQ.0 ) THEN + RCOND = ONE + RETURN + ELSE IF( ANORM.EQ.ZERO ) THEN + RETURN + END IF +* + SMLNUM = DLAMCH( 'Safe minimum' ) +* +* Estimate the norm of inv(A). +* + AINVNM = ZERO + NORMIN = 'N' + IF( ONENRM ) THEN + KASE1 = 1 + ELSE + KASE1 = 2 + END IF + KASE = 0 + 10 CONTINUE + CALL ZLACON( N, WORK( N+1 ), WORK, AINVNM, KASE ) + IF( KASE.NE.0 ) THEN + IF( KASE.EQ.KASE1 ) THEN +* +* Multiply by inv(L). +* + CALL ZLATRS( 'Lower', 'No transpose', 'Unit', NORMIN, N, A, + $ LDA, WORK, SL, RWORK, INFO ) +* +* Multiply by inv(U). +* + CALL ZLATRS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N, + $ A, LDA, WORK, SU, RWORK( N+1 ), INFO ) + ELSE +* +* Multiply by inv(U'). +* + CALL ZLATRS( 'Upper', 'Conjugate transpose', 'Non-unit', + $ NORMIN, N, A, LDA, WORK, SU, RWORK( N+1 ), + $ INFO ) +* +* Multiply by inv(L'). +* + CALL ZLATRS( 'Lower', 'Conjugate transpose', 'Unit', NORMIN, + $ N, A, LDA, WORK, SL, RWORK, INFO ) + END IF +* +* Divide X by 1/(SL*SU) if doing so will not cause overflow. +* + SCALE = SL*SU + NORMIN = 'Y' + IF( SCALE.NE.ONE ) THEN + IX = IZAMAX( N, WORK, 1 ) + IF( SCALE.LT.CABS1( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO ) + $ GO TO 20 + CALL ZDRSCL( N, SCALE, WORK, 1 ) + END IF + GO TO 10 + END IF +* +* Compute the estimate of the reciprocal condition number. +* + IF( AINVNM.NE.ZERO ) + $ RCOND = ( ONE / AINVNM ) / ANORM +* + 20 CONTINUE + RETURN +* +* End of ZGECON +* + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/lapack/zgetri.f Tue Feb 18 20:08:20 2003 +0000 @@ -0,0 +1,194 @@ + SUBROUTINE ZGETRI( N, A, LDA, IPIV, WORK, LWORK, INFO ) +* +* -- LAPACK routine (version 3.0) -- +* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., +* Courant Institute, Argonne National Lab, and Rice University +* June 30, 1999 +* +* .. Scalar Arguments .. + INTEGER INFO, LDA, LWORK, N +* .. +* .. Array Arguments .. + INTEGER IPIV( * ) + COMPLEX*16 A( LDA, * ), WORK( * ) +* .. +* +* Purpose +* ======= +* +* ZGETRI computes the inverse of a matrix using the LU factorization +* computed by ZGETRF. +* +* This method inverts U and then computes inv(A) by solving the system +* inv(A)*L = inv(U) for inv(A). +* +* Arguments +* ========= +* +* N (input) INTEGER +* The order of the matrix A. N >= 0. +* +* A (input/output) COMPLEX*16 array, dimension (LDA,N) +* On entry, the factors L and U from the factorization +* A = P*L*U as computed by ZGETRF. +* On exit, if INFO = 0, the inverse of the original matrix A. +* +* LDA (input) INTEGER +* The leading dimension of the array A. LDA >= max(1,N). +* +* IPIV (input) INTEGER array, dimension (N) +* The pivot indices from ZGETRF; for 1<=i<=N, row i of the +* matrix was interchanged with row IPIV(i). +* +* WORK (workspace/output) COMPLEX*16 array, dimension (LWORK) +* On exit, if INFO=0, then WORK(1) returns the optimal LWORK. +* +* LWORK (input) INTEGER +* The dimension of the array WORK. LWORK >= max(1,N). +* For optimal performance LWORK >= N*NB, where NB is +* the optimal blocksize returned by ILAENV. +* +* 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. +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument had an illegal value +* > 0: if INFO = i, U(i,i) is exactly zero; the matrix is +* singular and its inverse could not be computed. +* +* ===================================================================== +* +* .. Parameters .. + COMPLEX*16 ZERO, ONE + PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ), + $ ONE = ( 1.0D+0, 0.0D+0 ) ) +* .. +* .. Local Scalars .. + LOGICAL LQUERY + INTEGER I, IWS, J, JB, JJ, JP, LDWORK, LWKOPT, NB, + $ NBMIN, NN +* .. +* .. External Functions .. + INTEGER ILAENV + EXTERNAL ILAENV +* .. +* .. External Subroutines .. + EXTERNAL XERBLA, ZGEMM, ZGEMV, ZSWAP, ZTRSM, ZTRTRI +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. Executable Statements .. +* +* Test the input parameters. +* + INFO = 0 + NB = ILAENV( 1, 'ZGETRI', ' ', N, -1, -1, -1 ) + LWKOPT = N*NB + WORK( 1 ) = LWKOPT + LQUERY = ( LWORK.EQ.-1 ) + IF( N.LT.0 ) THEN + INFO = -1 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -3 + ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN + INFO = -6 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'ZGETRI', -INFO ) + RETURN + ELSE IF( LQUERY ) THEN + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 ) + $ RETURN +* +* Form inv(U). If INFO > 0 from ZTRTRI, then U is singular, +* and the inverse is not computed. +* + CALL ZTRTRI( 'Upper', 'Non-unit', N, A, LDA, INFO ) + IF( INFO.GT.0 ) + $ RETURN +* + NBMIN = 2 + LDWORK = N + IF( NB.GT.1 .AND. NB.LT.N ) THEN + IWS = MAX( LDWORK*NB, 1 ) + IF( LWORK.LT.IWS ) THEN + NB = LWORK / LDWORK + NBMIN = MAX( 2, ILAENV( 2, 'ZGETRI', ' ', N, -1, -1, -1 ) ) + END IF + ELSE + IWS = N + END IF +* +* Solve the equation inv(A)*L = inv(U) for inv(A). +* + IF( NB.LT.NBMIN .OR. NB.GE.N ) THEN +* +* Use unblocked code. +* + DO 20 J = N, 1, -1 +* +* Copy current column of L to WORK and replace with zeros. +* + DO 10 I = J + 1, N + WORK( I ) = A( I, J ) + A( I, J ) = ZERO + 10 CONTINUE +* +* Compute current column of inv(A). +* + IF( J.LT.N ) + $ CALL ZGEMV( 'No transpose', N, N-J, -ONE, A( 1, J+1 ), + $ LDA, WORK( J+1 ), 1, ONE, A( 1, J ), 1 ) + 20 CONTINUE + ELSE +* +* Use blocked code. +* + NN = ( ( N-1 ) / NB )*NB + 1 + DO 50 J = NN, 1, -NB + JB = MIN( NB, N-J+1 ) +* +* Copy current block column of L to WORK and replace with +* zeros. +* + DO 40 JJ = J, J + JB - 1 + DO 30 I = JJ + 1, N + WORK( I+( JJ-J )*LDWORK ) = A( I, JJ ) + A( I, JJ ) = ZERO + 30 CONTINUE + 40 CONTINUE +* +* Compute current block column of inv(A). +* + IF( J+JB.LE.N ) + $ CALL ZGEMM( 'No transpose', 'No transpose', N, JB, + $ N-J-JB+1, -ONE, A( 1, J+JB ), LDA, + $ WORK( J+JB ), LDWORK, ONE, A( 1, J ), LDA ) + CALL ZTRSM( 'Right', 'Lower', 'No transpose', 'Unit', N, JB, + $ ONE, WORK( J ), LDWORK, A( 1, J ), LDA ) + 50 CONTINUE + END IF +* +* Apply column interchanges. +* + DO 60 J = N - 1, 1, -1 + JP = IPIV( J ) + IF( JP.NE.J ) + $ CALL ZSWAP( N, A( 1, J ), 1, A( 1, JP ), 1 ) + 60 CONTINUE +* + WORK( 1 ) = IWS + RETURN +* +* End of ZGETRI +* + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/lapack/ztrti2.f Tue Feb 18 20:08:20 2003 +0000 @@ -0,0 +1,147 @@ + SUBROUTINE ZTRTI2( UPLO, DIAG, 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 +* September 30, 1994 +* +* .. Scalar Arguments .. + CHARACTER DIAG, UPLO + INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. + COMPLEX*16 A( LDA, * ) +* .. +* +* Purpose +* ======= +* +* ZTRTI2 computes the inverse of a complex upper or lower triangular +* matrix. +* +* This is the Level 2 BLAS version of the algorithm. +* +* Arguments +* ========= +* +* UPLO (input) CHARACTER*1 +* Specifies whether the matrix A is upper or lower triangular. +* = 'U': Upper triangular +* = 'L': Lower triangular +* +* DIAG (input) CHARACTER*1 +* Specifies whether or not the matrix A is unit triangular. +* = 'N': Non-unit triangular +* = 'U': Unit triangular +* +* N (input) INTEGER +* The order of the matrix A. N >= 0. +* +* A (input/output) COMPLEX*16 array, dimension (LDA,N) +* On entry, the triangular matrix A. If UPLO = 'U', the +* leading n by n upper triangular part of the array A contains +* the upper triangular matrix, and the strictly lower +* triangular part of A is not referenced. If UPLO = 'L', the +* leading n by n lower triangular part of the array A contains +* the lower triangular matrix, and the strictly upper +* triangular part of A is not referenced. If DIAG = 'U', the +* diagonal elements of A are also not referenced and are +* assumed to be 1. +* +* On exit, the (triangular) inverse of the original matrix, in +* the same storage format. +* +* 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 NOUNIT, UPPER + INTEGER J + COMPLEX*16 AJJ +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. External Subroutines .. + EXTERNAL XERBLA, ZSCAL, ZTRMV +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX +* .. +* .. Executable Statements .. +* +* Test the input parameters. +* + INFO = 0 + UPPER = LSAME( UPLO, 'U' ) + NOUNIT = LSAME( DIAG, 'N' ) + IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN + INFO = -1 + ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN + INFO = -2 + ELSE IF( N.LT.0 ) THEN + INFO = -3 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -5 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'ZTRTI2', -INFO ) + RETURN + END IF +* + IF( UPPER ) THEN +* +* Compute inverse of upper triangular matrix. +* + DO 10 J = 1, N + IF( NOUNIT ) THEN + A( J, J ) = ONE / A( J, J ) + AJJ = -A( J, J ) + ELSE + AJJ = -ONE + END IF +* +* Compute elements 1:j-1 of j-th column. +* + CALL ZTRMV( 'Upper', 'No transpose', DIAG, J-1, A, LDA, + $ A( 1, J ), 1 ) + CALL ZSCAL( J-1, AJJ, A( 1, J ), 1 ) + 10 CONTINUE + ELSE +* +* Compute inverse of lower triangular matrix. +* + DO 20 J = N, 1, -1 + IF( NOUNIT ) THEN + A( J, J ) = ONE / A( J, J ) + AJJ = -A( J, J ) + ELSE + AJJ = -ONE + END IF + IF( J.LT.N ) THEN +* +* Compute elements j+1:n of j-th column. +* + CALL ZTRMV( 'Lower', 'No transpose', DIAG, N-J, + $ A( J+1, J+1 ), LDA, A( J+1, J ), 1 ) + CALL ZSCAL( N-J, AJJ, A( J+1, J ), 1 ) + END IF + 20 CONTINUE + END IF +* + RETURN +* +* End of ZTRTI2 +* + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/lapack/ztrtri.f Tue Feb 18 20:08:20 2003 +0000 @@ -0,0 +1,178 @@ + SUBROUTINE ZTRTRI( UPLO, DIAG, 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 +* September 30, 1994 +* +* .. Scalar Arguments .. + CHARACTER DIAG, UPLO + INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. + COMPLEX*16 A( LDA, * ) +* .. +* +* Purpose +* ======= +* +* ZTRTRI computes the inverse of a complex upper or lower triangular +* matrix A. +* +* This is the Level 3 BLAS version of the algorithm. +* +* Arguments +* ========= +* +* UPLO (input) CHARACTER*1 +* = 'U': A is upper triangular; +* = 'L': A is lower triangular. +* +* DIAG (input) CHARACTER*1 +* = 'N': A is non-unit triangular; +* = 'U': A is unit triangular. +* +* N (input) INTEGER +* The order of the matrix A. N >= 0. +* +* A (input/output) COMPLEX*16 array, dimension (LDA,N) +* On entry, the triangular matrix A. If UPLO = 'U', the +* leading N-by-N upper triangular part of the array A contains +* the upper triangular matrix, and the strictly lower +* triangular part of A is not referenced. If UPLO = 'L', the +* leading N-by-N lower triangular part of the array A contains +* the lower triangular matrix, and the strictly upper +* triangular part of A is not referenced. If DIAG = 'U', the +* diagonal elements of A are also not referenced and are +* assumed to be 1. +* On exit, the (triangular) inverse of the original matrix, in +* the same storage format. +* +* 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, A(i,i) is exactly zero. The triangular +* matrix is singular and its inverse can not be computed. +* +* ===================================================================== +* +* .. Parameters .. + COMPLEX*16 ONE, ZERO + PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ), + $ ZERO = ( 0.0D+0, 0.0D+0 ) ) +* .. +* .. Local Scalars .. + LOGICAL NOUNIT, UPPER + INTEGER J, JB, NB, NN +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER ILAENV + EXTERNAL LSAME, ILAENV +* .. +* .. External Subroutines .. + EXTERNAL XERBLA, ZTRMM, ZTRSM, ZTRTI2 +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. Executable Statements .. +* +* Test the input parameters. +* + INFO = 0 + UPPER = LSAME( UPLO, 'U' ) + NOUNIT = LSAME( DIAG, 'N' ) + IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN + INFO = -1 + ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN + INFO = -2 + ELSE IF( N.LT.0 ) THEN + INFO = -3 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -5 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'ZTRTRI', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 ) + $ RETURN +* +* Check for singularity if non-unit. +* + IF( NOUNIT ) THEN + DO 10 INFO = 1, N + IF( A( INFO, INFO ).EQ.ZERO ) + $ RETURN + 10 CONTINUE + INFO = 0 + END IF +* +* Determine the block size for this environment. +* + NB = ILAENV( 1, 'ZTRTRI', UPLO // DIAG, N, -1, -1, -1 ) + IF( NB.LE.1 .OR. NB.GE.N ) THEN +* +* Use unblocked code +* + CALL ZTRTI2( UPLO, DIAG, N, A, LDA, INFO ) + ELSE +* +* Use blocked code +* + IF( UPPER ) THEN +* +* Compute inverse of upper triangular matrix +* + DO 20 J = 1, N, NB + JB = MIN( NB, N-J+1 ) +* +* Compute rows 1:j-1 of current block column +* + CALL ZTRMM( 'Left', 'Upper', 'No transpose', DIAG, J-1, + $ JB, ONE, A, LDA, A( 1, J ), LDA ) + CALL ZTRSM( 'Right', 'Upper', 'No transpose', DIAG, J-1, + $ JB, -ONE, A( J, J ), LDA, A( 1, J ), LDA ) +* +* Compute inverse of current diagonal block +* + CALL ZTRTI2( 'Upper', DIAG, JB, A( J, J ), LDA, INFO ) + 20 CONTINUE + ELSE +* +* Compute inverse of lower triangular matrix +* + NN = ( ( N-1 ) / NB )*NB + 1 + DO 30 J = NN, 1, -NB + JB = MIN( NB, N-J+1 ) + IF( J+JB.LE.N ) THEN +* +* Compute rows j+jb:n of current block column +* + CALL ZTRMM( 'Left', 'Lower', 'No transpose', DIAG, + $ N-J-JB+1, JB, ONE, A( J+JB, J+JB ), LDA, + $ A( J+JB, J ), LDA ) + CALL ZTRSM( 'Right', 'Lower', 'No transpose', DIAG, + $ N-J-JB+1, JB, -ONE, A( J, J ), LDA, + $ A( J+JB, J ), LDA ) + END IF +* +* Compute inverse of current diagonal block +* + CALL ZTRTI2( 'Lower', DIAG, JB, A( J, J ), LDA, INFO ) + 30 CONTINUE + END IF + END IF +* + RETURN +* +* End of ZTRTRI +* + END
--- a/libcruft/linpack/Makefile.in Sun Feb 16 04:29:00 2003 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,19 +0,0 @@ -# -# Makefile for octave's libcruft/linpack directory -# -# John W. Eaton -# jwe@bevo.che.wisc.edu -# University of Wisconsin-Madison -# Department of Chemical Engineering - -TOPDIR = ../.. - -srcdir = @srcdir@ -top_srcdir = @top_srcdir@ -VPATH = @srcdir@ - -EXTERNAL_DISTFILES = $(DISTFILES) - -include $(TOPDIR)/Makeconf - -include ../Makerules
--- a/libcruft/linpack/dgbfa.f Sun Feb 16 04:29:00 2003 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,174 +0,0 @@ - SUBROUTINE DGBFA(ABD,LDA,N,ML,MU,IPVT,INFO) - INTEGER LDA,N,ML,MU,IPVT(1),INFO - DOUBLE PRECISION ABD(LDA,1) -C -C DGBFA FACTORS A DOUBLE PRECISION BAND MATRIX BY ELIMINATION. -C -C DGBFA IS USUALLY CALLED BY DGBCO, BUT IT CAN BE CALLED -C DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED. -C -C ON ENTRY -C -C ABD DOUBLE PRECISION(LDA, N) -C CONTAINS THE MATRIX IN BAND STORAGE. THE COLUMNS -C OF THE MATRIX ARE STORED IN THE COLUMNS OF ABD AND -C THE DIAGONALS OF THE MATRIX ARE STORED IN ROWS -C ML+1 THROUGH 2*ML+MU+1 OF ABD . -C SEE THE COMMENTS BELOW FOR DETAILS. -C -C LDA INTEGER -C THE LEADING DIMENSION OF THE ARRAY ABD . -C LDA MUST BE .GE. 2*ML + MU + 1 . -C -C N INTEGER -C THE ORDER OF THE ORIGINAL MATRIX. -C -C ML INTEGER -C NUMBER OF DIAGONALS BELOW THE MAIN DIAGONAL. -C 0 .LE. ML .LT. N . -C -C MU INTEGER -C NUMBER OF DIAGONALS ABOVE THE MAIN DIAGONAL. -C 0 .LE. MU .LT. N . -C MORE EFFICIENT IF ML .LE. MU . -C ON RETURN -C -C ABD AN UPPER TRIANGULAR MATRIX IN BAND STORAGE AND -C THE MULTIPLIERS WHICH WERE USED TO OBTAIN IT. -C THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE -C L IS A PRODUCT OF PERMUTATION AND UNIT LOWER -C TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR. -C -C IPVT INTEGER(N) -C AN INTEGER VECTOR OF PIVOT INDICES. -C -C INFO INTEGER -C = 0 NORMAL VALUE. -C = K IF U(K,K) .EQ. 0.0 . THIS IS NOT AN ERROR -C CONDITION FOR THIS SUBROUTINE, BUT IT DOES -C INDICATE THAT DGBSL WILL DIVIDE BY ZERO IF -C CALLED. USE RCOND IN DGBCO FOR A RELIABLE -C INDICATION OF SINGULARITY. -C -C BAND STORAGE -C -C IF A IS A BAND MATRIX, THE FOLLOWING PROGRAM SEGMENT -C WILL SET UP THE INPUT. -C -C ML = (BAND WIDTH BELOW THE DIAGONAL) -C MU = (BAND WIDTH ABOVE THE DIAGONAL) -C M = ML + MU + 1 -C DO 20 J = 1, N -C I1 = MAX0(1, J-MU) -C I2 = MIN0(N, J+ML) -C DO 10 I = I1, I2 -C K = I - J + M -C ABD(K,J) = A(I,J) -C 10 CONTINUE -C 20 CONTINUE -C -C THIS USES ROWS ML+1 THROUGH 2*ML+MU+1 OF ABD . -C IN ADDITION, THE FIRST ML ROWS IN ABD ARE USED FOR -C ELEMENTS GENERATED DURING THE TRIANGULARIZATION. -C THE TOTAL NUMBER OF ROWS NEEDED IN ABD IS 2*ML+MU+1 . -C THE ML+MU BY ML+MU UPPER LEFT TRIANGLE AND THE -C ML BY ML LOWER RIGHT TRIANGLE ARE NOT REFERENCED. -C -C LINPACK. THIS VERSION DATED 08/14/78 . -C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. -C -C SUBROUTINES AND FUNCTIONS -C -C BLAS DAXPY,DSCAL,IDAMAX -C FORTRAN MAX0,MIN0 -C -C INTERNAL VARIABLES -C - DOUBLE PRECISION T - INTEGER I,IDAMAX,I0,J,JU,JZ,J0,J1,K,KP1,L,LM,M,MM,NM1 -C -C - M = ML + MU + 1 - INFO = 0 -C -C ZERO INITIAL FILL-IN COLUMNS -C - J0 = MU + 2 - J1 = MIN0(N,M) - 1 - IF (J1 .LT. J0) GO TO 30 - DO 20 JZ = J0, J1 - I0 = M + 1 - JZ - DO 10 I = I0, ML - ABD(I,JZ) = 0.0D0 - 10 CONTINUE - 20 CONTINUE - 30 CONTINUE - JZ = J1 - JU = 0 -C -C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING -C - NM1 = N - 1 - IF (NM1 .LT. 1) GO TO 130 - DO 120 K = 1, NM1 - KP1 = K + 1 -C -C ZERO NEXT FILL-IN COLUMN -C - JZ = JZ + 1 - IF (JZ .GT. N) GO TO 50 - IF (ML .LT. 1) GO TO 50 - DO 40 I = 1, ML - ABD(I,JZ) = 0.0D0 - 40 CONTINUE - 50 CONTINUE -C -C FIND L = PIVOT INDEX -C - LM = MIN0(ML,N-K) - L = IDAMAX(LM+1,ABD(M,K),1) + M - 1 - IPVT(K) = L + K - M -C -C ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED -C - IF (ABD(L,K) .EQ. 0.0D0) GO TO 100 -C -C INTERCHANGE IF NECESSARY -C - IF (L .EQ. M) GO TO 60 - T = ABD(L,K) - ABD(L,K) = ABD(M,K) - ABD(M,K) = T - 60 CONTINUE -C -C COMPUTE MULTIPLIERS -C - T = -1.0D0/ABD(M,K) - CALL DSCAL(LM,T,ABD(M+1,K),1) -C -C ROW ELIMINATION WITH COLUMN INDEXING -C - JU = MIN0(MAX0(JU,MU+IPVT(K)),N) - MM = M - IF (JU .LT. KP1) GO TO 90 - DO 80 J = KP1, JU - L = L - 1 - MM = MM - 1 - T = ABD(L,J) - IF (L .EQ. MM) GO TO 70 - ABD(L,J) = ABD(MM,J) - ABD(MM,J) = T - 70 CONTINUE - CALL DAXPY(LM,T,ABD(M+1,K),1,ABD(MM+1,J),1) - 80 CONTINUE - 90 CONTINUE - GO TO 110 - 100 CONTINUE - INFO = K - 110 CONTINUE - 120 CONTINUE - 130 CONTINUE - IPVT(N) = N - IF (ABD(M,N) .EQ. 0.0D0) INFO = N - RETURN - END
--- a/libcruft/linpack/dgbsl.f Sun Feb 16 04:29:00 2003 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,135 +0,0 @@ - SUBROUTINE DGBSL(ABD,LDA,N,ML,MU,IPVT,B,JOB) - INTEGER LDA,N,ML,MU,IPVT(1),JOB - DOUBLE PRECISION ABD(LDA,1),B(1) -C -C DGBSL SOLVES THE DOUBLE PRECISION BAND SYSTEM -C A * X = B OR TRANS(A) * X = B -C USING THE FACTORS COMPUTED BY DGBCO OR DGBFA. -C -C ON ENTRY -C -C ABD DOUBLE PRECISION(LDA, N) -C THE OUTPUT FROM DGBCO OR DGBFA. -C -C LDA INTEGER -C THE LEADING DIMENSION OF THE ARRAY ABD . -C -C N INTEGER -C THE ORDER OF THE ORIGINAL MATRIX. -C -C ML INTEGER -C NUMBER OF DIAGONALS BELOW THE MAIN DIAGONAL. -C -C MU INTEGER -C NUMBER OF DIAGONALS ABOVE THE MAIN DIAGONAL. -C -C IPVT INTEGER(N) -C THE PIVOT VECTOR FROM DGBCO OR DGBFA. -C -C B DOUBLE PRECISION(N) -C THE RIGHT HAND SIDE VECTOR. -C -C JOB INTEGER -C = 0 TO SOLVE A*X = B , -C = NONZERO TO SOLVE TRANS(A)*X = B , WHERE -C TRANS(A) IS THE TRANSPOSE. -C -C ON RETURN -C -C B THE SOLUTION VECTOR X . -C -C ERROR CONDITION -C -C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS A -C ZERO ON THE DIAGONAL. TECHNICALLY THIS INDICATES SINGULARITY -C BUT IT IS OFTEN CAUSED BY IMPROPER ARGUMENTS OR IMPROPER -C SETTING OF LDA . IT WILL NOT OCCUR IF THE SUBROUTINES ARE -C CALLED CORRECTLY AND IF DGBCO HAS SET RCOND .GT. 0.0 -C OR DGBFA HAS SET INFO .EQ. 0 . -C -C TO COMPUTE INVERSE(A) * C WHERE C IS A MATRIX -C WITH P COLUMNS -C CALL DGBCO(ABD,LDA,N,ML,MU,IPVT,RCOND,Z) -C IF (RCOND IS TOO SMALL) GO TO ... -C DO 10 J = 1, P -C CALL DGBSL(ABD,LDA,N,ML,MU,IPVT,C(1,J),0) -C 10 CONTINUE -C -C LINPACK. THIS VERSION DATED 08/14/78 . -C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. -C -C SUBROUTINES AND FUNCTIONS -C -C BLAS DAXPY,DDOT -C FORTRAN MIN0 -C -C INTERNAL VARIABLES -C - DOUBLE PRECISION DDOT,T - INTEGER K,KB,L,LA,LB,LM,M,NM1 -C - M = MU + ML + 1 - NM1 = N - 1 - IF (JOB .NE. 0) GO TO 50 -C -C JOB = 0 , SOLVE A * X = B -C FIRST SOLVE L*Y = B -C - IF (ML .EQ. 0) GO TO 30 - IF (NM1 .LT. 1) GO TO 30 - DO 20 K = 1, NM1 - LM = MIN0(ML,N-K) - L = IPVT(K) - T = B(L) - IF (L .EQ. K) GO TO 10 - B(L) = B(K) - B(K) = T - 10 CONTINUE - CALL DAXPY(LM,T,ABD(M+1,K),1,B(K+1),1) - 20 CONTINUE - 30 CONTINUE -C -C NOW SOLVE U*X = Y -C - DO 40 KB = 1, N - K = N + 1 - KB - B(K) = B(K)/ABD(M,K) - LM = MIN0(K,M) - 1 - LA = M - LM - LB = K - LM - T = -B(K) - CALL DAXPY(LM,T,ABD(LA,K),1,B(LB),1) - 40 CONTINUE - GO TO 100 - 50 CONTINUE -C -C JOB = NONZERO, SOLVE TRANS(A) * X = B -C FIRST SOLVE TRANS(U)*Y = B -C - DO 60 K = 1, N - LM = MIN0(K,M) - 1 - LA = M - LM - LB = K - LM - T = DDOT(LM,ABD(LA,K),1,B(LB),1) - B(K) = (B(K) - T)/ABD(M,K) - 60 CONTINUE -C -C NOW SOLVE TRANS(L)*X = Y -C - IF (ML .EQ. 0) GO TO 90 - IF (NM1 .LT. 1) GO TO 90 - DO 80 KB = 1, NM1 - K = N - KB - LM = MIN0(ML,N-K) - B(K) = B(K) + DDOT(LM,ABD(M+1,K),1,B(K+1),1) - L = IPVT(K) - IF (L .EQ. K) GO TO 70 - T = B(L) - B(L) = B(K) - B(K) = T - 70 CONTINUE - 80 CONTINUE - 90 CONTINUE - 100 CONTINUE - RETURN - END
--- a/libcruft/linpack/dgeco.f Sun Feb 16 04:29:00 2003 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,193 +0,0 @@ - SUBROUTINE DGECO(A,LDA,N,IPVT,RCOND,Z) - INTEGER LDA,N,IPVT(1) - DOUBLE PRECISION A(LDA,1),Z(1) - DOUBLE PRECISION RCOND -C -C DGECO FACTORS A DOUBLE PRECISION MATRIX BY GAUSSIAN ELIMINATION -C AND ESTIMATES THE CONDITION OF THE MATRIX. -C -C IF RCOND IS NOT NEEDED, DGEFA IS SLIGHTLY FASTER. -C TO SOLVE A*X = B , FOLLOW DGECO BY DGESL. -C TO COMPUTE INVERSE(A)*C , FOLLOW DGECO BY DGESL. -C TO COMPUTE DETERMINANT(A) , FOLLOW DGECO BY DGEDI. -C TO COMPUTE INVERSE(A) , FOLLOW DGECO BY DGEDI. -C -C ON ENTRY -C -C A DOUBLE PRECISION(LDA, N) -C THE MATRIX TO BE FACTORED. -C -C LDA INTEGER -C THE LEADING DIMENSION OF THE ARRAY A . -C -C N INTEGER -C THE ORDER OF THE MATRIX A . -C -C ON RETURN -C -C A AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS -C WHICH WERE USED TO OBTAIN IT. -C THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE -C L IS A PRODUCT OF PERMUTATION AND UNIT LOWER -C TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR. -C -C IPVT INTEGER(N) -C AN INTEGER VECTOR OF PIVOT INDICES. -C -C RCOND DOUBLE PRECISION -C AN ESTIMATE OF THE RECIPROCAL CONDITION OF A . -C FOR THE SYSTEM A*X = B , RELATIVE PERTURBATIONS -C IN A AND B OF SIZE EPSILON MAY CAUSE -C RELATIVE PERTURBATIONS IN X OF SIZE EPSILON/RCOND . -C IF RCOND IS SO SMALL THAT THE LOGICAL EXPRESSION -C 1.0 + RCOND .EQ. 1.0 -C IS TRUE, THEN A MAY BE SINGULAR TO WORKING -C PRECISION. IN PARTICULAR, RCOND IS ZERO IF -C EXACT SINGULARITY IS DETECTED OR THE ESTIMATE -C UNDERFLOWS. -C -C Z DOUBLE PRECISION(N) -C A WORK VECTOR WHOSE CONTENTS ARE USUALLY UNIMPORTANT. -C IF A IS CLOSE TO A SINGULAR MATRIX, THEN Z IS -C AN APPROXIMATE NULL VECTOR IN THE SENSE THAT -C NORM(A*Z) = RCOND*NORM(A)*NORM(Z) . -C -C LINPACK. THIS VERSION DATED 08/14/78 . -C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. -C -C SUBROUTINES AND FUNCTIONS -C -C LINPACK DGEFA -C BLAS DAXPY,DDOT,DSCAL,DASUM -C FORTRAN DABS,DMAX1,DSIGN -C -C INTERNAL VARIABLES -C - DOUBLE PRECISION DDOT,EK,T,WK,WKM - DOUBLE PRECISION ANORM,S,DASUM,SM,YNORM - INTEGER INFO,J,K,KB,KP1,L -C -C -C COMPUTE 1-NORM OF A -C - ANORM = 0.0D0 - DO 10 J = 1, N - ANORM = DMAX1(ANORM,DASUM(N,A(1,J),1)) - 10 CONTINUE -C -C FACTOR -C - CALL DGEFA(A,LDA,N,IPVT,INFO) -C -C RCOND = 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A)))) . -C ESTIMATE = NORM(Z)/NORM(Y) WHERE A*Z = Y AND TRANS(A)*Y = E . -C TRANS(A) IS THE TRANSPOSE OF A . THE COMPONENTS OF E ARE -C CHOSEN TO CAUSE MAXIMUM LOCAL GROWTH IN THE ELEMENTS OF W WHERE -C TRANS(U)*W = E . THE VECTORS ARE FREQUENTLY RESCALED TO AVOID -C OVERFLOW. -C -C SOLVE TRANS(U)*W = E -C - EK = 1.0D0 - DO 20 J = 1, N - Z(J) = 0.0D0 - 20 CONTINUE - DO 100 K = 1, N - IF (Z(K) .NE. 0.0D0) EK = DSIGN(EK,-Z(K)) - IF (DABS(EK-Z(K)) .LE. DABS(A(K,K))) GO TO 30 - S = DABS(A(K,K))/DABS(EK-Z(K)) - CALL DSCAL(N,S,Z,1) - EK = S*EK - 30 CONTINUE - WK = EK - Z(K) - WKM = -EK - Z(K) - S = DABS(WK) - SM = DABS(WKM) - IF (A(K,K) .EQ. 0.0D0) GO TO 40 - WK = WK/A(K,K) - WKM = WKM/A(K,K) - GO TO 50 - 40 CONTINUE - WK = 1.0D0 - WKM = 1.0D0 - 50 CONTINUE - KP1 = K + 1 - IF (KP1 .GT. N) GO TO 90 - DO 60 J = KP1, N - SM = SM + DABS(Z(J)+WKM*A(K,J)) - Z(J) = Z(J) + WK*A(K,J) - S = S + DABS(Z(J)) - 60 CONTINUE - IF (S .GE. SM) GO TO 80 - T = WKM - WK - WK = WKM - DO 70 J = KP1, N - Z(J) = Z(J) + T*A(K,J) - 70 CONTINUE - 80 CONTINUE - 90 CONTINUE - Z(K) = WK - 100 CONTINUE - S = 1.0D0/DASUM(N,Z,1) - CALL DSCAL(N,S,Z,1) -C -C SOLVE TRANS(L)*Y = W -C - DO 120 KB = 1, N - K = N + 1 - KB - IF (K .LT. N) Z(K) = Z(K) + DDOT(N-K,A(K+1,K),1,Z(K+1),1) - IF (DABS(Z(K)) .LE. 1.0D0) GO TO 110 - S = 1.0D0/DABS(Z(K)) - CALL DSCAL(N,S,Z,1) - 110 CONTINUE - L = IPVT(K) - T = Z(L) - Z(L) = Z(K) - Z(K) = T - 120 CONTINUE - S = 1.0D0/DASUM(N,Z,1) - CALL DSCAL(N,S,Z,1) -C - YNORM = 1.0D0 -C -C SOLVE L*V = Y -C - DO 140 K = 1, N - L = IPVT(K) - T = Z(L) - Z(L) = Z(K) - Z(K) = T - IF (K .LT. N) CALL DAXPY(N-K,T,A(K+1,K),1,Z(K+1),1) - IF (DABS(Z(K)) .LE. 1.0D0) GO TO 130 - S = 1.0D0/DABS(Z(K)) - CALL DSCAL(N,S,Z,1) - YNORM = S*YNORM - 130 CONTINUE - 140 CONTINUE - S = 1.0D0/DASUM(N,Z,1) - CALL DSCAL(N,S,Z,1) - YNORM = S*YNORM -C -C SOLVE U*Z = V -C - DO 160 KB = 1, N - K = N + 1 - KB - IF (DABS(Z(K)) .LE. DABS(A(K,K))) GO TO 150 - S = DABS(A(K,K))/DABS(Z(K)) - CALL DSCAL(N,S,Z,1) - YNORM = S*YNORM - 150 CONTINUE - IF (A(K,K) .NE. 0.0D0) Z(K) = Z(K)/A(K,K) - IF (A(K,K) .EQ. 0.0D0) Z(K) = 1.0D0 - T = -Z(K) - CALL DAXPY(K-1,T,A(1,K),1,Z(1),1) - 160 CONTINUE -C MAKE ZNORM = 1.0 - S = 1.0D0/DASUM(N,Z,1) - CALL DSCAL(N,S,Z,1) - YNORM = S*YNORM -C - IF (ANORM .NE. 0.0D0) RCOND = YNORM/ANORM - IF (ANORM .EQ. 0.0D0) RCOND = 0.0D0 - RETURN - END
--- a/libcruft/linpack/dgedi.f Sun Feb 16 04:29:00 2003 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,128 +0,0 @@ - SUBROUTINE DGEDI(A,LDA,N,IPVT,DET,WORK,JOB) - INTEGER LDA,N,IPVT(1),JOB - DOUBLE PRECISION A(LDA,1),DET(2),WORK(1) -C -C DGEDI COMPUTES THE DETERMINANT AND INVERSE OF A MATRIX -C USING THE FACTORS COMPUTED BY DGECO OR DGEFA. -C -C ON ENTRY -C -C A DOUBLE PRECISION(LDA, N) -C THE OUTPUT FROM DGECO OR DGEFA. -C -C LDA INTEGER -C THE LEADING DIMENSION OF THE ARRAY A . -C -C N INTEGER -C THE ORDER OF THE MATRIX A . -C -C IPVT INTEGER(N) -C THE PIVOT VECTOR FROM DGECO OR DGEFA. -C -C WORK DOUBLE PRECISION(N) -C WORK VECTOR. CONTENTS DESTROYED. -C -C JOB INTEGER -C = 11 BOTH DETERMINANT AND INVERSE. -C = 01 INVERSE ONLY. -C = 10 DETERMINANT ONLY. -C -C ON RETURN -C -C A INVERSE OF ORIGINAL MATRIX IF REQUESTED. -C OTHERWISE UNCHANGED. -C -C DET DOUBLE PRECISION(2) -C DETERMINANT OF ORIGINAL MATRIX IF REQUESTED. -C OTHERWISE NOT REFERENCED. -C DETERMINANT = DET(1) * 10.0**DET(2) -C WITH 1.0 .LE. DABS(DET(1)) .LT. 10.0 -C OR DET(1) .EQ. 0.0 . -C -C ERROR CONDITION -C -C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS -C A ZERO ON THE DIAGONAL AND THE INVERSE IS REQUESTED. -C IT WILL NOT OCCUR IF THE SUBROUTINES ARE CALLED CORRECTLY -C AND IF DGECO HAS SET RCOND .GT. 0.0 OR DGEFA HAS SET -C INFO .EQ. 0 . -C -C LINPACK. THIS VERSION DATED 08/14/78 . -C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. -C -C SUBROUTINES AND FUNCTIONS -C -C BLAS DAXPY,DSCAL,DSWAP -C FORTRAN DABS,MOD -C -C INTERNAL VARIABLES -C - DOUBLE PRECISION T - DOUBLE PRECISION TEN - INTEGER I,J,K,KB,KP1,L,NM1 -C -C -C COMPUTE DETERMINANT -C - IF (JOB/10 .EQ. 0) GO TO 70 - DET(1) = 1.0D0 - DET(2) = 0.0D0 - TEN = 10.0D0 - DO 50 I = 1, N - IF (IPVT(I) .NE. I) DET(1) = -DET(1) - DET(1) = A(I,I)*DET(1) -C ...EXIT - IF (DET(1) .EQ. 0.0D0) GO TO 60 - 10 IF (DABS(DET(1)) .GE. 1.0D0) GO TO 20 - DET(1) = TEN*DET(1) - DET(2) = DET(2) - 1.0D0 - GO TO 10 - 20 CONTINUE - 30 IF (DABS(DET(1)) .LT. TEN) GO TO 40 - DET(1) = DET(1)/TEN - DET(2) = DET(2) + 1.0D0 - GO TO 30 - 40 CONTINUE - 50 CONTINUE - 60 CONTINUE - 70 CONTINUE -C -C COMPUTE INVERSE(U) -C - IF (MOD(JOB,10) .EQ. 0) GO TO 150 - DO 100 K = 1, N - A(K,K) = 1.0D0/A(K,K) - T = -A(K,K) - CALL DSCAL(K-1,T,A(1,K),1) - KP1 = K + 1 - IF (N .LT. KP1) GO TO 90 - DO 80 J = KP1, N - T = A(K,J) - A(K,J) = 0.0D0 - CALL DAXPY(K,T,A(1,K),1,A(1,J),1) - 80 CONTINUE - 90 CONTINUE - 100 CONTINUE -C -C FORM INVERSE(U)*INVERSE(L) -C - NM1 = N - 1 - IF (NM1 .LT. 1) GO TO 140 - DO 130 KB = 1, NM1 - K = N - KB - KP1 = K + 1 - DO 110 I = KP1, N - WORK(I) = A(I,K) - A(I,K) = 0.0D0 - 110 CONTINUE - DO 120 J = KP1, N - T = WORK(J) - CALL DAXPY(N,T,A(1,J),1,A(1,K),1) - 120 CONTINUE - L = IPVT(K) - IF (L .NE. K) CALL DSWAP(N,A(1,K),1,A(1,L),1) - 130 CONTINUE - 140 CONTINUE - 150 CONTINUE - RETURN - END
--- a/libcruft/linpack/dgefa.f Sun Feb 16 04:29:00 2003 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,103 +0,0 @@ - SUBROUTINE DGEFA(A,LDA,N,IPVT,INFO) - INTEGER LDA,N,IPVT(1),INFO - DOUBLE PRECISION A(LDA,1) -C -C DGEFA FACTORS A DOUBLE PRECISION MATRIX BY GAUSSIAN ELIMINATION. -C -C DGEFA IS USUALLY CALLED BY DGECO, BUT IT CAN BE CALLED -C DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED. -C (TIME FOR DGECO) = (1 + 9/N)*(TIME FOR DGEFA) . -C -C ON ENTRY -C -C A DOUBLE PRECISION(LDA, N) -C THE MATRIX TO BE FACTORED. -C -C LDA INTEGER -C THE LEADING DIMENSION OF THE ARRAY A . -C -C N INTEGER -C THE ORDER OF THE MATRIX A . -C -C ON RETURN -C -C A AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS -C WHICH WERE USED TO OBTAIN IT. -C THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE -C L IS A PRODUCT OF PERMUTATION AND UNIT LOWER -C TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR. -C -C IPVT INTEGER(N) -C AN INTEGER VECTOR OF PIVOT INDICES. -C -C INFO INTEGER -C = 0 NORMAL VALUE. -C = K IF U(K,K) .EQ. 0.0 . THIS IS NOT AN ERROR -C CONDITION FOR THIS SUBROUTINE, BUT IT DOES -C INDICATE THAT DGESL OR DGEDI WILL DIVIDE BY ZERO -C IF CALLED. USE RCOND IN DGECO FOR A RELIABLE -C INDICATION OF SINGULARITY. -C -C LINPACK. THIS VERSION DATED 08/14/78 . -C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. -C -C SUBROUTINES AND FUNCTIONS -C -C BLAS DAXPY,DSCAL,IDAMAX -C -C INTERNAL VARIABLES -C - DOUBLE PRECISION T - INTEGER IDAMAX,J,K,KP1,L,NM1 -C -C -C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING -C - INFO = 0 - NM1 = N - 1 - IF (NM1 .LT. 1) GO TO 70 - DO 60 K = 1, NM1 - KP1 = K + 1 -C -C FIND L = PIVOT INDEX -C - L = IDAMAX(N-K+1,A(K,K),1) + K - 1 - IPVT(K) = L -C -C ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED -C - IF (A(L,K) .EQ. 0.0D0) GO TO 40 -C -C INTERCHANGE IF NECESSARY -C - IF (L .EQ. K) GO TO 10 - T = A(L,K) - A(L,K) = A(K,K) - A(K,K) = T - 10 CONTINUE -C -C COMPUTE MULTIPLIERS -C - T = -1.0D0/A(K,K) - CALL DSCAL(N-K,T,A(K+1,K),1) -C -C ROW ELIMINATION WITH COLUMN INDEXING -C - DO 30 J = KP1, N - T = A(L,J) - IF (L .EQ. K) GO TO 20 - A(L,J) = A(K,J) - A(K,J) = T - 20 CONTINUE - CALL DAXPY(N-K,T,A(K+1,K),1,A(K+1,J),1) - 30 CONTINUE - GO TO 50 - 40 CONTINUE - INFO = K - 50 CONTINUE - 60 CONTINUE - 70 CONTINUE - IPVT(N) = N - IF (A(N,N) .EQ. 0.0D0) INFO = N - RETURN - END
--- a/libcruft/linpack/dgesl.f Sun Feb 16 04:29:00 2003 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,117 +0,0 @@ - SUBROUTINE DGESL(A,LDA,N,IPVT,B,JOB) - INTEGER LDA,N,IPVT(1),JOB - DOUBLE PRECISION A(LDA,1),B(1) -C -C DGESL SOLVES THE DOUBLE PRECISION SYSTEM -C A * X = B OR TRANS(A) * X = B -C USING THE FACTORS COMPUTED BY DGECO OR DGEFA. -C -C ON ENTRY -C -C A DOUBLE PRECISION(LDA, N) -C THE OUTPUT FROM DGECO OR DGEFA. -C -C LDA INTEGER -C THE LEADING DIMENSION OF THE ARRAY A . -C -C N INTEGER -C THE ORDER OF THE MATRIX A . -C -C IPVT INTEGER(N) -C THE PIVOT VECTOR FROM DGECO OR DGEFA. -C -C B DOUBLE PRECISION(N) -C THE RIGHT HAND SIDE VECTOR. -C -C JOB INTEGER -C = 0 TO SOLVE A*X = B , -C = NONZERO TO SOLVE TRANS(A)*X = B WHERE -C TRANS(A) IS THE TRANSPOSE. -C -C ON RETURN -C -C B THE SOLUTION VECTOR X . -C -C ERROR CONDITION -C -C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS A -C ZERO ON THE DIAGONAL. TECHNICALLY THIS INDICATES SINGULARITY -C BUT IT IS OFTEN CAUSED BY IMPROPER ARGUMENTS OR IMPROPER -C SETTING OF LDA . IT WILL NOT OCCUR IF THE SUBROUTINES ARE -C CALLED CORRECTLY AND IF DGECO HAS SET RCOND .GT. 0.0 -C OR DGEFA HAS SET INFO .EQ. 0 . -C -C TO COMPUTE INVERSE(A) * C WHERE C IS A MATRIX -C WITH P COLUMNS -C CALL DGECO(A,LDA,N,IPVT,RCOND,Z) -C IF (RCOND IS TOO SMALL) GO TO ... -C DO 10 J = 1, P -C CALL DGESL(A,LDA,N,IPVT,C(1,J),0) -C 10 CONTINUE -C -C LINPACK. THIS VERSION DATED 08/14/78 . -C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. -C -C SUBROUTINES AND FUNCTIONS -C -C BLAS DAXPY,DDOT -C -C INTERNAL VARIABLES -C - DOUBLE PRECISION DDOT,T - INTEGER K,KB,L,NM1 -C - NM1 = N - 1 - IF (JOB .NE. 0) GO TO 50 -C -C JOB = 0 , SOLVE A * X = B -C FIRST SOLVE L*Y = B -C - IF (NM1 .LT. 1) GO TO 30 - DO 20 K = 1, NM1 - L = IPVT(K) - T = B(L) - IF (L .EQ. K) GO TO 10 - B(L) = B(K) - B(K) = T - 10 CONTINUE - CALL DAXPY(N-K,T,A(K+1,K),1,B(K+1),1) - 20 CONTINUE - 30 CONTINUE -C -C NOW SOLVE U*X = Y -C - DO 40 KB = 1, N - K = N + 1 - KB - B(K) = B(K)/A(K,K) - T = -B(K) - CALL DAXPY(K-1,T,A(1,K),1,B(1),1) - 40 CONTINUE - GO TO 100 - 50 CONTINUE -C -C JOB = NONZERO, SOLVE TRANS(A) * X = B -C FIRST SOLVE TRANS(U)*Y = B -C - DO 60 K = 1, N - T = DDOT(K-1,A(1,K),1,B(1),1) - B(K) = (B(K) - T)/A(K,K) - 60 CONTINUE -C -C NOW SOLVE TRANS(L)*X = Y -C - IF (NM1 .LT. 1) GO TO 90 - DO 80 KB = 1, NM1 - K = N - KB - B(K) = B(K) + DDOT(N-K,A(K+1,K),1,B(K+1),1) - L = IPVT(K) - IF (L .EQ. K) GO TO 70 - T = B(L) - B(L) = B(K) - B(K) = T - 70 CONTINUE - 80 CONTINUE - 90 CONTINUE - 100 CONTINUE - RETURN - END
--- a/libcruft/linpack/spofa.f Sun Feb 16 04:29:00 2003 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,71 +0,0 @@ - SUBROUTINE SPOFA(A,LDA,N,INFO) - INTEGER LDA,N,INFO - REAL A(LDA,1) -C -C SPOFA FACTORS A REAL SYMMETRIC POSITIVE DEFINITE MATRIX. -C -C SPOFA IS USUALLY CALLED BY SPOCO, BUT IT CAN BE CALLED -C DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED. -C (TIME FOR SPOCO) = (1 + 18/N)*(TIME FOR SPOFA) . -C -C ON ENTRY -C -C A REAL(LDA, N) -C THE SYMMETRIC MATRIX TO BE FACTORED. ONLY THE -C DIAGONAL AND UPPER TRIANGLE ARE USED. -C -C LDA INTEGER -C THE LEADING DIMENSION OF THE ARRAY A . -C -C N INTEGER -C THE ORDER OF THE MATRIX A . -C -C ON RETURN -C -C A AN UPPER TRIANGULAR MATRIX R SO THAT A = TRANS(R)*R -C WHERE TRANS(R) IS THE TRANSPOSE. -C THE STRICT LOWER TRIANGLE IS UNALTERED. -C IF INFO .NE. 0 , THE FACTORIZATION IS NOT COMPLETE. -C -C INFO INTEGER -C = 0 FOR NORMAL RETURN. -C = K SIGNALS AN ERROR CONDITION. THE LEADING MINOR -C OF ORDER K IS NOT POSITIVE DEFINITE. -C -C LINPACK. THIS VERSION DATED 08/14/78 . -C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. -C -C SUBROUTINES AND FUNCTIONS -C -C BLAS SDOT -C FORTRAN SQRT -C -C INTERNAL VARIABLES -C - REAL SDOT,T - REAL S - INTEGER J,JM1,K -C BEGIN BLOCK WITH ...EXITS TO 40 -C -C - DO 30 J = 1, N - INFO = J - S = 0.0E0 - JM1 = J - 1 - IF (JM1 .LT. 1) GO TO 20 - DO 10 K = 1, JM1 - T = A(K,J) - SDOT(K-1,A(1,K),1,A(1,J),1) - T = T/A(K,K) - A(K,J) = T - S = S + T*T - 10 CONTINUE - 20 CONTINUE - S = A(J,J) - S -C ......EXIT - IF (S .LE. 0.0E0) GO TO 40 - A(J,J) = SQRT(S) - 30 CONTINUE - INFO = 0 - 40 CONTINUE - RETURN - END
--- a/libcruft/linpack/zgeco.f Sun Feb 16 04:29:00 2003 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,201 +0,0 @@ - subroutine zgeco(a,lda,n,ipvt,rcond,z) - integer lda,n,ipvt(1) - complex*16 a(lda,1),z(1) - double precision rcond -c -c zgeco factors a complex*16 matrix by gaussian elimination -c and estimates the condition of the matrix. -c -c if rcond is not needed, zgefa is slightly faster. -c to solve a*x = b , follow zgeco by zgesl. -c to compute inverse(a)*c , follow zgeco by zgesl. -c to compute determinant(a) , follow zgeco by zgedi. -c to compute inverse(a) , follow zgeco by zgedi. -c -c on entry -c -c a complex*16(lda, n) -c the matrix to be factored. -c -c lda integer -c the leading dimension of the array a . -c -c n integer -c the order of the matrix a . -c -c on return -c -c a an upper triangular matrix and the multipliers -c which were used to obtain it. -c the factorization can be written a = l*u where -c l is a product of permutation and unit lower -c triangular matrices and u is upper triangular. -c -c ipvt integer(n) -c an integer vector of pivot indices. -c -c rcond double precision -c an estimate of the reciprocal condition of a . -c for the system a*x = b , relative perturbations -c in a and b of size epsilon may cause -c relative perturbations in x of size epsilon/rcond . -c if rcond is so small that the logical expression -c 1.0 + rcond .eq. 1.0 -c is true, then a may be singular to working -c precision. in particular, rcond is zero if -c exact singularity is detected or the estimate -c underflows. -c -c z complex*16(n) -c a work vector whose contents are usually unimportant. -c if a is close to a singular matrix, then z is -c an approximate null vector in the sense that -c norm(a*z) = rcond*norm(a)*norm(z) . -c -c linpack. this version dated 08/14/78 . -c cleve moler, university of new mexico, argonne national lab. -c -c subroutines and functions -c -c linpack zgefa -c blas zaxpy,zdotc,zdscal,dzasum -c fortran dabs,dmax1,dcmplx,dconjg -c -c internal variables -c - complex*16 zdotc,ek,t,wk,wkm - double precision anorm,s,dzasum,sm,ynorm - integer info,j,k,kb,kp1,l -c - complex*16 zdum,zdum1,zdum2,csign1 - double precision cabs1 - double precision dreal,dimag - complex*16 zdumr,zdumi - dreal(zdumr) = zdumr - dimag(zdumi) = (0.0d0,-1.0d0)*zdumi - cabs1(zdum) = dabs(dreal(zdum)) + dabs(dimag(zdum)) - csign1(zdum1,zdum2) = cabs1(zdum1)*(zdum2/cabs1(zdum2)) -c -c compute 1-norm of a -c - anorm = 0.0d0 - do 10 j = 1, n - anorm = dmax1(anorm,dzasum(n,a(1,j),1)) - 10 continue -c -c factor -c - call zgefa(a,lda,n,ipvt,info) -c -c rcond = 1/(norm(a)*(estimate of norm(inverse(a)))) . -c estimate = norm(z)/norm(y) where a*z = y and ctrans(a)*y = e . -c ctrans(a) is the conjugate transpose of a . -c the components of e are chosen to cause maximum local -c growth in the elements of w where ctrans(u)*w = e . -c the vectors are frequently rescaled to avoid overflow. -c -c solve ctrans(u)*w = e -c - ek = (1.0d0,0.0d0) - do 20 j = 1, n - z(j) = (0.0d0,0.0d0) - 20 continue - do 100 k = 1, n - if (cabs1(z(k)) .ne. 0.0d0) ek = csign1(ek,-z(k)) - if (cabs1(ek-z(k)) .le. cabs1(a(k,k))) go to 30 - s = cabs1(a(k,k))/cabs1(ek-z(k)) - call zdscal(n,s,z,1) - ek = dcmplx(s,0.0d0)*ek - 30 continue - wk = ek - z(k) - wkm = -ek - z(k) - s = cabs1(wk) - sm = cabs1(wkm) - if (cabs1(a(k,k)) .eq. 0.0d0) go to 40 - wk = wk/dconjg(a(k,k)) - wkm = wkm/dconjg(a(k,k)) - go to 50 - 40 continue - wk = (1.0d0,0.0d0) - wkm = (1.0d0,0.0d0) - 50 continue - kp1 = k + 1 - if (kp1 .gt. n) go to 90 - do 60 j = kp1, n - sm = sm + cabs1(z(j)+wkm*dconjg(a(k,j))) - z(j) = z(j) + wk*dconjg(a(k,j)) - s = s + cabs1(z(j)) - 60 continue - if (s .ge. sm) go to 80 - t = wkm - wk - wk = wkm - do 70 j = kp1, n - z(j) = z(j) + t*dconjg(a(k,j)) - 70 continue - 80 continue - 90 continue - z(k) = wk - 100 continue - s = 1.0d0/dzasum(n,z,1) - call zdscal(n,s,z,1) -c -c solve ctrans(l)*y = w -c - do 120 kb = 1, n - k = n + 1 - kb - if (k .lt. n) z(k) = z(k) + zdotc(n-k,a(k+1,k),1,z(k+1),1) - if (cabs1(z(k)) .le. 1.0d0) go to 110 - s = 1.0d0/cabs1(z(k)) - call zdscal(n,s,z,1) - 110 continue - l = ipvt(k) - t = z(l) - z(l) = z(k) - z(k) = t - 120 continue - s = 1.0d0/dzasum(n,z,1) - call zdscal(n,s,z,1) -c - ynorm = 1.0d0 -c -c solve l*v = y -c - do 140 k = 1, n - l = ipvt(k) - t = z(l) - z(l) = z(k) - z(k) = t - if (k .lt. n) call zaxpy(n-k,t,a(k+1,k),1,z(k+1),1) - if (cabs1(z(k)) .le. 1.0d0) go to 130 - s = 1.0d0/cabs1(z(k)) - call zdscal(n,s,z,1) - ynorm = s*ynorm - 130 continue - 140 continue - s = 1.0d0/dzasum(n,z,1) - call zdscal(n,s,z,1) - ynorm = s*ynorm -c -c solve u*z = v -c - do 160 kb = 1, n - k = n + 1 - kb - if (cabs1(z(k)) .le. cabs1(a(k,k))) go to 150 - s = cabs1(a(k,k))/cabs1(z(k)) - call zdscal(n,s,z,1) - ynorm = s*ynorm - 150 continue - if (cabs1(a(k,k)) .ne. 0.0d0) z(k) = z(k)/a(k,k) - if (cabs1(a(k,k)) .eq. 0.0d0) z(k) = (1.0d0,0.0d0) - t = -z(k) - call zaxpy(k-1,t,a(1,k),1,z(1),1) - 160 continue -c make znorm = 1.0 - s = 1.0d0/dzasum(n,z,1) - call zdscal(n,s,z,1) - ynorm = s*ynorm -c - if (anorm .ne. 0.0d0) rcond = ynorm/anorm - if (anorm .eq. 0.0d0) rcond = 0.0d0 - return - end
--- a/libcruft/linpack/zgedi.f Sun Feb 16 04:29:00 2003 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,135 +0,0 @@ - subroutine zgedi(a,lda,n,ipvt,det,work,job) - integer lda,n,ipvt(1),job - complex*16 a(lda,1),det(2),work(1) -c -c zgedi computes the determinant and inverse of a matrix -c using the factors computed by zgeco or zgefa. -c -c on entry -c -c a complex*16(lda, n) -c the output from zgeco or zgefa. -c -c lda integer -c the leading dimension of the array a . -c -c n integer -c the order of the matrix a . -c -c ipvt integer(n) -c the pivot vector from zgeco or zgefa. -c -c work complex*16(n) -c work vector. contents destroyed. -c -c job integer -c = 11 both determinant and inverse. -c = 01 inverse only. -c = 10 determinant only. -c -c on return -c -c a inverse of original matrix if requested. -c otherwise unchanged. -c -c det complex*16(2) -c determinant of original matrix if requested. -c otherwise not referenced. -c determinant = det(1) * 10.0**det(2) -c with 1.0 .le. cabs1(det(1)) .lt. 10.0 -c or det(1) .eq. 0.0 . -c -c error condition -c -c a division by zero will occur if the input factor contains -c a zero on the diagonal and the inverse is requested. -c it will not occur if the subroutines are called correctly -c and if zgeco has set rcond .gt. 0.0 or zgefa has set -c info .eq. 0 . -c -c linpack. this version dated 08/14/78 . -c cleve moler, university of new mexico, argonne national lab. -c -c subroutines and functions -c -c blas zaxpy,zscal,zswap -c fortran dabs,dcmplx,mod -c -c internal variables -c - complex*16 t - double precision ten - integer i,j,k,kb,kp1,l,nm1 -c - complex*16 zdum - double precision cabs1 - double precision dreal,dimag - complex*16 zdumr,zdumi - dreal(zdumr) = zdumr - dimag(zdumi) = (0.0d0,-1.0d0)*zdumi - cabs1(zdum) = dabs(dreal(zdum)) + dabs(dimag(zdum)) -c -c compute determinant -c - if (job/10 .eq. 0) go to 70 - det(1) = (1.0d0,0.0d0) - det(2) = (0.0d0,0.0d0) - ten = 10.0d0 - do 50 i = 1, n - if (ipvt(i) .ne. i) det(1) = -det(1) - det(1) = a(i,i)*det(1) -c ...exit - if (cabs1(det(1)) .eq. 0.0d0) go to 60 - 10 if (cabs1(det(1)) .ge. 1.0d0) go to 20 - det(1) = dcmplx(ten,0.0d0)*det(1) - det(2) = det(2) - (1.0d0,0.0d0) - go to 10 - 20 continue - 30 if (cabs1(det(1)) .lt. ten) go to 40 - det(1) = det(1)/dcmplx(ten,0.0d0) - det(2) = det(2) + (1.0d0,0.0d0) - go to 30 - 40 continue - 50 continue - 60 continue - 70 continue -c -c compute inverse(u) -c - if (mod(job,10) .eq. 0) go to 150 - do 100 k = 1, n - a(k,k) = (1.0d0,0.0d0)/a(k,k) - t = -a(k,k) - call zscal(k-1,t,a(1,k),1) - kp1 = k + 1 - if (n .lt. kp1) go to 90 - do 80 j = kp1, n - t = a(k,j) - a(k,j) = (0.0d0,0.0d0) - call zaxpy(k,t,a(1,k),1,a(1,j),1) - 80 continue - 90 continue - 100 continue -c -c form inverse(u)*inverse(l) -c - nm1 = n - 1 - if (nm1 .lt. 1) go to 140 - do 130 kb = 1, nm1 - k = n - kb - kp1 = k + 1 - do 110 i = kp1, n - work(i) = a(i,k) - a(i,k) = (0.0d0,0.0d0) - 110 continue - do 120 j = kp1, n - t = work(j) - call zaxpy(n,t,a(1,j),1,a(1,k),1) - 120 continue - l = ipvt(k) - if (l .ne. k) call zswap(n,a(1,k),1,a(1,l),1) - 130 continue - 140 continue - 150 continue - return - end
--- a/libcruft/linpack/zgefa.f Sun Feb 16 04:29:00 2003 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,111 +0,0 @@ - subroutine zgefa(a,lda,n,ipvt,info) - integer lda,n,ipvt(1),info - complex*16 a(lda,1) -c -c zgefa factors a complex*16 matrix by gaussian elimination. -c -c zgefa is usually called by zgeco, but it can be called -c directly with a saving in time if rcond is not needed. -c (time for zgeco) = (1 + 9/n)*(time for zgefa) . -c -c on entry -c -c a complex*16(lda, n) -c the matrix to be factored. -c -c lda integer -c the leading dimension of the array a . -c -c n integer -c the order of the matrix a . -c -c on return -c -c a an upper triangular matrix and the multipliers -c which were used to obtain it. -c the factorization can be written a = l*u where -c l is a product of permutation and unit lower -c triangular matrices and u is upper triangular. -c -c ipvt integer(n) -c an integer vector of pivot indices. -c -c info integer -c = 0 normal value. -c = k if u(k,k) .eq. 0.0 . this is not an error -c condition for this subroutine, but it does -c indicate that zgesl or zgedi will divide by zero -c if called. use rcond in zgeco for a reliable -c indication of singularity. -c -c linpack. this version dated 08/14/78 . -c cleve moler, university of new mexico, argonne national lab. -c -c subroutines and functions -c -c blas zaxpy,zscal,izamax -c fortran dabs -c -c internal variables -c - complex*16 t - integer izamax,j,k,kp1,l,nm1 -c - complex*16 zdum - double precision cabs1 - double precision dreal,dimag - complex*16 zdumr,zdumi - dreal(zdumr) = zdumr - dimag(zdumi) = (0.0d0,-1.0d0)*zdumi - cabs1(zdum) = dabs(dreal(zdum)) + dabs(dimag(zdum)) -c -c gaussian elimination with partial pivoting -c - info = 0 - nm1 = n - 1 - if (nm1 .lt. 1) go to 70 - do 60 k = 1, nm1 - kp1 = k + 1 -c -c find l = pivot index -c - l = izamax(n-k+1,a(k,k),1) + k - 1 - ipvt(k) = l -c -c zero pivot implies this column already triangularized -c - if (cabs1(a(l,k)) .eq. 0.0d0) go to 40 -c -c interchange if necessary -c - if (l .eq. k) go to 10 - t = a(l,k) - a(l,k) = a(k,k) - a(k,k) = t - 10 continue -c -c compute multipliers -c - t = -(1.0d0,0.0d0)/a(k,k) - call zscal(n-k,t,a(k+1,k),1) -c -c row elimination with column indexing -c - do 30 j = kp1, n - t = a(l,j) - if (l .eq. k) go to 20 - a(l,j) = a(k,j) - a(k,j) = t - 20 continue - call zaxpy(n-k,t,a(k+1,k),1,a(k+1,j),1) - 30 continue - go to 50 - 40 continue - info = k - 50 continue - 60 continue - 70 continue - ipvt(n) = n - if (cabs1(a(n,n)) .eq. 0.0d0) info = n - return - end
--- a/libcruft/linpack/zgesl.f Sun Feb 16 04:29:00 2003 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,122 +0,0 @@ - subroutine zgesl(a,lda,n,ipvt,b,job) - integer lda,n,ipvt(1),job - complex*16 a(lda,1),b(1) -c -c zgesl solves the complex*16 system -c a * x = b or ctrans(a) * x = b -c using the factors computed by zgeco or zgefa. -c -c on entry -c -c a complex*16(lda, n) -c the output from zgeco or zgefa. -c -c lda integer -c the leading dimension of the array a . -c -c n integer -c the order of the matrix a . -c -c ipvt integer(n) -c the pivot vector from zgeco or zgefa. -c -c b complex*16(n) -c the right hand side vector. -c -c job integer -c = 0 to solve a*x = b , -c = nonzero to solve ctrans(a)*x = b where -c ctrans(a) is the conjugate transpose. -c -c on return -c -c b the solution vector x . -c -c error condition -c -c a division by zero will occur if the input factor contains a -c zero on the diagonal. technically this indicates singularity -c but it is often caused by improper arguments or improper -c setting of lda . it will not occur if the subroutines are -c called correctly and if zgeco has set rcond .gt. 0.0 -c or zgefa has set info .eq. 0 . -c -c to compute inverse(a) * c where c is a matrix -c with p columns -c call zgeco(a,lda,n,ipvt,rcond,z) -c if (rcond is too small) go to ... -c do 10 j = 1, p -c call zgesl(a,lda,n,ipvt,c(1,j),0) -c 10 continue -c -c linpack. this version dated 08/14/78 . -c cleve moler, university of new mexico, argonne national lab. -c -c subroutines and functions -c -c blas zaxpy,zdotc -c fortran dconjg -c -c internal variables -c - complex*16 zdotc,t - integer k,kb,l,nm1 -c double precision dreal,dimag -c complex*16 zdumr,zdumi -c dreal(zdumr) = zdumr -c dimag(zdumi) = (0.0d0,-1.0d0)*zdumi -c - nm1 = n - 1 - if (job .ne. 0) go to 50 -c -c job = 0 , solve a * x = b -c first solve l*y = b -c - if (nm1 .lt. 1) go to 30 - do 20 k = 1, nm1 - l = ipvt(k) - t = b(l) - if (l .eq. k) go to 10 - b(l) = b(k) - b(k) = t - 10 continue - call zaxpy(n-k,t,a(k+1,k),1,b(k+1),1) - 20 continue - 30 continue -c -c now solve u*x = y -c - do 40 kb = 1, n - k = n + 1 - kb - b(k) = b(k)/a(k,k) - t = -b(k) - call zaxpy(k-1,t,a(1,k),1,b(1),1) - 40 continue - go to 100 - 50 continue -c -c job = nonzero, solve ctrans(a) * x = b -c first solve ctrans(u)*y = b -c - do 60 k = 1, n - t = zdotc(k-1,a(1,k),1,b(1),1) - b(k) = (b(k) - t)/dconjg(a(k,k)) - 60 continue -c -c now solve ctrans(l)*x = y -c - if (nm1 .lt. 1) go to 90 - do 80 kb = 1, nm1 - k = n - kb - b(k) = b(k) + zdotc(n-k,a(k+1,k),1,b(k+1),1) - l = ipvt(k) - if (l .eq. k) go to 70 - t = b(l) - b(l) = b(k) - b(k) = t - 70 continue - 80 continue - 90 continue - 100 continue - return - end
--- a/libcruft/odepack/lsode.f Sun Feb 16 04:29:00 2003 +0000 +++ b/libcruft/odepack/lsode.f Tue Feb 18 20:08:20 2003 +0000 @@ -921,9 +921,9 @@ C VNORM COMPUTES THE WEIGHTED R.M.S. NORM OF A VECTOR. C SRCOM IS A USER-CALLABLE ROUTINE TO SAVE AND RESTORE C THE CONTENTS OF THE INTERNAL COMMON BLOCKS. -C DGEFA AND DGESL ARE ROUTINES FROM LINPACK FOR SOLVING FULL +C DGETRF AND DGETRS ARE ROUTINES FROM LAPACK FOR SOLVING FULL C SYSTEMS OF LINEAR ALGEBRAIC EQUATIONS. -C DGBFA AND DGBSL ARE ROUTINES FROM LINPACK FOR SOLVING BANDED +C DGBTRF AND DGBTRS ARE ROUTINES FROM LAPACK FOR SOLVING BANDED C LINEAR SYSTEMS. C DAXPY, DSCAL, IDAMAX, AND DDOT ARE BASIC LINEAR ALGEBRA MODULES C (BLAS) USED BY THE ABOVE LINPACK ROUTINES.
--- a/libcruft/odepack/prepj.f Sun Feb 16 04:29:00 2003 +0000 +++ b/libcruft/odepack/prepj.f Tue Feb 18 20:08:20 2003 +0000 @@ -29,7 +29,7 @@ C J IS STORED IN WM AND REPLACED BY P. IF MITER .NE. 3, P IS THEN C SUBJECTED TO LU DECOMPOSITION IN PREPARATION FOR LATER SOLUTION C OF LINEAR SYSTEMS WITH P AS COEFFICIENT MATRIX. THIS IS DONE -C BY DGEFA IF MITER = 1 OR 2, AND BY DGBFA IF MITER = 4 OR 5. +C BY DGETRF IF MITER = 1 OR 2, AND BY DGBTRF IF MITER = 4 OR 5. C C IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION C WITH PREPJ USES THE FOLLOWING.. @@ -95,7 +95,7 @@ WM(J) = WM(J) + 1.0D0 250 J = J + NP1 C DO LU DECOMPOSITION ON P. -------------------------------------------- - CALL DGEFA (WM(3), N, N, IWM(21), IER) + CALL DGETRF ( N, N, WM(3), N, IWM(21), IER) IF (IER .NE. 0) IERPJ = 1 RETURN C IF MITER = 3, CONSTRUCT A DIAGONAL APPROXIMATION TO J AND P. --------- @@ -170,7 +170,7 @@ WM(II) = WM(II) + 1.0D0 580 II = II + MEBAND C DO LU DECOMPOSITION OF P. -------------------------------------------- - CALL DGBFA (WM(3), MEBAND, N, ML, MU, IWM(21), IER) + CALL DGBTRF ( N, N, ML, MU, WM(3), MEBAND, IWM(21), IER) IF (IER .NE. 0) IERPJ = 1 RETURN C----------------------- END OF SUBROUTINE PREPJ -----------------------
--- a/libcruft/odepack/solsy.f Sun Feb 16 04:29:00 2003 +0000 +++ b/libcruft/odepack/solsy.f Tue Feb 18 20:08:20 2003 +0000 @@ -18,10 +18,10 @@ C----------------------------------------------------------------------- C THIS ROUTINE MANAGES THE SOLUTION OF THE LINEAR SYSTEM ARISING FROM C A CHORD ITERATION. IT IS CALLED IF MITER .NE. 0. -C IF MITER IS 1 OR 2, IT CALLS DGESL TO ACCOMPLISH THIS. +C IF MITER IS 1 OR 2, IT CALLS DGETRS TO ACCOMPLISH THIS. C IF MITER = 3 IT UPDATES THE COEFFICIENT H*EL0 IN THE DIAGONAL C MATRIX, AND THEN COMPUTES THE SOLUTION. -C IF MITER IS 4 OR 5, IT CALLS DGBSL. +C IF MITER IS 4 OR 5, IT CALLS DGBTRS. C COMMUNICATION WITH SOLSY USES THE FOLLOWING VARIABLES.. C WM = REAL WORK SPACE CONTAINING THE INVERSE DIAGONAL MATRIX IF C MITER = 3 AND THE LU DECOMPOSITION OF THE MATRIX OTHERWISE. @@ -41,7 +41,7 @@ C----------------------------------------------------------------------- IERSL = 0 GO TO (100, 100, 300, 400, 400), MITER - 100 CALL DGESL (WM(3), N, N, IWM(21), X, 0) + 100 CALL DGETRS ( 'N', N, 1, WM(3), N, IWM(21), X, N, INLPCK) RETURN C 300 PHL0 = WM(2) @@ -62,7 +62,8 @@ 400 ML = IWM(1) MU = IWM(2) MEBAND = 2*ML + MU + 1 - CALL DGBSL (WM(3), MEBAND, N, ML, MU, IWM(21), X, 0) + CALL DGBTRS ( 'N', N, ML, MU, 1, WM(3), MEBAND, IWM(21), X, N, + * INLPCK) RETURN C----------------------- END OF SUBROUTINE SOLSY ----------------------- END
--- a/libcruft/odessa/odessa.f Sun Feb 16 04:29:00 2003 +0000 +++ b/libcruft/odessa/odessa.f Tue Feb 18 20:08:20 2003 +0000 @@ -1172,9 +1172,9 @@ C ODESSA_VNORM COMPUTES THE WEIGHTED R.M.S. NORM OF A VECTOR. C ODESSA_SVCOM AND ODESSA_RSCOM ARE USER-CALLABLE ROUTINES TO SAVE AND RESTORE, C RESPECTIVELY, THE CONTENTS OF THE INTERNAL COMMON BLOCKS. -C DGEFA AND DGESL ARE ROUTINES FROM LINPACK FOR SOLVING FULL +C DGETRF AND DGETRS ARE ROUTINES FROM LAPACK FOR SOLVING FULL C SYSTEMS OF LINEAR ALGEBRAIC EQUATIONS. -C DGBFA AND DGBSL ARE ROUTINES FROM LINPACK FOR SOLVING BANDED +C DGBTRF AND DGBTRS ARE ROUTINES FROM LAPACK FOR SOLVING BANDED C LINEAR SYSTEMS. C DAXPY, DSCAL, IDAMAX, AND DDOT ARE BASIC LINEAR ALGEBRA MODULES C (BLAS) USED BY THE ABOVE LINPACK ROUTINES.
--- a/libcruft/odessa/odessa_prepj.f Sun Feb 16 04:29:00 2003 +0000 +++ b/libcruft/odessa/odessa_prepj.f Tue Feb 18 20:08:20 2003 +0000 @@ -20,7 +20,7 @@ C J IS STORED IN WM AND REPLACED BY P. IF MITER .NE. 3, P IS THEN C SUBJECTED TO LU DECOMPOSITION (JOPT = 0) IN PREPARATION FOR LATER C SOLUTION OF LINEAR SYSTEMS WITH P AS COEFFICIENT MATRIX. THIS IS -C DONE BY DGEFA IF MITER = 1 OR 2, AND BY DGBFA IF MITER = 4 OR 5. +C DONE BY DGETRF IF MITER = 1 OR 2, AND BY DGBTRF IF MITER = 4 OR 5. C C IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION C WITH ODESSA_PREPJ USES THE FOLLOWING.. @@ -86,7 +86,7 @@ WM(J) = WM(J) + ONE 250 J = J + (N + 1) C DO LU DECOMPOSITION ON P. -------------------------------------------- - CALL DGEFA (WM(3), N, N, IWM(21), IER) + CALL DGETRF ( N, N, WM(3), N, IWM(21), IER) IF (IER .NE. 0) IERPJ = 1 RETURN C IF MITER = 3, CONSTRUCT A DIAGONAL APPROXIMATION TO J AND P. --------- @@ -159,7 +159,7 @@ WM(II) = WM(II) + ONE 580 II = II + MEBAND C DO LU DECOMPOSITION OF P. -------------------------------------------- - CALL DGBFA (WM(3), MEBAND, N, ML, MU, IWM(21), IER) + CALL DGBTRF ( N, N, ML, MU, WM(3), MEBAND, IWM(21), IER) IF (IER .NE. 0) IERPJ = 1 RETURN C---------------- END OF SUBROUTINE ODESSA_PREPJ -----------------------
--- a/libcruft/odessa/odessa_solsy.f Sun Feb 16 04:29:00 2003 +0000 +++ b/libcruft/odessa/odessa_solsy.f Tue Feb 18 20:08:20 2003 +0000 @@ -10,10 +10,10 @@ C----------------------------------------------------------------------- C THIS ROUTINE MANAGES THE SOLUTION OF THE LINEAR SYSTEM ARISING FROM C A CHORD ITERATION. IT IS CALLED IF MITER .NE. 0. -C IF MITER IS 1 OR 2, IT CALLS DGESL TO ACCOMPLISH THIS. +C IF MITER IS 1 OR 2, IT CALLS DGETRS TO ACCOMPLISH THIS. C IF MITER = 3 IT UPDATES THE COEFFICIENT H*EL0 IN THE DIAGONAL C MATRIX, AND THEN COMPUTES THE SOLUTION. -C IF MITER IS 4 OR 5, IT CALLS DGBSL. +C IF MITER IS 4 OR 5, IT CALLS DGBTRS. C COMMUNICATION WITH ODESSA_SOLSY USES THE FOLLOWING VARIABLES.. C WM = REAL WORK SPACE CONTAINING THE INVERSE DIAGONAL MATRIX IF C MITER = 3 AND THE LU DECOMPOSITION OF THE MATRIX OTHERWISE. @@ -33,7 +33,7 @@ C----------------------------------------------------------------------- IERSL = 0 GO TO (100, 100, 300, 400, 400), MITER - 100 CALL DGESL (WM(3), N, N, IWM(21), X, 0) + 100 CALL DGETRS ( 'N', N, 1, WM(3), N, IWM(21), X, N, INLPCK) RETURN C 300 PHL0 = WM(2) @@ -54,7 +54,8 @@ 400 ML = IWM(1) MU = IWM(2) MEBAND = 2*ML + MU + 1 - CALL DGBSL (WM(3), MEBAND, N, ML, MU, IWM(21), X, 0) + CALL DGBTRS ( 'N', N, ML, MU, 1, WM(3), MEBAND, IWM(21), X, N, + * INLPCK) RETURN C----------------------- END OF SUBROUTINE ODESSA_SOLSY ----------------------- END
--- a/libcruft/ranlib/setgmn.f Sun Feb 16 04:29:00 2003 +0000 +++ b/libcruft/ranlib/setgmn.f Tue Feb 18 20:08:20 2003 +0000 @@ -1,10 +1,10 @@ SUBROUTINE setgmn(meanv,covm,ldcovm,p,parm) C SUBROUTINE setgmn(meanv,covm,p,parm) C JJV changed this routine to take leading dimension of COVM -C JJV argument and pass it to SPOFA, making it easier to use +C JJV argument and pass it to SPOTRF, making it easier to use C JJV if the COVM which is used is contained in a larger matrix -C JJV and to make the routine more consistent with LINPACK. -C JJV Changes are in comments, declarations, and the call to SPOFA. +C JJV and to make the routine more consistent with LAPACK. +C JJV Changes are in comments, declarations, and the call to SPOTRF. C********************************************************************** C C SUBROUTINE SETGMN( MEANV, COVM, LDCOVM, P, PARM) @@ -58,7 +58,7 @@ INTEGER i,icount,info,j C .. C .. External Subroutines .. - EXTERNAL spofa + EXTERNAL spotrf C .. C .. Executable Statements .. C @@ -81,7 +81,8 @@ C Cholesky decomposition to find A s.t. trans(A)*(A) = COVM C C CALL spofa(covm,p,p,info) - CALL spofa(covm,ldcovm,p,info) +C CALL spofa(covm,ldcovm,p,info) + CALL spotrf ( 'Upper', p, covm, ldcovm, info) IF (.NOT. (info.NE.0)) GO TO 30 WRITE (*,*) ' COVM not positive definite in SETGMN' STOP ' COVM not positive definite in SETGMN'
--- a/liboctave/CMatrix.cc Sun Feb 16 04:29:00 2003 +0000 +++ b/liboctave/CMatrix.cc Tue Feb 18 20:08:20 2003 +0000 @@ -79,14 +79,19 @@ const Complex&, Complex*, const int&, long, long); - int F77_FUNC (zgeco, ZGECO) (Complex*, const int&, const int&, int*, - double&, Complex*); - - int F77_FUNC (zgedi, ZGEDI) (Complex*, const int&, const int&, int*, - Complex*, Complex*, const int&); - - int F77_FUNC (zgesl, ZGESL) (Complex*, const int&, const int&, int*, - Complex*, const int&); + int F77_FUNC (zgetrf, ZGETRF) (const int&, const int&, Complex*, const int&, + int*, int&); + + int F77_FUNC (zgetrs, ZGETRS) (const char*, const int&, const int&, + Complex*, const int&, + const int*, Complex*, const int&, int&); + + int F77_FUNC (zgetri, ZGETRI) (const int&, Complex*, const int&, const int*, + Complex*, const int&, int&); + + int F77_FUNC (zgecon, ZGECON) (const char*, const int&, Complex*, + const int&, const double&, double&, + Complex*, double*, int&); int F77_FUNC (zgelss, ZGELSS) (const int&, const int&, const int&, Complex*, const int&, Complex*, @@ -925,18 +930,19 @@ { int info; double rcond; - return inverse (info, rcond); + return inverse (info, rcond, 0, 0); } ComplexMatrix ComplexMatrix::inverse (int& info) const { double rcond; - return inverse (info, rcond); + return inverse (info, rcond, 0, 0); } ComplexMatrix -ComplexMatrix::inverse (int& info, double& rcond, int force) const +ComplexMatrix::inverse (int& info, double& rcond, int force, + int calc_cond) const { ComplexMatrix retval; @@ -947,44 +953,83 @@ (*current_liboctave_error_handler) ("inverse requires square matrix"); else { - info = 0; - Array<int> ipvt (nr); int *pipvt = ipvt.fortran_vec (); - Array<Complex> z (nr); - Complex *pz = z.fortran_vec (); - retval = *this; Complex *tmp_data = retval.fortran_vec (); - F77_XFCN (zgeco, ZGECO, (tmp_data, nr, nc, pipvt, rcond, pz)); + Array<Complex> z(1); + int lwork = 1; + + // Query the optimum work array size + + F77_XFCN (zgetri, ZGETRI, (nc, tmp_data, nr, pipvt, + z.fortran_vec (), lwork, info)); + + if (f77_exception_encountered) + { + (*current_liboctave_error_handler) + ("unrecoverable error in zgetri"); + return retval; + } + + lwork = static_cast<int> (real(z(0))); + lwork = (lwork < 2 *nc ? 2*nc : lwork); + z.resize (lwork); + Complex *pz = z.fortran_vec (); + + info = 0; + + /* Calculate the norm of the matrix, for later use */ + double anorm; + if (calc_cond) + anorm = retval.abs().sum().row(0).max(); + + F77_XFCN (zgetrf, ZGETRF, (nc, nc, tmp_data, nr, pipvt, info)); if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in zgeco"); + (*current_liboctave_error_handler) ("unrecoverable error in zgetrf"); else { - volatile double rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) + /* Throw-away extra info LAPACK gives so as to not change output */ + rcond = 0.; + if ( info != 0) info = -1; + else if (calc_cond) + { + /* Now calculate the condition number for non-singular matrix */ + char job = '1'; + Array<double> rz (2 * nc); + double *prz = rz.fortran_vec (); + F77_XFCN (zgecon, ZGECON, ( &job, nc, tmp_data, nr, anorm, + rcond, pz, prz, info)); + + if (f77_exception_encountered) + (*current_liboctave_error_handler) + ("unrecoverable error in zgecon"); + + if ( info != 0) + info = -1; + } if (info == -1 && ! force) retval = *this; // Restore contents. else { - Complex *dummy = 0; - - F77_XFCN (zgedi, ZGEDI, (tmp_data, nr, nc, pipvt, dummy, - pz, 1)); + F77_XFCN (zgetri, ZGETRI, (nc, tmp_data, nr, pipvt, + pz, lwork, info)); if (f77_exception_encountered) (*current_liboctave_error_handler) - ("unrecoverable error in zgedi"); + ("unrecoverable error in zgetri"); + + if ( info != 0) + info = -1; } } } - + return retval; } @@ -1356,18 +1401,18 @@ { int info; double rcond; - return determinant (info, rcond); + return determinant (info, rcond, 0); } ComplexDET ComplexMatrix::determinant (int& info) const { double rcond; - return determinant (info, rcond); + return determinant (info, rcond, 0); } ComplexDET -ComplexMatrix::determinant (int& info, double& rcond) const +ComplexMatrix::determinant (int& info, double& rcond, int calc_cond) const { ComplexDET retval; @@ -1383,45 +1428,81 @@ } else { - info = 0; - Array<int> ipvt (nr); int *pipvt = ipvt.fortran_vec (); - Array<Complex> z (nr); - Complex *pz = z.fortran_vec (); - ComplexMatrix atmp = *this; Complex *tmp_data = atmp.fortran_vec (); - F77_XFCN (zgeco, ZGECO, (tmp_data, nr, nr, pipvt, rcond, pz)); + info = 0; + + /* Calculate the norm of the matrix, for later use */ + double anorm = 0; + if (calc_cond) + anorm = atmp.abs().sum().row(0).max(); + + F77_XFCN (zgetrf, ZGETRF, (nr, nc, tmp_data, nr, pipvt, info)); if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in zgeco"); + (*current_liboctave_error_handler) ("unrecoverable error in zgetrf"); else { - volatile double rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) + /* Throw-away extra info LAPACK gives so as to not change output */ + rcond = 0.; + if ( info != 0) { info = -1; retval = ComplexDET (); - } - else + } + else { - Complex d[2]; - - F77_XFCN (zgedi, ZGEDI, (tmp_data, nr, nr, pipvt, d, pz, 10)); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dgedi"); - else - retval = ComplexDET (d); + if (calc_cond) + { + /* Now calc the condition number for non-singular matrix */ + char job = '1'; + Array<Complex> z (2*nr); + Complex *pz = z.fortran_vec (); + Array<double> rz (2*nr); + double *prz = rz.fortran_vec (); + + F77_XFCN (zgecon, ZGECON, ( &job, nc, tmp_data, nr, anorm, + rcond, pz, prz, info)); + + if (f77_exception_encountered) + (*current_liboctave_error_handler) + ("unrecoverable error in zgecon"); + } + + if ( info != 0) + { + info = -1; + retval = ComplexDET (); + } + else + { + Complex d[2] = { 1., 0.}; + for (int i=0; i<nc; i++) + { + if (ipvt(i) != (i+1)) d[0] = -d[0]; + d[0] = d[0] * atmp(i,i); + if (d[0] == 0.) break; + while (::abs(d[0]) < 1.) + { + d[0] = 10. * d[0]; + d[1] = d[1] - 1.; + } + while (::abs(d[0]) >= 10.) + { + d[0] = 0.1 * d[0]; + d[1] = d[1] + 1.; + } + } + retval = ComplexDET (d); + } } } } - + return retval; } @@ -1494,55 +1575,82 @@ Array<int> ipvt (nr); int *pipvt = ipvt.fortran_vec (); - Array<Complex> z (nr); - Complex *pz = z.fortran_vec (); - ComplexMatrix atmp = *this; Complex *tmp_data = atmp.fortran_vec (); - F77_XFCN (zgeco, ZGECO, (tmp_data, nr, nr, pipvt, rcond, pz)); + Array<Complex> z (2 * nc); + Complex *pz = z.fortran_vec (); + Array<double> rz (2 * nc); + double *prz = rz.fortran_vec (); + + /* Calculate the norm of the matrix, for later use */ + double anorm = atmp.abs().sum().row(0).max(); + + F77_XFCN (zgetrf, ZGETRF, (nr, nr, tmp_data, nr, pipvt, info)); if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in zgeco"); + (*current_liboctave_error_handler) ("unrecoverable error in zgetrf"); else { - volatile double rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) - { + /* Throw-away extra info LAPACK gives so as to not change output */ + rcond = 0.; + if ( info != 0) + { info = -2; if (sing_handler) sing_handler (rcond); else (*current_liboctave_error_handler) - ("matrix singular to machine precision, rcond = %g", - rcond); - } - else + ("matrix singular to machine precision"); + + } + else { - retval = b; - Complex *result = retval.fortran_vec (); - - int b_nc = b.cols (); - - for (volatile int j = 0; j < b_nc; j++) + /* Now calculate the condition number for non-singular matrix */ + char job = '1'; + F77_XFCN (zgecon, ZGECON, ( &job, nc, tmp_data, nr, anorm, + rcond, pz, prz, info)); + + if (f77_exception_encountered) + (*current_liboctave_error_handler) + ("unrecoverable error in zgecon"); + + if ( info != 0) + info = -2; + + volatile double rcond_plus_one = rcond + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcond)) { - F77_XFCN (zgesl, ZGESL, (tmp_data, nr, nr, pipvt, - &result[nr*j], 0)); + info = -2; + + if (sing_handler) + sing_handler (rcond); + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision, rcond = %g", + rcond); + } + else + { + retval = b; + Complex *result = retval.fortran_vec (); + + int b_nc = b.cols (); + + char job = 'N'; + F77_XFCN (zgetrs, ZGETRS, (&job, nr, b_nc, tmp_data, nr, + pipvt, result, b.rows(), info)); if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("unrecoverable error in dgesl"); - - break; - } + (*current_liboctave_error_handler) + ("unrecoverable error in zgetrs"); } } } } - + return retval; } @@ -1616,23 +1724,27 @@ Array<int> ipvt (nr); int *pipvt = ipvt.fortran_vec (); - Array<Complex> z (nr); - Complex *pz = z.fortran_vec (); - ComplexMatrix atmp = *this; Complex *tmp_data = atmp.fortran_vec (); - F77_XFCN (zgeco, ZGECO, (tmp_data, nr, nr, pipvt, rcond, pz)); + Array<Complex> z (2 * nc); + Complex *pz = z.fortran_vec (); + Array<double> rz (2 * nc); + double *prz = rz.fortran_vec (); + + /* Calculate the norm of the matrix, for later use */ + double anorm = atmp.abs().sum().row(0).max(); + + F77_XFCN (zgetrf, ZGETRF, (nr, nr, tmp_data, nr, pipvt, info)); if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zgeco"); + (*current_liboctave_error_handler) ("unrecoverable error in zgetrf"); else { - volatile double rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) - { + /* Throw-away extra info LAPACK gives so as to not change output */ + rcond = 0.; + if ( info != 0) + { info = -2; if (sing_handler) @@ -1641,21 +1753,51 @@ (*current_liboctave_error_handler) ("matrix singular to machine precision, rcond = %g", rcond); - } - else + } + else { - retval = b; - Complex *result = retval.fortran_vec (); - - F77_XFCN (zgesl, ZGESL, (tmp_data, nr, nr, pipvt, result, 0)); + /* Now calculate the condition number for non-singular matrix */ + char job = '1'; + F77_XFCN (zgecon, ZGECON, ( &job, nc, tmp_data, nr, anorm, + rcond, pz, prz, info)); if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dgesl"); + (*current_liboctave_error_handler) + ("unrecoverable error in zgecon"); + + if ( info != 0) + info = -2; + + volatile double rcond_plus_one = rcond + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcond)) + { + info = -2; + + if (sing_handler) + sing_handler (rcond); + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision, rcond = %g", + rcond); + } + else + { + retval = b; + Complex *result = retval.fortran_vec (); + + char job = 'N'; + F77_XFCN (zgetrs, ZGETRS, (&job, nr, 1, tmp_data, nr, pipvt, + result, b.length(), info)); + + if (f77_exception_encountered) + (*current_liboctave_error_handler) + ("unrecoverable error in zgetrs"); + + } } } } - return retval; } @@ -2488,6 +2630,20 @@ #undef COL_EXPR } +Matrix ComplexMatrix::abs (void) const +{ + int nr = rows (); + int nc = cols (); + + Matrix retval (nr, nc); + + for (int j = 0; j < nc; j++) + for (int i = 0; i < nr; i++) + retval (i, j) = ::abs (elem (i, j)); + + return retval; +} + ComplexColumnVector ComplexMatrix::diag (void) const { @@ -2600,7 +2756,7 @@ bool real_only = row_is_real_only (i); - double abs_min = real_only ? real (tmp_min) : abs (tmp_min); + double abs_min = real_only ? real (tmp_min) : ::abs (tmp_min); if (xisnan (tmp_min)) idx_j = -1; @@ -2610,7 +2766,7 @@ { Complex tmp = elem (i, j); - double abs_tmp = real_only ? real (tmp) : abs (tmp); + double abs_tmp = real_only ? real (tmp) : ::abs (tmp); if (xisnan (tmp)) { @@ -2662,7 +2818,7 @@ bool real_only = row_is_real_only (i); - double abs_max = real_only ? real (tmp_max) : abs (tmp_max); + double abs_max = real_only ? real (tmp_max) : ::abs (tmp_max); if (xisnan (tmp_max)) idx_j = -1; @@ -2672,7 +2828,7 @@ { Complex tmp = elem (i, j); - double abs_tmp = real_only ? real (tmp) : abs (tmp); + double abs_tmp = real_only ? real (tmp) : ::abs (tmp); if (xisnan (tmp)) { @@ -2724,7 +2880,7 @@ bool real_only = column_is_real_only (j); - double abs_min = real_only ? real (tmp_min) : abs (tmp_min); + double abs_min = real_only ? real (tmp_min) : ::abs (tmp_min); if (xisnan (tmp_min)) idx_i = -1; @@ -2734,7 +2890,7 @@ { Complex tmp = elem (i, j); - double abs_tmp = real_only ? real (tmp) : abs (tmp); + double abs_tmp = real_only ? real (tmp) : ::abs (tmp); if (xisnan (tmp)) { @@ -2786,7 +2942,7 @@ bool real_only = column_is_real_only (j); - double abs_max = real_only ? real (tmp_max) : abs (tmp_max); + double abs_max = real_only ? real (tmp_max) : ::abs (tmp_max); if (xisnan (tmp_max)) idx_i = -1; @@ -2796,7 +2952,7 @@ { Complex tmp = elem (i, j); - double abs_tmp = real_only ? real (tmp) : abs (tmp); + double abs_tmp = real_only ? real (tmp) : ::abs (tmp); if (xisnan (tmp)) {
--- a/liboctave/CMatrix.h Sun Feb 16 04:29:00 2003 +0000 +++ b/liboctave/CMatrix.h Tue Feb 18 20:08:20 2003 +0000 @@ -140,7 +140,8 @@ ComplexMatrix inverse (void) const; ComplexMatrix inverse (int& info) const; - ComplexMatrix inverse (int& info, double& rcond, int force = 0) const; + ComplexMatrix inverse (int& info, double& rcond, int force = 0, + int calc_cond = 1) const; ComplexMatrix pseudo_inverse (double tol = 0.0); @@ -152,7 +153,7 @@ ComplexDET determinant (void) const; ComplexDET determinant (int& info) const; - ComplexDET determinant (int& info, double& rcond) const; + ComplexDET determinant (int& info, double& rcond, int calc_cond = 1) const; ComplexMatrix solve (const Matrix& b) const; ComplexMatrix solve (const Matrix& b, int& info) const; @@ -251,6 +252,7 @@ ComplexMatrix prod (int dim = -1) const; ComplexMatrix sum (int dim = -1) const; ComplexMatrix sumsq (int dim = -1) const; + Matrix abs (void) const; ComplexColumnVector diag (void) const; ComplexColumnVector diag (int k) const;
--- a/liboctave/ChangeLog Sun Feb 16 04:29:00 2003 +0000 +++ b/liboctave/ChangeLog Tue Feb 18 20:08:20 2003 +0000 @@ -1,3 +1,23 @@ +2003-02-18 David Bateman <dbateman@free.fr> + + * dMatrix.cc (Matrix::inverse, Matrix::determinant, Matrix::solve): + Use Lapack instead of Linpack. + * Cmatrix.cc (ComplexMatrix::inverse, ComplexMatrix::determinant, + ComplexMatrix::solve): Likewise. + + * dMatrix.cc (Matrix::determinant, Matrix::inverse): New arg, + calc_cond. If 0, skip condition number calculation. + * CMatrix.cc (ComplexMatrix::determinant, ComplexMatrix::inverse): + Likewise. + + * CmplxLU.cc (ComplexLU::ComplexLU): Allow non-square matrices. + * dbleLU.cc (LU::LU): Likewise. + * base-lu.cc (base_lu::L), base_lu::U, base_lu::P): Likewise. + +2002-10-31 John W. Eaton <jwe@bevo.che.wisc.edu> + + * octave.test/arith/prod-4.m, octave.test/arith/sum-4.m: + 2003-02-14 John W. Eaton <jwe@bevo.che.wisc.edu> * Array2-idx.h (Array2<T>::index): Fix thinko.
--- a/liboctave/CmplxLU.cc Sun Feb 16 04:29:00 2003 +0000 +++ b/liboctave/CmplxLU.cc Tue Feb 18 20:08:20 2003 +0000 @@ -51,16 +51,9 @@ { int a_nr = a.rows (); int a_nc = a.cols (); + int mn = (a_nr < a_nc ? a_nr : a_nc); - if (a_nr == 0 || a_nc == 0 || a_nr != a_nc) - { - (*current_liboctave_error_handler) ("ComplexLU requires square matrix"); - return; - } - - int n = a_nr; - - ipvt.resize (n); + ipvt.resize (mn); int *pipvt = ipvt.fortran_vec (); a_fact = a; @@ -68,10 +61,10 @@ int info = 0; - F77_XFCN (zgetrf, ZGETRF, (n, n, tmp_data, n, pipvt, info)); + F77_XFCN (zgetrf, ZGETRF, (a_nr, a_nc, tmp_data, a_nr, pipvt, info)); if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in zgesv"); + (*current_liboctave_error_handler) ("unrecoverable error in zgetrf"); else ipvt -= 1; }
--- a/liboctave/base-lu.cc Sun Feb 16 04:29:00 2003 +0000 +++ b/liboctave/base-lu.cc Tue Feb 18 20:08:20 2003 +0000 @@ -34,14 +34,16 @@ lu_type base_lu <lu_type, lu_elt_type, p_type, p_elt_type> :: L (void) const { - int n = ipvt.length (); + int a_nr = a_fact.rows (); + int a_nc = a_fact.cols (); + int mn = (a_nr < a_nc ? a_nr : a_nc); - lu_type l (n, n, lu_elt_type (0.0)); + lu_type l (a_nr, mn, lu_elt_type (0.0)); - for (int i = 0; i < n; i++) + for (int i = 0; i < a_nr; i++) { l.xelem (i, i) = 1.0; - for (int j = 0; j < i; j++) + for (int j = 0; j < (i < a_nc ? i : a_nc); j++) l.xelem (i, j) = a_fact.xelem (i, j); } @@ -52,13 +54,15 @@ lu_type base_lu <lu_type, lu_elt_type, p_type, p_elt_type> :: U (void) const { - int n = ipvt.length (); - - lu_type u (n, n, lu_elt_type (0.0)); + int a_nr = a_fact.rows (); + int a_nc = a_fact.cols (); + int mn = (a_nr < a_nc ? a_nr : a_nc); - for (int i = 0; i < n; i++) + lu_type u (mn, a_nc, lu_elt_type (0.0)); + + for (int i = 0; i < mn; i++) { - for (int j = i; j < n; j++) + for (int j = i; j < a_nc; j++) u.xelem (i, j) = a_fact.xelem (i, j); } @@ -69,14 +73,14 @@ p_type base_lu <lu_type, lu_elt_type, p_type, p_elt_type> :: P (void) const { - int n = ipvt.length (); + int a_nr = a_fact.rows (); - Array<int> pvt (n); + Array<int> pvt (a_nr); - for (int i = 0; i < n; i++) + for (int i = 0; i < a_nr; i++) pvt.xelem (i) = i; - for (int i = 0; i < n - 1; i++) + for (int i = 0; i < ipvt.length(); i++) { int k = ipvt.xelem (i); @@ -88,9 +92,9 @@ } } - p_type p (n, n, p_elt_type (0.0)); + p_type p (a_nr, a_nr, p_elt_type (0.0)); - for (int i = 0; i < n; i++) + for (int i = 0; i < a_nr; i++) p.xelem (i, pvt.xelem (i)) = 1.0; return p;
--- a/liboctave/dMatrix.cc Sun Feb 16 04:29:00 2003 +0000 +++ b/liboctave/dMatrix.cc Tue Feb 18 20:08:20 2003 +0000 @@ -75,15 +75,19 @@ const double&, double*, const int&, long, long); - int F77_FUNC (dgeco, DGECO) (double*, const int&, const int&, int*, - double&, double*); - - int F77_FUNC (dgesl, DGESL) (const double*, const int&, const int&, - const int*, double*, const int&); - - int F77_FUNC (dgedi, DGEDI) (double*, const int&, const int&, - const int*, double*, double*, - const int&); + int F77_FUNC (dgetrf, DGETRF) (const int&, const int&, double*, const int&, + int*, int&); + + int F77_FUNC (dgetrs, DGETRS) (const char*, const int&, const int&, + const double*, const int&, + const int*, double*, const int&, int&); + + int F77_FUNC (dgetri, DGETRI) (const int&, double*, const int&, const int*, + double*, const int&, int&); + + int F77_FUNC (dgecon, DGECON) (const char*, const int&, double*, + const int&, const double&, double&, + double*, int*, int&); int F77_FUNC (dgelss, DGELSS) (const int&, const int&, const int&, double*, const int&, double*, @@ -595,18 +599,18 @@ { int info; double rcond; - return inverse (info, rcond); + return inverse (info, rcond, 0, 0); } Matrix Matrix::inverse (int& info) const { double rcond; - return inverse (info, rcond); + return inverse (info, rcond, 0, 0); } Matrix -Matrix::inverse (int& info, double& rcond, int force) const +Matrix::inverse (int& info, double& rcond, int force, int calc_cond) const { Matrix retval; @@ -617,40 +621,78 @@ (*current_liboctave_error_handler) ("inverse requires square matrix"); else { - info = 0; - Array<int> ipvt (nr); int *pipvt = ipvt.fortran_vec (); - Array<double> z (nr); - double *pz = z.fortran_vec (); - retval = *this; double *tmp_data = retval.fortran_vec (); - F77_XFCN (dgeco, DGECO, (tmp_data, nr, nc, pipvt, rcond, pz)); + Array<double> z(1); + int lwork = -1; + + // Query the optimum work array size + F77_XFCN (dgetri, DGETRI, (nc, tmp_data, nr, pipvt, + z.fortran_vec (), lwork, info)); + + if (f77_exception_encountered) + { + (*current_liboctave_error_handler) + ("unrecoverable error in dgetri"); + return retval; + } + + lwork = static_cast<int> (z(0)); + lwork = (lwork < 2 *nc ? 2*nc : lwork); + z.resize (lwork); + double *pz = z.fortran_vec (); + + info = 0; + + /* Calculate the norm of the matrix, for later use */ + double anorm = 0; + if (calc_cond) + anorm = retval.abs().sum().row(0).max(); + + F77_XFCN (dgetrf, DGETRF, (nc, nc, tmp_data, nr, pipvt, info)); if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in dgeco"); + (*current_liboctave_error_handler) ("unrecoverable error in dgetrf"); else { - volatile double rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) + /* Throw-away extra info LAPACK gives so as to not change output */ + rcond = 0.; + if ( info != 0) info = -1; + else if (calc_cond) + { + /* Now calculate the condition number for non-singular matrix */ + char job = '1'; + Array<int> iz (nc); + int *piz = iz.fortran_vec (); + F77_XFCN (dgecon, DGECON, ( &job, nc, tmp_data, nr, anorm, + rcond, pz, piz, info)); + + if (f77_exception_encountered) + (*current_liboctave_error_handler) + ("unrecoverable error in dgecon"); + + if ( info != 0) + info = -1; + } if (info == -1 && ! force) retval = *this; // Restore matrix contents. else { - double *dummy = 0; - - F77_XFCN (dgedi, DGEDI, (tmp_data, nr, nc, pipvt, dummy, - pz, 1)); + F77_XFCN (dgetri, DGETRI, (nc, tmp_data, nr, pipvt, + pz, lwork, info)); if (f77_exception_encountered) (*current_liboctave_error_handler) - ("unrecoverable error in dgedi"); + ("unrecoverable error in dgetri"); + + if ( info != 0) + info = -1; } } } @@ -1024,18 +1066,18 @@ { int info; double rcond; - return determinant (info, rcond); + return determinant (info, rcond, 0); } DET Matrix::determinant (int& info) const { double rcond; - return determinant (info, rcond); + return determinant (info, rcond, 0); } DET -Matrix::determinant (int& info, double& rcond) const +Matrix::determinant (int& info, double& rcond, int calc_cond) const { DET retval; @@ -1051,41 +1093,77 @@ } else { - info = 0; - Array<int> ipvt (nr); int *pipvt = ipvt.fortran_vec (); - Array<double> z (nr); - double *pz = z.fortran_vec (); - Matrix atmp = *this; double *tmp_data = atmp.fortran_vec (); - F77_XFCN (dgeco, DGECO, (tmp_data, nr, nr, pipvt, rcond, pz)); + info = 0; + + /* Calculate the norm of the matrix, for later use */ + double anorm = 0; + if (calc_cond) + anorm = atmp.abs().sum().row(0).max(); + + F77_XFCN (dgetrf, DGETRF, (nr, nr, tmp_data, nr, pipvt, info)); if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in dgeco"); + (*current_liboctave_error_handler) ("unrecoverable error in dgetrf"); else { - volatile double rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) + /* Throw-away extra info LAPACK gives so as to not change output */ + rcond = 0.; + if ( info != 0) { - info = -1; - retval = DET (); - } - else + info = -1; + retval = DET (); + } + else { - double d[2]; - - F77_XFCN (dgedi, DGEDI, (tmp_data, nr, nr, pipvt, d, pz, 10)); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dgedi"); - else - retval = DET (d); + if (calc_cond) + { + /* Now calc the condition number for non-singular matrix */ + char job = '1'; + Array<double> z (4 * nc); + double *pz = z.fortran_vec (); + Array<int> iz (nc); + int *piz = iz.fortran_vec (); + + F77_XFCN (dgecon, DGECON, ( &job, nc, tmp_data, nr, anorm, + rcond, pz, piz, info)); + + if (f77_exception_encountered) + (*current_liboctave_error_handler) + ("unrecoverable error in dgecon"); + } + + if ( info != 0) + { + info = -1; + retval = DET (); + } + else + { + double d[2] = { 1., 0.}; + for (int i=0; i<nc; i++) + { + if (ipvt(i) != (i+1)) d[0] = -d[0]; + d[0] *= atmp(i,i); + if (d[0] == 0.) break; + while (fabs(d[0]) < 1.) + { + d[0] = 10. * d[0]; + d[1] = d[1] - 1.; + } + while (fabs(d[0]) >= 10.) + { + d[0] = 0.1 * d[0]; + d[1] = d[1] + 1.; + } + } + retval = DET (d); + } } } } @@ -1098,14 +1176,14 @@ { int info; double rcond; - return solve (b, info, rcond); + return solve (b, info, rcond, 0); } Matrix Matrix::solve (const Matrix& b, int& info) const { double rcond; - return solve (b, info, rcond); + return solve (b, info, rcond, 0); } Matrix @@ -1133,21 +1211,26 @@ Array<int> ipvt (nr); int *pipvt = ipvt.fortran_vec (); - Array<double> z (nr); - double *pz = z.fortran_vec (); - Matrix atmp = *this; double *tmp_data = atmp.fortran_vec (); - F77_XFCN (dgeco, DGECO, (tmp_data, nr, nr, pipvt, rcond, pz)); + Array<double> z (4 * nc); + double *pz = z.fortran_vec (); + Array<int> iz (nc); + int *piz = iz.fortran_vec (); + + /* Calculate the norm of the matrix, for later use */ + double anorm = atmp.abs().sum().row(0).max(); + + F77_XFCN (dgetrf, DGETRF, (nr, nr, tmp_data, nr, pipvt, info)); if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in dgeco"); + (*current_liboctave_error_handler) ("unrecoverable error in dgetrf"); else { - volatile double rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) + /* Throw-away extra info LAPACK gives so as to not change output */ + rcond = 0.; + if ( info != 0) { info = -2; @@ -1155,28 +1238,50 @@ sing_handler (rcond); else (*current_liboctave_error_handler) - ("matrix singular to machine precision, rcond = %g", - rcond); - } - else + ("matrix singular to machine precision"); + + } + else { - retval = b; - double *result = retval.fortran_vec (); - - int b_nc = b.cols (); - - for (volatile int j = 0; j < b_nc; j++) + /* Now calculate the condition number for non-singular matrix */ + char job = '1'; + F77_XFCN (dgecon, DGECON, ( &job, nc, tmp_data, nr, anorm, + rcond, pz, piz, info)); + + if (f77_exception_encountered) + (*current_liboctave_error_handler) + ("unrecoverable error in dgecon"); + + if ( info != 0) + info = -2; + + volatile double rcond_plus_one = rcond + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcond)) { - F77_XFCN (dgesl, DGESL, (tmp_data, nr, nr, pipvt, - &result[nr*j], 0)); - + info = -2; + + if (sing_handler) + sing_handler (rcond); + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision, rcond = %g", + rcond); + } + else + { + retval = b; + double *result = retval.fortran_vec (); + + int b_nc = b.cols (); + + char job = 'N'; + F77_XFCN (dgetrs, DGETRS, (&job, nr, b_nc, tmp_data, nr, + pipvt, result, b.rows(), info)); + if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("unrecoverable error in dgesl"); - - break; - } + (*current_liboctave_error_handler) + ("unrecoverable error in dgetrs"); } } } @@ -1253,22 +1358,26 @@ Array<int> ipvt (nr); int *pipvt = ipvt.fortran_vec (); - Array<double> z (nr); - double *pz = z.fortran_vec (); - Matrix atmp = *this; double *tmp_data = atmp.fortran_vec (); - F77_XFCN (dgeco, DGECO, (tmp_data, nr, nr, pipvt, rcond, pz)); + Array<double> z (4 * nc); + double *pz = z.fortran_vec (); + Array<int> iz (nc); + int *piz = iz.fortran_vec (); + + /* Calculate the norm of the matrix, for later use */ + double anorm = atmp.abs().sum().row(0).max(); + + F77_XFCN (dgetrf, DGETRF, (nr, nr, tmp_data, nr, pipvt, info)); if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dgeco"); + (*current_liboctave_error_handler) ("unrecoverable error in dgetrf"); else { - volatile double rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) + /* Throw-away extra info LAPACK gives so as to not change output */ + rcond = 0.; + if ( info > 0) { info = -2; @@ -1276,23 +1385,53 @@ sing_handler (rcond); else (*current_liboctave_error_handler) - ("matrix singular to machine precision, rcond = %g", - rcond); - } - else + ("matrix singular to machine precision"); + + } + else { - retval = b; - double *result = retval.fortran_vec (); - - F77_XFCN (dgesl, DGESL, (tmp_data, nr, nr, pipvt, result, 0)); - + /* Now calculate the condition number for non-singular matrix */ + char job = '1'; + F77_XFCN (dgecon, DGECON, ( &job, nc, tmp_data, nr, anorm, + rcond, pz, piz, info)); + if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dgesl"); + (*current_liboctave_error_handler) + ("unrecoverable error in dgecon"); + + if ( info != 0) + info = -2; + + volatile double rcond_plus_one = rcond + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcond)) + { + info = -2; + + if (sing_handler) + sing_handler (rcond); + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision, rcond = %g", + rcond); + } + else + { + retval = b; + double *result = retval.fortran_vec (); + + char job = 'N'; + F77_XFCN (dgetrs, DGETRS, (&job, nr, 1, tmp_data, nr, pipvt, + result, b.length(), info)); + + if (f77_exception_encountered) + (*current_liboctave_error_handler) + ("unrecoverable error in dgetrs"); + } } } } - + return retval; }
--- a/liboctave/dMatrix.h Sun Feb 16 04:29:00 2003 +0000 +++ b/liboctave/dMatrix.h Tue Feb 18 20:08:20 2003 +0000 @@ -115,7 +115,8 @@ Matrix inverse (void) const; Matrix inverse (int& info) const; - Matrix inverse (int& info, double& rcond, int force = 0) const; + Matrix inverse (int& info, double& rcond, int force = 0, + int calc_cond = 1) const; Matrix pseudo_inverse (double tol = 0.0); @@ -127,7 +128,7 @@ DET determinant (void) const; DET determinant (int& info) const; - DET determinant (int& info, double& rcond) const; + DET determinant (int& info, double& rcond, int calc_cond = 1) const; Matrix solve (const Matrix& b) const; Matrix solve (const Matrix& b, int& info) const;
--- a/liboctave/dbleLU.cc Sun Feb 16 04:29:00 2003 +0000 +++ b/liboctave/dbleLU.cc Tue Feb 18 20:08:20 2003 +0000 @@ -51,16 +51,9 @@ { int a_nr = a.rows (); int a_nc = a.cols (); + int mn = (a_nr < a_nc ? a_nr : a_nc); - if (a_nr == 0 || a_nc == 0 || a_nr != a_nc) - { - (*current_liboctave_error_handler) ("LU requires square matrix"); - return; - } - - int n = a_nr; - - ipvt.resize (n); + ipvt.resize (mn); int *pipvt = ipvt.fortran_vec (); a_fact = a; @@ -68,10 +61,10 @@ int info = 0; - F77_XFCN (dgetrf, DGETRF, (n, n, tmp_data, n, pipvt, info)); + F77_XFCN (dgetrf, DGETRF, (a_nr, a_nc, tmp_data, a_nr, pipvt, info)); if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in dgesv"); + (*current_liboctave_error_handler) ("unrecoverable error in dgetrf"); else ipvt -= 1; }
--- a/scripts/ChangeLog Sun Feb 16 04:29:00 2003 +0000 +++ b/scripts/ChangeLog Tue Feb 18 20:08:20 2003 +0000 @@ -1,3 +1,7 @@ +2003-02-18 David Bateman <dbateman@free.fr> + + * mkpkgadd: Scan C++ files as well + 2003-02-13 Schloegl Alois <alois.schloegl@tugraz.at> * strings/findstr.m: Return empty set for zero-length target.
--- a/scripts/mkpkgadd Sun Feb 16 04:29:00 2003 +0000 +++ b/scripts/mkpkgadd Tue Feb 18 20:08:20 2003 +0000 @@ -9,8 +9,14 @@ cd $dir -files=`ls *.m` +m_files=`ls *.m` +cxx_files=`ls *.cc` -if [ -n "$files" ]; then - sed -n 's/^[#%][#%]* *PKG_ADD: *//p' $files +if [ -n "$m_files" ]; then + sed -n 's/^[#%][#%]* *PKG_ADD: *//p' $m_files fi + +if [ -n "$cxx_files" ]; then + sed -n -e 's,^//* *PKG_ADD: *,,p' \ + -e 's,^/\** *PKG_ADD: *\(.*\) \*/$,\1,p' $cxx_files +fi
--- a/src/ChangeLog Sun Feb 16 04:29:00 2003 +0000 +++ b/src/ChangeLog Tue Feb 18 20:08:20 2003 +0000 @@ -1,3 +1,13 @@ +2003-02-18 David Bateman <dbateman@free.fr> + + * DLD-FUNCTIONS/lu.cc (Flu): Allow non-square matrices. + +2003-02-17 John W. Eaton <jwe@bevo.che.wisc.edu> + + * load-save.cc (read_binary_file_header, do_load, do_save, + write_header): No longer static. + (load_save_format): Move enum decl to load-save.h. + 2003-02-15 John W. Eaton <jwe@bevo.che.wisc.edu> * oct-stdstrm.h, oct-stdstrm.cc (octave_base_stdiostream,
--- a/src/DLD-FUNCTIONS/det.cc Sun Feb 16 04:29:00 2003 +0000 +++ b/src/DLD-FUNCTIONS/det.cc Tue Feb 18 20:08:20 2003 +0000 @@ -36,7 +36,7 @@ DEFUN_DLD (det, args, nargout, "-*- texinfo -*-\n\ @deftypefn {Loadable Function} {[@var{d}, @var{rcond}] = } det (@var{a})\n\ -Compute the determinant of @var{a} using @sc{Linpack}. Return an estimate\n\ +Compute the determinant of @var{a} using @sc{Lapack}. Return an estimate\n\ of the reciprocal condition number if requested.\n\ @end deftypefn") {
--- a/src/DLD-FUNCTIONS/lu.cc Sun Feb 16 04:29:00 2003 +0000 +++ b/src/DLD-FUNCTIONS/lu.cc Tue Feb 18 20:08:20 2003 +0000 @@ -65,6 +65,8 @@ 0 1\n\ 1 0\n\ @end example\n\ +\n\ +The matrix is not required to be square..\n\ @end deftypefn") { octave_value_list retval; @@ -89,12 +91,6 @@ else if (arg_is_empty > 0) return octave_value_list (3, Matrix ()); - if (nr != nc) - { - gripe_square_matrix_required ("lu"); - return retval; - } - if (arg.is_real_type ()) { Matrix m = arg.matrix_value ();
--- a/src/load-save.cc Sun Feb 16 04:29:00 2003 +0000 +++ b/src/load-save.cc Tue Feb 18 20:08:20 2003 +0000 @@ -94,20 +94,7 @@ #define OCT_RBV DBL_MAX / 100.0 #endif -enum load_save_format - { - LS_ASCII, - LS_BINARY, - LS_MAT_ASCII, - LS_MAT_BINARY, - LS_MAT5_BINARY, -#ifdef HAVE_HDF5 - LS_HDF5, -#endif /* HAVE_HDF5 */ - LS_UNKNOWN - }; - - enum arrayclasstype +enum arrayclasstype { mxCELL_CLASS=1, // cell array mxSTRUCT_CLASS, // structure @@ -2835,10 +2822,9 @@ return false; } -static int +int read_binary_file_header (std::istream& is, bool& swap, - oct_mach_info::float_format& flt_fmt, - bool quiet = false) + oct_mach_info::float_format& flt_fmt, bool quiet) { const int magic_len = 10; char magic[magic_len+1]; @@ -2953,7 +2939,7 @@ return retval; } -static octave_value +octave_value do_load (std::istream& stream, const std::string& orig_fname, bool force, load_save_format format, oct_mach_info::float_format flt_fmt, bool list_only, bool swap, bool verbose, bool import, @@ -4636,7 +4622,7 @@ // Save the info from sr on stream os in the format specified by fmt. -static void +void do_save (std::ostream& os, symbol_record *sr, load_save_format fmt, int save_as_floats, bool& infnan_warned) { @@ -4750,7 +4736,7 @@ return retval; } -static void +void write_header (std::ostream& os, load_save_format format) { switch (format)
--- a/src/load-save.h Sun Feb 16 04:29:00 2003 +0000 +++ b/src/load-save.h Tue Feb 18 20:08:20 2003 +0000 @@ -29,14 +29,47 @@ class octave_value; +enum load_save_format + { + LS_ASCII, + LS_BINARY, + LS_MAT_ASCII, + LS_MAT_BINARY, + LS_MAT5_BINARY, +#ifdef HAVE_HDF5 + LS_HDF5, +#endif /* HAVE_HDF5 */ + LS_UNKNOWN + }; + extern bool save_ascii_data_for_plotting (std::ostream& os, const octave_value& t, const std::string& name = std::string ()); -extern bool save_three_d (std::ostream& os, const octave_value& t, - bool parametric = false); +extern bool +save_three_d (std::ostream& os, const octave_value& t, + bool parametric = false); + +extern void +save_user_variables (void); + +extern int +read_binary_file_header (std::istream& is, bool& swap, + oct_mach_info::float_format& flt_fmt, + bool quiet = false); -extern void save_user_variables (void); +extern octave_value +do_load (std::istream& stream, const std::string& orig_fname, bool force, + load_save_format format, oct_mach_info::float_format flt_fmt, + bool list_only, bool swap, bool verbose, bool import, + const string_vector& argv, int argv_idx, int argc, int nargout); + +extern void +do_save (std::ostream& os, symbol_record *sr, load_save_format fmt, + int save_as_floats, bool& infnan_warned); + +extern void +write_header (std::ostream& os, load_save_format format); #endif
--- a/src/variables.cc Sun Feb 16 04:29:00 2003 +0000 +++ b/src/variables.cc Tue Feb 18 20:08:20 2003 +0000 @@ -833,7 +833,7 @@ { symbol_record *sr = fbi_sym_tab->lookup (name); - // It is a prorgramming error to look for builtins that aren't. + // It is a programming error to look for builtins that aren't. // Use != here to avoid possible conversion to int of smaller type // than the sr pointer. @@ -860,7 +860,7 @@ int status = 0; symbol_record *sr = fbi_sym_tab->lookup (name); - // It is a prorgramming error to look for builtins that aren't. + // It is a programming error to look for builtins that aren't. // Use != here to avoid possible conversion to int of smaller type // than the sr pointer. @@ -885,7 +885,7 @@ { symbol_record *sr = fbi_sym_tab->lookup (name); - // It is a prorgramming error to look for builtins that aren't. + // It is a programming error to look for builtins that aren't. // Use != here to avoid possible conversion to int of smaller type // than the sr pointer.
--- a/test/octave.test/linalg/linalg.exp Sun Feb 16 04:29:00 2003 +0000 +++ b/test/octave.test/linalg/linalg.exp Tue Feb 18 20:08:20 2003 +0000 @@ -159,7 +159,7 @@ do_test lu-5.m set test lu-6 -set prog_output "error:.*" +set prog_output "ans = 1" do_test lu-6.m set test qr-1
--- a/test/octave.test/linalg/lu-6.m Sun Feb 16 04:29:00 2003 +0000 +++ b/test/octave.test/linalg/lu-6.m Tue Feb 18 20:08:20 2003 +0000 @@ -1,1 +1,4 @@ -lu ([1, 2; 3, 4; 5, 6]) +[l u p] = lu ([1, 2; 3, 4; 5, 6]); +(abs (l - [1, 0; 1/5, 1; 3/5, 1/2]) < sqrt (eps) + && abs (u - [5, 6; 0, 4/5]) < sqrt (eps) + && abs (p - [0, 0, 1; 1, 0, 0; 0 1 0]) < sqrt (eps))