src/gmxlib/gmx_lapack/sgeqrf.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012, by the GROMACS development team, led by
   5  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   6  * others, as listed in the AUTHORS file in the top-level source
   7  * directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 #include "gmx_lapack.h"
  36 #include "lapack_limits.h"
  37
  38 void
  39 F77_FUNC(sgeqrf,SGEQRF)(int *m,
  40         int *n,
  41         float *a,
  42         int *lda,
  43         float *tau,
  44         float *work,
  45         int *lwork,
  46         int *info)
  47 {
  48     int a_dim1, a_offset, i__1, i__2, i__3, i__4;
  49
  50     int i__, k, ib, nb, nx, iws, nbmin, iinfo;
  51     int ldwork, lwkopt;
  52
  53     a_dim1 = *lda;
  54     a_offset = 1 + a_dim1;
  55     a -= a_offset;
  56     --tau;
  57     --work;
  58
  59     *info = 0;
  60     nb = DGEQRF_BLOCKSIZE;
  61     lwkopt = *n * nb;
  62     work[1] = (float) lwkopt;
  63         if (*lwork==-1)
  64         return;
  65
  66
  67     k = (*m < *n) ? *m : *n;
  68     if (k == 0) {
  69         work[1] = 1.;
  70         return;
  71     }
  72
  73     nbmin = 2;
  74     nx = 0;
  75     iws = *n;
  76     if (nb > 1 && nb < k) {
  77
  78       nx = DGEQRF_CROSSOVER;
  79         if (nx < k) {
  80
  81             ldwork = *n;
  82             iws = ldwork * nb;
  83             if (*lwork < iws) {
  84
  85                 nb = *lwork / ldwork;
  86                 nbmin = DGEQRF_MINBLOCKSIZE;
  87             }
  88         }
  89     }
  90
  91     if (nb >= nbmin && nb < k && nx < k) {
  92         i__1 = k - nx;
  93         i__2 = nb;
  94         for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
  95
  96             i__3 = k - i__ + 1;
  97             ib = (i__3 < nb) ? i__3 : nb;
  98
  99             i__3 = *m - i__ + 1;
 100             F77_FUNC(sgeqr2,SGEQR2)(&i__3, &ib, &a[i__ + i__ * a_dim1], lda, &tau[i__], &work[
 101                     1], &iinfo);
 102             if (i__ + ib <= *n) {
 103
 104                 i__3 = *m - i__ + 1;
 105                 F77_FUNC(slarft,SLARFT)("Forward", "Columnwise", &i__3, &ib, &a[i__ + i__ *
 106                         a_dim1], lda, &tau[i__], &work[1], &ldwork);
 107
 108                 i__3 = *m - i__ + 1;
 109                 i__4 = *n - i__ - ib + 1;
 110                 F77_FUNC(slarfb,SLARFB)("Left", "Transpose", "Forward", "Columnwise", &i__3, &
 111                         i__4, &ib, &a[i__ + i__ * a_dim1], lda, &work[1], &
 112                         ldwork, &a[i__ + (i__ + ib) * a_dim1], lda, &work[ib
 113                         + 1], &ldwork);
 114             }
 115         }
 116     } else {
 117         i__ = 1;
 118     }
 119
 120     if (i__ <= k) {
 121         i__2 = *m - i__ + 1;
 122         i__1 = *n - i__ + 1;
 123         F77_FUNC(sgeqr2,SGEQR2)(&i__2, &i__1, &a[i__ + i__ * a_dim1], lda, &tau[i__], &work[1]
 124                 , &iinfo);
 125     }
 126
 127     work[1] = (float) iws;
 128     return;
 129
 130 }
 131