src/gmxlib/gmx_lapack/dgebrd.c

   1 #include "gmx_lapack.h"
   2 #include "gmx_blas.h"
   3 #include "lapack_limits.h"
   4
   5
   6 void
   7 F77_FUNC(dgebrd,DGEBRD)(int *m,
   8         int *n,
   9         double *a,
  10         int *lda,
  11         double *d__,
  12         double *e,
  13         double *tauq,
  14         double *taup,
  15         double *work,
  16         int *lwork,
  17         int *info)
  18 {
  19     /* System generated locals */
  20     int a_dim1, a_offset, i_1, i_2, i_3, i_4;
  21
  22     /* Local variables */
  23     int i_, j, nx,nb;
  24     double ws;
  25     int nbmin, iinfo, minmn;
  26     int ldwrkx, ldwrky;
  27     double one = 1.0;
  28     double minusone = -1.0;
  29
  30     a_dim1 = *lda;
  31     a_offset = 1 + a_dim1;
  32     a -= a_offset;
  33     --d__;
  34     --e;
  35     --tauq;
  36     --taup;
  37     --work;
  38
  39     nb = DGEBRD_BLOCKSIZE;
  40     *info = 0;
  41     if (*lwork==-1) {
  42       work[1] = (double) ( (*m + *n) * nb);
  43       return;
  44     }
  45     minmn = (*m < *n) ? *m : *n;
  46     if (minmn == 0) {
  47       work[1] = 1.;
  48       return;
  49     }
  50
  51     ws = (*m > *n) ? *m : *n;
  52     ldwrkx = *m;
  53     ldwrky = *n;
  54
  55     if (nb > 1 && nb < minmn) {
  56         nx = DGEBRD_CROSSOVER;
  57         if (nx < minmn) {
  58             ws = (double) ((*m + *n) * nb);
  59             if ((double) (*lwork) < ws) {
  60               nbmin = DGEBRD_MINBLOCKSIZE;
  61                 if (*lwork >= (*m + *n) * nbmin) {
  62                     nb = *lwork / (*m + *n);
  63                 } else {
  64                     nb = 1;
  65                     nx = minmn;
  66                 }
  67             }
  68         }
  69     } else {
  70         nx = minmn;
  71     }
  72
  73     i_1 = minmn - nx;
  74     i_2 = nb;
  75     for (i_ = 1; i_2 < 0 ? i_ >= i_1 : i_ <= i_1; i_ += i_2) {
  76
  77         i_3 = *m - i_ + 1;
  78         i_4 = *n - i_ + 1;
  79         F77_FUNC(dlabrd,DLABRD)(&i_3, &i_4, &nb, &a[i_ + i_ * a_dim1], lda, &d__[i_],
  80                 &e[i_], &tauq[i_], &taup[i_], &work[1], &ldwrkx,
  81                 &work[ldwrkx * nb + 1], &ldwrky);
  82
  83         i_3 = *m - i_ - nb + 1;
  84         i_4 = *n - i_ - nb + 1;
  85         F77_FUNC(dgemm,DGEMM)("N", "T", &i_3, &i_4, &nb, &minusone,
  86                &a[i_ + nb + i_ * a_dim1], lda, &work[ldwrkx * nb + nb + 1],
  87                &ldwrky, &one, &a[i_ + nb + (i_ + nb) * a_dim1], lda);
  88         i_3 = *m - i_ - nb + 1;
  89         i_4 = *n - i_ - nb + 1;
  90         F77_FUNC(dgemm,DGEMM)("N", "N", &i_3, &i_4, &nb, &minusone, &work[nb + 1], &ldwrkx,
  91                &a[i_ + (i_ + nb) * a_dim1], lda, &one,
  92                &a[i_ + nb + (i_ + nb) * a_dim1], lda);
  93
  94         if (*m >= *n) {
  95             i_3 = i_ + nb - 1;
  96             for (j = i_; j <= i_3; ++j) {
  97                 a[j + j * a_dim1] = d__[j];
  98                 a[j + (j + 1) * a_dim1] = e[j];
  99             }
 100         } else {
 101             i_3 = i_ + nb - 1;
 102             for (j = i_; j <= i_3; ++j) {
 103                 a[j + j * a_dim1] = d__[j];
 104                 a[j + 1 + j * a_dim1] = e[j];
 105             }
 106         }
 107     }
 108
 109     i_2 = *m - i_ + 1;
 110     i_1 = *n - i_ + 1;
 111     F77_FUNC(dgebd2,DGEBD2)(&i_2, &i_1, &a[i_ + i_ * a_dim1], lda, &d__[i_], &e[i_], &
 112             tauq[i_], &taup[i_], &work[1], &iinfo);
 113     work[1] = ws;
 114     return;
 115
 116 }