src/gromacs/linearalgebra/gmx_lapack/dorgqr.c

   1 #include "../gmx_lapack.h"
   2 #include "lapack_limits.h"
   3
   4
   5 void
   6 F77_FUNC(dorgqr,DORGQR)(int *m,
   7         int *n,
   8         int *k,
   9         double *a,
  10         int *lda,
  11         double *tau,
  12         double *work,
  13         int *lwork,
  14         int *info)
  15 {
  16     int a_dim1, a_offset, i__1, i__2, i__3;
  17
  18     int i__, j, l, ib, nb, ki, kk, nx, iws, nbmin, iinfo;
  19     int ldwork, lwkopt;
  20     int lquery;
  21
  22     a_dim1 = *lda;
  23     a_offset = 1 + a_dim1;
  24     a -= a_offset;
  25     --tau;
  26     --work;
  27
  28     ki = 0;
  29     *info = 0;
  30     nb = DORGQR_BLOCKSIZE;
  31     lwkopt = (*n) * nb;
  32     work[1] = (double) lwkopt;
  33     lquery = *lwork == -1;
  34     if (*m < 0) {
  35         *info = -1;
  36     } else if (*n < 0 || *n > *m) {
  37         *info = -2;
  38     } else if (*k < 0 || *k > *n) {
  39         *info = -3;
  40     } else if (*lda < (*m)) {
  41         *info = -5;
  42     } else if (*lwork < (*n) && ! lquery) {
  43         *info = -8;
  44     }
  45     if (*info != 0) {
  46         i__1 = -(*info);
  47         return;
  48     } else if (lquery) {
  49         return;
  50     }
  51
  52     if (*n <= 0) {
  53         work[1] = 1.;
  54         return;
  55     }
  56
  57     nbmin = 2;
  58     nx = 0;
  59     iws = *n;
  60     if (nb > 1 && nb < *k) {
  61
  62         nx = DORGQR_CROSSOVER;
  63         if (nx < *k) {
  64
  65             ldwork = *n;
  66             iws = ldwork * nb;
  67             if (*lwork < iws) {
  68
  69                 nb = *lwork / ldwork;
  70                 nbmin = DORGQR_MINBLOCKSIZE;
  71             }
  72         }
  73     }
  74
  75     if (nb >= nbmin && nb < *k && nx < *k) {
  76
  77         ki = (*k - nx - 1) / nb * nb;
  78         i__1 = *k, i__2 = ki + nb;
  79         kk = (i__1<i__2) ? i__1 : i__2;
  80
  81         i__1 = *n;
  82         for (j = kk + 1; j <= i__1; ++j) {
  83             i__2 = kk;
  84             for (i__ = 1; i__ <= i__2; ++i__) {
  85                 a[i__ + j * a_dim1] = 0.;
  86             }
  87         }
  88     } else {
  89         kk = 0;
  90     }
  91
  92     if (kk < *n) {
  93         i__1 = *m - kk;
  94         i__2 = *n - kk;
  95         i__3 = *k - kk;
  96         F77_FUNC(dorg2r,DORG2R)(&i__1, &i__2, &i__3, &a[kk + 1 + (kk + 1) * a_dim1], lda, &
  97                 tau[kk + 1], &work[1], &iinfo);
  98     }
  99
 100     if (kk > 0) {
 101
 102         i__1 = -nb;
 103         for (i__ = ki + 1; i__1 < 0 ? i__ >= 1 : i__ <= 1; i__ += i__1) {
 104             i__2 = nb, i__3 = *k - i__ + 1;
 105             ib = (i__2<i__3) ? i__2 : i__3;
 106             if (i__ + ib <= *n) {
 107
 108                 i__2 = *m - i__ + 1;
 109                 F77_FUNC(dlarft,DLARFT)("Forward", "Columnwise", &i__2, &ib, &a[i__ + i__ *
 110                         a_dim1], lda, &tau[i__], &work[1], &ldwork);
 111
 112                 i__2 = *m - i__ + 1;
 113                 i__3 = *n - i__ - ib + 1;
 114                 F77_FUNC(dlarfb,DLARFB)("Left", "No transpose", "Forward", "Columnwise", &
 115                         i__2, &i__3, &ib, &a[i__ + i__ * a_dim1], lda, &work[
 116                         1], &ldwork, &a[i__ + (i__ + ib) * a_dim1], lda, &
 117                         work[ib + 1], &ldwork);
 118             }
 119
 120             i__2 = *m - i__ + 1;
 121             F77_FUNC(dorg2r,DORG2R)(&i__2, &ib, &ib, &a[i__ + i__ * a_dim1], lda, &tau[i__], &
 122                     work[1], &iinfo);
 123
 124             i__2 = i__ + ib - 1;
 125             for (j = i__; j <= i__2; ++j) {
 126                 i__3 = i__ - 1;
 127                 for (l = 1; l <= i__3; ++l) {
 128                     a[l + j * a_dim1] = 0.;
 129                 }
 130             }
 131         }
 132     }
 133
 134     work[1] = (double) iws;
 135     return;
 136
 137 }