src/gromacs/linearalgebra/gmx_blas/dsyr2k.c

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