00001 /* dpftrs.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 doublereal c_b10 = 1.; 00019 00020 /* Subroutine */ int dpftrs_(char *transr, char *uplo, integer *n, integer * 00021 nrhs, doublereal *a, doublereal *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 dtfsm_(char *, char *, char *, char *, char *, 00030 integer *, integer *, doublereal *, doublereal *, doublereal *, 00031 integer *); 00032 logical lower; 00033 extern /* Subroutine */ int xerbla_(char *, integer *); 00034 00035 00036 /* -- LAPACK routine (version 3.2) -- */ 00037 00038 /* -- Contributed by Fred Gustavson of the IBM Watson Research Center -- */ 00039 /* -- November 2008 -- */ 00040 00041 /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */ 00042 /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */ 00043 00044 /* .. Scalar Arguments .. */ 00045 /* .. */ 00046 /* .. Array Arguments .. */ 00047 /* .. */ 00048 00049 /* Purpose */ 00050 /* ======= */ 00051 00052 /* DPFTRS solves a system of linear equations A*X = B with a symmetric */ 00053 /* positive definite matrix A using the Cholesky factorization */ 00054 /* A = U**T*U or A = L*L**T computed by DPFTRF. */ 00055 00056 /* Arguments */ 00057 /* ========= */ 00058 00059 /* TRANSR (input) CHARACTER */ 00060 /* = 'N': The Normal TRANSR of RFP A is stored; */ 00061 /* = 'T': The Transpose TRANSR of RFP A is stored. */ 00062 00063 /* UPLO (input) CHARACTER */ 00064 /* = 'U': Upper triangle of RFP A is stored; */ 00065 /* = 'L': Lower triangle of RFP A is stored. */ 00066 00067 /* N (input) INTEGER */ 00068 /* The order of the matrix A. N >= 0. */ 00069 00070 /* NRHS (input) INTEGER */ 00071 /* The number of right hand sides, i.e., the number of columns */ 00072 /* of the matrix B. NRHS >= 0. */ 00073 00074 /* A (input) DOUBLE PRECISION array, dimension ( N*(N+1)/2 ). */ 00075 /* The triangular factor U or L from the Cholesky factorization */ 00076 /* of RFP A = U**T*U or RFP A = L*L**T, as computed by DPFTRF. */ 00077 /* See note below for more details about RFP A. */ 00078 00079 /* B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) */ 00080 /* On entry, the right hand side matrix B. */ 00081 /* On exit, the solution matrix X. */ 00082 00083 /* LDB (input) INTEGER */ 00084 /* The leading dimension of the array B. LDB >= max(1,N). */ 00085 00086 /* INFO (output) INTEGER */ 00087 /* = 0: successful exit */ 00088 /* < 0: if INFO = -i, the i-th argument had an illegal value */ 00089 00090 /* Notes */ 00091 /* ===== */ 00092 00093 /* We first consider Rectangular Full Packed (RFP) Format when N is */ 00094 /* even. We give an example where N = 6. */ 00095 00096 /* AP is Upper AP is Lower */ 00097 00098 /* 00 01 02 03 04 05 00 */ 00099 /* 11 12 13 14 15 10 11 */ 00100 /* 22 23 24 25 20 21 22 */ 00101 /* 33 34 35 30 31 32 33 */ 00102 /* 44 45 40 41 42 43 44 */ 00103 /* 55 50 51 52 53 54 55 */ 00104 00105 00106 /* Let TRANSR = 'N'. RFP holds AP as follows: */ 00107 /* For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last */ 00108 /* three columns of AP upper. The lower triangle A(4:6,0:2) consists of */ 00109 /* the transpose of the first three columns of AP upper. */ 00110 /* For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first */ 00111 /* three columns of AP lower. The upper triangle A(0:2,0:2) consists of */ 00112 /* the transpose of the last three columns of AP lower. */ 00113 /* This covers the case N even and TRANSR = 'N'. */ 00114 00115 /* RFP A RFP A */ 00116 00117 /* 03 04 05 33 43 53 */ 00118 /* 13 14 15 00 44 54 */ 00119 /* 23 24 25 10 11 55 */ 00120 /* 33 34 35 20 21 22 */ 00121 /* 00 44 45 30 31 32 */ 00122 /* 01 11 55 40 41 42 */ 00123 /* 02 12 22 50 51 52 */ 00124 00125 /* Now let TRANSR = 'T'. RFP A in both UPLO cases is just the */ 00126 /* transpose of RFP A above. One therefore gets: */ 00127 00128 00129 /* RFP A RFP A */ 00130 00131 /* 03 13 23 33 00 01 02 33 00 10 20 30 40 50 */ 00132 /* 04 14 24 34 44 11 12 43 44 11 21 31 41 51 */ 00133 /* 05 15 25 35 45 55 22 53 54 55 22 32 42 52 */ 00134 00135 00136 /* We first consider Rectangular Full Packed (RFP) Format when N is */ 00137 /* odd. We give an example where N = 5. */ 00138 00139 /* AP is Upper AP is Lower */ 00140 00141 /* 00 01 02 03 04 00 */ 00142 /* 11 12 13 14 10 11 */ 00143 /* 22 23 24 20 21 22 */ 00144 /* 33 34 30 31 32 33 */ 00145 /* 44 40 41 42 43 44 */ 00146 00147 00148 /* Let TRANSR = 'N'. RFP holds AP as follows: */ 00149 /* For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last */ 00150 /* three columns of AP upper. The lower triangle A(3:4,0:1) consists of */ 00151 /* the transpose of the first two columns of AP upper. */ 00152 /* For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first */ 00153 /* three columns of AP lower. The upper triangle A(0:1,1:2) consists of */ 00154 /* the transpose of the last two columns of AP lower. */ 00155 /* This covers the case N odd and TRANSR = 'N'. */ 00156 00157 /* RFP A RFP A */ 00158 00159 /* 02 03 04 00 33 43 */ 00160 /* 12 13 14 10 11 44 */ 00161 /* 22 23 24 20 21 22 */ 00162 /* 00 33 34 30 31 32 */ 00163 /* 01 11 44 40 41 42 */ 00164 00165 /* Now let TRANSR = 'T'. RFP A in both UPLO cases is just the */ 00166 /* transpose of RFP A above. One therefore gets: */ 00167 00168 /* RFP A RFP A */ 00169 00170 /* 02 12 22 00 01 00 10 20 30 40 50 */ 00171 /* 03 13 23 33 11 33 11 21 31 41 51 */ 00172 /* 04 14 24 34 44 43 44 22 32 42 52 */ 00173 00174 /* ===================================================================== */ 00175 00176 /* .. Parameters .. */ 00177 /* .. */ 00178 /* .. Local Scalars .. */ 00179 /* .. */ 00180 /* .. External Functions .. */ 00181 /* .. */ 00182 /* .. External Subroutines .. */ 00183 /* .. */ 00184 /* .. Intrinsic Functions .. */ 00185 /* .. */ 00186 /* .. Executable Statements .. */ 00187 00188 /* Test the input parameters. */ 00189 00190 /* Parameter adjustments */ 00191 b_dim1 = *ldb; 00192 b_offset = 1 + b_dim1; 00193 b -= b_offset; 00194 00195 /* Function Body */ 00196 *info = 0; 00197 normaltransr = lsame_(transr, "N"); 00198 lower = lsame_(uplo, "L"); 00199 if (! normaltransr && ! lsame_(transr, "T")) { 00200 *info = -1; 00201 } else if (! lower && ! lsame_(uplo, "U")) { 00202 *info = -2; 00203 } else if (*n < 0) { 00204 *info = -3; 00205 } else if (*nrhs < 0) { 00206 *info = -4; 00207 } else if (*ldb < max(1,*n)) { 00208 *info = -7; 00209 } 00210 if (*info != 0) { 00211 i__1 = -(*info); 00212 xerbla_("DPFTRS", &i__1); 00213 return 0; 00214 } 00215 00216 /* Quick return if possible */ 00217 00218 if (*n == 0 || *nrhs == 0) { 00219 return 0; 00220 } 00221 00222 /* start execution: there are two triangular solves */ 00223 00224 if (lower) { 00225 dtfsm_(transr, "L", uplo, "N", "N", n, nrhs, &c_b10, a, &b[b_offset], 00226 ldb); 00227 dtfsm_(transr, "L", uplo, "T", "N", n, nrhs, &c_b10, a, &b[b_offset], 00228 ldb); 00229 } else { 00230 dtfsm_(transr, "L", uplo, "T", "N", n, nrhs, &c_b10, a, &b[b_offset], 00231 ldb); 00232 dtfsm_(transr, "L", uplo, "N", "N", n, nrhs, &c_b10, a, &b[b_offset], 00233 ldb); 00234 } 00235 00236 return 0; 00237 00238 /* End of DPFTRS */ 00239 00240 } /* dpftrs_ */