00001 #include "dsdplapack.h"
00002
00003 typedef long int integer;
00004 typedef long int logical;
00005
00006 #define max(a,b) ((a) >= (b) ? (a) : (b))
00007 #define dmax(a,b) (double)max(a,b)
00008
00009 static int lsame2(const char c1[], const char c2[]){
00010 if (c1[0]==c2[0]) return 1;
00011 else return 0;
00012 }
00013
00014 int dtrsm2(char *side, char *uplo, char *transa, char *diag,
00015 integer *m, integer *n, double *alpha, double *a, integer *
00016 lda, double *b, integer *ldb)
00017 {
00018
00019 integer a_dim1, a_offset, b_dim1, b_offset, i1, i2, i3;
00020
00021 static integer info;
00022 static double temp,alpha2,aref;
00023 static integer i, j, k;
00024 static logical lside;
00025 static integer nrowa;
00026 static logical upper;
00027 static logical nounit;
00028 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
00029 #define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1]
00030
00031 a_dim1 = *lda;
00032 a_offset = 1 + a_dim1 * 1;
00033 a -= a_offset;
00034 b_dim1 = *ldb;
00035 b_offset = 1 + b_dim1 * 1;
00036 b -= b_offset;
00037
00038 lside = lsame2(side, "L");
00039 if (lside) {
00040 nrowa = *m;
00041 } else {
00042 nrowa = *n;
00043 }
00044 nounit = lsame2(diag, "N");
00045 upper = lsame2(uplo, "U");
00046 info = 0;
00047 if (! lside && ! lsame2(side, "R")) {
00048 info = 1;
00049 } else if (! upper && ! lsame2(uplo, "L")) {
00050 info = 2;
00051 } else if (! lsame2(transa, "N") && ! lsame2(transa,
00052 "T") && ! lsame2(transa, "C")) {
00053 info = 3;
00054 } else if (! lsame2(diag, "U") && ! lsame2(diag,
00055 "N")) {
00056 info = 4;
00057 } else if (*m < 0) {
00058 info = 5;
00059 } else if (*n < 0) {
00060 info = 6;
00061 } else if (*lda < max(1,nrowa)) {
00062 info = 9;
00063 } else if (*ldb < max(1,*m)) {
00064 info = 11;
00065 }
00066 if (info != 0) {
00067 return info;
00068 }
00069
00070 if (*n == 0) {
00071 return 0;
00072 }
00073
00074 if (*alpha == 0.) {
00075 i1 = *n;
00076 for (j = 1; j <= i1; ++j) {
00077 i2 = *m;
00078 for (i = 1; i <= i2; ++i) {
00079 b_ref(i, j) = 0.;
00080
00081 }
00082
00083 }
00084 return 0;
00085 }
00086
00087 if (lside) {
00088 if (lsame2(transa, "N")) {
00089
00090 if (upper) {
00091 i1 = *n;
00092 for (j = 1; j <= i1; ++j) {
00093 if (*alpha != 1.) {
00094 i2 = *m;
00095 alpha2=*alpha;
00096 for (i = 1; i <= i2; ++i) {
00097 b_ref(i, j) *= alpha2;
00098
00099 }
00100 }
00101 for (k = *m; k >= 1; --k) {
00102 if (b_ref(k, j) != 0.) {
00103 if (nounit) {
00104 b_ref(k, j) /= a_ref(k, k);
00105 }
00106 i2 = k - 1;
00107 aref=-b_ref(k, j);
00108 for (i = 1; i <= i2; ++i) {
00109 b_ref(i, j) += aref * a_ref(i, k);
00110
00111 }
00112 }
00113
00114 }
00115
00116 }
00117 } else {
00118 i1 = *n;
00119 for (j = 1; j <= i1; ++j) {
00120 if (*alpha != 1.) {
00121 i2 = *m;
00122 alpha2=*alpha;
00123 for (i = 1; i <= i2; ++i) {
00124 b_ref(i, j) *= alpha2;
00125
00126 }
00127 }
00128 i2 = *m;
00129 for (k = 1; k <= i2; ++k) {
00130 if (b_ref(k, j) != 0.) {
00131 if (nounit) {
00132 b_ref(k, j) /= a_ref(k, k);
00133 }
00134 i3 = *m;
00135 aref= -b_ref(k, j);
00136 for (i = k + 1; i <= i3; ++i) {
00137 b_ref(i, j) += aref * a_ref(i, k);
00138
00139 }
00140 }
00141
00142 }
00143
00144 }
00145 }
00146 } else {
00147
00148 if (upper) {
00149 i1 = *n;
00150 for (j = 1; j <= i1; ++j) {
00151 i2 = *m;
00152 alpha2=*alpha;
00153 for (i = 1; i <= i2; ++i) {
00154 temp = alpha2 * b_ref(i, j);
00155 i3 = i - 1;
00156 for (k = 1; k <= i3; ++k) {
00157 temp -= a_ref(k, i) * b_ref(k, j);
00158
00159 }
00160 if (nounit) {
00161 temp /= a_ref(i, i);
00162 }
00163 b_ref(i, j) = temp;
00164
00165 }
00166
00167 }
00168 } else {
00169 i1 = *n;
00170 for (j = 1; j <= i1; ++j) {
00171 for (i = *m; i >= 1; --i) {
00172 temp = alpha2 * b_ref(i, j);
00173 i2 = *m;
00174 for (k = i + 1; k <= i2; ++k) {
00175 temp -= a_ref(k, i) * b_ref(k, j);
00176
00177 }
00178 if (nounit) {
00179 temp /= a_ref(i, i);
00180 }
00181 b_ref(i, j) = temp;
00182
00183 }
00184
00185 }
00186 }
00187 }
00188 } else {
00189 if (lsame2(transa, "N")) {
00190
00191 if (upper) {
00192 i1 = *n;
00193 for (j = 1; j <= i1; ++j) {
00194 if (*alpha != 1.) {
00195 i2 = *m;
00196 for (i = 1; i <= i2; ++i) {
00197 b_ref(i, j) *= alpha2;
00198
00199 }
00200 }
00201 i2 = j - 1;
00202 for (k = 1; k <= i2; ++k) {
00203 if (a_ref(k, j) != 0.) {
00204 i3 = *m;
00205 aref=-a_ref(k, j);
00206 for (i = 1; i <= i3; ++i) {
00207 b_ref(i, j) += aref * b_ref(i, k);
00208
00209 }
00210 }
00211
00212 }
00213 if (nounit) {
00214 temp = 1. / a_ref(j, j);
00215 i2 = *m;
00216 for (i = 1; i <= i2; ++i) {
00217 b_ref(i, j) *= temp;
00218
00219 }
00220 }
00221
00222 }
00223 } else {
00224 for (j = *n; j >= 1; --j) {
00225 if (*alpha != 1.) {
00226 i1 = *m;
00227 for (i = 1; i <= i1; ++i) {
00228 b_ref(i, j) *= alpha2;
00229
00230 }
00231 }
00232 i1 = *n;
00233 for (k = j + 1; k <= i1; ++k) {
00234 if (a_ref(k, j) != 0.) {
00235 i2 = *m;
00236 aref=-a_ref(k, j);
00237 for (i = 1; i <= i2; ++i) {
00238 b_ref(i, j) += aref * b_ref(i, k);
00239
00240 }
00241 }
00242
00243 }
00244 if (nounit) {
00245 temp = 1. / a_ref(j, j);
00246 i1 = *m;
00247 for (i = 1; i <= i1; ++i) {
00248 b_ref(i, j) *= temp;
00249
00250 }
00251 }
00252
00253 }
00254 }
00255 } else {
00256
00257 if (upper) {
00258 for (k = *n; k >= 1; --k) {
00259 if (nounit) {
00260 temp = 1. / a_ref(k, k);
00261 i1 = *m;
00262 for (i = 1; i <= i1; ++i) {
00263 b_ref(i, k) *= temp;
00264
00265 }
00266 }
00267 i1 = k - 1;
00268 for (j = 1; j <= i1; ++j) {
00269 if (a_ref(j, k) != 0.) {
00270 temp = a_ref(j, k);
00271 i2 = *m;
00272 for (i = 1; i <= i2; ++i) {
00273 b_ref(i, j) -= temp * b_ref(i, k);
00274
00275 }
00276 }
00277
00278 }
00279 if (*alpha != 1.) {
00280 i1 = *m;
00281 for (i = 1; i <= i1; ++i) {
00282 b_ref(i, k) *= alpha2;
00283
00284 }
00285 }
00286
00287 }
00288 } else {
00289 i1 = *n;
00290 for (k = 1; k <= i1; ++k) {
00291 if (nounit) {
00292 temp = 1. / a_ref(k, k);
00293 i2 = *m;
00294 for (i = 1; i <= i2; ++i) {
00295 b_ref(i, k) *= temp;
00296
00297 }
00298 }
00299 i2 = *n;
00300 for (j = k + 1; j <= i2; ++j) {
00301 if (a_ref(j, k) != 0.) {
00302 temp = a_ref(j, k);
00303 i3 = *m;
00304 for (i = 1; i <= i3; ++i) {
00305 b_ref(i, j) -= temp * b_ref(i, k);
00306
00307 }
00308 }
00309
00310 }
00311 if (*alpha != 1.) {
00312 i2 = *m;
00313 for (i = 1; i <= i2; ++i) {
00314 b_ref(i, k) *= alpha2;
00315
00316 }
00317 }
00318
00319 }
00320 }
00321 }
00322 }
00323 return 0;
00324
00325 }
00326 #undef b_ref
00327 #undef a_ref
00328