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