src/gmxlib/gmx_lapack/sgebd2.c

   1 #include "gmx_lapack.h"
   2
   3 void
   4 F77_FUNC(sgebd2,SGEBD2)(int *m,
   5         int *n,
   6         float *a,
   7         int *lda,
   8         float *d,
   9         float *e,
  10         float *tauq,
  11         float *taup,
  12         float *work,
  13         int *info)
  14 {
  15
  16   int i,i1,i2,i3;
  17
  18     *info = 0;
  19
  20   if(*m>=*n) {
  21     /* reduce to upper bidiag. form */
  22     for(i=0;i<*n;i++) {
  23       i1 = *m - i;
  24       i2 = ( (i+1) < (*m-1)) ? (i+1) : (*m-1);
  25       i3 = 1;
  26       F77_FUNC(slarfg,SLARFG)(&i1,&(a[i*(*lda)+i]),&(a[i*(*lda)+i2]),&i3,&(tauq[i]));
  27       d[i] = a[i*(*lda)+i];
  28       a[i*(*lda)+i] = 1.0;
  29       i2 = *n - i - 1;
  30       F77_FUNC(slarf,SLARF)("L",&i1,&i2,&(a[i*(*lda)+i]),&i3,&(tauq[i]),&(a[(i+1)*(*lda)+i]),lda,work);
  31       a[i*(*lda)+i] = d[i];
  32
  33       if(i<(*n-1)) {
  34
  35         i1 = *n - i -1;
  36         i2 = ( (i+2) < (*n-1)) ? (i+2) : (*n-1);
  37         F77_FUNC(slarfg,SLARFG)(&i1,&(a[(i+1)*(*lda)+i]),&(a[i2*(*lda)+i]),lda,&(taup[i]));
  38
  39         e[i] = a[(i+1)*(*lda)+i];
  40         a[(i+1)*(*lda)+i] = 1.0;
  41
  42         i1 = *m - i - 1;
  43         i2 = *n - i - 1;
  44         F77_FUNC(slarf,SLARF)("R",&i1,&i2,&(a[(i+1)*(*lda)+i]),lda,&(taup[i]),&(a[(i+1)*(*lda)+i+1]),lda,work);
  45         a[(i+1)*(*lda)+i] = e[i];
  46       } else
  47         taup[i] = 0.0;
  48     }
  49   } else {
  50     /* reduce to lower bidiag. form */
  51     for(i=0;i<*m;i++) {
  52       i1 = *n - i;
  53       i2 = ( (i+1) < (*n-1)) ? (i+1) : (*n-1);
  54       i3 = 1;
  55       F77_FUNC(slarfg,SLARFG)(&i1,&(a[i*(*lda)+i]),&(a[i2*(*lda)+i]),lda,&(taup[i]));
  56       d[i] = a[i*(*lda)+i];
  57       a[i*(*lda)+i] = 1.0;
  58
  59       i2 = *m - i - 1;
  60       i3 = ( (i+1) < (*m-1)) ? (i+1) : (*m-1);
  61       F77_FUNC(slarf,SLARF)("R",&i2,&i1,&(a[i*(*lda)+i]),lda,&(taup[i]),&(a[(i)*(*lda)+i3]),lda,work);
  62       a[i*(*lda)+i] = d[i];
  63
  64       if(i<(*m-1)) {
  65
  66         i1 = *m - i - 1;
  67         i2 = ( (i+2) < (*m-1)) ? (i+2) : (*m-1);
  68         i3 = 1;
  69         F77_FUNC(slarfg,SLARFG)(&i1,&(a[(i)*(*lda)+i+1]),&(a[i*(*lda)+i2]),&i3,&(tauq[i]));
  70
  71         e[i] = a[(i)*(*lda)+i+1];
  72         a[(i)*(*lda)+i+1] = 1.0;
  73
  74         i1 = *m - i - 1;
  75         i2 = *n - i - 1;
  76         i3 = 1;
  77         F77_FUNC(slarf,SLARF)("L",&i1,&i2,&(a[(i)*(*lda)+i+1]),&i3,&(tauq[i]),&(a[(i+1)*(*lda)+i+1]),lda,work);
  78         a[(i)*(*lda)+i+1] = e[i];
  79       } else
  80         tauq[i] = 0.0;
  81     }
  82   }
  83   return;
  84 }