src/gmxlib/gmx_lapack/slasq5.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 "gmx_lapack.h"
  37
  38 void
  39 F77_FUNC(slasq5,SLASQ5)(int *i0,
  40         int *n0,
  41         float *z__,
  42         int *pp,
  43         float *tau,
  44         float *dmin__,
  45         float *dmin1,
  46         float *dmin2,
  47         float *dn,
  48         float *dnm1,
  49         float *dnm2,
  50         int *ieee)
  51 {
  52     int i__1;
  53     float d__1, d__2;
  54
  55     static float d__;
  56     static int j4, j4p2;
  57     static float emin, temp;
  58
  59     --z__;
  60
  61     if (*n0 - *i0 - 1 <= 0) {
  62         return;
  63     }
  64
  65     j4 = (*i0 << 2) + *pp - 3;
  66     emin = z__[j4 + 4];
  67     d__ = z__[j4] - *tau;
  68     *dmin__ = d__;
  69     *dmin1 = -z__[j4];
  70
  71     if (*ieee) {
  72
  73         if (*pp == 0) {
  74             i__1 = 4*(*n0 - 3);
  75             for (j4 = *i0 << 2; j4 <= i__1; j4 += 4) {
  76                 z__[j4 - 2] = d__ + z__[j4 - 1];
  77                 temp = z__[j4 + 1] / z__[j4 - 2];
  78                 d__ = d__ * temp - *tau;
  79                 if(d__<*dmin__)
  80                   *dmin__ = d__;
  81                 z__[j4] = z__[j4 - 1] * temp;
  82                 d__1 = z__[j4];
  83                 if(d__1<emin)
  84                   emin = d__1;
  85             }
  86         } else {
  87             i__1 = 4*(*n0 - 3);
  88             for (j4 = *i0 << 2; j4 <= i__1; j4 += 4) {
  89                 z__[j4 - 3] = d__ + z__[j4];
  90                 temp = z__[j4 + 2] / z__[j4 - 3];
  91                 d__ = d__ * temp - *tau;
  92                 if(d__<*dmin__)
  93                   *dmin__ = d__;
  94                 z__[j4 - 1] = z__[j4] * temp;
  95                 d__1 = z__[j4 - 1];
  96                 if(d__1<emin)
  97                   emin = d__1;
  98             }
  99         }
 100
 101         *dnm2 = d__;
 102         *dmin2 = *dmin__;
 103         j4 = 4*(*n0 - 2) - *pp;
 104         j4p2 = j4 + (*pp << 1) - 1;
 105         z__[j4 - 2] = *dnm2 + z__[j4p2];
 106         z__[j4] = z__[j4p2 + 2] * (z__[j4p2] / z__[j4 - 2]);
 107         *dnm1 = z__[j4p2 + 2] * (*dnm2 / z__[j4 - 2]) - *tau;
 108         if(*dnm1<*dmin__)
 109           *dmin__ = *dnm1;
 110
 111         *dmin1 = *dmin__;
 112         j4 += 4;
 113         j4p2 = j4 + (*pp << 1) - 1;
 114         z__[j4 - 2] = *dnm1 + z__[j4p2];
 115         z__[j4] = z__[j4p2 + 2] * (z__[j4p2] / z__[j4 - 2]);
 116         *dn = z__[j4p2 + 2] * (*dnm1 / z__[j4 - 2]) - *tau;
 117         if(*dn<*dmin__)
 118           *dmin__ = *dn;
 119
 120     } else {
 121
 122         if (*pp == 0) {
 123             i__1 = 4*(*n0 - 3);
 124             for (j4 = *i0 << 2; j4 <= i__1; j4 += 4) {
 125                 z__[j4 - 2] = d__ + z__[j4 - 1];
 126                 if (d__ < 0.) {
 127                     return;
 128                 } else {
 129                     z__[j4] = z__[j4 + 1] * (z__[j4 - 1] / z__[j4 - 2]);
 130                     d__ = z__[j4 + 1] * (d__ / z__[j4 - 2]) - *tau;
 131                 }
 132                 if(d__<*dmin__)
 133                   *dmin__ = d__;
 134                 d__1 = emin, d__2 = z__[j4];
 135                 emin = (d__1<d__2) ? d__1 : d__2;
 136             }
 137         } else {
 138             i__1 = 4*(*n0 - 3);
 139             for (j4 = *i0 << 2; j4 <= i__1; j4 += 4) {
 140                 z__[j4 - 3] = d__ + z__[j4];
 141                 if (d__ < 0.) {
 142                     return;
 143                 } else {
 144                     z__[j4 - 1] = z__[j4 + 2] * (z__[j4] / z__[j4 - 3]);
 145                     d__ = z__[j4 + 2] * (d__ / z__[j4 - 3]) - *tau;
 146                 }
 147                 if(d__<*dmin__)
 148                   *dmin__ = d__;
 149                 d__1 = emin, d__2 = z__[j4 - 1];
 150                 emin = (d__1<d__2) ? d__1 : d__2;
 151             }
 152         }
 153
 154         *dnm2 = d__;
 155         *dmin2 = *dmin__;
 156         j4 = 4*(*n0 - 2) - *pp;
 157         j4p2 = j4 + (*pp << 1) - 1;
 158         z__[j4 - 2] = *dnm2 + z__[j4p2];
 159         if (*dnm2 < 0.) {
 160             return;
 161         } else {
 162             z__[j4] = z__[j4p2 + 2] * (z__[j4p2] / z__[j4 - 2]);
 163             *dnm1 = z__[j4p2 + 2] * (*dnm2 / z__[j4 - 2]) - *tau;
 164         }
 165         if(*dnm1<*dmin__)
 166           *dmin__ = *dnm1;
 167
 168         *dmin1 = *dmin__;
 169         j4 += 4;
 170         j4p2 = j4 + (*pp << 1) - 1;
 171         z__[j4 - 2] = *dnm1 + z__[j4p2];
 172         if (*dnm1 < 0.) {
 173             return;
 174         } else {
 175             z__[j4] = z__[j4p2 + 2] * (z__[j4p2] / z__[j4 - 2]);
 176             *dn = z__[j4p2 + 2] * (*dnm1 / z__[j4 - 2]) - *tau;
 177         }
 178         if(*dn<*dmin__)
 179           *dmin__ = *dn;
 180
 181     }
 182
 183     z__[j4 + 2] = *dn;
 184     z__[(*n0 << 2) - *pp] = emin;
 185     return;
 186
 187 }
 188