src/gmxlib/gmx_lapack/sormqr.c

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