template_blas_gemv.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_GEMV_HEADER
00036 #define TEMPLATE_BLAS_GEMV_HEADER
00037
00038 #include "template_blas_common.h"
00039
00040 template<class Treal>
00041 int template_blas_gemv(const char *trans, const integer *m, const integer *n, const Treal *
00042 alpha, const Treal *a, const integer *lda, const Treal *x, const integer *incx,
00043 const Treal *beta, Treal *y, const integer *incy)
00044 {
00045
00046 integer a_dim1, a_offset, i__1, i__2;
00047
00048 integer info;
00049 Treal temp;
00050 integer lenx, leny, i__, j;
00051 integer ix, iy, jx, jy, kx, ky;
00052 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
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
00121
00122 a_dim1 = *lda;
00123 a_offset = 1 + a_dim1 * 1;
00124 a -= a_offset;
00125 --x;
00126 --y;
00127
00128 info = 0;
00129 if (! template_blas_lsame(trans, "N") && ! template_blas_lsame(trans, "T") && ! template_blas_lsame(trans, "C")
00130 ) {
00131 info = 1;
00132 } else if (*m < 0) {
00133 info = 2;
00134 } else if (*n < 0) {
00135 info = 3;
00136 } else if (*lda < maxMACRO(1,*m)) {
00137 info = 6;
00138 } else if (*incx == 0) {
00139 info = 8;
00140 } else if (*incy == 0) {
00141 info = 11;
00142 }
00143 if (info != 0) {
00144 template_blas_erbla("GEMV ", &info);
00145 return 0;
00146 }
00147
00148 if (*m == 0 || *n == 0 || (*alpha == 0. && *beta == 1.) ) {
00149 return 0;
00150 }
00151
00152
00153 if (template_blas_lsame(trans, "N")) {
00154 lenx = *n;
00155 leny = *m;
00156 } else {
00157 lenx = *m;
00158 leny = *n;
00159 }
00160 if (*incx > 0) {
00161 kx = 1;
00162 } else {
00163 kx = 1 - (lenx - 1) * *incx;
00164 }
00165 if (*incy > 0) {
00166 ky = 1;
00167 } else {
00168 ky = 1 - (leny - 1) * *incy;
00169 }
00170
00171
00172
00173 if (*beta != 1.) {
00174 if (*incy == 1) {
00175 if (*beta == 0.) {
00176 i__1 = leny;
00177 for (i__ = 1; i__ <= i__1; ++i__) {
00178 y[i__] = 0.;
00179
00180 }
00181 } else {
00182 i__1 = leny;
00183 for (i__ = 1; i__ <= i__1; ++i__) {
00184 y[i__] = *beta * y[i__];
00185
00186 }
00187 }
00188 } else {
00189 iy = ky;
00190 if (*beta == 0.) {
00191 i__1 = leny;
00192 for (i__ = 1; i__ <= i__1; ++i__) {
00193 y[iy] = 0.;
00194 iy += *incy;
00195
00196 }
00197 } else {
00198 i__1 = leny;
00199 for (i__ = 1; i__ <= i__1; ++i__) {
00200 y[iy] = *beta * y[iy];
00201 iy += *incy;
00202
00203 }
00204 }
00205 }
00206 }
00207 if (*alpha == 0.) {
00208 return 0;
00209 }
00210 if (template_blas_lsame(trans, "N")) {
00211
00212 jx = kx;
00213 if (*incy == 1) {
00214 i__1 = *n;
00215 for (j = 1; j <= i__1; ++j) {
00216 if (x[jx] != 0.) {
00217 temp = *alpha * x[jx];
00218 i__2 = *m;
00219 for (i__ = 1; i__ <= i__2; ++i__) {
00220 y[i__] += temp * a_ref(i__, j);
00221
00222 }
00223 }
00224 jx += *incx;
00225
00226 }
00227 } else {
00228 i__1 = *n;
00229 for (j = 1; j <= i__1; ++j) {
00230 if (x[jx] != 0.) {
00231 temp = *alpha * x[jx];
00232 iy = ky;
00233 i__2 = *m;
00234 for (i__ = 1; i__ <= i__2; ++i__) {
00235 y[iy] += temp * a_ref(i__, j);
00236 iy += *incy;
00237
00238 }
00239 }
00240 jx += *incx;
00241
00242 }
00243 }
00244 } else {
00245
00246 jy = ky;
00247 if (*incx == 1) {
00248 i__1 = *n;
00249 for (j = 1; j <= i__1; ++j) {
00250 temp = 0.;
00251 i__2 = *m;
00252 for (i__ = 1; i__ <= i__2; ++i__) {
00253 temp += a_ref(i__, j) * x[i__];
00254
00255 }
00256 y[jy] += *alpha * temp;
00257 jy += *incy;
00258
00259 }
00260 } else {
00261 i__1 = *n;
00262 for (j = 1; j <= i__1; ++j) {
00263 temp = 0.;
00264 ix = kx;
00265 i__2 = *m;
00266 for (i__ = 1; i__ <= i__2; ++i__) {
00267 temp += a_ref(i__, j) * x[ix];
00268 ix += *incx;
00269
00270 }
00271 y[jy] += *alpha * temp;
00272 jy += *incy;
00273
00274 }
00275 }
00276 }
00277 return 0;
00278
00279 }
00280 #undef a_ref
00281
00282 #endif