dtrsm2.c

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 

Generated on Wed May 21 19:31:19 2008 for DSDP by  doxygen 1.5.5