00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037 #ifndef TEMPLATE_BLAS_TRSV_HEADER
00038 #define TEMPLATE_BLAS_TRSV_HEADER
00039
00040
00041 template<class Treal>
00042 int template_blas_trsv(const char *uplo, const char *trans, const char *diag, const integer *n,
00043 const Treal *a, const integer *lda, Treal *x, const integer *incx)
00044 {
00045
00046 integer a_dim1, a_offset, i__1, i__2;
00047
00048 integer info;
00049 Treal temp;
00050 integer i__, j;
00051 integer ix, jx, kx;
00052 logical nounit;
00053 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122 a_dim1 = *lda;
00123 a_offset = 1 + a_dim1 * 1;
00124 a -= a_offset;
00125 --x;
00126
00127 kx = 0;
00128
00129 info = 0;
00130 if (! template_blas_lsame(uplo, "U") && ! template_blas_lsame(uplo, "L")) {
00131 info = 1;
00132 } else if (! template_blas_lsame(trans, "N") && ! template_blas_lsame(trans,
00133 "T") && ! template_blas_lsame(trans, "C")) {
00134 info = 2;
00135 } else if (! template_blas_lsame(diag, "U") && ! template_blas_lsame(diag,
00136 "N")) {
00137 info = 3;
00138 } else if (*n < 0) {
00139 info = 4;
00140 } else if (*lda < maxMACRO(1,*n)) {
00141 info = 6;
00142 } else if (*incx == 0) {
00143 info = 8;
00144 }
00145 if (info != 0) {
00146 template_blas_erbla("TRSV ", &info);
00147 return 0;
00148 }
00149
00150 if (*n == 0) {
00151 return 0;
00152 }
00153 nounit = template_blas_lsame(diag, "N");
00154
00155
00156 if (*incx <= 0) {
00157 kx = 1 - (*n - 1) * *incx;
00158 } else if (*incx != 1) {
00159 kx = 1;
00160 }
00161
00162
00163 if (template_blas_lsame(trans, "N")) {
00164
00165 if (template_blas_lsame(uplo, "U")) {
00166 if (*incx == 1) {
00167 for (j = *n; j >= 1; --j) {
00168 if (x[j] != 0.) {
00169 if (nounit) {
00170 x[j] /= a_ref(j, j);
00171 }
00172 temp = x[j];
00173 for (i__ = j - 1; i__ >= 1; --i__) {
00174 x[i__] -= temp * a_ref(i__, j);
00175
00176 }
00177 }
00178
00179 }
00180 } else {
00181 jx = kx + (*n - 1) * *incx;
00182 for (j = *n; j >= 1; --j) {
00183 if (x[jx] != 0.) {
00184 if (nounit) {
00185 x[jx] /= a_ref(j, j);
00186 }
00187 temp = x[jx];
00188 ix = jx;
00189 for (i__ = j - 1; i__ >= 1; --i__) {
00190 ix -= *incx;
00191 x[ix] -= temp * a_ref(i__, j);
00192
00193 }
00194 }
00195 jx -= *incx;
00196
00197 }
00198 }
00199 } else {
00200 if (*incx == 1) {
00201 i__1 = *n;
00202 for (j = 1; j <= i__1; ++j) {
00203 if (x[j] != 0.) {
00204 if (nounit) {
00205 x[j] /= a_ref(j, j);
00206 }
00207 temp = x[j];
00208 i__2 = *n;
00209 for (i__ = j + 1; i__ <= i__2; ++i__) {
00210 x[i__] -= temp * a_ref(i__, j);
00211
00212 }
00213 }
00214
00215 }
00216 } else {
00217 jx = kx;
00218 i__1 = *n;
00219 for (j = 1; j <= i__1; ++j) {
00220 if (x[jx] != 0.) {
00221 if (nounit) {
00222 x[jx] /= a_ref(j, j);
00223 }
00224 temp = x[jx];
00225 ix = jx;
00226 i__2 = *n;
00227 for (i__ = j + 1; i__ <= i__2; ++i__) {
00228 ix += *incx;
00229 x[ix] -= temp * a_ref(i__, j);
00230
00231 }
00232 }
00233 jx += *incx;
00234
00235 }
00236 }
00237 }
00238 } else {
00239
00240 if (template_blas_lsame(uplo, "U")) {
00241 if (*incx == 1) {
00242 i__1 = *n;
00243 for (j = 1; j <= i__1; ++j) {
00244 temp = x[j];
00245 i__2 = j - 1;
00246 for (i__ = 1; i__ <= i__2; ++i__) {
00247 temp -= a_ref(i__, j) * x[i__];
00248
00249 }
00250 if (nounit) {
00251 temp /= a_ref(j, j);
00252 }
00253 x[j] = temp;
00254
00255 }
00256 } else {
00257 jx = kx;
00258 i__1 = *n;
00259 for (j = 1; j <= i__1; ++j) {
00260 temp = x[jx];
00261 ix = kx;
00262 i__2 = j - 1;
00263 for (i__ = 1; i__ <= i__2; ++i__) {
00264 temp -= a_ref(i__, j) * x[ix];
00265 ix += *incx;
00266
00267 }
00268 if (nounit) {
00269 temp /= a_ref(j, j);
00270 }
00271 x[jx] = temp;
00272 jx += *incx;
00273
00274 }
00275 }
00276 } else {
00277 if (*incx == 1) {
00278 for (j = *n; j >= 1; --j) {
00279 temp = x[j];
00280 i__1 = j + 1;
00281 for (i__ = *n; i__ >= i__1; --i__) {
00282 temp -= a_ref(i__, j) * x[i__];
00283
00284 }
00285 if (nounit) {
00286 temp /= a_ref(j, j);
00287 }
00288 x[j] = temp;
00289
00290 }
00291 } else {
00292 kx += (*n - 1) * *incx;
00293 jx = kx;
00294 for (j = *n; j >= 1; --j) {
00295 temp = x[jx];
00296 ix = kx;
00297 i__1 = j + 1;
00298 for (i__ = *n; i__ >= i__1; --i__) {
00299 temp -= a_ref(i__, j) * x[ix];
00300 ix -= *incx;
00301
00302 }
00303 if (nounit) {
00304 temp /= a_ref(j, j);
00305 }
00306 x[jx] = temp;
00307 jx -= *incx;
00308
00309 }
00310 }
00311 }
00312 }
00313 return 0;
00314
00315 }
00316 #undef a_ref
00317
00318 #endif