src/gmxlib/gmx_lapack/dormbr.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(dormbr,DORMBR)(const char *vect,
  40         const char *side,
  41         const char *trans,
  42         int *m,
  43         int *n,
  44         int *k,
  45         double *a,
  46         int *lda,
  47         double *tau,
  48         double *c__,
  49         int *ldc,
  50         double *work,
  51         int *lwork,
  52         int *info)
  53 {
  54     int a_dim1, a_offset, c_dim1, c_offset, i__1;
  55
  56
  57     int i1, i2, nb, mi, ni, nq, nw;
  58     int left;
  59     int iinfo;
  60     int notran;
  61     int applyq;
  62     char transt[1];
  63     int lwkopt;
  64     int lquery;
  65
  66     a_dim1 = *lda;
  67     a_offset = 1 + a_dim1;
  68     a -= a_offset;
  69     --tau;
  70     c_dim1 = *ldc;
  71     c_offset = 1 + c_dim1;
  72     c__ -= c_offset;
  73     --work;
  74     *info = 0;
  75     applyq = (*vect=='Q' || *vect=='q');
  76     left = (*side=='L' || *side=='l');
  77     notran = (*trans=='N' || *trans=='n');
  78     lquery = *lwork == -1;
  79
  80     if (left) {
  81         nq = *m;
  82         nw = *n;
  83     } else {
  84         nq = *n;
  85         nw = *m;
  86     }
  87
  88     nb = DORMQR_BLOCKSIZE;
  89     lwkopt = nw * nb;
  90     work[1] = (double) lwkopt;
  91
  92     if (*info != 0) {
  93         i__1 = -(*info);
  94         return;
  95     } else if (lquery) {
  96         return;
  97     }
  98
  99     work[1] = 1.;
 100     if (*m == 0 || *n == 0) {
 101         return;
 102     }
 103
 104     if (applyq) {
 105
 106         if (nq >= *k) {
 107
 108             F77_FUNC(dormqr,DORMQR)(side, trans, m, n, k, &a[a_offset], lda, &tau[1], &c__[
 109                     c_offset], ldc, &work[1], lwork, &iinfo);
 110         } else if (nq > 1) {
 111
 112             if (left) {
 113                 mi = *m - 1;
 114                 ni = *n;
 115                 i1 = 2;
 116                 i2 = 1;
 117             } else {
 118                 mi = *m;
 119                 ni = *n - 1;
 120                 i1 = 1;
 121                 i2 = 2;
 122             }
 123             i__1 = nq - 1;
 124             F77_FUNC(dormqr,DORMQR)(side, trans, &mi, &ni, &i__1, &a[a_dim1 + 2], lda, &tau[1]
 125                     , &c__[i1 + i2 * c_dim1], ldc, &work[1], lwork, &iinfo);
 126         }
 127     } else {
 128
 129         if (notran) {
 130             *(unsigned char *)transt = 'T';
 131         } else {
 132             *(unsigned char *)transt = 'N';
 133         }
 134         if (nq > *k) {
 135
 136             F77_FUNC(dormlq,DORMLQ)(side, transt, m, n, k, &a[a_offset], lda, &tau[1], &c__[
 137                     c_offset], ldc, &work[1], lwork, &iinfo);
 138         } else if (nq > 1) {
 139
 140             if (left) {
 141                 mi = *m - 1;
 142                 ni = *n;
 143                 i1 = 2;
 144                 i2 = 1;
 145             } else {
 146                 mi = *m;
 147                 ni = *n - 1;
 148                 i1 = 1;
 149                 i2 = 2;
 150             }
 151             i__1 = nq - 1;
 152             F77_FUNC(dormlq,DORMLQ)(side, transt, &mi, &ni, &i__1, &a[(a_dim1 << 1) + 1], lda,
 153                      &tau[1], &c__[i1 + i2 * c_dim1], ldc, &work[1], lwork, &
 154                     iinfo);
 155         }
 156     }
 157     work[1] = (double) lwkopt;
 158     return;
 159
 160
 161 }
 162
 163