src/gmxlib/gmx_lapack/slagts.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 void
  43 F77_FUNC(slagts,SLAGTS)(int *job,
  44         int *n,
  45         float *a,
  46         float *b,
  47         float *c__,
  48         float *d__,
  49         int *in,
  50         float *y,
  51         float *tol,
  52         int *info)
  53 {
  54     int i__1;
  55     float d__1, d__2, d__4, d__5;
  56
  57     int k;
  58     float ak, eps, temp, pert, absak, sfmin;
  59     float bignum,minval;
  60     --y;
  61     --in;
  62     --d__;
  63     --c__;
  64     --b;
  65     --a;
  66
  67     *info = 0;
  68     if (fabs(*job) > 2 || *job == 0) {
  69         *info = -1;
  70     } else if (*n < 0) {
  71         *info = -2;
  72     }
  73     if (*info != 0) {
  74         i__1 = -(*info);
  75         return;
  76     }
  77
  78     if (*n == 0) {
  79         return;
  80     }
  81     eps = GMX_FLOAT_EPS;
  82     minval = GMX_FLOAT_MIN;
  83     sfmin = minval / eps;
  84
  85     bignum = 1. / sfmin;
  86
  87     if (*job < 0) {
  88         if (*tol <= 0.) {
  89             *tol = fabs(a[1]);
  90             if (*n > 1) {
  91                 d__1 = *tol;
  92                 d__2 = fabs(a[2]);
  93                 d__1 = (d__1>d__2) ? d__1 : d__2;
  94                 d__2 = fabs(b[1]);
  95                 *tol = (d__1>d__2) ? d__1 : d__2;
  96             }
  97             i__1 = *n;
  98             for (k = 3; k <= i__1; ++k) {
  99               d__4 = *tol;
 100               d__5 = fabs(a[k]);
 101               d__4 = (d__4>d__5) ? d__4 : d__5;
 102               d__5 = fabs(b[k - 1]);
 103               d__4 = (d__4>d__5) ? d__4 : d__5;
 104               d__5 = fabs(d__[k - 2]);
 105               *tol = (d__4>d__5) ? d__4 : d__5;
 106             }
 107             *tol *= eps;
 108             if (fabs(*tol)<GMX_FLOAT_MIN) {
 109                 *tol = eps;
 110             }
 111         }
 112     }
 113
 114     if (fabs(fabs(*job)-1.0)<GMX_FLOAT_MIN) {
 115         i__1 = *n;
 116         for (k = 2; k <= i__1; ++k) {
 117             if (in[k - 1] == 0) {
 118                 y[k] -= c__[k - 1] * y[k - 1];
 119             } else {
 120                 temp = y[k - 1];
 121                 y[k - 1] = y[k];
 122                 y[k] = temp - c__[k - 1] * y[k];
 123             }
 124         }
 125         if (*job == 1) {
 126             for (k = *n; k >= 1; --k) {
 127                 if (k <= *n - 2) {
 128                     temp = y[k] - b[k] * y[k + 1] - d__[k] * y[k + 2];
 129                 } else if (k == *n - 1) {
 130                     temp = y[k] - b[k] * y[k + 1];
 131                 } else {
 132                     temp = y[k];
 133                 }
 134                 ak = a[k];
 135                 absak = fabs(ak);
 136                 if (absak < 1.) {
 137                     if (absak < sfmin) {
 138                         if (fabs(absak)<GMX_FLOAT_MIN || fabs(temp) * sfmin > absak) {
 139                             *info = k;
 140                             return;
 141                         } else {
 142                             temp *= bignum;
 143                             ak *= bignum;
 144                         }
 145                     } else if (fabs(temp) > absak * bignum) {
 146                         *info = k;
 147                         return;
 148                     }
 149                 }
 150                 y[k] = temp / ak;
 151             }
 152         } else {
 153             for (k = *n; k >= 1; --k) {
 154                 if (k <= *n - 2) {
 155                     temp = y[k] - b[k] * y[k + 1] - d__[k] * y[k + 2];
 156                 } else if (k == *n - 1) {
 157                     temp = y[k] - b[k] * y[k + 1];
 158                 } else {
 159                     temp = y[k];
 160                 }
 161                 ak = a[k];
 162
 163                 pert = *tol;
 164                 if(ak<0)
 165                   pert *= -1.0;
 166 L40:
 167                 absak = fabs(ak);
 168                 if (absak < 1.) {
 169                     if (absak < sfmin) {
 170                         if (fabs(absak)<GMX_FLOAT_MIN || fabs(temp) * sfmin > absak) {
 171                             ak += pert;
 172                             pert *= 2;
 173                             goto L40;
 174                         } else {
 175                             temp *= bignum;
 176                             ak *= bignum;
 177                         }
 178                     } else if (fabs(temp) > absak * bignum) {
 179                         ak += pert;
 180                         pert *= 2;
 181                         goto L40;
 182                     }
 183                 }
 184                 y[k] = temp / ak;
 185             }
 186         }
 187     } else {
 188
 189         if (*job == 2) {
 190             i__1 = *n;
 191             for (k = 1; k <= i__1; ++k) {
 192                 if (k >= 3) {
 193                     temp = y[k] - b[k - 1] * y[k - 1] - d__[k - 2] * y[k - 2];
 194                 } else if (k == 2) {
 195                     temp = y[k] - b[k - 1] * y[k - 1];
 196                 } else {
 197                     temp = y[k];
 198                 }
 199                 ak = a[k];
 200                 absak = fabs(ak);
 201                 if (absak < 1.) {
 202                     if (absak < sfmin) {
 203                         if (fabs(absak)<GMX_FLOAT_MIN || fabs(temp) * sfmin > absak) {
 204                             *info = k;
 205                             return;
 206                         } else {
 207                             temp *= bignum;
 208                             ak *= bignum;
 209                         }
 210                     } else if (fabs(temp) > absak * bignum) {
 211                         *info = k;
 212                         return;
 213                     }
 214                 }
 215                 y[k] = temp / ak;
 216             }
 217         } else {
 218             i__1 = *n;
 219             for (k = 1; k <= i__1; ++k) {
 220                 if (k >= 3) {
 221                     temp = y[k] - b[k - 1] * y[k - 1] - d__[k - 2] * y[k - 2];
 222                 } else if (k == 2) {
 223                     temp = y[k] - b[k - 1] * y[k - 1];
 224                 } else {
 225                     temp = y[k];
 226                 }
 227                 ak = a[k];
 228
 229                 pert = *tol;
 230                 if(ak<0)
 231                   pert *= -1.0;
 232
 233 L70:
 234                 absak = fabs(ak);
 235                 if (absak < 1.) {
 236                     if (absak < sfmin) {
 237                         if (fabs(absak)<GMX_FLOAT_MIN || fabs(temp) * sfmin > absak) {
 238                             ak += pert;
 239                             pert *= 2;
 240                             goto L70;
 241                         } else {
 242                             temp *= bignum;
 243                             ak *= bignum;
 244                         }
 245                     } else if (fabs(temp) > absak * bignum) {
 246                         ak += pert;
 247                         pert *= 2;
 248                         goto L70;
 249                     }
 250                 }
 251                 y[k] = temp / ak;
 252             }
 253         }
 254
 255         for (k = *n; k >= 2; --k) {
 256             if (in[k - 1] == 0) {
 257                 y[k - 1] -= c__[k - 1] * y[k];
 258             } else {
 259                 temp = y[k - 1];
 260                 y[k - 1] = y[k];
 261                 y[k] = temp - c__[k - 1] * y[k];
 262             }
 263         }
 264     }
 265
 266     return;
 267 }
 268
 269