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