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