DSDP
|
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 /* Subroutine */ 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 /* System generated locals */ 00019 integer a_dim1, a_offset, b_dim1, b_offset, i1, i2, i3; 00020 /* Local variables */ 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 /* Function Body */ 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 /* Quick return if possible. */ 00070 if (*n == 0) { 00071 return 0; 00072 } 00073 /* And when alpha.eq.zero. */ 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 /* L10: */ 00081 } 00082 /* L20: */ 00083 } 00084 return 0; 00085 } 00086 /* Start the operations. */ 00087 if (lside) { 00088 if (lsame2(transa, "N")) { 00089 /* Form B := alpha*inv( A )*B. */ 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 /* L30: */ 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 /* L40: */ 00111 } 00112 } 00113 /* L50: */ 00114 } 00115 /* L60: */ 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 /* L70: */ 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 /* L80: */ 00139 } 00140 } 00141 /* L90: */ 00142 } 00143 /* L100: */ 00144 } 00145 } 00146 } else { 00147 /* Form B := alpha*inv( A' )*B. */ 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 /* L110: */ 00159 } 00160 if (nounit) { 00161 temp /= a_ref(i, i); 00162 } 00163 b_ref(i, j) = temp; 00164 /* L120: */ 00165 } 00166 /* L130: */ 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 /* L140: */ 00177 } 00178 if (nounit) { 00179 temp /= a_ref(i, i); 00180 } 00181 b_ref(i, j) = temp; 00182 /* L150: */ 00183 } 00184 /* L160: */ 00185 } 00186 } 00187 } 00188 } else { 00189 if (lsame2(transa, "N")) { 00190 /* Form B := alpha*B*inv( A ). */ 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 /* L170: */ 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 /* L180: */ 00209 } 00210 } 00211 /* L190: */ 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 /* L200: */ 00219 } 00220 } 00221 /* L210: */ 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 /* L220: */ 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 /* L230: */ 00240 } 00241 } 00242 /* L240: */ 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 /* L250: */ 00250 } 00251 } 00252 /* L260: */ 00253 } 00254 } 00255 } else { 00256 /* Form B := alpha*B*inv( A' ). */ 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 /* L270: */ 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 /* L280: */ 00275 } 00276 } 00277 /* L290: */ 00278 } 00279 if (*alpha != 1.) { 00280 i1 = *m; 00281 for (i = 1; i <= i1; ++i) { 00282 b_ref(i, k) *= alpha2; 00283 /* L300: */ 00284 } 00285 } 00286 /* L310: */ 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 /* L320: */ 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 /* L330: */ 00307 } 00308 } 00309 /* L340: */ 00310 } 00311 if (*alpha != 1.) { 00312 i2 = *m; 00313 for (i = 1; i <= i2; ++i) { 00314 b_ref(i, k) *= alpha2; 00315 /* L350: */ 00316 } 00317 } 00318 /* L360: */ 00319 } 00320 } 00321 } 00322 } 00323 return 0; 00324 /* End of DTRSM . */ 00325 } /* dtrsm_ */ 00326 #undef b_ref 00327 #undef a_ref 00328