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