src/gmxlib/gmx_blas/ssyr2k.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012, by the GROMACS development team, led by
   5  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   6  * others, as listed in the AUTHORS file in the top-level source
   7  * directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 #include <ctype.h>
  36 #include <math.h>
  37
  38 #include <types/simple.h>
  39 #include "gmx_blas.h"
  40
  41 void
  42 F77_FUNC(ssyr2k,SSYR2K)(const char *uplo,
  43                         const char *trans,
  44                         int *n__,
  45                         int *k__,
  46                         float *alpha__,
  47                         float *a,
  48                         int *lda__,
  49                         float *b,
  50                         int *ldb__,
  51                         float *beta__,
  52                         float *c,
  53                         int *ldc__)
  54 {
  55   char ch1,ch2;
  56   int nrowa;
  57   int i,j,l;
  58   float temp1,temp2;
  59
  60   int n = *n__;
  61   int k = *k__;
  62   int lda = *lda__;
  63   int ldb = *ldb__;
  64   int ldc = *ldc__;
  65
  66   float alpha = *alpha__;
  67   float beta  = *beta__;
  68
  69   ch1 = toupper(*uplo);
  70   ch2 = toupper(*trans);
  71
  72   if(ch2 == 'N')
  73     nrowa = n;
  74   else
  75     nrowa = k;
  76
  77   if(n==0 || ( ( fabs(alpha)<GMX_FLOAT_MIN || k==0 ) && fabs(beta-1.0)<GMX_FLOAT_EPS))
  78     return;
  79
  80   if(fabs(alpha)<GMX_FLOAT_MIN) {
  81     if(ch1=='U') {
  82       if(fabs(beta)<GMX_FLOAT_MIN)
  83         for(j=1;j<=n;j++)
  84           for(i=1;i<=j;i++)
  85             c[(j-1)*(ldc)+(i-1)] = 0.0;
  86       else
  87         for(j=1;j<=n;j++)
  88           for(i=1;i<=j;i++)
  89             c[(j-1)*(ldc)+(i-1)] *= beta;
  90     } else {
  91       /* lower */
  92       if(fabs(beta)<GMX_FLOAT_MIN)
  93         for(j=1;j<=n;j++)
  94           for(i=j;i<=n;i++)
  95             c[(j-1)*(ldc)+(i-1)] = 0.0;
  96       else
  97         for(j=1;j<=n;j++)
  98           for(i=j;i<=n;i++)
  99             c[(j-1)*(ldc)+(i-1)] *= beta;
 100     }
 101     return;
 102   }
 103
 104   if(ch2=='N') {
 105     if(ch1=='U') {
 106       for(j=1;j<=n;j++) {
 107         if(fabs(beta)<GMX_FLOAT_MIN)
 108           for(i=1;i<=j;i++)
 109              c[(j-1)*(ldc)+(i-1)] = 0.0;
 110         else if(fabs(beta-1.0)>GMX_FLOAT_EPS)
 111           for(i=1;i<=j;i++)
 112             c[(j-1)*(ldc)+(i-1)] *= beta;
 113         for(l=1;l<=k;l++) {
 114           if( fabs(a[(l-1)*(lda)+(j-1)])>GMX_FLOAT_MIN ||
 115               fabs(b[(l-1)*(ldb)+(j-1)])>GMX_FLOAT_MIN) {
 116             temp1 = alpha * b[(l-1)*(ldb)+(j-1)];
 117             temp2 = alpha * a[(l-1)*(lda)+(j-1)];
 118             for(i=1;i<=j;i++)
 119               c[(j-1)*(ldc)+(i-1)] +=
 120                 a[(l-1)*(lda)+(i-1)] * temp1 +
 121                 b[(l-1)*(ldb)+(i-1)] * temp2;
 122           }
 123         }
 124       }
 125     } else {
 126       /* lower */
 127       for(j=1;j<=n;j++) {
 128         if(fabs(beta)<GMX_FLOAT_MIN)
 129           for(i=j;i<=n;i++)
 130             c[(j-1)*(ldc)+(i-1)] = 0.0;
 131         else if(fabs(beta-1.0)>GMX_FLOAT_EPS)
 132           for(i=j;i<=n;i++)
 133             c[(j-1)*(ldc)+(i-1)] *= beta;
 134         for(l=1;l<=k;l++) {
 135           if( fabs(a[(l-1)*(lda)+(j-1)])>GMX_FLOAT_MIN ||
 136               fabs(b[(l-1)*(ldb)+(j-1)])>GMX_FLOAT_MIN) {
 137             temp1 = alpha * b[(l-1)*(ldb)+(j-1)];
 138             temp2 = alpha * a[(l-1)*(lda)+(j-1)];
 139             for(i=j;i<=n;i++)
 140               c[(j-1)*(ldc)+(i-1)] +=
 141                 a[(l-1)*(lda)+(i-1)] * temp1 +
 142                 b[(l-1)*(ldb)+(i-1)] * temp2;
 143           }
 144         }
 145       }
 146     }
 147   } else {
 148     /* transpose */
 149     if(ch1=='U') {
 150       for(j=1;j<=n;j++)
 151         for(i=1;i<=j;i++) {
 152           temp1 = 0.0;
 153           temp2 = 0.0;
 154           for (l=1;l<=k;l++) {
 155              temp1 += a[(i-1)*(lda)+(l-1)] * b[(j-1)*(ldb)+(l-1)];
 156              temp2 += b[(i-1)*(ldb)+(l-1)] * a[(j-1)*(lda)+(l-1)];
 157           }
 158           if(fabs(beta)<GMX_FLOAT_MIN)
 159             c[(j-1)*(ldc)+(i-1)] = alpha * (temp1 + temp2);
 160           else
 161             c[(j-1)*(ldc)+(i-1)] = beta * c[(j-1)*(ldc)+(i-1)] +
 162               alpha * (temp1 + temp2);
 163         }
 164     } else {
 165       /* lower */
 166       for(j=1;j<=n;j++)
 167         for(i=j;i<=n;i++) {
 168           temp1 = 0.0;
 169           temp2 = 0.0;
 170           for (l=1;l<=k;l++) {
 171              temp1 += a[(i-1)*(lda)+(l-1)] * b[(j-1)*(ldb)+(l-1)];
 172              temp2 += b[(i-1)*(ldb)+(l-1)] * a[(j-1)*(lda)+(l-1)];
 173           }
 174           if(fabs(beta)<GMX_FLOAT_MIN)
 175             c[(j-1)*(ldc)+(i-1)] = alpha * (temp1 + temp2);
 176           else
 177             c[(j-1)*(ldc)+(i-1)] = beta * c[(j-1)*(ldc)+(i-1)] +
 178               alpha * (temp1 + temp2);
 179         }
 180     }
 181   }
 182   return;
 183 }