src/gmxlib/gmx_lapack/sgelqf.c

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