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