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_TRMV_HEADER
00038 #define TEMPLATE_BLAS_TRMV_HEADER
00039
00040
00041 template<class Treal>
00042 int template_blas_trmv(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 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("TRMV ", &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 i__1 = *n;
00166 for (j = 1; j <= i__1; ++j) {
00167 if (x[j] != 0.) {
00168 temp = x[j];
00169 i__2 = j - 1;
00170 for (i__ = 1; i__ <= i__2; ++i__) {
00171 x[i__] += temp * a_ref(i__, j);
00172
00173 }
00174 if (nounit) {
00175 x[j] *= a_ref(j, j);
00176 }
00177 }
00178
00179 }
00180 } else {
00181 jx = kx;
00182 i__1 = *n;
00183 for (j = 1; j <= i__1; ++j) {
00184 if (x[jx] != 0.) {
00185 temp = x[jx];
00186 ix = kx;
00187 i__2 = j - 1;
00188 for (i__ = 1; i__ <= i__2; ++i__) {
00189 x[ix] += temp * a_ref(i__, j);
00190 ix += *incx;
00191
00192 }
00193 if (nounit) {
00194 x[jx] *= a_ref(j, j);
00195 }
00196 }
00197 jx += *incx;
00198
00199 }
00200 }
00201 } else {
00202 if (*incx == 1) {
00203 for (j = *n; j >= 1; --j) {
00204 if (x[j] != 0.) {
00205 temp = x[j];
00206 i__1 = j + 1;
00207 for (i__ = *n; i__ >= i__1; --i__) {
00208 x[i__] += temp * a_ref(i__, j);
00209
00210 }
00211 if (nounit) {
00212 x[j] *= a_ref(j, j);
00213 }
00214 }
00215
00216 }
00217 } else {
00218 kx += (*n - 1) * *incx;
00219 jx = kx;
00220 for (j = *n; j >= 1; --j) {
00221 if (x[jx] != 0.) {
00222 temp = x[jx];
00223 ix = kx;
00224 i__1 = j + 1;
00225 for (i__ = *n; i__ >= i__1; --i__) {
00226 x[ix] += temp * a_ref(i__, j);
00227 ix -= *incx;
00228
00229 }
00230 if (nounit) {
00231 x[jx] *= a_ref(j, j);
00232 }
00233 }
00234 jx -= *incx;
00235
00236 }
00237 }
00238 }
00239 } else {
00240
00241 if (template_blas_lsame(uplo, "U")) {
00242 if (*incx == 1) {
00243 for (j = *n; j >= 1; --j) {
00244 temp = x[j];
00245 if (nounit) {
00246 temp *= a_ref(j, j);
00247 }
00248 for (i__ = j - 1; i__ >= 1; --i__) {
00249 temp += a_ref(i__, j) * x[i__];
00250
00251 }
00252 x[j] = temp;
00253
00254 }
00255 } else {
00256 jx = kx + (*n - 1) * *incx;
00257 for (j = *n; j >= 1; --j) {
00258 temp = x[jx];
00259 ix = jx;
00260 if (nounit) {
00261 temp *= a_ref(j, j);
00262 }
00263 for (i__ = j - 1; i__ >= 1; --i__) {
00264 ix -= *incx;
00265 temp += a_ref(i__, j) * x[ix];
00266
00267 }
00268 x[jx] = temp;
00269 jx -= *incx;
00270
00271 }
00272 }
00273 } else {
00274 if (*incx == 1) {
00275 i__1 = *n;
00276 for (j = 1; j <= i__1; ++j) {
00277 temp = x[j];
00278 if (nounit) {
00279 temp *= a_ref(j, j);
00280 }
00281 i__2 = *n;
00282 for (i__ = j + 1; i__ <= i__2; ++i__) {
00283 temp += a_ref(i__, j) * x[i__];
00284
00285 }
00286 x[j] = temp;
00287
00288 }
00289 } else {
00290 jx = kx;
00291 i__1 = *n;
00292 for (j = 1; j <= i__1; ++j) {
00293 temp = x[jx];
00294 ix = jx;
00295 if (nounit) {
00296 temp *= a_ref(j, j);
00297 }
00298 i__2 = *n;
00299 for (i__ = j + 1; i__ <= i__2; ++i__) {
00300 ix += *incx;
00301 temp += a_ref(i__, j) * x[ix];
00302
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