src/gmxlib/gmx_lapack/dgebd2.c

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