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