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