src/gmxlib/gmx_lapack/slarft.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 <math.h>
  36 #include <types/simple.h>
  37
  38 #include "gmx_blas.h"
  39 #include "gmx_lapack.h"
  40
  41 void
  42 F77_FUNC(slarft,SLARFT)(const char *direct,
  43         const char *storev,
  44         int *n,
  45         int *k,
  46         float *v,
  47         int *ldv,
  48         float *tau,
  49         float *t,
  50         int *ldt)
  51 {
  52     /* System generated locals */
  53     int t_dim1, t_offset, v_dim1, v_offset, i__1, i__2, i__3;
  54     float d__1;
  55
  56     /* Local variables */
  57     int i__, j;
  58     float vii;
  59     int c__1 = 1;
  60     float zero = 0.0;
  61
  62     v_dim1 = *ldv;
  63     v_offset = 1 + v_dim1;
  64     v -= v_offset;
  65     --tau;
  66     t_dim1 = *ldt;
  67     t_offset = 1 + t_dim1;
  68     t -= t_offset;
  69
  70     if (*n == 0) {
  71         return;
  72     }
  73
  74     if (*direct=='F' || *direct=='f') {
  75         i__1 = *k;
  76         for (i__ = 1; i__ <= i__1; ++i__) {
  77             if (fabs(tau[i__])<GMX_FLOAT_MIN) {
  78
  79                 i__2 = i__;
  80                 for (j = 1; j <= i__2; ++j) {
  81                     t[j + i__ * t_dim1] = 0.;
  82                 }
  83             } else {
  84
  85                 vii = v[i__ + i__ * v_dim1];
  86                 v[i__ + i__ * v_dim1] = 1.;
  87                 if (*storev=='C' || *storev=='c') {
  88
  89                     i__2 = *n - i__ + 1;
  90                     i__3 = i__ - 1;
  91                     d__1 = -tau[i__];
  92                     F77_FUNC(sgemv,SGEMV)("Transpose", &i__2, &i__3, &d__1, &v[i__ + v_dim1],
  93                              ldv, &v[i__ + i__ * v_dim1], &c__1, &zero, &t[
  94                             i__ * t_dim1 + 1], &c__1);
  95                 } else {
  96
  97                     i__2 = i__ - 1;
  98                     i__3 = *n - i__ + 1;
  99                     d__1 = -tau[i__];
 100                     F77_FUNC(sgemv,SGEMV)("No transpose", &i__2, &i__3, &d__1, &v[i__ *
 101                             v_dim1 + 1], ldv, &v[i__ + i__ * v_dim1], ldv, &
 102                             zero, &t[i__ * t_dim1 + 1], &c__1);
 103                 }
 104                 v[i__ + i__ * v_dim1] = vii;
 105
 106
 107                 i__2 = i__ - 1;
 108                 F77_FUNC(strmv,STRMV)("Upper", "No transpose", "Non-unit", &i__2, &t[
 109                         t_offset], ldt, &t[i__ * t_dim1 + 1], &c__1);
 110                 t[i__ + i__ * t_dim1] = tau[i__];
 111             }
 112         }
 113     } else {
 114         for (i__ = *k; i__ >= 1; --i__) {
 115             if (fabs(tau[i__])<GMX_FLOAT_MIN) {
 116
 117                 i__1 = *k;
 118                 for (j = i__; j <= i__1; ++j) {
 119                     t[j + i__ * t_dim1] = 0.;
 120                 }
 121             } else {
 122
 123                 if (i__ < *k) {
 124                     if (*storev=='C' || *storev=='c') {
 125                         vii = v[*n - *k + i__ + i__ * v_dim1];
 126                         v[*n - *k + i__ + i__ * v_dim1] = 1.;
 127
 128                         i__1 = *n - *k + i__;
 129                         i__2 = *k - i__;
 130                         d__1 = -tau[i__];
 131                         F77_FUNC(sgemv,SGEMV)("Transpose", &i__1, &i__2, &d__1, &v[(i__ + 1)
 132                                 * v_dim1 + 1], ldv, &v[i__ * v_dim1 + 1], &
 133                                 c__1, &zero, &t[i__ + 1 + i__ * t_dim1], &
 134                                 c__1);
 135                         v[*n - *k + i__ + i__ * v_dim1] = vii;
 136                     } else {
 137                         vii = v[i__ + (*n - *k + i__) * v_dim1];
 138                         v[i__ + (*n - *k + i__) * v_dim1] = 1.;
 139
 140                         i__1 = *k - i__;
 141                         i__2 = *n - *k + i__;
 142                         d__1 = -tau[i__];
 143                         F77_FUNC(sgemv,SGEMV)("No transpose", &i__1, &i__2, &d__1, &v[i__ +
 144                                 1 + v_dim1], ldv, &v[i__ + v_dim1], ldv, &
 145                                 zero, &t[i__ + 1 + i__ * t_dim1], &c__1);
 146                         v[i__ + (*n - *k + i__) * v_dim1] = vii;
 147                     }
 148
 149                     i__1 = *k - i__;
 150                     F77_FUNC(strmv,STRMV)("Lower", "No transpose", "Non-unit", &i__1, &t[i__
 151                             + 1 + (i__ + 1) * t_dim1], ldt, &t[i__ + 1 + i__ *
 152                              t_dim1], &c__1)
 153                             ;
 154                 }
 155                 t[i__ + i__ * t_dim1] = tau[i__];
 156             }
 157         }
 158     }
 159     return;
 160
 161
 162 }