src/gromacs/linearalgebra/gmx_lapack/ssytd2.c

   1 #include <ctype.h>
   2 #include <math.h>
   3
   4 #include "gromacs/utility/real.h"
   5
   6 #include "../gmx_blas.h"
   7 #include "../gmx_lapack.h"
   8
   9 void
  10 F77_FUNC(ssytd2,SSYTD2)(const char *    uplo,
  11         int *     n,
  12         float *  a,
  13         int *     lda,
  14         float *  d,
  15         float *  e,
  16         float *  tau,
  17     int gmx_unused *     info)
  18 {
  19   float minusone,zero;
  20   float taui,alpha,tmp;
  21   int ti1,ti2,ti3,i;
  22   const char ch=toupper(*uplo);
  23
  24   zero = 0.0;
  25   minusone = -1.0;
  26
  27   if(*n<=0)
  28     return;
  29
  30   if(ch=='U') {
  31     for(i=*n-1;i>=1;i--) {
  32
  33       ti1 = 1;
  34       F77_FUNC(slarfg,SLARFG)(&i,&(a[i*(*lda)+(i-1)]),&(a[i*(*lda)+0]),&ti1,&taui);
  35       e[i-1] = a[i*(*lda) + (i-1)];
  36       if(fabs(taui)>GMX_FLOAT_MIN) {
  37         a[i*(*lda)+(i-1)] = 1.0;
  38
  39         ti1 = 1;
  40         F77_FUNC(ssymv,SSYMV)("U",&i,&taui,a,lda,&(a[i*(*lda)+0]),&ti1,&zero,tau,&ti1);
  41
  42         tmp = F77_FUNC(sdot,SDOT)(&i,tau,&ti1,&(a[i*(*lda)+0]),&ti1);
  43
  44         alpha = -0.5*taui*tmp;
  45
  46         F77_FUNC(saxpy,SAXPY)(&i,&alpha,&(a[i*(*lda)+0]),&ti1,tau,&ti1);
  47
  48         F77_FUNC(ssyr2,SSYR2)("U",&i,&minusone,&(a[i*(*lda)+0]),&ti1,tau,&ti1,a,lda);
  49
  50         a[i*(*lda)+(i-1)] = e[i-1];
  51
  52       }
  53       d[i] = a[i*(*lda)+i];
  54       tau[i-1] = taui;
  55     }
  56     d[0] = a[0];
  57
  58   } else {
  59     /* lower */
  60
  61     for(i=1;i<*n;i++) {
  62
  63       ti1 = *n - i;
  64       ti2 = ( *n < i+2) ? *n : i+2;
  65       ti3 = 1;
  66       F77_FUNC(slarfg,SLARFG)(&ti1,&(a[(i-1)*(*lda)+(i)]),&(a[(i-1)*(*lda)+ti2-1]),&ti3,&taui);
  67
  68       e[i-1] = a[(i-1)*(*lda) + (i)];
  69
  70       if(fabs(taui)>GMX_FLOAT_MIN) {
  71         a[(i-1)*(*lda)+(i)] = 1.0;
  72
  73         ti1 = *n - i;
  74         ti2 = 1;
  75         F77_FUNC(ssymv,SSYMV)(uplo,&ti1,&taui,&(a[i*(*lda)+i]),lda,&(a[(i-1)*(*lda)+i]),
  76                &ti2,&zero,&(tau[i-1]),&ti2);
  77
  78         tmp = F77_FUNC(sdot,SDOT)(&ti1,&(tau[i-1]),&ti2,&(a[(i-1)*(*lda)+i]),&ti2);
  79
  80         alpha = -0.5*taui*tmp;
  81
  82         F77_FUNC(saxpy,SAXPY)(&ti1,&alpha,&(a[(i-1)*(*lda)+i]),&ti2,&(tau[i-1]),&ti2);
  83
  84         F77_FUNC(ssyr2,SSYR2)(uplo,&ti1,&minusone,&(a[(i-1)*(*lda)+i]),&ti2,&(tau[i-1]),&ti2,
  85                &(a[(i)*(*lda)+i]),lda);
  86
  87         a[(i-1)*(*lda)+(i)] = e[i-1];
  88
  89       }
  90       d[i-1] = a[(i-1)*(*lda)+i-1];
  91       tau[i-1] = taui;
  92     }
  93     d[*n-1] = a[(*n-1)*(*lda)+(*n-1)];
  94
  95   }
  96   return;
  97 }