src/gmxlib/gmx_lapack/slagtf.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 <math.h>
  36 #include <types/simple.h>
  37
  38 #include "gmx_lapack.h"
  39 #include "lapack_limits.h"
  40
  41
  42
  43 void
  44 F77_FUNC(slagtf,SLAGTF)(int *n,
  45         float *a,
  46         float *lambda,
  47         float *b,
  48         float *c__,
  49         float *tol,
  50         float *d__,
  51         int *in,
  52         int *info)
  53 {
  54     int i__1;
  55
  56     int k;
  57     float tl, eps, piv1, piv2, temp, mult, scale1, scale2;
  58
  59     --in;
  60     --d__;
  61     --c__;
  62     --b;
  63     --a;
  64
  65     *info = 0;
  66     if (*n < 0) {
  67         *info = -1;
  68         i__1 = -(*info);
  69         return;
  70     }
  71
  72     if (*n == 0)
  73         return;
  74
  75     a[1] -= *lambda;
  76     in[*n] = 0;
  77     if (*n == 1) {
  78         if (fabs(a[1])<GMX_FLOAT_MIN) {
  79             in[1] = 1;
  80         }
  81         return;
  82     }
  83
  84     eps = GMX_FLOAT_EPS;
  85
  86     tl = (*tol>eps) ? *tol : eps;
  87     scale1 = fabs(a[1]) + fabs(b[1]);
  88     i__1 = *n - 1;
  89     for (k = 1; k <= i__1; ++k) {
  90         a[k + 1] -= *lambda;
  91         scale2 = fabs(c__[k]) + fabs(a[k + 1]);
  92         if (k < *n - 1) {
  93             scale2 += fabs(b[k + 1]);
  94         }
  95         if (fabs(a[k])<GMX_FLOAT_MIN) {
  96             piv1 = 0.;
  97         } else {
  98             piv1 = fabs(a[k]) / scale1;
  99         }
 100         if (fabs(c__[k])<GMX_FLOAT_MIN) {
 101             in[k] = 0;
 102             piv2 = 0.;
 103             scale1 = scale2;
 104             if (k < *n - 1) {
 105                 d__[k] = 0.;
 106             }
 107         } else {
 108             piv2 = fabs(c__[k]) / scale2;
 109             if (piv2 <= piv1) {
 110                 in[k] = 0;
 111                 scale1 = scale2;
 112                 c__[k] /= a[k];
 113                 a[k + 1] -= c__[k] * b[k];
 114                 if (k < *n - 1) {
 115                     d__[k] = 0.;
 116                 }
 117             } else {
 118                 in[k] = 1;
 119                 mult = a[k] / c__[k];
 120                 a[k] = c__[k];
 121                 temp = a[k + 1];
 122                 a[k + 1] = b[k] - mult * temp;
 123                 if (k < *n - 1) {
 124                     d__[k] = b[k + 1];
 125                     b[k + 1] = -mult * d__[k];
 126                 }
 127                 b[k] = temp;
 128                 c__[k] = mult;
 129             }
 130         }
 131         if (((piv1>piv2) ? piv1 : piv2) <= tl && in[*n] == 0) {
 132             in[*n] = k;
 133         }
 134     }
 135     if (fabs(a[*n]) <= scale1 * tl && in[*n] == 0) {
 136         in[*n] = *n;
 137     }
 138
 139     return;
 140
 141 }
 142
 143