src/gromacs/linearalgebra/gmx_blas/ssyr2k.c

   1 #include <ctype.h>
   2 #include <math.h>
   3
   4 #include "gromacs/utility/real.h"
   5 #include "../gmx_blas.h"
   6
   7 void
   8 F77_FUNC(ssyr2k,SSYR2K)(const char *uplo,
   9                         const char *trans,
  10                         int *n__,
  11                         int *k__,
  12                         float *alpha__,
  13                         float *a,
  14                         int *lda__,
  15                         float *b,
  16                         int *ldb__,
  17                         float *beta__,
  18                         float *c,
  19                         int *ldc__)
  20 {
  21   char ch1,ch2;
  22   int nrowa;
  23   int i,j,l;
  24   float temp1,temp2;
  25
  26   int n = *n__;
  27   int k = *k__;
  28   int lda = *lda__;
  29   int ldb = *ldb__;
  30   int ldc = *ldc__;
  31
  32   float alpha = *alpha__;
  33   float beta  = *beta__;
  34
  35   ch1 = toupper(*uplo);
  36   ch2 = toupper(*trans);
  37
  38   if(ch2 == 'N')
  39     nrowa = n;
  40   else
  41     nrowa = k;
  42
  43   if(n==0 || ( ( fabs(alpha)<GMX_FLOAT_MIN || k==0 ) && fabs(beta-1.0)<GMX_FLOAT_EPS))
  44     return;
  45
  46   if(fabs(alpha)<GMX_FLOAT_MIN) {
  47     if(ch1=='U') {
  48       if(fabs(beta)<GMX_FLOAT_MIN)
  49         for(j=1;j<=n;j++)
  50           for(i=1;i<=j;i++)
  51             c[(j-1)*(ldc)+(i-1)] = 0.0;
  52       else
  53         for(j=1;j<=n;j++)
  54           for(i=1;i<=j;i++)
  55             c[(j-1)*(ldc)+(i-1)] *= beta;
  56     } else {
  57       /* lower */
  58       if(fabs(beta)<GMX_FLOAT_MIN)
  59         for(j=1;j<=n;j++)
  60           for(i=j;i<=n;i++)
  61             c[(j-1)*(ldc)+(i-1)] = 0.0;
  62       else
  63         for(j=1;j<=n;j++)
  64           for(i=j;i<=n;i++)
  65             c[(j-1)*(ldc)+(i-1)] *= beta;
  66     }
  67     return;
  68   }
  69
  70   if(ch2=='N') {
  71     if(ch1=='U') {
  72       for(j=1;j<=n;j++) {
  73         if(fabs(beta)<GMX_FLOAT_MIN)
  74           for(i=1;i<=j;i++)
  75              c[(j-1)*(ldc)+(i-1)] = 0.0;
  76         else if(fabs(beta-1.0)>GMX_FLOAT_EPS)
  77           for(i=1;i<=j;i++)
  78             c[(j-1)*(ldc)+(i-1)] *= beta;
  79         for(l=1;l<=k;l++) {
  80           if( fabs(a[(l-1)*(lda)+(j-1)])>GMX_FLOAT_MIN ||
  81               fabs(b[(l-1)*(ldb)+(j-1)])>GMX_FLOAT_MIN) {
  82             temp1 = alpha * b[(l-1)*(ldb)+(j-1)];
  83             temp2 = alpha * a[(l-1)*(lda)+(j-1)];
  84             for(i=1;i<=j;i++)
  85               c[(j-1)*(ldc)+(i-1)] +=
  86                 a[(l-1)*(lda)+(i-1)] * temp1 +
  87                 b[(l-1)*(ldb)+(i-1)] * temp2;
  88           }
  89         }
  90       }
  91     } else {
  92       /* lower */
  93       for(j=1;j<=n;j++) {
  94         if(fabs(beta)<GMX_FLOAT_MIN)
  95           for(i=j;i<=n;i++)
  96             c[(j-1)*(ldc)+(i-1)] = 0.0;
  97         else if(fabs(beta-1.0)>GMX_FLOAT_EPS)
  98           for(i=j;i<=n;i++)
  99             c[(j-1)*(ldc)+(i-1)] *= beta;
 100         for(l=1;l<=k;l++) {
 101           if( fabs(a[(l-1)*(lda)+(j-1)])>GMX_FLOAT_MIN ||
 102               fabs(b[(l-1)*(ldb)+(j-1)])>GMX_FLOAT_MIN) {
 103             temp1 = alpha * b[(l-1)*(ldb)+(j-1)];
 104             temp2 = alpha * a[(l-1)*(lda)+(j-1)];
 105             for(i=j;i<=n;i++)
 106               c[(j-1)*(ldc)+(i-1)] +=
 107                 a[(l-1)*(lda)+(i-1)] * temp1 +
 108                 b[(l-1)*(ldb)+(i-1)] * temp2;
 109           }
 110         }
 111       }
 112     }
 113   } else {
 114     /* transpose */
 115     if(ch1=='U') {
 116       for(j=1;j<=n;j++)
 117         for(i=1;i<=j;i++) {
 118           temp1 = 0.0;
 119           temp2 = 0.0;
 120           for (l=1;l<=k;l++) {
 121              temp1 += a[(i-1)*(lda)+(l-1)] * b[(j-1)*(ldb)+(l-1)];
 122              temp2 += b[(i-1)*(ldb)+(l-1)] * a[(j-1)*(lda)+(l-1)];
 123           }
 124           if(fabs(beta)<GMX_FLOAT_MIN)
 125             c[(j-1)*(ldc)+(i-1)] = alpha * (temp1 + temp2);
 126           else
 127             c[(j-1)*(ldc)+(i-1)] = beta * c[(j-1)*(ldc)+(i-1)] +
 128               alpha * (temp1 + temp2);
 129         }
 130     } else {
 131       /* lower */
 132       for(j=1;j<=n;j++)
 133         for(i=j;i<=n;i++) {
 134           temp1 = 0.0;
 135           temp2 = 0.0;
 136           for (l=1;l<=k;l++) {
 137              temp1 += a[(i-1)*(lda)+(l-1)] * b[(j-1)*(ldb)+(l-1)];
 138              temp2 += b[(i-1)*(ldb)+(l-1)] * a[(j-1)*(lda)+(l-1)];
 139           }
 140           if(fabs(beta)<GMX_FLOAT_MIN)
 141             c[(j-1)*(ldc)+(i-1)] = alpha * (temp1 + temp2);
 142           else
 143             c[(j-1)*(ldc)+(i-1)] = beta * c[(j-1)*(ldc)+(i-1)] +
 144               alpha * (temp1 + temp2);
 145         }
 146     }
 147   }
 148   return;
 149 }