src/gromacs/linearalgebra/gmx_blas/strsm.c

   1 #include <ctype.h>
   2 #include <math.h>
   3
   4 #include "types/simple.h"
   5 #include "../gmx_blas.h"
   6
   7 void
   8 F77_FUNC(strsm,STRSM)(const char * side,
   9                       const char * uplo,
  10                       const char * transa,
  11                       const char * diag,
  12                       int *  m__,
  13                       int *  n__,
  14                       float *alpha__,
  15                       float *a,
  16                       int *  lda__,
  17                       float *b,
  18                       int *  ldb__)
  19 {
  20   const char xside  = toupper(*side);
  21   const char xuplo  = toupper(*uplo);
  22   const char xtrans = toupper(*transa);
  23   const char xdiag  = toupper(*diag);
  24   int i,j,k;
  25   float temp;
  26
  27   int m = *m__;
  28   int n = *n__;
  29   int lda = *lda__;
  30   int ldb = *ldb__;
  31   float alpha = *alpha__;
  32
  33   if(n<=0)
  34     return;
  35
  36
  37   if(fabs(alpha)<GMX_FLOAT_MIN) {
  38     for(j=0;j<n;j++)
  39       for(i=0;i<m;i++)
  40         b[j*(ldb)+i] = 0.0;
  41     return;
  42   }
  43
  44   if(xside=='L') {
  45     /* left side */
  46     if(xtrans=='N') {
  47       /* No transpose */
  48       if(xuplo=='U') {
  49         /* upper */
  50         for(j=0;j<n;j++) {
  51           if(fabs(alpha-1.0)>GMX_FLOAT_EPS) {
  52             for(i=0;i<m;i++)
  53               b[j*(ldb)+i] *= alpha;
  54           }
  55           for(k=m-1;k>=0;k--) {
  56             if( fabs(b[j*(ldb)+k])>GMX_FLOAT_MIN) {
  57               if(xdiag=='N')
  58                 b[j*(ldb)+k] /= a[k*(lda)+k];
  59               for(i=0;i<k;i++)
  60                 b[j*(ldb)+i] -= b[j*(ldb)+k]*a[k*(lda)+i];
  61             }
  62           }
  63         }
  64       } else {
  65         /* lower */
  66         for(j=0;j<n;j++) {
  67           if(fabs(alpha-1.0)>GMX_FLOAT_EPS)
  68             for(i=0;i<m;i++)
  69               b[j*(ldb)+i] *= alpha;
  70           for(k=0;k<m;k++) {
  71             if( fabs(b[j*(ldb)+k])>GMX_FLOAT_MIN) {
  72               if(xdiag=='N')
  73                 b[j*(ldb)+k] /= a[k*(lda)+k];
  74               for(i=k+1;i<m;i++)
  75                 b[j*(ldb)+i] -= b[j*(ldb)+k]*a[k*(lda)+i];
  76             }
  77           }
  78         }
  79       }
  80     } else {
  81       /* Transpose */
  82       if(xuplo=='U') {
  83         /* upper */
  84         for(j=0;j<n;j++) {
  85           for(i=0;i<m;i++) {
  86             temp = alpha * b[j*(ldb)+i];
  87             for(k=0;k<i;k++)
  88               temp -= a[i*(lda)+k] * b[j*(ldb)+k];
  89             if(xdiag=='N')
  90                 temp /= a[i*(lda)+i];
  91             b[j*(ldb)+i] = temp;
  92           }
  93         }
  94       } else {
  95         /* lower */
  96         for(j=0;j<n;j++) {
  97           for(i=m-1;i>=0;i--) {
  98             temp = alpha * b[j*(ldb)+i];
  99             for(k=i+1;k<m;k++)
 100               temp -= a[i*(lda)+k] * b[j*(ldb)+k];
 101             if(xdiag=='N')
 102                 temp /= a[i*(lda)+i];
 103             b[j*(ldb)+i] = temp;
 104           }
 105         }
 106       }
 107     }
 108   } else {
 109     /* right side */
 110     if(xtrans=='N') {
 111       /* No transpose */
 112       if(xuplo=='U') {
 113         /* upper */
 114         for(j=0;j<n;j++) {
 115           if(fabs(alpha-1.0)>GMX_FLOAT_EPS)
 116             for(i=0;i<m;i++)
 117               b[j*(ldb)+i] *= alpha;
 118           for(k=0;k<j;k++) {
 119             if( fabs(a[j*(lda)+k])>GMX_FLOAT_MIN) {
 120               for(i=0;i<m;i++)
 121                 b[j*(ldb)+i] -= a[j*(lda)+k]*b[k*(ldb)+i];
 122             }
 123           }
 124           if(xdiag=='N') {
 125             temp = 1.0/a[j*(lda)+j];
 126             for(i=0;i<m;i++)
 127               b[j*(ldb)+i] *= temp;
 128           }
 129         }
 130       } else {
 131         /* lower */
 132         for(j=n-1;j>=0;j--) {
 133           if(fabs(alpha-1.0)>GMX_FLOAT_EPS)
 134             for(i=0;i<m;i++)
 135               b[j*(ldb)+i] *= alpha;
 136           for(k=j+1;k<n;k++) {
 137             if( fabs(a[j*(lda)+k])>GMX_FLOAT_MIN ) {
 138               for(i=0;i<m;i++)
 139                 b[j*(ldb)+i] -= a[j*(lda)+k]*b[k*(ldb)+i];
 140             }
 141           }
 142           if(xdiag=='N') {
 143             temp = 1.0/a[j*(lda)+j];
 144             for(i=0;i<m;i++)
 145               b[j*(ldb)+i] *= temp;
 146           }
 147         }
 148       }
 149     } else {
 150       /* Transpose */
 151       if(xuplo=='U') {
 152         /* upper */
 153         for(k=n-1;k>=0;k--) {
 154           if(xdiag=='N') {
 155             temp = 1.0/a[k*(lda)+k];
 156             for(i=0;i<m;i++)
 157               b[k*(ldb)+i] *= temp;
 158           }
 159           for(j=0;j<k;j++) {
 160             if( fabs(a[k*(lda)+j])>GMX_FLOAT_MIN) {
 161               temp = a[k*(lda)+j];
 162               for(i=0;i<m;i++)
 163                 b[j*(ldb)+i] -= temp * b[k*(ldb)+i];
 164             }
 165           }
 166           if(fabs(alpha-1.0)>GMX_FLOAT_EPS)
 167             for(i=0;i<m;i++)
 168               b[k*(ldb)+i] *= alpha;
 169         }
 170       } else {
 171         /* lower */
 172         for(k=0;k<n;k++) {
 173           if(xdiag=='N') {
 174             temp = 1.0/a[k*(lda)+k];
 175             for(i=0;i<m;i++)
 176               b[k*(ldb)+i] *= temp;
 177           }
 178           for(j=k+1;j<n;j++) {
 179             if( fabs(a[k*(lda)+j])>GMX_FLOAT_MIN) {
 180               temp = a[k*(lda)+j];
 181               for(i=0;i<m;i++)
 182                 b[j*(ldb)+i] -= temp * b[k*(ldb)+i];
 183             }
 184           }
 185           if(fabs(alpha-1.0)>GMX_FLOAT_EPS)
 186             for(i=0;i<m;i++)
 187               b[k*(ldb)+i] *= alpha;
 188         }
 189       }
 190     }
 191   }
 192 }