00001 /* cpftrs.f -- translated by f2c (version 20061008). 00002 You must link the resulting object file with libf2c: 00003 on Microsoft Windows system, link with libf2c.lib; 00004 on Linux or Unix systems, link with .../path/to/libf2c.a -lm 00005 or, if you install libf2c.a in a standard place, with -lf2c -lm 00006 -- in that order, at the end of the command line, as in 00007 cc *.o -lf2c -lm 00008 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., 00009 00010 http://www.netlib.org/f2c/libf2c.zip 00011 */ 00012 00013 #include "f2c.h" 00014 #include "blaswrap.h" 00015 00016 /* Table of constant values */ 00017 00018 static complex c_b1 = {1.f,0.f}; 00019 00020 /* Subroutine */ int cpftrs_(char *transr, char *uplo, integer *n, integer * 00021 nrhs, complex *a, complex *b, integer *ldb, integer *info) 00022 { 00023 /* System generated locals */ 00024 integer b_dim1, b_offset, i__1; 00025 00026 /* Local variables */ 00027 logical normaltransr; 00028 extern logical lsame_(char *, char *); 00029 extern /* Subroutine */ int ctfsm_(char *, char *, char *, char *, char *, 00030 integer *, integer *, complex *, complex *, complex *, integer *); 00031 logical lower; 00032 extern /* Subroutine */ int xerbla_(char *, integer *); 00033 00034 00035 /* -- LAPACK routine (version 3.2) -- */ 00036 00037 /* -- Contributed by Fred Gustavson of the IBM Watson Research Center -- */ 00038 /* -- November 2008 -- */ 00039 00040 /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */ 00041 /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */ 00042 00043 /* .. Scalar Arguments .. */ 00044 /* .. */ 00045 /* .. Array Arguments .. */ 00046 /* .. */ 00047 00048 /* Purpose */ 00049 /* ======= */ 00050 00051 /* CPFTRS solves a system of linear equations A*X = B with a Hermitian */ 00052 /* positive definite matrix A using the Cholesky factorization */ 00053 /* A = U**H*U or A = L*L**H computed by CPFTRF. */ 00054 00055 /* Arguments */ 00056 /* ========= */ 00057 00058 /* TRANSR (input) CHARACTER */ 00059 /* = 'N': The Normal TRANSR of RFP A is stored; */ 00060 /* = 'C': The Conjugate-transpose TRANSR of RFP A is stored. */ 00061 00062 /* UPLO (input) CHARACTER */ 00063 /* = 'U': Upper triangle of RFP A is stored; */ 00064 /* = 'L': Lower triangle of RFP A is stored. */ 00065 00066 /* N (input) INTEGER */ 00067 /* The order of the matrix A. N >= 0. */ 00068 00069 /* NRHS (input) INTEGER */ 00070 /* The number of right hand sides, i.e., the number of columns */ 00071 /* of the matrix B. NRHS >= 0. */ 00072 00073 /* A (input) COMPLEX array, dimension ( N*(N+1)/2 ); */ 00074 /* The triangular factor U or L from the Cholesky factorization */ 00075 /* of RFP A = U**H*U or RFP A = L*L**H, as computed by CPFTRF. */ 00076 /* See note below for more details about RFP A. */ 00077 00078 /* B (input/output) COMPLEX array, dimension (LDB,NRHS) */ 00079 /* On entry, the right hand side matrix B. */ 00080 /* On exit, the solution matrix X. */ 00081 00082 /* LDB (input) INTEGER */ 00083 /* The leading dimension of the array B. LDB >= max(1,N). */ 00084 00085 /* INFO (output) INTEGER */ 00086 /* = 0: successful exit */ 00087 /* < 0: if INFO = -i, the i-th argument had an illegal value */ 00088 00089 /* Note: */ 00090 /* ===== */ 00091 00092 /* We first consider Standard Packed Format when N is even. */ 00093 /* We give an example where N = 6. */ 00094 00095 /* AP is Upper AP is Lower */ 00096 00097 /* 00 01 02 03 04 05 00 */ 00098 /* 11 12 13 14 15 10 11 */ 00099 /* 22 23 24 25 20 21 22 */ 00100 /* 33 34 35 30 31 32 33 */ 00101 /* 44 45 40 41 42 43 44 */ 00102 /* 55 50 51 52 53 54 55 */ 00103 00104 00105 /* Let TRANSR = 'N'. RFP holds AP as follows: */ 00106 /* For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last */ 00107 /* three columns of AP upper. The lower triangle A(4:6,0:2) consists of */ 00108 /* conjugate-transpose of the first three columns of AP upper. */ 00109 /* For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first */ 00110 /* three columns of AP lower. The upper triangle A(0:2,0:2) consists of */ 00111 /* conjugate-transpose of the last three columns of AP lower. */ 00112 /* To denote conjugate we place -- above the element. This covers the */ 00113 /* case N even and TRANSR = 'N'. */ 00114 00115 /* RFP A RFP A */ 00116 00117 /* -- -- -- */ 00118 /* 03 04 05 33 43 53 */ 00119 /* -- -- */ 00120 /* 13 14 15 00 44 54 */ 00121 /* -- */ 00122 /* 23 24 25 10 11 55 */ 00123 00124 /* 33 34 35 20 21 22 */ 00125 /* -- */ 00126 /* 00 44 45 30 31 32 */ 00127 /* -- -- */ 00128 /* 01 11 55 40 41 42 */ 00129 /* -- -- -- */ 00130 /* 02 12 22 50 51 52 */ 00131 00132 /* Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate- */ 00133 /* transpose of RFP A above. One therefore gets: */ 00134 00135 00136 /* RFP A RFP A */ 00137 00138 /* -- -- -- -- -- -- -- -- -- -- */ 00139 /* 03 13 23 33 00 01 02 33 00 10 20 30 40 50 */ 00140 /* -- -- -- -- -- -- -- -- -- -- */ 00141 /* 04 14 24 34 44 11 12 43 44 11 21 31 41 51 */ 00142 /* -- -- -- -- -- -- -- -- -- -- */ 00143 /* 05 15 25 35 45 55 22 53 54 55 22 32 42 52 */ 00144 00145 00146 /* We next consider Standard Packed Format when N is odd. */ 00147 /* We give an example where N = 5. */ 00148 00149 /* AP is Upper AP is Lower */ 00150 00151 /* 00 01 02 03 04 00 */ 00152 /* 11 12 13 14 10 11 */ 00153 /* 22 23 24 20 21 22 */ 00154 /* 33 34 30 31 32 33 */ 00155 /* 44 40 41 42 43 44 */ 00156 00157 00158 /* Let TRANSR = 'N'. RFP holds AP as follows: */ 00159 /* For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last */ 00160 /* three columns of AP upper. The lower triangle A(3:4,0:1) consists of */ 00161 /* conjugate-transpose of the first two columns of AP upper. */ 00162 /* For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first */ 00163 /* three columns of AP lower. The upper triangle A(0:1,1:2) consists of */ 00164 /* conjugate-transpose of the last two columns of AP lower. */ 00165 /* To denote conjugate we place -- above the element. This covers the */ 00166 /* case N odd and TRANSR = 'N'. */ 00167 00168 /* RFP A RFP A */ 00169 00170 /* -- -- */ 00171 /* 02 03 04 00 33 43 */ 00172 /* -- */ 00173 /* 12 13 14 10 11 44 */ 00174 00175 /* 22 23 24 20 21 22 */ 00176 /* -- */ 00177 /* 00 33 34 30 31 32 */ 00178 /* -- -- */ 00179 /* 01 11 44 40 41 42 */ 00180 00181 /* Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate- */ 00182 /* transpose of RFP A above. One therefore gets: */ 00183 00184 00185 /* RFP A RFP A */ 00186 00187 /* -- -- -- -- -- -- -- -- -- */ 00188 /* 02 12 22 00 01 00 10 20 30 40 50 */ 00189 /* -- -- -- -- -- -- -- -- -- */ 00190 /* 03 13 23 33 11 33 11 21 31 41 51 */ 00191 /* -- -- -- -- -- -- -- -- -- */ 00192 /* 04 14 24 34 44 43 44 22 32 42 52 */ 00193 00194 /* ===================================================================== */ 00195 00196 /* .. Parameters .. */ 00197 /* .. */ 00198 /* .. Local Scalars .. */ 00199 /* .. */ 00200 /* .. External Functions .. */ 00201 /* .. */ 00202 /* .. External Subroutines .. */ 00203 /* .. */ 00204 /* .. Intrinsic Functions .. */ 00205 /* .. */ 00206 /* .. Executable Statements .. */ 00207 00208 /* Test the input parameters. */ 00209 00210 /* Parameter adjustments */ 00211 b_dim1 = *ldb; 00212 b_offset = 1 + b_dim1; 00213 b -= b_offset; 00214 00215 /* Function Body */ 00216 *info = 0; 00217 normaltransr = lsame_(transr, "N"); 00218 lower = lsame_(uplo, "L"); 00219 if (! normaltransr && ! lsame_(transr, "C")) { 00220 *info = -1; 00221 } else if (! lower && ! lsame_(uplo, "U")) { 00222 *info = -2; 00223 } else if (*n < 0) { 00224 *info = -3; 00225 } else if (*nrhs < 0) { 00226 *info = -4; 00227 } else if (*ldb < max(1,*n)) { 00228 *info = -7; 00229 } 00230 if (*info != 0) { 00231 i__1 = -(*info); 00232 xerbla_("CPFTRS", &i__1); 00233 return 0; 00234 } 00235 00236 /* Quick return if possible */ 00237 00238 if (*n == 0 || *nrhs == 0) { 00239 return 0; 00240 } 00241 00242 /* start execution: there are two triangular solves */ 00243 00244 if (lower) { 00245 ctfsm_(transr, "L", uplo, "N", "N", n, nrhs, &c_b1, a, &b[b_offset], 00246 ldb); 00247 ctfsm_(transr, "L", uplo, "C", "N", n, nrhs, &c_b1, a, &b[b_offset], 00248 ldb); 00249 } else { 00250 ctfsm_(transr, "L", uplo, "C", "N", n, nrhs, &c_b1, a, &b[b_offset], 00251 ldb); 00252 ctfsm_(transr, "L", uplo, "N", "N", n, nrhs, &c_b1, a, &b[b_offset], 00253 ldb); 00254 } 00255 00256 return 0; 00257 00258 /* End of CPFTRS */ 00259 00260 } /* cpftrs_ */