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_TRMM_HEADER
00038 #define TEMPLATE_BLAS_TRMM_HEADER
00039
00040 #include "template_blas_common.h"
00041
00042 template<class Treal>
00043 int template_blas_trmm(const char *side, const char *uplo, const char *transa, const char *diag,
00044 const integer *m, const integer *n, const Treal *alpha, const Treal *a, const integer *
00045 lda, Treal *b, const integer *ldb)
00046 {
00047
00048 integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3;
00049
00050 integer info;
00051 Treal temp;
00052 integer i__, j, k;
00053 logical lside;
00054 integer nrowa;
00055 logical upper;
00056 logical nounit;
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
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 a_dim1 = *lda;
00144 a_offset = 1 + a_dim1 * 1;
00145 a -= a_offset;
00146 b_dim1 = *ldb;
00147 b_offset = 1 + b_dim1 * 1;
00148 b -= b_offset;
00149
00150 lside = template_blas_lsame(side, "L");
00151 if (lside) {
00152 nrowa = *m;
00153 } else {
00154 nrowa = *n;
00155 }
00156 nounit = template_blas_lsame(diag, "N");
00157 upper = template_blas_lsame(uplo, "U");
00158 info = 0;
00159 if (! lside && ! template_blas_lsame(side, "R")) {
00160 info = 1;
00161 } else if (! upper && ! template_blas_lsame(uplo, "L")) {
00162 info = 2;
00163 } else if (! template_blas_lsame(transa, "N") && ! template_blas_lsame(transa,
00164 "T") && ! template_blas_lsame(transa, "C")) {
00165 info = 3;
00166 } else if (! template_blas_lsame(diag, "U") && ! template_blas_lsame(diag,
00167 "N")) {
00168 info = 4;
00169 } else if (*m < 0) {
00170 info = 5;
00171 } else if (*n < 0) {
00172 info = 6;
00173 } else if (*lda < maxMACRO(1,nrowa)) {
00174 info = 9;
00175 } else if (*ldb < maxMACRO(1,*m)) {
00176 info = 11;
00177 }
00178 if (info != 0) {
00179 template_blas_erbla("TRMM ", &info);
00180 return 0;
00181 }
00182
00183 if (*n == 0) {
00184 return 0;
00185 }
00186
00187 if (*alpha == 0.) {
00188 i__1 = *n;
00189 for (j = 1; j <= i__1; ++j) {
00190 i__2 = *m;
00191 for (i__ = 1; i__ <= i__2; ++i__) {
00192 b_ref(i__, j) = 0.;
00193
00194 }
00195
00196 }
00197 return 0;
00198 }
00199
00200 if (lside) {
00201 if (template_blas_lsame(transa, "N")) {
00202
00203 if (upper) {
00204 i__1 = *n;
00205 for (j = 1; j <= i__1; ++j) {
00206 i__2 = *m;
00207 for (k = 1; k <= i__2; ++k) {
00208 if (b_ref(k, j) != 0.) {
00209 temp = *alpha * b_ref(k, j);
00210 i__3 = k - 1;
00211 for (i__ = 1; i__ <= i__3; ++i__) {
00212 b_ref(i__, j) = b_ref(i__, j) + temp * a_ref(
00213 i__, k);
00214
00215 }
00216 if (nounit) {
00217 temp *= a_ref(k, k);
00218 }
00219 b_ref(k, j) = temp;
00220 }
00221
00222 }
00223
00224 }
00225 } else {
00226 i__1 = *n;
00227 for (j = 1; j <= i__1; ++j) {
00228 for (k = *m; k >= 1; --k) {
00229 if (b_ref(k, j) != 0.) {
00230 temp = *alpha * b_ref(k, j);
00231 b_ref(k, j) = temp;
00232 if (nounit) {
00233 b_ref(k, j) = b_ref(k, j) * a_ref(k, k);
00234 }
00235 i__2 = *m;
00236 for (i__ = k + 1; i__ <= i__2; ++i__) {
00237 b_ref(i__, j) = b_ref(i__, j) + temp * a_ref(
00238 i__, k);
00239
00240 }
00241 }
00242
00243 }
00244
00245 }
00246 }
00247 } else {
00248
00249 if (upper) {
00250 i__1 = *n;
00251 for (j = 1; j <= i__1; ++j) {
00252 for (i__ = *m; i__ >= 1; --i__) {
00253 temp = b_ref(i__, j);
00254 if (nounit) {
00255 temp *= a_ref(i__, i__);
00256 }
00257 i__2 = i__ - 1;
00258 for (k = 1; k <= i__2; ++k) {
00259 temp += a_ref(k, i__) * b_ref(k, j);
00260
00261 }
00262 b_ref(i__, j) = *alpha * temp;
00263
00264 }
00265
00266 }
00267 } else {
00268 i__1 = *n;
00269 for (j = 1; j <= i__1; ++j) {
00270 i__2 = *m;
00271 for (i__ = 1; i__ <= i__2; ++i__) {
00272 temp = b_ref(i__, j);
00273 if (nounit) {
00274 temp *= a_ref(i__, i__);
00275 }
00276 i__3 = *m;
00277 for (k = i__ + 1; k <= i__3; ++k) {
00278 temp += a_ref(k, i__) * b_ref(k, j);
00279
00280 }
00281 b_ref(i__, j) = *alpha * temp;
00282
00283 }
00284
00285 }
00286 }
00287 }
00288 } else {
00289 if (template_blas_lsame(transa, "N")) {
00290
00291 if (upper) {
00292 for (j = *n; j >= 1; --j) {
00293 temp = *alpha;
00294 if (nounit) {
00295 temp *= a_ref(j, j);
00296 }
00297 i__1 = *m;
00298 for (i__ = 1; i__ <= i__1; ++i__) {
00299 b_ref(i__, j) = temp * b_ref(i__, j);
00300
00301 }
00302 i__1 = j - 1;
00303 for (k = 1; k <= i__1; ++k) {
00304 if (a_ref(k, j) != 0.) {
00305 temp = *alpha * a_ref(k, j);
00306 i__2 = *m;
00307 for (i__ = 1; i__ <= i__2; ++i__) {
00308 b_ref(i__, j) = b_ref(i__, j) + temp * b_ref(
00309 i__, k);
00310
00311 }
00312 }
00313
00314 }
00315
00316 }
00317 } else {
00318 i__1 = *n;
00319 for (j = 1; j <= i__1; ++j) {
00320 temp = *alpha;
00321 if (nounit) {
00322 temp *= a_ref(j, j);
00323 }
00324 i__2 = *m;
00325 for (i__ = 1; i__ <= i__2; ++i__) {
00326 b_ref(i__, j) = temp * b_ref(i__, j);
00327
00328 }
00329 i__2 = *n;
00330 for (k = j + 1; k <= i__2; ++k) {
00331 if (a_ref(k, j) != 0.) {
00332 temp = *alpha * a_ref(k, j);
00333 i__3 = *m;
00334 for (i__ = 1; i__ <= i__3; ++i__) {
00335 b_ref(i__, j) = b_ref(i__, j) + temp * b_ref(
00336 i__, k);
00337
00338 }
00339 }
00340
00341 }
00342
00343 }
00344 }
00345 } else {
00346
00347 if (upper) {
00348 i__1 = *n;
00349 for (k = 1; k <= i__1; ++k) {
00350 i__2 = k - 1;
00351 for (j = 1; j <= i__2; ++j) {
00352 if (a_ref(j, k) != 0.) {
00353 temp = *alpha * a_ref(j, k);
00354 i__3 = *m;
00355 for (i__ = 1; i__ <= i__3; ++i__) {
00356 b_ref(i__, j) = b_ref(i__, j) + temp * b_ref(
00357 i__, k);
00358
00359 }
00360 }
00361
00362 }
00363 temp = *alpha;
00364 if (nounit) {
00365 temp *= a_ref(k, k);
00366 }
00367 if (temp != 1.) {
00368 i__2 = *m;
00369 for (i__ = 1; i__ <= i__2; ++i__) {
00370 b_ref(i__, k) = temp * b_ref(i__, k);
00371
00372 }
00373 }
00374
00375 }
00376 } else {
00377 for (k = *n; k >= 1; --k) {
00378 i__1 = *n;
00379 for (j = k + 1; j <= i__1; ++j) {
00380 if (a_ref(j, k) != 0.) {
00381 temp = *alpha * a_ref(j, k);
00382 i__2 = *m;
00383 for (i__ = 1; i__ <= i__2; ++i__) {
00384 b_ref(i__, j) = b_ref(i__, j) + temp * b_ref(
00385 i__, k);
00386
00387 }
00388 }
00389
00390 }
00391 temp = *alpha;
00392 if (nounit) {
00393 temp *= a_ref(k, k);
00394 }
00395 if (temp != 1.) {
00396 i__1 = *m;
00397 for (i__ = 1; i__ <= i__1; ++i__) {
00398 b_ref(i__, k) = temp * b_ref(i__, k);
00399
00400 }
00401 }
00402
00403 }
00404 }
00405 }
00406 }
00407 return 0;
00408
00409 }
00410 #undef b_ref
00411 #undef a_ref
00412
00413 #endif