src/gromacs/linearalgebra/gmx_lapack/slasq4.c

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