src/gmxlib/gmx_lapack/sormtr.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(sormtr,SORMTR)(const char *side,
  41         const char *uplo,
  42         const char *trans,
  43         int *m,
  44         int *n,
  45         float *a,
  46         int *lda,
  47         float *tau,
  48         float *c__,
  49         int *ldc,
  50         float *work,
  51         int *lwork,
  52         int *info)
  53 {
  54     int a_dim1, a_offset, c_dim1, c_offset, i__2;
  55
  56     int i1, i2, nb, mi, ni, nq, nw;
  57     int left;
  58     int iinfo;
  59     int upper;
  60     int lwkopt;
  61     int lquery;
  62
  63
  64     a_dim1 = *lda;
  65     a_offset = 1 + a_dim1;
  66     a -= a_offset;
  67     --tau;
  68     c_dim1 = *ldc;
  69     c_offset = 1 + c_dim1;
  70     c__ -= c_offset;
  71     --work;
  72
  73     *info = 0;
  74     left = (*side=='L' || *side=='l');
  75     upper = (*uplo=='U' || *uplo=='u');
  76     lquery = *lwork == -1;
  77
  78     if (left) {
  79         nq = *m;
  80         nw = *n;
  81     } else {
  82         nq = *n;
  83         nw = *m;
  84     }
  85
  86
  87     nb = DORMQL_BLOCKSIZE;
  88     lwkopt = nw * nb;
  89     work[1] = (float) lwkopt;
  90
  91     if (*info != 0) {
  92         i__2 = -(*info);
  93         return;
  94     } else if (lquery) {
  95         return;
  96     }
  97
  98     if (*m == 0 || *n == 0 || nq == 1) {
  99         work[1] = 1.;
 100         return;
 101     }
 102
 103     if (left) {
 104         mi = *m - 1;
 105         ni = *n;
 106     } else {
 107         mi = *m;
 108         ni = *n - 1;
 109     }
 110
 111     if (upper) {
 112         i__2 = nq - 1;
 113         F77_FUNC(sormql,SORMQL)(side, trans, &mi, &ni, &i__2, &a[(a_dim1 << 1) + 1], lda, &
 114                 tau[1], &c__[c_offset], ldc, &work[1], lwork, &iinfo);
 115     } else {
 116         if (left) {
 117             i1 = 2;
 118             i2 = 1;
 119         } else {
 120             i1 = 1;
 121             i2 = 2;
 122         }
 123         i__2 = nq - 1;
 124         F77_FUNC(sormqr,SORMQR)(side, trans, &mi, &ni, &i__2, &a[a_dim1 + 2], lda, &tau[1], &
 125                 c__[i1 + i2 * c_dim1], ldc, &work[1], lwork, &iinfo);
 126     }
 127     work[1] = (float) lwkopt;
 128     return;
 129
 130 }
 131
 132