template_blas_trsv.h
Go to the documentation of this file.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 #ifndef TEMPLATE_BLAS_TRSV_HEADER
00036 #define TEMPLATE_BLAS_TRSV_HEADER
00037
00038
00039 template<class Treal>
00040 int template_blas_trsv(const char *uplo, const char *trans, const char *diag, const integer *n,
00041 const Treal *a, const integer *lda, Treal *x, const integer *incx)
00042 {
00043
00044 integer a_dim1, a_offset, i__1, i__2;
00045
00046 integer info;
00047 Treal temp;
00048 integer i__, j;
00049 integer ix, jx, kx;
00050 logical nounit;
00051 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
00052
00053
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 a_dim1 = *lda;
00121 a_offset = 1 + a_dim1 * 1;
00122 a -= a_offset;
00123 --x;
00124
00125 kx = 0;
00126
00127 info = 0;
00128 if (! template_blas_lsame(uplo, "U") && ! template_blas_lsame(uplo, "L")) {
00129 info = 1;
00130 } else if (! template_blas_lsame(trans, "N") && ! template_blas_lsame(trans,
00131 "T") && ! template_blas_lsame(trans, "C")) {
00132 info = 2;
00133 } else if (! template_blas_lsame(diag, "U") && ! template_blas_lsame(diag,
00134 "N")) {
00135 info = 3;
00136 } else if (*n < 0) {
00137 info = 4;
00138 } else if (*lda < maxMACRO(1,*n)) {
00139 info = 6;
00140 } else if (*incx == 0) {
00141 info = 8;
00142 }
00143 if (info != 0) {
00144 template_blas_erbla("TRSV ", &info);
00145 return 0;
00146 }
00147
00148 if (*n == 0) {
00149 return 0;
00150 }
00151 nounit = template_blas_lsame(diag, "N");
00152
00153
00154 if (*incx <= 0) {
00155 kx = 1 - (*n - 1) * *incx;
00156 } else if (*incx != 1) {
00157 kx = 1;
00158 }
00159
00160
00161 if (template_blas_lsame(trans, "N")) {
00162
00163 if (template_blas_lsame(uplo, "U")) {
00164 if (*incx == 1) {
00165 for (j = *n; j >= 1; --j) {
00166 if (x[j] != 0.) {
00167 if (nounit) {
00168 x[j] /= a_ref(j, j);
00169 }
00170 temp = x[j];
00171 for (i__ = j - 1; i__ >= 1; --i__) {
00172 x[i__] -= temp * a_ref(i__, j);
00173
00174 }
00175 }
00176
00177 }
00178 } else {
00179 jx = kx + (*n - 1) * *incx;
00180 for (j = *n; j >= 1; --j) {
00181 if (x[jx] != 0.) {
00182 if (nounit) {
00183 x[jx] /= a_ref(j, j);
00184 }
00185 temp = x[jx];
00186 ix = jx;
00187 for (i__ = j - 1; i__ >= 1; --i__) {
00188 ix -= *incx;
00189 x[ix] -= temp * a_ref(i__, j);
00190
00191 }
00192 }
00193 jx -= *incx;
00194
00195 }
00196 }
00197 } else {
00198 if (*incx == 1) {
00199 i__1 = *n;
00200 for (j = 1; j <= i__1; ++j) {
00201 if (x[j] != 0.) {
00202 if (nounit) {
00203 x[j] /= a_ref(j, j);
00204 }
00205 temp = x[j];
00206 i__2 = *n;
00207 for (i__ = j + 1; i__ <= i__2; ++i__) {
00208 x[i__] -= temp * a_ref(i__, j);
00209
00210 }
00211 }
00212
00213 }
00214 } else {
00215 jx = kx;
00216 i__1 = *n;
00217 for (j = 1; j <= i__1; ++j) {
00218 if (x[jx] != 0.) {
00219 if (nounit) {
00220 x[jx] /= a_ref(j, j);
00221 }
00222 temp = x[jx];
00223 ix = jx;
00224 i__2 = *n;
00225 for (i__ = j + 1; i__ <= i__2; ++i__) {
00226 ix += *incx;
00227 x[ix] -= temp * a_ref(i__, j);
00228
00229 }
00230 }
00231 jx += *incx;
00232
00233 }
00234 }
00235 }
00236 } else {
00237
00238 if (template_blas_lsame(uplo, "U")) {
00239 if (*incx == 1) {
00240 i__1 = *n;
00241 for (j = 1; j <= i__1; ++j) {
00242 temp = x[j];
00243 i__2 = j - 1;
00244 for (i__ = 1; i__ <= i__2; ++i__) {
00245 temp -= a_ref(i__, j) * x[i__];
00246
00247 }
00248 if (nounit) {
00249 temp /= a_ref(j, j);
00250 }
00251 x[j] = temp;
00252
00253 }
00254 } else {
00255 jx = kx;
00256 i__1 = *n;
00257 for (j = 1; j <= i__1; ++j) {
00258 temp = x[jx];
00259 ix = kx;
00260 i__2 = j - 1;
00261 for (i__ = 1; i__ <= i__2; ++i__) {
00262 temp -= a_ref(i__, j) * x[ix];
00263 ix += *incx;
00264
00265 }
00266 if (nounit) {
00267 temp /= a_ref(j, j);
00268 }
00269 x[jx] = temp;
00270 jx += *incx;
00271
00272 }
00273 }
00274 } else {
00275 if (*incx == 1) {
00276 for (j = *n; j >= 1; --j) {
00277 temp = x[j];
00278 i__1 = j + 1;
00279 for (i__ = *n; i__ >= i__1; --i__) {
00280 temp -= a_ref(i__, j) * x[i__];
00281
00282 }
00283 if (nounit) {
00284 temp /= a_ref(j, j);
00285 }
00286 x[j] = temp;
00287
00288 }
00289 } else {
00290 kx += (*n - 1) * *incx;
00291 jx = kx;
00292 for (j = *n; j >= 1; --j) {
00293 temp = x[jx];
00294 ix = kx;
00295 i__1 = j + 1;
00296 for (i__ = *n; i__ >= i__1; --i__) {
00297 temp -= a_ref(i__, j) * x[ix];
00298 ix -= *incx;
00299
00300 }
00301 if (nounit) {
00302 temp /= a_ref(j, j);
00303 }
00304 x[jx] = temp;
00305 jx -= *incx;
00306
00307 }
00308 }
00309 }
00310 }
00311 return 0;
00312
00313 }
00314 #undef a_ref
00315
00316 #endif