00001 SUBROUTINE CPTTS2( IUPLO, N, NRHS, D, E, B, LDB ) 00002 * 00003 * -- LAPACK 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 IUPLO, LDB, N, NRHS 00010 * .. 00011 * .. Array Arguments .. 00012 REAL D( * ) 00013 COMPLEX B( LDB, * ), E( * ) 00014 * .. 00015 * 00016 * Purpose 00017 * ======= 00018 * 00019 * CPTTS2 solves a tridiagonal system of the form 00020 * A * X = B 00021 * using the factorization A = U'*D*U or A = L*D*L' computed by CPTTRF. 00022 * D is a diagonal matrix specified in the vector D, U (or L) is a unit 00023 * bidiagonal matrix whose superdiagonal (subdiagonal) is specified in 00024 * the vector E, and X and B are N by NRHS matrices. 00025 * 00026 * Arguments 00027 * ========= 00028 * 00029 * IUPLO (input) INTEGER 00030 * Specifies the form of the factorization and whether the 00031 * vector E is the superdiagonal of the upper bidiagonal factor 00032 * U or the subdiagonal of the lower bidiagonal factor L. 00033 * = 1: A = U'*D*U, E is the superdiagonal of U 00034 * = 0: A = L*D*L', E is the subdiagonal of L 00035 * 00036 * N (input) INTEGER 00037 * The order of the tridiagonal matrix A. N >= 0. 00038 * 00039 * NRHS (input) INTEGER 00040 * The number of right hand sides, i.e., the number of columns 00041 * of the matrix B. NRHS >= 0. 00042 * 00043 * D (input) REAL array, dimension (N) 00044 * The n diagonal elements of the diagonal matrix D from the 00045 * factorization A = U'*D*U or A = L*D*L'. 00046 * 00047 * E (input) COMPLEX array, dimension (N-1) 00048 * If IUPLO = 1, the (n-1) superdiagonal elements of the unit 00049 * bidiagonal factor U from the factorization A = U'*D*U. 00050 * If IUPLO = 0, the (n-1) subdiagonal elements of the unit 00051 * bidiagonal factor L from the factorization A = L*D*L'. 00052 * 00053 * B (input/output) REAL array, dimension (LDB,NRHS) 00054 * On entry, the right hand side vectors B for the system of 00055 * linear equations. 00056 * On exit, the solution vectors, X. 00057 * 00058 * LDB (input) INTEGER 00059 * The leading dimension of the array B. LDB >= max(1,N). 00060 * 00061 * ===================================================================== 00062 * 00063 * .. Local Scalars .. 00064 INTEGER I, J 00065 * .. 00066 * .. External Subroutines .. 00067 EXTERNAL CSSCAL 00068 * .. 00069 * .. Intrinsic Functions .. 00070 INTRINSIC CONJG 00071 * .. 00072 * .. Executable Statements .. 00073 * 00074 * Quick return if possible 00075 * 00076 IF( N.LE.1 ) THEN 00077 IF( N.EQ.1 ) 00078 $ CALL CSSCAL( NRHS, 1. / D( 1 ), B, LDB ) 00079 RETURN 00080 END IF 00081 * 00082 IF( IUPLO.EQ.1 ) THEN 00083 * 00084 * Solve A * X = B using the factorization A = U'*D*U, 00085 * overwriting each right hand side vector with its solution. 00086 * 00087 IF( NRHS.LE.2 ) THEN 00088 J = 1 00089 5 CONTINUE 00090 * 00091 * Solve U' * x = b. 00092 * 00093 DO 10 I = 2, N 00094 B( I, J ) = B( I, J ) - B( I-1, J )*CONJG( E( I-1 ) ) 00095 10 CONTINUE 00096 * 00097 * Solve D * U * x = b. 00098 * 00099 DO 20 I = 1, N 00100 B( I, J ) = B( I, J ) / D( I ) 00101 20 CONTINUE 00102 DO 30 I = N - 1, 1, -1 00103 B( I, J ) = B( I, J ) - B( I+1, J )*E( I ) 00104 30 CONTINUE 00105 IF( J.LT.NRHS ) THEN 00106 J = J + 1 00107 GO TO 5 00108 END IF 00109 ELSE 00110 DO 60 J = 1, NRHS 00111 * 00112 * Solve U' * x = b. 00113 * 00114 DO 40 I = 2, N 00115 B( I, J ) = B( I, J ) - B( I-1, J )*CONJG( E( I-1 ) ) 00116 40 CONTINUE 00117 * 00118 * Solve D * U * x = b. 00119 * 00120 B( N, J ) = B( N, J ) / D( N ) 00121 DO 50 I = N - 1, 1, -1 00122 B( I, J ) = B( I, J ) / D( I ) - B( I+1, J )*E( I ) 00123 50 CONTINUE 00124 60 CONTINUE 00125 END IF 00126 ELSE 00127 * 00128 * Solve A * X = B using the factorization A = L*D*L', 00129 * overwriting each right hand side vector with its solution. 00130 * 00131 IF( NRHS.LE.2 ) THEN 00132 J = 1 00133 65 CONTINUE 00134 * 00135 * Solve L * x = b. 00136 * 00137 DO 70 I = 2, N 00138 B( I, J ) = B( I, J ) - B( I-1, J )*E( I-1 ) 00139 70 CONTINUE 00140 * 00141 * Solve D * L' * x = b. 00142 * 00143 DO 80 I = 1, N 00144 B( I, J ) = B( I, J ) / D( I ) 00145 80 CONTINUE 00146 DO 90 I = N - 1, 1, -1 00147 B( I, J ) = B( I, J ) - B( I+1, J )*CONJG( E( I ) ) 00148 90 CONTINUE 00149 IF( J.LT.NRHS ) THEN 00150 J = J + 1 00151 GO TO 65 00152 END IF 00153 ELSE 00154 DO 120 J = 1, NRHS 00155 * 00156 * Solve L * x = b. 00157 * 00158 DO 100 I = 2, N 00159 B( I, J ) = B( I, J ) - B( I-1, J )*E( I-1 ) 00160 100 CONTINUE 00161 * 00162 * Solve D * L' * x = b. 00163 * 00164 B( N, J ) = B( N, J ) / D( N ) 00165 DO 110 I = N - 1, 1, -1 00166 B( I, J ) = B( I, J ) / D( I ) - 00167 $ B( I+1, J )*CONJG( E( I ) ) 00168 110 CONTINUE 00169 120 CONTINUE 00170 END IF 00171 END IF 00172 * 00173 RETURN 00174 * 00175 * End of CPTTS2 00176 * 00177 END
1.7.2