src/gmxlib/gmx_lapack/sgeqrf.c

   1 #include "gmx_lapack.h"
   2 #include "lapack_limits.h"
   3
   4 void
   5 F77_FUNC(sgeqrf,SGEQRF)(int *m,
   6         int *n,
   7         float *a,
   8         int *lda,
   9         float *tau,
  10         float *work,
  11         int *lwork,
  12         int *info)
  13 {
  14     int a_dim1, a_offset, i__1, i__2, i__3, i__4;
  15
  16     int i__, k, ib, nb, nx, iws, nbmin, iinfo;
  17     int ldwork, lwkopt;
  18
  19     a_dim1 = *lda;
  20     a_offset = 1 + a_dim1;
  21     a -= a_offset;
  22     --tau;
  23     --work;
  24
  25     *info = 0;
  26     nb = DGEQRF_BLOCKSIZE;
  27     lwkopt = *n * nb;
  28     work[1] = (float) lwkopt;
  29         if (*lwork==-1)
  30         return;
  31
  32
  33     k = (*m < *n) ? *m : *n;
  34     if (k == 0) {
  35         work[1] = 1.;
  36         return;
  37     }
  38
  39     nbmin = 2;
  40     nx = 0;
  41     iws = *n;
  42     if (nb > 1 && nb < k) {
  43
  44       nx = DGEQRF_CROSSOVER;
  45         if (nx < k) {
  46
  47             ldwork = *n;
  48             iws = ldwork * nb;
  49             if (*lwork < iws) {
  50
  51                 nb = *lwork / ldwork;
  52                 nbmin = DGEQRF_MINBLOCKSIZE;
  53             }
  54         }
  55     }
  56
  57     if (nb >= nbmin && nb < k && nx < k) {
  58         i__1 = k - nx;
  59         i__2 = nb;
  60         for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
  61
  62             i__3 = k - i__ + 1;
  63             ib = (i__3 < nb) ? i__3 : nb;
  64
  65             i__3 = *m - i__ + 1;
  66             F77_FUNC(sgeqr2,SGEQR2)(&i__3, &ib, &a[i__ + i__ * a_dim1], lda, &tau[i__], &work[
  67                     1], &iinfo);
  68             if (i__ + ib <= *n) {
  69
  70                 i__3 = *m - i__ + 1;
  71                 F77_FUNC(slarft,SLARFT)("Forward", "Columnwise", &i__3, &ib, &a[i__ + i__ *
  72                         a_dim1], lda, &tau[i__], &work[1], &ldwork);
  73
  74                 i__3 = *m - i__ + 1;
  75                 i__4 = *n - i__ - ib + 1;
  76                 F77_FUNC(slarfb,SLARFB)("Left", "Transpose", "Forward", "Columnwise", &i__3, &
  77                         i__4, &ib, &a[i__ + i__ * a_dim1], lda, &work[1], &
  78                         ldwork, &a[i__ + (i__ + ib) * a_dim1], lda, &work[ib
  79                         + 1], &ldwork);
  80             }
  81         }
  82     } else {
  83         i__ = 1;
  84     }
  85
  86     if (i__ <= k) {
  87         i__2 = *m - i__ + 1;
  88         i__1 = *n - i__ + 1;
  89         F77_FUNC(sgeqr2,SGEQR2)(&i__2, &i__1, &a[i__ + i__ * a_dim1], lda, &tau[i__], &work[1]
  90                 , &iinfo);
  91     }
  92
  93     work[1] = (float) iws;
  94     return;
  95
  96 }
  97