 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DGTTS2( ITRANS, N, NRHS, DL, D, DU, DU2, IPIV, B, LDB ) 00002 * 00003 * -- LAPACK auxiliary routine (version 3.2) -- 00004 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00005 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00006 * November 2006 00007 * 00008 * .. Scalar Arguments .. 00009 INTEGER ITRANS, LDB, N, NRHS 00010 * .. 00011 * .. Array Arguments .. 00012 INTEGER IPIV( * ) 00013 DOUBLE PRECISION B( LDB, * ), D( * ), DL( * ), DU( * ), DU2( * ) 00014 * .. 00015 * 00016 * Purpose 00017 * ======= 00018 * 00019 * DGTTS2 solves one of the systems of equations 00020 * A*X = B or A**T*X = B, 00021 * with a tridiagonal matrix A using the LU factorization computed 00022 * by DGTTRF. 00023 * 00024 * Arguments 00025 * ========= 00026 * 00027 * ITRANS (input) INTEGER 00028 * Specifies the form of the system of equations. 00029 * = 0: A * X = B (No transpose) 00030 * = 1: A**T* X = B (Transpose) 00031 * = 2: A**T* X = B (Conjugate transpose = Transpose) 00032 * 00033 * N (input) INTEGER 00034 * The order of the matrix A. 00035 * 00036 * NRHS (input) INTEGER 00037 * The number of right hand sides, i.e., the number of columns 00038 * of the matrix B. NRHS >= 0. 00039 * 00040 * DL (input) DOUBLE PRECISION array, dimension (N-1) 00041 * The (n-1) multipliers that define the matrix L from the 00042 * LU factorization of A. 00043 * 00044 * D (input) DOUBLE PRECISION array, dimension (N) 00045 * The n diagonal elements of the upper triangular matrix U from 00046 * the LU factorization of A. 00047 * 00048 * DU (input) DOUBLE PRECISION array, dimension (N-1) 00049 * The (n-1) elements of the first super-diagonal of U. 00050 * 00051 * DU2 (input) DOUBLE PRECISION array, dimension (N-2) 00052 * The (n-2) elements of the second super-diagonal of U. 00053 * 00054 * IPIV (input) INTEGER array, dimension (N) 00055 * The pivot indices; for 1 <= i <= n, row i of the matrix was 00056 * interchanged with row IPIV(i). IPIV(i) will always be either 00057 * i or i+1; IPIV(i) = i indicates a row interchange was not 00058 * required. 00059 * 00060 * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) 00061 * On entry, the matrix of right hand side vectors B. 00062 * On exit, B is overwritten by the solution vectors X. 00063 * 00064 * LDB (input) INTEGER 00065 * The leading dimension of the array B. LDB >= max(1,N). 00066 * 00067 * ===================================================================== 00068 * 00069 * .. Local Scalars .. 00070 INTEGER I, IP, J 00071 DOUBLE PRECISION TEMP 00072 * .. 00073 * .. Executable Statements .. 00074 * 00075 * Quick return if possible 00076 * 00077 IF( N.EQ.0 .OR. NRHS.EQ.0 ) 00078 $ RETURN 00079 * 00080 IF( ITRANS.EQ.0 ) THEN 00081 * 00082 * Solve A*X = B using the LU factorization of A, 00083 * overwriting each right hand side vector with its solution. 00084 * 00085 IF( NRHS.LE.1 ) THEN 00086 J = 1 00087 10 CONTINUE 00088 * 00089 * Solve L*x = b. 00090 * 00091 DO 20 I = 1, N - 1 00092 IP = IPIV( I ) 00093 TEMP = B( I+1-IP+I, J ) - DL( I )*B( IP, J ) 00094 B( I, J ) = B( IP, J ) 00095 B( I+1, J ) = TEMP 00096 20 CONTINUE 00097 * 00098 * Solve U*x = b. 00099 * 00100 B( N, J ) = B( N, J ) / D( N ) 00101 IF( N.GT.1 ) 00102 $ B( N-1, J ) = ( B( N-1, J )-DU( N-1 )*B( N, J ) ) / 00103 $ D( N-1 ) 00104 DO 30 I = N - 2, 1, -1 00105 B( I, J ) = ( B( I, J )-DU( I )*B( I+1, J )-DU2( I )* 00106 $ B( I+2, J ) ) / D( I ) 00107 30 CONTINUE 00108 IF( J.LT.NRHS ) THEN 00109 J = J + 1 00110 GO TO 10 00111 END IF 00112 ELSE 00113 DO 60 J = 1, NRHS 00114 * 00115 * Solve L*x = b. 00116 * 00117 DO 40 I = 1, N - 1 00118 IF( IPIV( I ).EQ.I ) THEN 00119 B( I+1, J ) = B( I+1, J ) - DL( I )*B( I, J ) 00120 ELSE 00121 TEMP = B( I, J ) 00122 B( I, J ) = B( I+1, J ) 00123 B( I+1, J ) = TEMP - DL( I )*B( I, J ) 00124 END IF 00125 40 CONTINUE 00126 * 00127 * Solve U*x = b. 00128 * 00129 B( N, J ) = B( N, J ) / D( N ) 00130 IF( N.GT.1 ) 00131 $ B( N-1, J ) = ( B( N-1, J )-DU( N-1 )*B( N, J ) ) / 00132 $ D( N-1 ) 00133 DO 50 I = N - 2, 1, -1 00134 B( I, J ) = ( B( I, J )-DU( I )*B( I+1, J )-DU2( I )* 00135 $ B( I+2, J ) ) / D( I ) 00136 50 CONTINUE 00137 60 CONTINUE 00138 END IF 00139 ELSE 00140 * 00141 * Solve A**T * X = B. 00142 * 00143 IF( NRHS.LE.1 ) THEN 00144 * 00145 * Solve U**T*x = b. 00146 * 00147 J = 1 00148 70 CONTINUE 00149 B( 1, J ) = B( 1, J ) / D( 1 ) 00150 IF( N.GT.1 ) 00151 $ B( 2, J ) = ( B( 2, J )-DU( 1 )*B( 1, J ) ) / D( 2 ) 00152 DO 80 I = 3, N 00153 B( I, J ) = ( B( I, J )-DU( I-1 )*B( I-1, J )-DU2( I-2 )* 00154 $ B( I-2, J ) ) / D( I ) 00155 80 CONTINUE 00156 * 00157 * Solve L**T*x = b. 00158 * 00159 DO 90 I = N - 1, 1, -1 00160 IP = IPIV( I ) 00161 TEMP = B( I, J ) - DL( I )*B( I+1, J ) 00162 B( I, J ) = B( IP, J ) 00163 B( IP, J ) = TEMP 00164 90 CONTINUE 00165 IF( J.LT.NRHS ) THEN 00166 J = J + 1 00167 GO TO 70 00168 END IF 00169 * 00170 ELSE 00171 DO 120 J = 1, NRHS 00172 * 00173 * Solve U**T*x = b. 00174 * 00175 B( 1, J ) = B( 1, J ) / D( 1 ) 00176 IF( N.GT.1 ) 00177 $ B( 2, J ) = ( B( 2, J )-DU( 1 )*B( 1, J ) ) / D( 2 ) 00178 DO 100 I = 3, N 00179 B( I, J ) = ( B( I, J )-DU( I-1 )*B( I-1, J )- 00180 $ DU2( I-2 )*B( I-2, J ) ) / D( I ) 00181 100 CONTINUE 00182 DO 110 I = N - 1, 1, -1 00183 IF( IPIV( I ).EQ.I ) THEN 00184 B( I, J ) = B( I, J ) - DL( I )*B( I+1, J ) 00185 ELSE 00186 TEMP = B( I+1, J ) 00187 B( I+1, J ) = B( I, J ) - DL( I )*TEMP 00188 B( I, J ) = TEMP 00189 END IF 00190 110 CONTINUE 00191 120 CONTINUE 00192 END IF 00193 END IF 00194 * 00195 * End of DGTTS2 00196 * 00197 END
1.7.3 
