src/gmxlib/gmx_lapack/slaebz.c

   1 #include <math.h>
   2 #include "gmx_lapack.h"
   3
   4 void
   5 F77_FUNC(slaebz,SLAEBZ)(int *ijob,
   6         int *nitmax,
   7         int *n,
   8         int *mmax,
   9         int *minp,
  10         int *nbmin,
  11         float *abstol,
  12         float *reltol,
  13         float *pivmin,
  14         float *d__,
  15         float *e,
  16         float *e2,
  17         int *nval,
  18         float *ab,
  19         float *c__,
  20         int *mout,
  21         int *nab,
  22         float *work,
  23         int *iwork,
  24         int *info)
  25 {
  26     int nab_dim1, nab_offset, ab_dim1, ab_offset, i__1, i__2, i__3, i__4,
  27             i__5, i__6;
  28     float d__1, d__2, d__3, d__4;
  29
  30     int j, kf, ji, kl, jp, jit;
  31     float tmp1, tmp2;
  32     int itmp1, itmp2, kfnew, klnew;
  33
  34     nab_dim1 = *mmax;
  35     nab_offset = 1 + nab_dim1;
  36     nab -= nab_offset;
  37     ab_dim1 = *mmax;
  38     ab_offset = 1 + ab_dim1;
  39     ab -= ab_offset;
  40     --d__;
  41     --e;
  42     --e2;
  43     --nval;
  44     --c__;
  45     --work;
  46     --iwork;
  47
  48     *info = 0;
  49     if (*ijob < 1 || *ijob > 3) {
  50         *info = -1;
  51         return;
  52     }
  53
  54     if (*ijob == 1) {
  55
  56         *mout = 0;
  57
  58         i__1 = *minp;
  59         for (ji = 1; ji <= i__1; ++ji) {
  60             for (jp = 1; jp <= 2; ++jp) {
  61                 tmp1 = d__[1] - ab[ji + jp * ab_dim1];
  62                 if (fabs(tmp1) < *pivmin) {
  63                     tmp1 = -(*pivmin);
  64                 }
  65                 nab[ji + jp * nab_dim1] = 0;
  66                 if (tmp1 <= 0.) {
  67                     nab[ji + jp * nab_dim1] = 1;
  68                 }
  69
  70                 i__2 = *n;
  71                 for (j = 2; j <= i__2; ++j) {
  72                     tmp1 = d__[j] - e2[j - 1] / tmp1 - ab[ji + jp * ab_dim1];
  73                     if (fabs(tmp1) < *pivmin) {
  74                         tmp1 = -(*pivmin);
  75                     }
  76                     if (tmp1 <= 0.) {
  77                         ++nab[ji + jp * nab_dim1];
  78                     }
  79                 }
  80             }
  81             *mout = *mout + nab[ji + (nab_dim1 << 1)] - nab[ji + nab_dim1];
  82         }
  83         return;
  84     }
  85
  86     kf = 1;
  87     kl = *minp;
  88
  89     if (*ijob == 2) {
  90         i__1 = *minp;
  91         for (ji = 1; ji <= i__1; ++ji) {
  92             c__[ji] = (ab[ji + ab_dim1] + ab[ji + (ab_dim1 << 1)]) * .5;
  93         }
  94     }
  95
  96     i__1 = *nitmax;
  97     for (jit = 1; jit <= i__1; ++jit) {
  98
  99         if (kl - kf + 1 >= *nbmin && *nbmin > 0) {
 100
 101             i__2 = kl;
 102             for (ji = kf; ji <= i__2; ++ji) {
 103
 104                 work[ji] = d__[1] - c__[ji];
 105                 iwork[ji] = 0;
 106                 if (work[ji] <= *pivmin) {
 107                     iwork[ji] = 1;
 108                     d__1 = work[ji], d__2 = -(*pivmin);
 109                     work[ji] = (d__1<d__2) ? d__1 : d__2;
 110                 }
 111
 112                 i__3 = *n;
 113                 for (j = 2; j <= i__3; ++j) {
 114                     work[ji] = d__[j] - e2[j - 1] / work[ji] - c__[ji];
 115                     if (work[ji] <= *pivmin) {
 116                         ++iwork[ji];
 117                         d__1 = work[ji], d__2 = -(*pivmin);
 118                         work[ji] = (d__1<d__2) ? d__1 : d__2;
 119                     }
 120                 }
 121             }
 122
 123             if (*ijob <= 2) {
 124
 125                 klnew = kl;
 126                 i__2 = kl;
 127                 for (ji = kf; ji <= i__2; ++ji) {
 128
 129                   i__5 = nab[ji + nab_dim1];
 130                   i__6 = iwork[ji];
 131                   i__3 = nab[ji + (nab_dim1 << 1)];
 132                   i__4 = (i__5>i__6) ? i__5 : i__6;
 133                     iwork[ji] = (i__3<i__4) ? i__3 : i__4;
 134
 135                     if (iwork[ji] == nab[ji + (nab_dim1 << 1)]) {
 136
 137                         ab[ji + (ab_dim1 << 1)] = c__[ji];
 138
 139                     } else if (iwork[ji] == nab[ji + nab_dim1]) {
 140
 141                         ab[ji + ab_dim1] = c__[ji];
 142                     } else {
 143                         ++klnew;
 144                         if (klnew <= *mmax) {
 145
 146                             ab[klnew + (ab_dim1 << 1)] = ab[ji + (ab_dim1 <<
 147                                     1)];
 148                             nab[klnew + (nab_dim1 << 1)] = nab[ji + (nab_dim1
 149                                     << 1)];
 150                             ab[klnew + ab_dim1] = c__[ji];
 151                             nab[klnew + nab_dim1] = iwork[ji];
 152                             ab[ji + (ab_dim1 << 1)] = c__[ji];
 153                             nab[ji + (nab_dim1 << 1)] = iwork[ji];
 154                         } else {
 155                             *info = *mmax + 1;
 156                         }
 157                     }
 158                 }
 159                 if (*info != 0) {
 160                     return;
 161                 }
 162                 kl = klnew;
 163             } else {
 164
 165                 i__2 = kl;
 166                 for (ji = kf; ji <= i__2; ++ji) {
 167                     if (iwork[ji] <= nval[ji]) {
 168                         ab[ji + ab_dim1] = c__[ji];
 169                         nab[ji + nab_dim1] = iwork[ji];
 170                     }
 171                     if (iwork[ji] >= nval[ji]) {
 172                         ab[ji + (ab_dim1 << 1)] = c__[ji];
 173                         nab[ji + (nab_dim1 << 1)] = iwork[ji];
 174                     }
 175                 }
 176             }
 177
 178         } else {
 179
 180             klnew = kl;
 181             i__2 = kl;
 182             for (ji = kf; ji <= i__2; ++ji) {
 183
 184                 tmp1 = c__[ji];
 185                 tmp2 = d__[1] - tmp1;
 186                 itmp1 = 0;
 187                 if (tmp2 <= *pivmin) {
 188                     itmp1 = 1;
 189                     d__1 = tmp2, d__2 = -(*pivmin);
 190                     tmp2 = (d__1<d__2) ? d__1 : d__2;
 191                 }
 192
 193                 i__3 = *n;
 194                 for (j = 2; j <= i__3; ++j) {
 195                     tmp2 = d__[j] - e2[j - 1] / tmp2 - tmp1;
 196                     if (tmp2 <= *pivmin) {
 197                         ++itmp1;
 198                         d__1 = tmp2, d__2 = -(*pivmin);
 199                         tmp2 = (d__1<d__2) ? d__1 : d__2;
 200                     }
 201                 }
 202
 203                 if (*ijob <= 2) {
 204
 205                     i__5 = nab[ji + nab_dim1];
 206                     i__3 = nab[ji + (nab_dim1 << 1)];
 207                     i__4 = (i__5>itmp1) ? i__5 : itmp1;
 208                     itmp1 = (i__3<i__4) ? i__3 : i__4;
 209
 210                     if (itmp1 == nab[ji + (nab_dim1 << 1)]) {
 211
 212                         ab[ji + (ab_dim1 << 1)] = tmp1;
 213
 214                     } else if (itmp1 == nab[ji + nab_dim1]) {
 215
 216                         ab[ji + ab_dim1] = tmp1;
 217                     } else if (klnew < *mmax) {
 218
 219                         ++klnew;
 220                         ab[klnew + (ab_dim1 << 1)] = ab[ji + (ab_dim1 << 1)];
 221                         nab[klnew + (nab_dim1 << 1)] = nab[ji + (nab_dim1 <<
 222                                 1)];
 223                         ab[klnew + ab_dim1] = tmp1;
 224                         nab[klnew + nab_dim1] = itmp1;
 225                         ab[ji + (ab_dim1 << 1)] = tmp1;
 226                         nab[ji + (nab_dim1 << 1)] = itmp1;
 227                     } else {
 228                         *info = *mmax + 1;
 229                         return;
 230                     }
 231                 } else {
 232
 233                     if (itmp1 <= nval[ji]) {
 234                         ab[ji + ab_dim1] = tmp1;
 235                         nab[ji + nab_dim1] = itmp1;
 236                     }
 237                     if (itmp1 >= nval[ji]) {
 238                         ab[ji + (ab_dim1 << 1)] = tmp1;
 239                         nab[ji + (nab_dim1 << 1)] = itmp1;
 240                     }
 241                 }
 242             }
 243             kl = klnew;
 244
 245         }
 246
 247         kfnew = kf;
 248         i__2 = kl;
 249         for (ji = kf; ji <= i__2; ++ji) {
 250             tmp1 = fabs(ab[ji + (ab_dim1 << 1)] - ab[ji + ab_dim1]);
 251             d__3 = fabs(ab[ji + (ab_dim1 << 1)]);
 252             d__4 = fabs(ab[ji + ab_dim1]);
 253             tmp2 = (d__3>d__4) ? d__3 : d__4;
 254             d__1 = (*abstol>*pivmin) ? *abstol : *pivmin;
 255             d__2 = *reltol * tmp2;
 256             if (tmp1 < ((d__1>d__2) ? d__1 : d__2) || nab[ji + nab_dim1] >= nab[ji + (
 257                     nab_dim1 << 1)]) {
 258
 259                 if (ji > kfnew) {
 260                     tmp1 = ab[ji + ab_dim1];
 261                     tmp2 = ab[ji + (ab_dim1 << 1)];
 262                     itmp1 = nab[ji + nab_dim1];
 263                     itmp2 = nab[ji + (nab_dim1 << 1)];
 264                     ab[ji + ab_dim1] = ab[kfnew + ab_dim1];
 265                     ab[ji + (ab_dim1 << 1)] = ab[kfnew + (ab_dim1 << 1)];
 266                     nab[ji + nab_dim1] = nab[kfnew + nab_dim1];
 267                     nab[ji + (nab_dim1 << 1)] = nab[kfnew + (nab_dim1 << 1)];
 268                     ab[kfnew + ab_dim1] = tmp1;
 269                     ab[kfnew + (ab_dim1 << 1)] = tmp2;
 270                     nab[kfnew + nab_dim1] = itmp1;
 271                     nab[kfnew + (nab_dim1 << 1)] = itmp2;
 272                     if (*ijob == 3) {
 273                         itmp1 = nval[ji];
 274                         nval[ji] = nval[kfnew];
 275                         nval[kfnew] = itmp1;
 276                     }
 277                 }
 278                 ++kfnew;
 279             }
 280         }
 281         kf = kfnew;
 282
 283         i__2 = kl;
 284         for (ji = kf; ji <= i__2; ++ji) {
 285             c__[ji] = (ab[ji + ab_dim1] + ab[ji + (ab_dim1 << 1)]) * .5;
 286         }
 287
 288         if (kf > kl) {
 289             break;
 290         }
 291     }
 292
 293     i__1 = kl + 1 - kf;
 294     if(i__1>0)
 295       *info = i__1;
 296
 297     *mout = kl;
 298
 299     return;
 300
 301 }
 302
 303