src/gmxlib/gmx_lapack/dormql.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012,2013, 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(dormql,DORMQL)(const char *side, const char *trans, int *m, int *n,
  40         int *k, double *a, int *lda, double *tau, double *
  41         c__, int *ldc, double *work, int *lwork, int *info)
  42 {
  43     int a_dim1, a_offset, c_dim1, c_offset, i__1, i__2, i__4, i__5;
  44     int c__65 = 65;
  45
  46     int i__;
  47     double t[4160];
  48     int i1, i2, i3, ib, nb, mi, ni, nq, nw, iws;
  49     int left;
  50     int nbmin, iinfo;
  51     int notran;
  52     int ldwork, lwkopt;
  53     int lquery;
  54
  55
  56     a_dim1 = *lda;
  57     a_offset = 1 + a_dim1;
  58     a -= a_offset;
  59     --tau;
  60     c_dim1 = *ldc;
  61     c_offset = 1 + c_dim1;
  62     c__ -= c_offset;
  63     --work;
  64
  65     *info = 0;
  66     left = (*side=='L' || *side=='l');
  67     notran = (*trans=='N' || *trans=='n');
  68     lquery = *lwork == -1;
  69
  70     if (left) {
  71         nq = *m;
  72         nw = *n;
  73     } else {
  74         nq = *n;
  75         nw = *m;
  76     }
  77
  78     nb = DORMQL_BLOCKSIZE;
  79     lwkopt = nw * nb;
  80     work[1] = (double) lwkopt;
  81
  82     if (*info != 0) {
  83         i__1 = -(*info);
  84         return;
  85     } else if (lquery) {
  86         return;
  87     }
  88
  89     if (*m == 0 || *n == 0 || *k == 0) {
  90         work[1] = 1.;
  91         return;
  92     }
  93
  94     nbmin = 2;
  95     ldwork = nw;
  96     if (nb > 1 && nb < *k) {
  97         iws = nw * nb;
  98         if (*lwork < iws) {
  99             nb = *lwork / ldwork;
 100             nbmin = DORMQL_MINBLOCKSIZE;
 101         }
 102     } else {
 103         iws = nw;
 104     }
 105
 106     if (nb < nbmin || nb >= *k) {
 107
 108         F77_FUNC(dorm2l,DORM2L)(side, trans, m, n, k, &a[a_offset], lda, &tau[1], &c__[
 109                 c_offset], ldc, &work[1], &iinfo);
 110     } else {
 111
 112         if ((left && notran) || (! left && ! notran)) {
 113             i1 = 1;
 114             i2 = *k;
 115             i3 = nb;
 116         } else {
 117             i1 = (*k - 1) / nb * nb + 1;
 118             i2 = 1;
 119             i3 = -nb;
 120         }
 121
 122         if (left) {
 123             ni = *n;
 124         } else {
 125             mi = *m;
 126         }
 127
 128         i__1 = i2;
 129         i__2 = i3;
 130         for (i__ = i1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
 131             i__4 = nb, i__5 = *k - i__ + 1;
 132             ib = (i__4<i__5) ? i__4 : i__5;
 133
 134             i__4 = nq - *k + i__ + ib - 1;
 135             F77_FUNC(dlarft,DLARFT)("Backward", "Columnwise", &i__4, &ib, &a[i__ * a_dim1 + 1]
 136                     , lda, &tau[i__], t, &c__65);
 137             if (left) {
 138
 139                 mi = *m - *k + i__ + ib - 1;
 140             } else {
 141
 142                 ni = *n - *k + i__ + ib - 1;
 143             }
 144
 145             F77_FUNC(dlarfb,DLARFB)(side, trans, "Backward", "Columnwise", &mi, &ni, &ib, &a[
 146                     i__ * a_dim1 + 1], lda, t, &c__65, &c__[c_offset], ldc, &
 147                     work[1], &ldwork);
 148         }
 149     }
 150     work[1] = (double) lwkopt;
 151     return;
 152
 153 }
 154
 155