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_GEMM_HEADER
00038 #define TEMPLATE_BLAS_GEMM_HEADER
00039
00040 #include "template_blas_common.h"
00041
00042 template<class Treal>
00043 int template_blas_gemm(const char *transa, const char *transb, const integer *m, const integer *
00044 n, const integer *k, const Treal *alpha, const Treal *a, const integer *lda,
00045 const Treal *b, const integer *ldb, const Treal *beta, Treal *c__,
00046 const integer *ldc)
00047 {
00048
00049 integer a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, i__1, i__2,
00050 i__3;
00051
00052 integer info;
00053 logical nota, notb;
00054 Treal temp;
00055 integer i__, j, l;
00056 integer nrowa, nrowb;
00057 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
00058 #define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1]
00059 #define c___ref(a_1,a_2) c__[(a_2)*c_dim1 + a_1]
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
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152 a_dim1 = *lda;
00153 a_offset = 1 + a_dim1 * 1;
00154 a -= a_offset;
00155 b_dim1 = *ldb;
00156 b_offset = 1 + b_dim1 * 1;
00157 b -= b_offset;
00158 c_dim1 = *ldc;
00159 c_offset = 1 + c_dim1 * 1;
00160 c__ -= c_offset;
00161
00162 nota = template_blas_lsame(transa, "N");
00163 notb = template_blas_lsame(transb, "N");
00164 if (nota) {
00165 nrowa = *m;
00166 } else {
00167 nrowa = *k;
00168 }
00169 if (notb) {
00170 nrowb = *k;
00171 } else {
00172 nrowb = *n;
00173 }
00174
00175 info = 0;
00176 if (! nota && ! template_blas_lsame(transa, "C") && ! template_blas_lsame(
00177 transa, "T")) {
00178 info = 1;
00179 } else if (! notb && ! template_blas_lsame(transb, "C") && !
00180 template_blas_lsame(transb, "T")) {
00181 info = 2;
00182 } else if (*m < 0) {
00183 info = 3;
00184 } else if (*n < 0) {
00185 info = 4;
00186 } else if (*k < 0) {
00187 info = 5;
00188 } else if (*lda < maxMACRO(1,nrowa)) {
00189 info = 8;
00190 } else if (*ldb < maxMACRO(1,nrowb)) {
00191 info = 10;
00192 } else if (*ldc < maxMACRO(1,*m)) {
00193 info = 13;
00194 }
00195 if (info != 0) {
00196 template_blas_erbla("DGEMM ", &info);
00197 return 0;
00198 }
00199
00200 if (*m == 0 || *n == 0 || ( (*alpha == 0. || *k == 0) && *beta == 1.) ) {
00201 return 0;
00202 }
00203
00204 if (*alpha == 0.) {
00205 if (*beta == 0.) {
00206 i__1 = *n;
00207 for (j = 1; j <= i__1; ++j) {
00208 i__2 = *m;
00209 for (i__ = 1; i__ <= i__2; ++i__) {
00210 c___ref(i__, j) = 0.;
00211
00212 }
00213
00214 }
00215 } else {
00216 i__1 = *n;
00217 for (j = 1; j <= i__1; ++j) {
00218 i__2 = *m;
00219 for (i__ = 1; i__ <= i__2; ++i__) {
00220 c___ref(i__, j) = *beta * c___ref(i__, j);
00221
00222 }
00223
00224 }
00225 }
00226 return 0;
00227 }
00228
00229 if (notb) {
00230 if (nota) {
00231
00232 i__1 = *n;
00233 for (j = 1; j <= i__1; ++j) {
00234 if (*beta == 0.) {
00235 i__2 = *m;
00236 for (i__ = 1; i__ <= i__2; ++i__) {
00237 c___ref(i__, j) = 0.;
00238
00239 }
00240 } else if (*beta != 1.) {
00241 i__2 = *m;
00242 for (i__ = 1; i__ <= i__2; ++i__) {
00243 c___ref(i__, j) = *beta * c___ref(i__, j);
00244
00245 }
00246 }
00247 i__2 = *k;
00248 for (l = 1; l <= i__2; ++l) {
00249 if (b_ref(l, j) != 0.) {
00250 temp = *alpha * b_ref(l, j);
00251 i__3 = *m;
00252 for (i__ = 1; i__ <= i__3; ++i__) {
00253 c___ref(i__, j) = c___ref(i__, j) + temp * a_ref(
00254 i__, l);
00255
00256 }
00257 }
00258
00259 }
00260
00261 }
00262 } else {
00263
00264 i__1 = *n;
00265 for (j = 1; j <= i__1; ++j) {
00266 i__2 = *m;
00267 for (i__ = 1; i__ <= i__2; ++i__) {
00268 temp = 0.;
00269 i__3 = *k;
00270 for (l = 1; l <= i__3; ++l) {
00271 temp += a_ref(l, i__) * b_ref(l, j);
00272
00273 }
00274 if (*beta == 0.) {
00275 c___ref(i__, j) = *alpha * temp;
00276 } else {
00277 c___ref(i__, j) = *alpha * temp + *beta * c___ref(i__,
00278 j);
00279 }
00280
00281 }
00282
00283 }
00284 }
00285 } else {
00286 if (nota) {
00287
00288 i__1 = *n;
00289 for (j = 1; j <= i__1; ++j) {
00290 if (*beta == 0.) {
00291 i__2 = *m;
00292 for (i__ = 1; i__ <= i__2; ++i__) {
00293 c___ref(i__, j) = 0.;
00294
00295 }
00296 } else if (*beta != 1.) {
00297 i__2 = *m;
00298 for (i__ = 1; i__ <= i__2; ++i__) {
00299 c___ref(i__, j) = *beta * c___ref(i__, j);
00300
00301 }
00302 }
00303 i__2 = *k;
00304 for (l = 1; l <= i__2; ++l) {
00305 if (b_ref(j, l) != 0.) {
00306 temp = *alpha * b_ref(j, l);
00307 i__3 = *m;
00308 for (i__ = 1; i__ <= i__3; ++i__) {
00309 c___ref(i__, j) = c___ref(i__, j) + temp * a_ref(
00310 i__, l);
00311
00312 }
00313 }
00314
00315 }
00316
00317 }
00318 } else {
00319
00320 i__1 = *n;
00321 for (j = 1; j <= i__1; ++j) {
00322 i__2 = *m;
00323 for (i__ = 1; i__ <= i__2; ++i__) {
00324 temp = 0.;
00325 i__3 = *k;
00326 for (l = 1; l <= i__3; ++l) {
00327 temp += a_ref(l, i__) * b_ref(j, l);
00328
00329 }
00330 if (*beta == 0.) {
00331 c___ref(i__, j) = *alpha * temp;
00332 } else {
00333 c___ref(i__, j) = *alpha * temp + *beta * c___ref(i__,
00334 j);
00335 }
00336
00337 }
00338
00339 }
00340 }
00341 }
00342 return 0;
00343
00344 }
00345 #undef c___ref
00346 #undef b_ref
00347 #undef a_ref
00348
00349 #endif