src/gromacs/linearalgebra/gmx_lapack/slarrvx.c

   1 #include <math.h>
   2
   3 #include "types/simple.h"
   4
   5 #include "../gmx_blas.h"
   6 #include "../gmx_lapack.h"
   7 #include "lapack_limits.h"
   8
   9
  10 void
  11 F77_FUNC(slarrvx,SLARRVX)(int *n,
  12         float *d__,
  13         float *l,
  14         int *isplit,
  15         int *m,
  16         float *w,
  17         int *iblock,
  18         int *indexw,
  19         float *gersch,
  20         float *tol,
  21         float *z__,
  22         int *ldz,
  23         int *isuppz,
  24         float *work,
  25         int *iwork,
  26         int *info)
  27 {
  28     int z_dim1, z_offset, i__1, i__2, i__3, i__4, i__5, i__6;
  29     float d__1, d__2;
  30     float c_b5 = 0.;
  31     int c__1 = 1;
  32     int c__2 = 2;
  33
  34     int i__, j, k, p, q;
  35     int im, in;
  36     float gap, eps, tmp;
  37     int zto;
  38     float ztz;
  39     int iend, jblk;
  40     int wend, iter, temp[1], ktot;
  41     int itmp1, itmp2;
  42     int indld;
  43     float sigma;
  44     int ndone, iinfo, iindr;
  45     float resid;
  46     int nomgs;
  47     int nclus;
  48     int zfrom, iindc1, iindc2;
  49     float lambda;
  50     int ibegin;
  51     int indgap, indlld;
  52     float mingma;
  53     int oldien, oldncl, wbegin;
  54     float relgap;
  55     int oldcls;
  56     int ndepth, inderr, iindwk;
  57     int newcls, oldfst;
  58     float minrgp=0.0;
  59     int indwrk, oldlst;
  60     float reltol;
  61     int newfrs, newftt, parity;
  62     float mgstol, nrminv, rqcorr;
  63     int newlst, newsiz;
  64
  65
  66     --d__;
  67     --l;
  68     --isplit;
  69     --w;
  70     --iblock;
  71     --indexw;
  72     --gersch;
  73     z_dim1 = *ldz;
  74     z_offset = 1 + z_dim1;
  75     z__ -= z_offset;
  76     --isuppz;
  77     --work;
  78     --iwork;
  79
  80     inderr = *n;
  81     indld = *n << 1;
  82     indlld = *n * 3;
  83     indgap = *n << 2;
  84     indwrk = *n * 5 + 1;
  85
  86     iindr = *n;
  87     iindc1 = *n << 1;
  88     iindc2 = *n * 3;
  89     iindwk = (*n << 2) + 1;
  90
  91     eps = GMX_FLOAT_EPS;
  92
  93     i__1 = *n << 1;
  94     for (i__ = 1; i__ <= i__1; ++i__) {
  95         iwork[i__] = 0;
  96     }
  97     F77_FUNC(slaset,SLASET)("Full", n, m, &c_b5, &c_b5, &z__[z_offset], ldz);
  98     mgstol = eps * 100.;
  99
 100     ibegin = 1;
 101     wbegin = 1;
 102     i__1 = iblock[*m];
 103     for (jblk = 1; jblk <= i__1; ++jblk) {
 104         iend = isplit[jblk];
 105
 106         wend = wbegin - 1;
 107 L171:
 108         if (wend < *m) {
 109             if (iblock[wend + 1] == jblk) {
 110                 ++wend;
 111                 goto L171;
 112             }
 113         }
 114         if (wend < wbegin) {
 115             ibegin = iend + 1;
 116             continue;
 117         }
 118
 119         if (ibegin == iend) {
 120             z__[ibegin + wbegin * z_dim1] = 1.;
 121             isuppz[(wbegin << 1) - 1] = ibegin;
 122             isuppz[wbegin * 2] = ibegin;
 123             ibegin = iend + 1;
 124             wbegin = wend + 1;
 125             continue;
 126         }
 127         oldien = ibegin - 1;
 128         in = iend - oldien;
 129         d__1 = .001, d__2 = 1. / (float) in;
 130         reltol = (d__1<d__2) ? d__1 : d__2;
 131         im = wend - wbegin + 1;
 132         F77_FUNC(scopy,SCOPY)(&im, &w[wbegin], &c__1, &work[1], &c__1);
 133         i__2 = im - 1;
 134         for (i__ = 1; i__ <= i__2; ++i__) {
 135             work[inderr + i__] = eps * fabs(work[i__]);
 136             work[indgap + i__] = work[i__ + 1] - work[i__];
 137         }
 138         work[inderr + im] = eps * fabs(work[im]);
 139         d__2 = fabs(work[im]);
 140         work[indgap + im] = (d__2>eps) ? d__2 : eps;
 141         ndone = 0;
 142
 143         ndepth = 0;
 144         parity = 1;
 145         nclus = 1;
 146         iwork[iindc1 + 1] = 1;
 147         iwork[iindc1 + 2] = im;
 148
 149 L40:
 150         if (ndone < im) {
 151             oldncl = nclus;
 152             nclus = 0;
 153             parity = 1 - parity;
 154             if (parity == 0) {
 155                 oldcls = iindc1;
 156                 newcls = iindc2;
 157             } else {
 158                 oldcls = iindc2;
 159                 newcls = iindc1;
 160             }
 161             i__2 = oldncl;
 162             for (i__ = 1; i__ <= i__2; ++i__) {
 163
 164                 j = oldcls + (i__ << 1);
 165                 oldfst = iwork[j - 1];
 166                 oldlst = iwork[j];
 167                 if (ndepth > 0) {
 168                     j = wbegin + oldfst - 1;
 169                     F77_FUNC(scopy,SCOPY)(&in, &z__[ibegin + j * z_dim1], &c__1, &d__[ibegin]
 170                             , &c__1);
 171                     i__3 = in - 1;
 172                     F77_FUNC(scopy,SCOPY)(&i__3, &z__[ibegin + (j + 1) * z_dim1], &c__1, &l[
 173                             ibegin], &c__1);
 174                     F77_FUNC(slaset,SLASET)("Full", &in, &c__2, &c_b5, &c_b5, &z__[ibegin + j
 175                             * z_dim1], ldz);
 176                 }
 177                 k = ibegin;
 178                 i__3 = in - 1;
 179                 for (j = 1; j <= i__3; ++j) {
 180                     tmp = d__[k] * l[k];
 181                     work[indld + j] = tmp;
 182                     work[indlld + j] = tmp * l[k];
 183                     ++k;
 184                 }
 185                 if (ndepth > 0) {
 186
 187                     p = indexw[wbegin - 1 + oldfst];
 188                     q = indexw[wbegin - 1 + oldlst];
 189                     d__1 = eps * 4.;
 190                     i__3 = p - oldfst;
 191                     F77_FUNC(slarrbx,SLARRBX)(&in, &d__[ibegin], &l[ibegin], &work[indld + 1], &
 192                             work[indlld + 1], &p, &q, &reltol, &d__1, &i__3, &
 193                             work[1], &work[indgap + 1], &work[inderr + 1], &
 194                             work[indwrk + in], &iwork[iindwk], &iinfo);
 195                 }
 196                 newfrs = oldfst;
 197                 i__3 = oldlst;
 198                 for (j = oldfst; j <= i__3; ++j) {
 199                     if (j == oldlst || work[indgap + j] >=
 200                         reltol * fabs(work[j])) {
 201                         newlst = j;
 202                     } else {
 203
 204                         relgap = work[indgap + j] / fabs(work[j]);
 205                         if (j == newfrs) {
 206                             minrgp = relgap;
 207                         } else {
 208                             minrgp = (minrgp<relgap) ? minrgp : relgap;
 209                         }
 210                         continue;
 211                     }
 212                     newsiz = newlst - newfrs + 1;
 213                     newftt = wbegin + newfrs - 1;
 214                     nomgs = newsiz == 1 || newsiz > 1 || minrgp < mgstol;
 215                     if (newsiz > 1 && nomgs) {
 216
 217                         F77_FUNC(slarrfx,SLARRFX)(&in, &d__[ibegin], &l[ibegin], &work[indld +
 218                                 1], &work[indlld + 1], &newfrs, &newlst, &
 219                                 work[1], &sigma, &z__[ibegin + newftt *
 220                                 z_dim1], &z__[ibegin + (newftt + 1) * z_dim1],
 221                                  &work[indwrk], info);
 222                         if (*info == 0) {
 223                             tmp = eps * fabs(sigma);
 224                             i__4 = newlst;
 225                             for (k = newfrs; k <= i__4; ++k) {
 226                                 work[k] -= sigma;
 227                                 d__1 = work[indgap + k];
 228                                 work[indgap + k] = (d__1>tmp) ? d__1 : tmp;
 229                                 work[inderr + k] += tmp;
 230                             }
 231                             ++nclus;
 232                             k = newcls + (nclus << 1);
 233                             iwork[k - 1] = newfrs;
 234                             iwork[k] = newlst;
 235                         } else {
 236                             *info = 0;
 237                             if (minrgp >= mgstol) {
 238                                 nomgs = 0;
 239                             } else {
 240
 241                                 work[indwrk] = d__[ibegin];
 242                                 i__4 = in - 1;
 243                                 for (k = 1; k <= i__4; ++k) {
 244                                     work[indwrk + k] = d__[ibegin + k] + work[
 245                                             indlld + k];
 246                                 }
 247                                 i__4 = newsiz;
 248                                 for (k = 1; k <= i__4; ++k) {
 249                                     iwork[iindwk + k - 1] = 1;
 250                                 }
 251                                 i__4 = newlst;
 252                                 for (k = newfrs; k <= i__4; ++k) {
 253                                     isuppz[2*(oldien + k) - 1] = 1;
 254                                     isuppz[(oldien + k) * 2] = in;
 255                                 }
 256                                 temp[0] = in;
 257                                 F77_FUNC(sstein,SSTEIN)(&in, &work[indwrk], &work[indld + 1],
 258                                         &newsiz, &work[newfrs], &iwork[iindwk]
 259                                         , temp, &z__[ibegin + newftt * z_dim1]
 260                                         , ldz, &work[indwrk + in], &iwork[
 261                                         iindwk + in], &iwork[iindwk + (in*2)], &iinfo);
 262                                 if (iinfo != 0) {
 263                                     *info = 2;
 264                                     return;
 265                                 }
 266                                 ndone += newsiz;
 267                             }
 268                         }
 269                     } else {
 270                         ktot = newftt;
 271                         i__4 = newlst;
 272                         for (k = newfrs; k <= i__4; ++k) {
 273                             iter = 0;
 274 L90:
 275                             lambda = work[k];
 276
 277                             F77_FUNC(slar1vx,SLAR1VX)(&in, &c__1, &in, &lambda, &d__[ibegin], &
 278                                     l[ibegin], &work[indld + 1], &work[indlld
 279                                     + 1], &w[wbegin + k - 1], &gersch[(oldien
 280                                     << 1) + 1], &z__[ibegin + ktot * z_dim1],
 281                                     &ztz, &mingma, &iwork[iindr + ktot], &
 282                                     isuppz[(ktot << 1) - 1], &work[indwrk]);
 283                             tmp = 1. / ztz;
 284                             nrminv = sqrt(tmp);
 285                             resid = fabs(mingma) * nrminv;
 286                             rqcorr = mingma * tmp;
 287                             if (k == in) {
 288                                 gap = work[indgap + k - 1];
 289                             } else if (k == 1) {
 290                                 gap = work[indgap + k];
 291                             } else {
 292                                 d__1 = work[indgap + k - 1], d__2 = work[
 293                                         indgap + k];
 294                                 gap = (d__1<d__2) ? d__1 : d__2;
 295                             }
 296                             ++iter;
 297                             if (resid > *tol * gap && fabs(rqcorr) > eps * 4. *
 298                                      fabs(lambda)) {
 299                                 work[k] = lambda + rqcorr;
 300                                 if (iter < 8) {
 301                                     goto L90;
 302                                 }
 303                             }
 304                             iwork[ktot] = 1;
 305                             if (newsiz == 1) {
 306                                 ++ndone;
 307                             }
 308                             zfrom = isuppz[(ktot << 1) - 1];
 309                             zto = isuppz[ktot * 2];
 310                             i__5 = zto - zfrom + 1;
 311                             F77_FUNC(sscal,SSCAL)(&i__5, &nrminv, &z__[ibegin + zfrom - 1 +
 312                                     ktot * z_dim1], &c__1);
 313                             ++ktot;
 314                         }
 315                         if (newsiz > 1) {
 316                             itmp1 = isuppz[(newftt << 1) - 1];
 317                             itmp2 = isuppz[newftt * 2];
 318                             ktot = oldien + newlst;
 319                             i__4 = ktot;
 320                             for (p = newftt + 1; p <= i__4; ++p) {
 321                                 i__5 = p - 1;
 322                                 for (q = newftt; q <= i__5; ++q) {
 323                                     tmp = -F77_FUNC(sdot,SDOT)(&in, &z__[ibegin + p *
 324                                             z_dim1], &c__1, &z__[ibegin + q *
 325                                             z_dim1], &c__1);
 326                                     F77_FUNC(saxpy,SAXPY)(&in, &tmp, &z__[ibegin + q *
 327                                             z_dim1], &c__1, &z__[ibegin + p *
 328                                             z_dim1], &c__1);
 329                                 }
 330                                 tmp = 1. / F77_FUNC(snrm2,SNRM2)(&in, &z__[ibegin + p *
 331                                         z_dim1], &c__1);
 332                                 F77_FUNC(sscal,SSCAL)(&in, &tmp, &z__[ibegin + p * z_dim1], &
 333                                         c__1);
 334                                 i__5 = itmp1, i__6 = isuppz[(p << 1) - 1];
 335                                 itmp1 = (i__5<i__6) ? i__5 : i__6;
 336                                 i__5 = itmp2, i__6 = isuppz[p * 2];
 337                                 itmp2 = (i__5>i__6) ? i__5 : i__6;
 338                             }
 339                             i__4 = ktot;
 340                             for (p = newftt; p <= i__4; ++p) {
 341                                 isuppz[(p << 1) - 1] = itmp1;
 342                                 isuppz[p * 2] = itmp2;
 343                             }
 344                             ndone += newsiz;
 345                         }
 346                     }
 347                     newfrs = j + 1;
 348                 }
 349             }
 350             ++ndepth;
 351             goto L40;
 352         }
 353         j = wbegin << 1;
 354         i__2 = wend;
 355         for (i__ = wbegin; i__ <= i__2; ++i__) {
 356             isuppz[j - 1] += oldien;
 357             isuppz[j] += oldien;
 358             j += 2;
 359
 360         }
 361         ibegin = iend + 1;
 362         wbegin = wend + 1;
 363     }
 364
 365     return;
 366
 367 }