src/gmxlib/gmx_lapack/dlasq4.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
  40 void
  41 F77_FUNC(dlasq4,DLASQ4)(int *i0,
  42         int *n0,
  43         double *z__,
  44         int *pp,
  45         int *n0in,
  46         double *dmin__,
  47         double *dmin1,
  48         double *dmin2,
  49         double *dn,
  50         double *dn1,
  51         double *dn2,
  52         double *tau,
  53         int *ttype)
  54 {
  55     double g = 0.;
  56     int i__1;
  57     double d__1, d__2;
  58
  59     double s, a2, b1, b2;
  60     int i4, nn, np;
  61     double gam, gap1, gap2;
  62
  63
  64     if (*dmin__ <= 0.) {
  65         *tau = -(*dmin__);
  66         *ttype = -1;
  67         return;
  68     }
  69
  70     s = 0.0;
  71
  72     nn = (*n0 << 2) + *pp;
  73     if (*n0in == *n0) {
  74
  75         if ( fabs(*dmin__ - *dn)<GMX_DOUBLE_EPS*fabs(*dmin__ + *dn) ||
  76          fabs(*dmin__ - *dn1)<GMX_DOUBLE_EPS*fabs(*dmin__ + *dn1)) {
  77
  78             b1 = sqrt(z__[nn - 3]) * sqrt(z__[nn - 5]);
  79             b2 = sqrt(z__[nn - 7]) * sqrt(z__[nn - 9]);
  80             a2 = z__[nn - 7] + z__[nn - 5];
  81
  82         if ( fabs(*dmin__ - *dn)<GMX_DOUBLE_EPS*fabs(*dmin__ + *dn) &&
  83              fabs(*dmin1 - *dn1)<GMX_DOUBLE_EPS*fabs(*dmin1 + *dn1)) {
  84
  85             gap2 = *dmin2 - a2 - *dmin2 * .25;
  86                 if (gap2 > 0. && gap2 > b2) {
  87                     gap1 = a2 - *dn - b2 / gap2 * b2;
  88                 } else {
  89                     gap1 = a2 - *dn - (b1 + b2);
  90                 }
  91                 if (gap1 > 0. && gap1 > b1) {
  92                     d__1 = *dn - b1 / gap1 * b1, d__2 = *dmin__ * .5;
  93                     s = (d__1>d__2) ? d__1 : d__2;
  94                     *ttype = -2;
  95                 } else {
  96                     s = 0.;
  97                     if (*dn > b1) {
  98                         s = *dn - b1;
  99                     }
 100                     if (a2 > b1 + b2) {
 101                         d__1 = s, d__2 = a2 - (b1 + b2);
 102                         s = (d__1<d__2) ? d__1 : d__2;
 103                     }
 104                     d__1 = s, d__2 = *dmin__ * .333;
 105                     s = (d__1>d__2) ? d__1 : d__2;
 106                     *ttype = -3;
 107                 }
 108             } else {
 109
 110
 111                 *ttype = -4;
 112                 s = *dmin__ * .25;
 113                 if (fabs(*dmin__ - *dn)<GMX_DOUBLE_EPS*fabs(*dmin__ + *dn)) {
 114                     gam = *dn;
 115                     a2 = 0.;
 116                     if (z__[nn - 5] > z__[nn - 7]) {
 117                         return;
 118                     }
 119                     b2 = z__[nn - 5] / z__[nn - 7];
 120                     np = nn - 9;
 121                 } else {
 122                     np = nn - (*pp << 1);
 123                     b2 = z__[np - 2];
 124                     gam = *dn1;
 125                     if (z__[np - 4] > z__[np - 2]) {
 126                         return;
 127                     }
 128                     a2 = z__[np - 4] / z__[np - 2];
 129                     if (z__[nn - 9] > z__[nn - 11]) {
 130                         return;
 131                     }
 132                     b2 = z__[nn - 9] / z__[nn - 11];
 133                     np = nn - 13;
 134                 }
 135
 136
 137                 a2 += b2;
 138                 i__1 = (*i0 << 2) - 1 + *pp;
 139                 for (i4 = np; i4 >= i__1; i4 += -4) {
 140                     if (fabs(b2)<GMX_DOUBLE_MIN) {
 141                         goto L20;
 142                     }
 143                     b1 = b2;
 144                     if (z__[i4] > z__[i4 - 2]) {
 145                         return;
 146                     }
 147                     b2 *= z__[i4] / z__[i4 - 2];
 148                     a2 += b2;
 149                     if (((b2>b1) ? b2 : b1) * 100. < a2 || .563 < a2) {
 150                         goto L20;
 151                     }
 152                 }
 153 L20:
 154                 a2 *= 1.05;
 155
 156
 157                 if (a2 < .563) {
 158                     s = gam * (1. - sqrt(a2)) / (a2 + 1.);
 159                 }
 160             }
 161         } else if (fabs(*dmin__ - *dn2)<GMX_DOUBLE_EPS*fabs(*dmin__ + *dn2)) {
 162
 163             *ttype = -5;
 164             s = *dmin__ * .25;
 165
 166             np = nn - (*pp << 1);
 167             b1 = z__[np - 2];
 168             b2 = z__[np - 6];
 169             gam = *dn2;
 170             if (z__[np - 8] > b2 || z__[np - 4] > b1) {
 171                 return;
 172             }
 173             a2 = z__[np - 8] / b2 * (z__[np - 4] / b1 + 1.);
 174
 175
 176             if (*n0 - *i0 > 2) {
 177                 b2 = z__[nn - 13] / z__[nn - 15];
 178                 a2 += b2;
 179                 i__1 = (*i0 << 2) - 1 + *pp;
 180                 for (i4 = nn - 17; i4 >= i__1; i4 += -4) {
 181                     if (fabs(b2)<GMX_DOUBLE_MIN) {
 182                         goto L40;
 183                     }
 184                     b1 = b2;
 185                     if (z__[i4] > z__[i4 - 2]) {
 186                         return;
 187                     }
 188                     b2 *= z__[i4] / z__[i4 - 2];
 189                     a2 += b2;
 190                     if (((b2>b1) ? b2 : b1) * 100. < a2 || .563 < a2) {
 191                         goto L40;
 192                     }
 193                 }
 194 L40:
 195                 a2 *= 1.05;
 196             }
 197
 198             if (a2 < .563) {
 199                 s = gam * (1. - sqrt(a2)) / (a2 + 1.);
 200             }
 201         } else {
 202
 203             if (*ttype == -6) {
 204                 g += (1. - g) * .333;
 205             } else if (*ttype == -18) {
 206                 g = .083250000000000005;
 207             } else {
 208                 g = .25;
 209             }
 210             s = g * *dmin__;
 211             *ttype = -6;
 212         }
 213
 214     } else if (*n0in == *n0 + 1) {
 215
 216         if ( fabs(*dmin1 - *dn1)<GMX_DOUBLE_EPS*fabs(*dmin1 + *dn1) &&
 217              fabs(*dmin2 - *dn2)<GMX_DOUBLE_EPS*fabs(*dmin2 + *dn2)) {
 218
 219             *ttype = -7;
 220             s = *dmin1 * .333;
 221             if (z__[nn - 5] > z__[nn - 7]) {
 222                 return;
 223             }
 224             b1 = z__[nn - 5] / z__[nn - 7];
 225             b2 = b1;
 226             if (fabs(b2)<GMX_DOUBLE_MIN) {
 227                 goto L60;
 228             }
 229             i__1 = (*i0 << 2) - 1 + *pp;
 230             for (i4 = (*n0 << 2) - 9 + *pp; i4 >= i__1; i4 += -4) {
 231                 a2 = b1;
 232                 if (z__[i4] > z__[i4 - 2]) {
 233                     return;
 234                 }
 235                 b1 *= z__[i4] / z__[i4 - 2];
 236                 b2 += b1;
 237                 if (((a2>b1) ? a2 : b1) * 100. < b2) {
 238                     goto L60;
 239                 }
 240             }
 241 L60:
 242             b2 = sqrt(b2 * 1.05);
 243             d__1 = b2;
 244             a2 = *dmin1 / (d__1 * d__1 + 1.);
 245             gap2 = *dmin2 * .5 - a2;
 246             if (gap2 > 0. && gap2 > b2 * a2) {
 247                 d__1 = s, d__2 = a2 * (1. - a2 * 1.01 * (b2 / gap2) * b2);
 248                 s = (d__1>d__2) ? d__1 : d__2;
 249             } else {
 250                 d__1 = s, d__2 = a2 * (1. - b2 * 1.01);
 251                 s = (d__1>d__2) ? d__1 : d__2;
 252                 *ttype = -8;
 253             }
 254         } else {
 255
 256             s = *dmin1 * .25;
 257             if (fabs(*dmin1 - *dn1)<GMX_DOUBLE_EPS*fabs(*dmin1 + *dn1)) {
 258                 s = *dmin1 * .5;
 259             }
 260             *ttype = -9;
 261         }
 262
 263     } else if (*n0in == *n0 + 2) {
 264
 265         if (fabs(*dmin2 - *dn2)<GMX_DOUBLE_EPS*fabs(*dmin2 + *dn2) &&
 266         z__[nn - 5] * 2. < z__[nn - 7]) {
 267             *ttype = -10;
 268             s = *dmin2 * .333;
 269             if (z__[nn - 5] > z__[nn - 7]) {
 270                 return;
 271             }
 272             b1 = z__[nn - 5] / z__[nn - 7];
 273             b2 = b1;
 274             if (fabs(b2)<GMX_DOUBLE_MIN) {
 275                 goto L80;
 276             }
 277             i__1 = (*i0 << 2) - 1 + *pp;
 278             for (i4 = (*n0 << 2) - 9 + *pp; i4 >= i__1; i4 += -4) {
 279                 if (z__[i4] > z__[i4 - 2]) {
 280                     return;
 281                 }
 282                 b1 *= z__[i4] / z__[i4 - 2];
 283                 b2 += b1;
 284                 if (b1 * 100. < b2) {
 285                     goto L80;
 286                 }
 287             }
 288 L80:
 289             b2 = sqrt(b2 * 1.05);
 290             d__1 = b2;
 291             a2 = *dmin2 / (d__1 * d__1 + 1.);
 292             gap2 = z__[nn - 7] + z__[nn - 9] - sqrt(z__[nn - 11]) * sqrt(z__[
 293                     nn - 9]) - a2;
 294             if (gap2 > 0. && gap2 > b2 * a2) {
 295                 d__1 = s, d__2 = a2 * (1. - a2 * 1.01 * (b2 / gap2) * b2);
 296                 s = (d__1>d__2) ? d__1 : d__2;
 297             } else {
 298                 d__1 = s, d__2 = a2 * (1. - b2 * 1.01);
 299                 s = (d__1>d__2) ? d__1 : d__2;
 300             }
 301         } else {
 302             s = *dmin2 * .25;
 303             *ttype = -11;
 304         }
 305     } else if (*n0in > *n0 + 2) {
 306
 307         s = 0.;
 308         *ttype = -12;
 309     }
 310
 311     *tau = s;
 312     return;
 313
 314 }
 315
 316