src/gmxlib/gmx_lapack/dormlq.c

   1 #include "gmx_lapack.h"
   2 #include "lapack_limits.h"
   3
   4
   5 void
   6 F77_FUNC(dormlq,DORMLQ)(const char *side,
   7         const char *trans,
   8         int *m,
   9         int *n,
  10         int *k,
  11         double *a,
  12         int *lda,
  13         double *tau,
  14         double *c__,
  15         int *ldc,
  16         double *work,
  17         int *lwork,
  18         int *info)
  19 {
  20     int a_dim1, a_offset, c_dim1, c_offset, i__1, i__2, i__4,
  21             i__5;
  22
  23
  24     int i__;
  25     double t[4160]      /* was [65][64] */;
  26     int i1, i2, i3, ib, ic, jc, nb, mi, ni, nq, nw, iws;
  27     int left;
  28     int nbmin, iinfo;
  29     int notran;
  30     int ldwork;
  31     char transt[1];
  32     int lwkopt;
  33     int lquery;
  34     int ldt = 65;
  35
  36     a_dim1 = *lda;
  37     a_offset = 1 + a_dim1;
  38     a -= a_offset;
  39     --tau;
  40     c_dim1 = *ldc;
  41     c_offset = 1 + c_dim1;
  42     c__ -= c_offset;
  43     --work;
  44
  45     ic = jc = 0;
  46
  47     *info = 0;
  48     left = (*side=='L' || *side=='l');
  49     notran = (*trans=='N' || *trans=='n');
  50     lquery = *lwork == -1;
  51
  52     if (left) {
  53         nq = *m;
  54         nw = *n;
  55     } else {
  56         nq = *n;
  57         nw = *m;
  58     }
  59
  60     nb = DORMLQ_BLOCKSIZE;
  61     lwkopt = nw * nb;
  62     work[1] = (double) lwkopt;
  63
  64     if (*info != 0) {
  65         i__1 = -(*info);
  66         return;
  67     } else if (lquery) {
  68         return;
  69     }
  70
  71     if (*m == 0 || *n == 0 || *k == 0) {
  72         work[1] = 1.;
  73         return;
  74     }
  75
  76     nbmin = 2;
  77     ldwork = nw;
  78     if (nb > 1 && nb < *k) {
  79         iws = nw * nb;
  80         if (*lwork < iws) {
  81             nb = *lwork / ldwork;
  82             nbmin = DORMLQ_MINBLOCKSIZE;
  83         }
  84     } else {
  85         iws = nw;
  86     }
  87
  88     if (nb < nbmin || nb >= *k) {
  89
  90
  91         F77_FUNC(dorml2,DORML2)(side, trans, m, n, k, &a[a_offset], lda, &tau[1], &c__[
  92                 c_offset], ldc, &work[1], &iinfo);
  93     } else {
  94
  95         if ((left && notran) || (!left && !notran)) {
  96             i1 = 1;
  97             i2 = *k;
  98             i3 = nb;
  99         } else {
 100             i1 = (*k - 1) / nb * nb + 1;
 101             i2 = 1;
 102             i3 = -nb;
 103         }
 104
 105         if (left) {
 106             ni = *n;
 107             jc = 1;
 108         } else {
 109             mi = *m;
 110             ic = 1;
 111         }
 112
 113         if (notran) {
 114             *(unsigned char *)transt = 'T';
 115         } else {
 116             *(unsigned char *)transt = 'N';
 117         }
 118
 119         i__1 = i2;
 120         i__2 = i3;
 121         for (i__ = i1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
 122             i__4 = nb, i__5 = *k - i__ + 1;
 123             ib = (i__4<i__5) ? i__4 : i__5;
 124
 125             i__4 = nq - i__ + 1;
 126             F77_FUNC(dlarft,DLARFT)("Forward", "Rowwise", &i__4, &ib, &a[i__ + i__ * a_dim1],
 127                     lda, &tau[i__], t, &ldt);
 128             if (left) {
 129
 130                 mi = *m - i__ + 1;
 131                 ic = i__;
 132             } else {
 133
 134                 ni = *n - i__ + 1;
 135                 jc = i__;
 136             }
 137
 138             F77_FUNC(dlarfb,DLARFB)(side, transt, "Forward", "Rowwise", &mi, &ni, &ib, &a[i__
 139                     + i__ * a_dim1], lda, t, &ldt, &c__[ic + jc * c_dim1],
 140                     ldc, &work[1], &ldwork);
 141         }
 142     }
 143     work[1] = (double) lwkopt;
 144     return;
 145
 146 }
 147
 148