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