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