src/gromacs/linearalgebra/gmx_blas/dsyr2k.cpp

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