src/gmxlib/gmx_lapack/dorgqr.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
  39 void
  40 F77_FUNC(dorgqr,DORGQR)(int *m,
  41         int *n,
  42         int *k,
  43         double *a,
  44         int *lda,
  45         double *tau,
  46         double *work,
  47         int *lwork,
  48         int *info)
  49 {
  50     int a_dim1, a_offset, i__1, i__2, i__3;
  51
  52     int i__, j, l, ib, nb, ki, kk, nx, iws, nbmin, iinfo;
  53     int ldwork, lwkopt;
  54     int lquery;
  55
  56     a_dim1 = *lda;
  57     a_offset = 1 + a_dim1;
  58     a -= a_offset;
  59     --tau;
  60     --work;
  61
  62     ki = 0;
  63     *info = 0;
  64     nb = DORGQR_BLOCKSIZE;
  65     lwkopt = (*n) * nb;
  66     work[1] = (double) lwkopt;
  67     lquery = *lwork == -1;
  68     if (*m < 0) {
  69         *info = -1;
  70     } else if (*n < 0 || *n > *m) {
  71         *info = -2;
  72     } else if (*k < 0 || *k > *n) {
  73         *info = -3;
  74     } else if (*lda < (*m)) {
  75         *info = -5;
  76     } else if (*lwork < (*n) && ! lquery) {
  77         *info = -8;
  78     }
  79     if (*info != 0) {
  80         i__1 = -(*info);
  81         return;
  82     } else if (lquery) {
  83         return;
  84     }
  85
  86     if (*n <= 0) {
  87         work[1] = 1.;
  88         return;
  89     }
  90
  91     nbmin = 2;
  92     nx = 0;
  93     iws = *n;
  94     if (nb > 1 && nb < *k) {
  95
  96         nx = DORGQR_CROSSOVER;
  97         if (nx < *k) {
  98
  99             ldwork = *n;
 100             iws = ldwork * nb;
 101             if (*lwork < iws) {
 102
 103                 nb = *lwork / ldwork;
 104                 nbmin = DORGQR_MINBLOCKSIZE;
 105             }
 106         }
 107     }
 108
 109     if (nb >= nbmin && nb < *k && nx < *k) {
 110
 111         ki = (*k - nx - 1) / nb * nb;
 112         i__1 = *k, i__2 = ki + nb;
 113         kk = (i__1<i__2) ? i__1 : i__2;
 114
 115         i__1 = *n;
 116         for (j = kk + 1; j <= i__1; ++j) {
 117             i__2 = kk;
 118             for (i__ = 1; i__ <= i__2; ++i__) {
 119                 a[i__ + j * a_dim1] = 0.;
 120             }
 121         }
 122     } else {
 123         kk = 0;
 124     }
 125
 126     if (kk < *n) {
 127         i__1 = *m - kk;
 128         i__2 = *n - kk;
 129         i__3 = *k - kk;
 130         F77_FUNC(dorg2r,DORG2R)(&i__1, &i__2, &i__3, &a[kk + 1 + (kk + 1) * a_dim1], lda, &
 131                 tau[kk + 1], &work[1], &iinfo);
 132     }
 133
 134     if (kk > 0) {
 135
 136         i__1 = -nb;
 137         for (i__ = ki + 1; i__1 < 0 ? i__ >= 1 : i__ <= 1; i__ += i__1) {
 138             i__2 = nb, i__3 = *k - i__ + 1;
 139             ib = (i__2<i__3) ? i__2 : i__3;
 140             if (i__ + ib <= *n) {
 141
 142                 i__2 = *m - i__ + 1;
 143                 F77_FUNC(dlarft,DLARFT)("Forward", "Columnwise", &i__2, &ib, &a[i__ + i__ *
 144                         a_dim1], lda, &tau[i__], &work[1], &ldwork);
 145
 146                 i__2 = *m - i__ + 1;
 147                 i__3 = *n - i__ - ib + 1;
 148                 F77_FUNC(dlarfb,DLARFB)("Left", "No transpose", "Forward", "Columnwise", &
 149                         i__2, &i__3, &ib, &a[i__ + i__ * a_dim1], lda, &work[
 150                         1], &ldwork, &a[i__ + (i__ + ib) * a_dim1], lda, &
 151                         work[ib + 1], &ldwork);
 152             }
 153
 154             i__2 = *m - i__ + 1;
 155             F77_FUNC(dorg2r,DORG2R)(&i__2, &ib, &ib, &a[i__ + i__ * a_dim1], lda, &tau[i__], &
 156                     work[1], &iinfo);
 157
 158             i__2 = i__ + ib - 1;
 159             for (j = i__; j <= i__2; ++j) {
 160                 i__3 = i__ - 1;
 161                 for (l = 1; l <= i__3; ++l) {
 162                     a[l + j * a_dim1] = 0.;
 163                 }
 164             }
 165         }
 166     }
 167
 168     work[1] = (double) iws;
 169     return;
 170
 171 }