src/gromacs/linearalgebra/gmx_blas/ssyr2.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(ssyr2,SSYR2)(const char *    uplo,
   9                       int *     n__,
  10                       float *  alpha__,
  11                       float *  x,
  12                       int *     incx__,
  13                       float *  y,
  14                       int *     incy__,
  15                       float *  a,
  16                       int *     lda__)
  17 {
  18     int kx,ky,ix,iy,jx,jy,j,i;
  19     float temp1,temp2;
  20     const char ch=toupper(*uplo);
  21
  22     int n = *n__;
  23     int lda = *lda__;
  24     int incx = *incx__;
  25     int incy = *incy__;
  26     float alpha = *alpha__;
  27
  28     if(n<=0 || fabs(alpha)<GMX_FLOAT_MIN || incx==0 || incy==0 ||
  29        (ch != 'U' && ch != 'L'))
  30         return;
  31
  32     jx = jy = kx = ky = 0;
  33
  34     /* init start points for non-unit increments */
  35     if(incx!=1 || incy!=1) {
  36         if(incx>0)
  37             kx = 1;
  38         else
  39             kx = 1 - (n - 1)*(incx);
  40         if(incy>0)
  41             ky = 1;
  42         else
  43             ky = 1 - (n - 1)*(incy);
  44
  45         jx = kx;
  46         jy = ky;
  47     }
  48
  49     if(ch == 'U') {
  50         /* Data in upper part of A */
  51         if(incx==1 && incy==1) {
  52             /* Unit increments for both x and y */
  53             for(j=1;j<=n;j++) {
  54                 if( fabs(x[j-1])>GMX_FLOAT_MIN  || fabs(y[j-1])>GMX_FLOAT_MIN ) {
  55                     temp1 = alpha * y[j-1];
  56                     temp2 = alpha * x[j-1];
  57                     for(i=1;i<=j;i++)
  58                         a[(j-1)*(lda)+(i-1)] += x[i-1]*temp1 + y[i-1]*temp2;
  59                 }
  60             }
  61         } else {
  62
  63             /* non-unit increments */
  64             for(j=1;j<=n;j++) {
  65
  66                 if( fabs(x[jx-1])>GMX_FLOAT_MIN || fabs(y[jy-1])>GMX_FLOAT_MIN ) {
  67                     temp1 = alpha * y[jy-1];
  68                     temp2 = alpha * x[jx-1];
  69                     ix = kx;
  70                     iy = ky;
  71                     for(i=1;i<=j;i++) {
  72                         a[(j-1)*(lda)+(i-1)] += x[ix-1]*temp1 + y[iy-1]*temp2;
  73                         ix += incx;
  74                         iy += incy;
  75                     }
  76                 }
  77                 jx += incx;
  78                 jy += incy;
  79             }
  80         }
  81     } else {
  82         /* Data in lower part of A */
  83         if(incx==1 && incy==1) {
  84             /* Unit increments for both x and y */
  85             for(j=1;j<=n;j++) {
  86                 if( fabs(x[j-1])>GMX_FLOAT_MIN  || fabs(y[j-1])>GMX_FLOAT_MIN ) {
  87                     temp1 = alpha * y[j-1];
  88                     temp2 = alpha * x[j-1];
  89                     for(i=j;i<=n;i++)
  90                         a[(j-1)*(lda)+(i-1)] += x[i-1]*temp1 + y[i-1]*temp2;
  91                 }
  92             }
  93         } else {
  94
  95             /* non-unit increments */
  96             for(j=1;j<=n;j++) {
  97
  98                 if( fabs(x[jx-1])>GMX_FLOAT_MIN || fabs(y[jy-1])>GMX_FLOAT_MIN ) {
  99                     temp1 = alpha * y[jy-1];
 100                     temp2 = alpha * x[jx-1];
 101                     ix = jx;
 102                     iy = jy;
 103                     for(i=j;i<=n;i++) {
 104                         a[(j-1)*(lda)+(i-1)] += x[ix-1]*temp1 + y[iy-1]*temp2;
 105                         ix += incx;
 106                         iy += incy;
 107                     }
 108                 }
 109                 jx += incx;
 110                 jy += incy;
 111             }
 112         }
 113     }
 114
 115     return;
 116 }