 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE CPFTRS( TRANSR, UPLO, N, NRHS, A, B, LDB, INFO ) 00002 * 00003 * -- LAPACK routine (version 3.3.1) -- 00004 * 00005 * -- Contributed by Fred Gustavson of the IBM Watson Research Center -- 00006 * -- April 2011 -- 00007 * 00008 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00009 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00010 * 00011 * .. Scalar Arguments .. 00012 CHARACTER TRANSR, UPLO 00013 INTEGER INFO, LDB, N, NRHS 00014 * .. 00015 * .. Array Arguments .. 00016 COMPLEX A( 0: * ), B( LDB, * ) 00017 * .. 00018 * 00019 * Purpose 00020 * ======= 00021 * 00022 * CPFTRS solves a system of linear equations A*X = B with a Hermitian 00023 * positive definite matrix A using the Cholesky factorization 00024 * A = U**H*U or A = L*L**H computed by CPFTRF. 00025 * 00026 * Arguments 00027 * ========= 00028 * 00029 * TRANSR (input) CHARACTER*1 00030 * = 'N': The Normal TRANSR of RFP A is stored; 00031 * = 'C': The Conjugate-transpose TRANSR of RFP A is stored. 00032 * 00033 * UPLO (input) CHARACTER*1 00034 * = 'U': Upper triangle of RFP A is stored; 00035 * = 'L': Lower triangle of RFP A is stored. 00036 * 00037 * N (input) INTEGER 00038 * The order of the matrix A. N >= 0. 00039 * 00040 * NRHS (input) INTEGER 00041 * The number of right hand sides, i.e., the number of columns 00042 * of the matrix B. NRHS >= 0. 00043 * 00044 * A (input) COMPLEX array, dimension ( N*(N+1)/2 ); 00045 * The triangular factor U or L from the Cholesky factorization 00046 * of RFP A = U**H*U or RFP A = L*L**H, as computed by CPFTRF. 00047 * See note below for more details about RFP A. 00048 * 00049 * B (input/output) COMPLEX array, dimension (LDB,NRHS) 00050 * On entry, the right hand side matrix B. 00051 * On exit, the solution matrix X. 00052 * 00053 * LDB (input) INTEGER 00054 * The leading dimension of the array B. LDB >= max(1,N). 00055 * 00056 * INFO (output) INTEGER 00057 * = 0: successful exit 00058 * < 0: if INFO = -i, the i-th argument had an illegal value 00059 * 00060 * Further Details 00061 * =============== 00062 * 00063 * We first consider Standard Packed Format when N is even. 00064 * We give an example where N = 6. 00065 * 00066 * AP is Upper AP is Lower 00067 * 00068 * 00 01 02 03 04 05 00 00069 * 11 12 13 14 15 10 11 00070 * 22 23 24 25 20 21 22 00071 * 33 34 35 30 31 32 33 00072 * 44 45 40 41 42 43 44 00073 * 55 50 51 52 53 54 55 00074 * 00075 * 00076 * Let TRANSR = 'N'. RFP holds AP as follows: 00077 * For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last 00078 * three columns of AP upper. The lower triangle A(4:6,0:2) consists of 00079 * conjugate-transpose of the first three columns of AP upper. 00080 * For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first 00081 * three columns of AP lower. The upper triangle A(0:2,0:2) consists of 00082 * conjugate-transpose of the last three columns of AP lower. 00083 * To denote conjugate we place -- above the element. This covers the 00084 * case N even and TRANSR = 'N'. 00085 * 00086 * RFP A RFP A 00087 * 00088 * -- -- -- 00089 * 03 04 05 33 43 53 00090 * -- -- 00091 * 13 14 15 00 44 54 00092 * -- 00093 * 23 24 25 10 11 55 00094 * 00095 * 33 34 35 20 21 22 00096 * -- 00097 * 00 44 45 30 31 32 00098 * -- -- 00099 * 01 11 55 40 41 42 00100 * -- -- -- 00101 * 02 12 22 50 51 52 00102 * 00103 * Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate- 00104 * transpose of RFP A above. One therefore gets: 00105 * 00106 * 00107 * RFP A RFP A 00108 * 00109 * -- -- -- -- -- -- -- -- -- -- 00110 * 03 13 23 33 00 01 02 33 00 10 20 30 40 50 00111 * -- -- -- -- -- -- -- -- -- -- 00112 * 04 14 24 34 44 11 12 43 44 11 21 31 41 51 00113 * -- -- -- -- -- -- -- -- -- -- 00114 * 05 15 25 35 45 55 22 53 54 55 22 32 42 52 00115 * 00116 * 00117 * We next consider Standard Packed Format when N is odd. 00118 * We give an example where N = 5. 00119 * 00120 * AP is Upper AP is Lower 00121 * 00122 * 00 01 02 03 04 00 00123 * 11 12 13 14 10 11 00124 * 22 23 24 20 21 22 00125 * 33 34 30 31 32 33 00126 * 44 40 41 42 43 44 00127 * 00128 * 00129 * Let TRANSR = 'N'. RFP holds AP as follows: 00130 * For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last 00131 * three columns of AP upper. The lower triangle A(3:4,0:1) consists of 00132 * conjugate-transpose of the first two columns of AP upper. 00133 * For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first 00134 * three columns of AP lower. The upper triangle A(0:1,1:2) consists of 00135 * conjugate-transpose of the last two columns of AP lower. 00136 * To denote conjugate we place -- above the element. This covers the 00137 * case N odd and TRANSR = 'N'. 00138 * 00139 * RFP A RFP A 00140 * 00141 * -- -- 00142 * 02 03 04 00 33 43 00143 * -- 00144 * 12 13 14 10 11 44 00145 * 00146 * 22 23 24 20 21 22 00147 * -- 00148 * 00 33 34 30 31 32 00149 * -- -- 00150 * 01 11 44 40 41 42 00151 * 00152 * Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate- 00153 * transpose of RFP A above. One therefore gets: 00154 * 00155 * 00156 * RFP A RFP A 00157 * 00158 * -- -- -- -- -- -- -- -- -- 00159 * 02 12 22 00 01 00 10 20 30 40 50 00160 * -- -- -- -- -- -- -- -- -- 00161 * 03 13 23 33 11 33 11 21 31 41 51 00162 * -- -- -- -- -- -- -- -- -- 00163 * 04 14 24 34 44 43 44 22 32 42 52 00164 * 00165 * ===================================================================== 00166 * 00167 * .. Parameters .. 00168 COMPLEX CONE 00169 PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ) ) 00170 * .. 00171 * .. Local Scalars .. 00172 LOGICAL LOWER, NORMALTRANSR 00173 * .. 00174 * .. External Functions .. 00175 LOGICAL LSAME 00176 EXTERNAL LSAME 00177 * .. 00178 * .. External Subroutines .. 00179 EXTERNAL XERBLA, CTFSM 00180 * .. 00181 * .. Intrinsic Functions .. 00182 INTRINSIC MAX 00183 * .. 00184 * .. Executable Statements .. 00185 * 00186 * Test the input parameters. 00187 * 00188 INFO = 0 00189 NORMALTRANSR = LSAME( TRANSR, 'N' ) 00190 LOWER = LSAME( UPLO, 'L' ) 00191 IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'C' ) ) THEN 00192 INFO = -1 00193 ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN 00194 INFO = -2 00195 ELSE IF( N.LT.0 ) THEN 00196 INFO = -3 00197 ELSE IF( NRHS.LT.0 ) THEN 00198 INFO = -4 00199 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00200 INFO = -7 00201 END IF 00202 IF( INFO.NE.0 ) THEN 00203 CALL XERBLA( 'CPFTRS', -INFO ) 00204 RETURN 00205 END IF 00206 * 00207 * Quick return if possible 00208 * 00209 IF( N.EQ.0 .OR. NRHS.EQ.0 ) 00210 $ RETURN 00211 * 00212 * start execution: there are two triangular solves 00213 * 00214 IF( LOWER ) THEN 00215 CALL CTFSM( TRANSR, 'L', UPLO, 'N', 'N', N, NRHS, CONE, A, B, 00216 $ LDB ) 00217 CALL CTFSM( TRANSR, 'L', UPLO, 'C', 'N', N, NRHS, CONE, A, B, 00218 $ LDB ) 00219 ELSE 00220 CALL CTFSM( TRANSR, 'L', UPLO, 'C', 'N', N, NRHS, CONE, A, B, 00221 $ LDB ) 00222 CALL CTFSM( TRANSR, 'L', UPLO, 'N', 'N', N, NRHS, CONE, A, B, 00223 $ LDB ) 00224 END IF 00225 * 00226 RETURN 00227 * 00228 * End of CPFTRS 00229 * 00230 END
1.7.3 
