Mercurial > octave-nkf
changeset 7053:570a382ce556
[project @ 2007-10-23 23:17:36 by jwe]
| author | jwe |
|---|---|
| date | Tue, 23 Oct 2007 23:17:36 +0000 |
| parents | ee70ac66041f |
| children | 1c5b2b5b58ba |
| files | libcruft/ChangeLog libcruft/lapack/Makefile.in libcruft/lapack/dgtts2.f libcruft/lapack/zgtts2.f |
| diffstat | 4 files changed, 474 insertions(+), 2 deletions(-) [+] |
line wrap: on
line diff
--- a/libcruft/ChangeLog Tue Oct 23 17:46:48 2007 +0000 +++ b/libcruft/ChangeLog Tue Oct 23 23:17:36 2007 +0000 @@ -1,3 +1,8 @@ +2007-10-23 John W. Eaton <jwe@octave.org> + + * lapack/dgtts2.f, lapack/zgtts2.f: New files. + * lapack/Makefile.in (FSRC): Add them to the list. + 2007-10-16 John W. Eaton <jwe@octave.org> * lapack/dlacn2.f, lapack/dlacn2.f, lapack/dlahr2.f,
--- a/libcruft/lapack/Makefile.in Tue Oct 23 17:46:48 2007 +0000 +++ b/libcruft/lapack/Makefile.in Tue Oct 23 23:17:36 2007 +0000 @@ -30,7 +30,7 @@ dgebd2.f dgebrd.f dgecon.f dgeesx.f dgeev.f dgehd2.f dgehrd.f \ dgelq2.f dgelqf.f dgelss.f dgelsy.f dgeqp3.f dgeqpf.f dgeqr2.f \ dgeqrf.f dgesv.f dgesvd.f dgetf2.f dgetrf.f dgetri.f dgetrs.f \ - dggbak.f dggbal.f dgghrd.f dgtsv.f dgttrf.f dgttrs.f dhgeqz.f \ + dggbak.f dggbal.f dgghrd.f dgtsv.f dgttrf.f dgttrs.f dgtts2.f dhgeqz.f \ dhseqr.f dlabad.f dlabrd.f dlacn2.f dlacon.f dlacpy.f dladiv.f \ dlae2.f dlaev2.f dlaexc.f dlag2.f dlahqr.f dlahr2.f dlahrd.f \ dlaic1.f dlaln2.f dlamc1.f dlamc2.f dlamc3.f dlamc4.f dlamc5.f \ @@ -53,7 +53,7 @@ zgeesx.f zgeev.f zgehd2.f zgehrd.f zgelq2.f zgelqf.f zgelss.f \ zgelsy.f zgeqp3.f zgeqpf.f zgeqr2.f zgeqrf.f zgesv.f zgesvd.f \ zgetf2.f zgetrf.f zgetri.f zgetrs.f zggbal.f zgtsv.f zgttrf.f \ - zgttrs.f zheev.f zhetd2.f zhetrd.f zhseqr.f zlabrd.f zlacgv.f \ + zgttrs.f zgtts2.f zheev.f zhetd2.f zhetrd.f zhseqr.f zlabrd.f zlacgv.f \ zlacn2.f zlacon.f zlacpy.f zladiv.f zlahqr.f zlahr2.f zlahrd.f \ zlaic1.f zlange.f zlanhe.f zlanhs.f zlantr.f zlaqp2.f zlaqps.f \ zlaqr0.f zlaqr1.f zlaqr2.f zlaqr3.f zlaqr4.f zlaqr5.f zlarf.f \
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/lapack/dgtts2.f Tue Oct 23 23:17:36 2007 +0000 @@ -0,0 +1,196 @@ + SUBROUTINE DGTTS2( ITRANS, N, NRHS, DL, D, DU, DU2, IPIV, B, LDB ) +* +* -- LAPACK auxiliary routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* .. Scalar Arguments .. + INTEGER ITRANS, LDB, N, NRHS +* .. +* .. Array Arguments .. + INTEGER IPIV( * ) + DOUBLE PRECISION B( LDB, * ), D( * ), DL( * ), DU( * ), DU2( * ) +* .. +* +* Purpose +* ======= +* +* DGTTS2 solves one of the systems of equations +* A*X = B or A'*X = B, +* with a tridiagonal matrix A using the LU factorization computed +* by DGTTRF. +* +* Arguments +* ========= +* +* ITRANS (input) INTEGER +* Specifies the form of the system of equations. +* = 0: A * X = B (No transpose) +* = 1: A'* X = B (Transpose) +* = 2: A'* X = B (Conjugate transpose = Transpose) +* +* N (input) INTEGER +* The order of the matrix A. +* +* NRHS (input) INTEGER +* The number of right hand sides, i.e., the number of columns +* of the matrix B. NRHS >= 0. +* +* DL (input) DOUBLE PRECISION array, dimension (N-1) +* The (n-1) multipliers that define the matrix L from the +* LU factorization of A. +* +* D (input) DOUBLE PRECISION array, dimension (N) +* The n diagonal elements of the upper triangular matrix U from +* the LU factorization of A. +* +* DU (input) DOUBLE PRECISION array, dimension (N-1) +* The (n-1) elements of the first super-diagonal of U. +* +* DU2 (input) DOUBLE PRECISION array, dimension (N-2) +* The (n-2) elements of the second super-diagonal of U. +* +* IPIV (input) INTEGER array, dimension (N) +* The pivot indices; for 1 <= i <= n, row i of the matrix was +* interchanged with row IPIV(i). IPIV(i) will always be either +* i or i+1; IPIV(i) = i indicates a row interchange was not +* required. +* +* B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) +* On entry, the matrix of right hand side vectors B. +* On exit, B is overwritten by the solution vectors X. +* +* LDB (input) INTEGER +* The leading dimension of the array B. LDB >= max(1,N). +* +* ===================================================================== +* +* .. Local Scalars .. + INTEGER I, IP, J + DOUBLE PRECISION TEMP +* .. +* .. Executable Statements .. +* +* Quick return if possible +* + IF( N.EQ.0 .OR. NRHS.EQ.0 ) + $ RETURN +* + IF( ITRANS.EQ.0 ) THEN +* +* Solve A*X = B using the LU factorization of A, +* overwriting each right hand side vector with its solution. +* + IF( NRHS.LE.1 ) THEN + J = 1 + 10 CONTINUE +* +* Solve L*x = b. +* + DO 20 I = 1, N - 1 + IP = IPIV( I ) + TEMP = B( I+1-IP+I, J ) - DL( I )*B( IP, J ) + B( I, J ) = B( IP, J ) + B( I+1, J ) = TEMP + 20 CONTINUE +* +* Solve U*x = b. +* + B( N, J ) = B( N, J ) / D( N ) + IF( N.GT.1 ) + $ B( N-1, J ) = ( B( N-1, J )-DU( N-1 )*B( N, J ) ) / + $ D( N-1 ) + DO 30 I = N - 2, 1, -1 + B( I, J ) = ( B( I, J )-DU( I )*B( I+1, J )-DU2( I )* + $ B( I+2, J ) ) / D( I ) + 30 CONTINUE + IF( J.LT.NRHS ) THEN + J = J + 1 + GO TO 10 + END IF + ELSE + DO 60 J = 1, NRHS +* +* Solve L*x = b. +* + DO 40 I = 1, N - 1 + IF( IPIV( I ).EQ.I ) THEN + B( I+1, J ) = B( I+1, J ) - DL( I )*B( I, J ) + ELSE + TEMP = B( I, J ) + B( I, J ) = B( I+1, J ) + B( I+1, J ) = TEMP - DL( I )*B( I, J ) + END IF + 40 CONTINUE +* +* Solve U*x = b. +* + B( N, J ) = B( N, J ) / D( N ) + IF( N.GT.1 ) + $ B( N-1, J ) = ( B( N-1, J )-DU( N-1 )*B( N, J ) ) / + $ D( N-1 ) + DO 50 I = N - 2, 1, -1 + B( I, J ) = ( B( I, J )-DU( I )*B( I+1, J )-DU2( I )* + $ B( I+2, J ) ) / D( I ) + 50 CONTINUE + 60 CONTINUE + END IF + ELSE +* +* Solve A' * X = B. +* + IF( NRHS.LE.1 ) THEN +* +* Solve U'*x = b. +* + J = 1 + 70 CONTINUE + B( 1, J ) = B( 1, J ) / D( 1 ) + IF( N.GT.1 ) + $ B( 2, J ) = ( B( 2, J )-DU( 1 )*B( 1, J ) ) / D( 2 ) + DO 80 I = 3, N + B( I, J ) = ( B( I, J )-DU( I-1 )*B( I-1, J )-DU2( I-2 )* + $ B( I-2, J ) ) / D( I ) + 80 CONTINUE +* +* Solve L'*x = b. +* + DO 90 I = N - 1, 1, -1 + IP = IPIV( I ) + TEMP = B( I, J ) - DL( I )*B( I+1, J ) + B( I, J ) = B( IP, J ) + B( IP, J ) = TEMP + 90 CONTINUE + IF( J.LT.NRHS ) THEN + J = J + 1 + GO TO 70 + END IF +* + ELSE + DO 120 J = 1, NRHS +* +* Solve U'*x = b. +* + B( 1, J ) = B( 1, J ) / D( 1 ) + IF( N.GT.1 ) + $ B( 2, J ) = ( B( 2, J )-DU( 1 )*B( 1, J ) ) / D( 2 ) + DO 100 I = 3, N + B( I, J ) = ( B( I, J )-DU( I-1 )*B( I-1, J )- + $ DU2( I-2 )*B( I-2, J ) ) / D( I ) + 100 CONTINUE + DO 110 I = N - 1, 1, -1 + IF( IPIV( I ).EQ.I ) THEN + B( I, J ) = B( I, J ) - DL( I )*B( I+1, J ) + ELSE + TEMP = B( I+1, J ) + B( I+1, J ) = B( I, J ) - DL( I )*TEMP + B( I, J ) = TEMP + END IF + 110 CONTINUE + 120 CONTINUE + END IF + END IF +* +* End of DGTTS2 +* + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/lapack/zgtts2.f Tue Oct 23 23:17:36 2007 +0000 @@ -0,0 +1,271 @@ + SUBROUTINE ZGTTS2( ITRANS, N, NRHS, DL, D, DU, DU2, IPIV, B, LDB ) +* +* -- LAPACK auxiliary routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* .. Scalar Arguments .. + INTEGER ITRANS, LDB, N, NRHS +* .. +* .. Array Arguments .. + INTEGER IPIV( * ) + COMPLEX*16 B( LDB, * ), D( * ), DL( * ), DU( * ), DU2( * ) +* .. +* +* Purpose +* ======= +* +* ZGTTS2 solves one of the systems of equations +* A * X = B, A**T * X = B, or A**H * X = B, +* with a tridiagonal matrix A using the LU factorization computed +* by ZGTTRF. +* +* Arguments +* ========= +* +* ITRANS (input) INTEGER +* Specifies the form of the system of equations. +* = 0: A * X = B (No transpose) +* = 1: A**T * X = B (Transpose) +* = 2: A**H * X = B (Conjugate transpose) +* +* N (input) INTEGER +* The order of the matrix A. +* +* NRHS (input) INTEGER +* The number of right hand sides, i.e., the number of columns +* of the matrix B. NRHS >= 0. +* +* DL (input) COMPLEX*16 array, dimension (N-1) +* The (n-1) multipliers that define the matrix L from the +* LU factorization of A. +* +* D (input) COMPLEX*16 array, dimension (N) +* The n diagonal elements of the upper triangular matrix U from +* the LU factorization of A. +* +* DU (input) COMPLEX*16 array, dimension (N-1) +* The (n-1) elements of the first super-diagonal of U. +* +* DU2 (input) COMPLEX*16 array, dimension (N-2) +* The (n-2) elements of the second super-diagonal of U. +* +* IPIV (input) INTEGER array, dimension (N) +* The pivot indices; for 1 <= i <= n, row i of the matrix was +* interchanged with row IPIV(i). IPIV(i) will always be either +* i or i+1; IPIV(i) = i indicates a row interchange was not +* required. +* +* B (input/output) COMPLEX*16 array, dimension (LDB,NRHS) +* On entry, the matrix of right hand side vectors B. +* On exit, B is overwritten by the solution vectors X. +* +* LDB (input) INTEGER +* The leading dimension of the array B. LDB >= max(1,N). +* +* ===================================================================== +* +* .. Local Scalars .. + INTEGER I, J + COMPLEX*16 TEMP +* .. +* .. Intrinsic Functions .. + INTRINSIC DCONJG +* .. +* .. Executable Statements .. +* +* Quick return if possible +* + IF( N.EQ.0 .OR. NRHS.EQ.0 ) + $ RETURN +* + IF( ITRANS.EQ.0 ) THEN +* +* Solve A*X = B using the LU factorization of A, +* overwriting each right hand side vector with its solution. +* + IF( NRHS.LE.1 ) THEN + J = 1 + 10 CONTINUE +* +* Solve L*x = b. +* + DO 20 I = 1, N - 1 + IF( IPIV( I ).EQ.I ) THEN + B( I+1, J ) = B( I+1, J ) - DL( I )*B( I, J ) + ELSE + TEMP = B( I, J ) + B( I, J ) = B( I+1, J ) + B( I+1, J ) = TEMP - DL( I )*B( I, J ) + END IF + 20 CONTINUE +* +* Solve U*x = b. +* + B( N, J ) = B( N, J ) / D( N ) + IF( N.GT.1 ) + $ B( N-1, J ) = ( B( N-1, J )-DU( N-1 )*B( N, J ) ) / + $ D( N-1 ) + DO 30 I = N - 2, 1, -1 + B( I, J ) = ( B( I, J )-DU( I )*B( I+1, J )-DU2( I )* + $ B( I+2, J ) ) / D( I ) + 30 CONTINUE + IF( J.LT.NRHS ) THEN + J = J + 1 + GO TO 10 + END IF + ELSE + DO 60 J = 1, NRHS +* +* Solve L*x = b. +* + DO 40 I = 1, N - 1 + IF( IPIV( I ).EQ.I ) THEN + B( I+1, J ) = B( I+1, J ) - DL( I )*B( I, J ) + ELSE + TEMP = B( I, J ) + B( I, J ) = B( I+1, J ) + B( I+1, J ) = TEMP - DL( I )*B( I, J ) + END IF + 40 CONTINUE +* +* Solve U*x = b. +* + B( N, J ) = B( N, J ) / D( N ) + IF( N.GT.1 ) + $ B( N-1, J ) = ( B( N-1, J )-DU( N-1 )*B( N, J ) ) / + $ D( N-1 ) + DO 50 I = N - 2, 1, -1 + B( I, J ) = ( B( I, J )-DU( I )*B( I+1, J )-DU2( I )* + $ B( I+2, J ) ) / D( I ) + 50 CONTINUE + 60 CONTINUE + END IF + ELSE IF( ITRANS.EQ.1 ) THEN +* +* Solve A**T * X = B. +* + IF( NRHS.LE.1 ) THEN + J = 1 + 70 CONTINUE +* +* Solve U**T * x = b. +* + B( 1, J ) = B( 1, J ) / D( 1 ) + IF( N.GT.1 ) + $ B( 2, J ) = ( B( 2, J )-DU( 1 )*B( 1, J ) ) / D( 2 ) + DO 80 I = 3, N + B( I, J ) = ( B( I, J )-DU( I-1 )*B( I-1, J )-DU2( I-2 )* + $ B( I-2, J ) ) / D( I ) + 80 CONTINUE +* +* Solve L**T * x = b. +* + DO 90 I = N - 1, 1, -1 + IF( IPIV( I ).EQ.I ) THEN + B( I, J ) = B( I, J ) - DL( I )*B( I+1, J ) + ELSE + TEMP = B( I+1, J ) + B( I+1, J ) = B( I, J ) - DL( I )*TEMP + B( I, J ) = TEMP + END IF + 90 CONTINUE + IF( J.LT.NRHS ) THEN + J = J + 1 + GO TO 70 + END IF + ELSE + DO 120 J = 1, NRHS +* +* Solve U**T * x = b. +* + B( 1, J ) = B( 1, J ) / D( 1 ) + IF( N.GT.1 ) + $ B( 2, J ) = ( B( 2, J )-DU( 1 )*B( 1, J ) ) / D( 2 ) + DO 100 I = 3, N + B( I, J ) = ( B( I, J )-DU( I-1 )*B( I-1, J )- + $ DU2( I-2 )*B( I-2, J ) ) / D( I ) + 100 CONTINUE +* +* Solve L**T * x = b. +* + DO 110 I = N - 1, 1, -1 + IF( IPIV( I ).EQ.I ) THEN + B( I, J ) = B( I, J ) - DL( I )*B( I+1, J ) + ELSE + TEMP = B( I+1, J ) + B( I+1, J ) = B( I, J ) - DL( I )*TEMP + B( I, J ) = TEMP + END IF + 110 CONTINUE + 120 CONTINUE + END IF + ELSE +* +* Solve A**H * X = B. +* + IF( NRHS.LE.1 ) THEN + J = 1 + 130 CONTINUE +* +* Solve U**H * x = b. +* + B( 1, J ) = B( 1, J ) / DCONJG( D( 1 ) ) + IF( N.GT.1 ) + $ B( 2, J ) = ( B( 2, J )-DCONJG( DU( 1 ) )*B( 1, J ) ) / + $ DCONJG( D( 2 ) ) + DO 140 I = 3, N + B( I, J ) = ( B( I, J )-DCONJG( DU( I-1 ) )*B( I-1, J )- + $ DCONJG( DU2( I-2 ) )*B( I-2, J ) ) / + $ DCONJG( D( I ) ) + 140 CONTINUE +* +* Solve L**H * x = b. +* + DO 150 I = N - 1, 1, -1 + IF( IPIV( I ).EQ.I ) THEN + B( I, J ) = B( I, J ) - DCONJG( DL( I ) )*B( I+1, J ) + ELSE + TEMP = B( I+1, J ) + B( I+1, J ) = B( I, J ) - DCONJG( DL( I ) )*TEMP + B( I, J ) = TEMP + END IF + 150 CONTINUE + IF( J.LT.NRHS ) THEN + J = J + 1 + GO TO 130 + END IF + ELSE + DO 180 J = 1, NRHS +* +* Solve U**H * x = b. +* + B( 1, J ) = B( 1, J ) / DCONJG( D( 1 ) ) + IF( N.GT.1 ) + $ B( 2, J ) = ( B( 2, J )-DCONJG( DU( 1 ) )*B( 1, J ) ) + $ / DCONJG( D( 2 ) ) + DO 160 I = 3, N + B( I, J ) = ( B( I, J )-DCONJG( DU( I-1 ) )* + $ B( I-1, J )-DCONJG( DU2( I-2 ) )* + $ B( I-2, J ) ) / DCONJG( D( I ) ) + 160 CONTINUE +* +* Solve L**H * x = b. +* + DO 170 I = N - 1, 1, -1 + IF( IPIV( I ).EQ.I ) THEN + B( I, J ) = B( I, J ) - DCONJG( DL( I ) )* + $ B( I+1, J ) + ELSE + TEMP = B( I+1, J ) + B( I+1, J ) = B( I, J ) - DCONJG( DL( I ) )*TEMP + B( I, J ) = TEMP + END IF + 170 CONTINUE + 180 CONTINUE + END IF + END IF +* +* End of ZGTTS2 +* + END