src/gromacs/linearalgebra/gmx_lapack/dstegr.c

   1 #include <math.h>
   2 #include "../gmx_blas.h"
   3 #include "../gmx_lapack.h"
   4 #include "lapack_limits.h"
   5
   6 #include "types/simple.h"
   7
   8 void
   9 F77_FUNC(dstegr,DSTEGR)(const char *jobz,
  10         const char *range,
  11         int *n,
  12         double *d__,
  13         double *e,
  14         double *vl,
  15         double *vu,
  16         int *il,
  17         int *iu,
  18         double *abstol,
  19         int *m,
  20         double *w,
  21         double *z__,
  22         int *ldz,
  23         int *isuppz,
  24         double *work,
  25         int *lwork,
  26         int *iwork,
  27         int *liwork,
  28         int *info)
  29 {
  30     int z_dim1, z_offset, i__1, i__2;
  31     double d__1, d__2;
  32     int c__1 = 1;
  33
  34     int i__, j;
  35     int jj;
  36     double eps, tol, tmp, rmin, rmax;
  37     int itmp;
  38     double tnrm;
  39     double scale;
  40     int iinfo, iindw;
  41     int lwmin;
  42     int wantz;
  43     int iindbl;
  44     int valeig,alleig,indeig;
  45     double safmin,minval;
  46     double bignum;
  47     int iindwk, indgrs;
  48     double thresh;
  49     int iinspl, indwrk, liwmin, nsplit;
  50     double smlnum;
  51     int lquery;
  52
  53
  54     --d__;
  55     --e;
  56     --w;
  57     z_dim1 = *ldz;
  58     z_offset = 1 + z_dim1;
  59     z__ -= z_offset;
  60     --isuppz;
  61     --work;
  62     --iwork;
  63
  64     wantz = (*jobz=='V' || *jobz=='v');
  65     alleig = (*range=='A' || *range=='a');
  66     valeig = (*range=='V' || *range=='v');
  67     indeig = (*range=='I' || *range=='i');
  68
  69     lquery = *lwork == -1 || *liwork == -1;
  70     lwmin = *n * 17;
  71     liwmin = *n * 10;
  72
  73     *info = 0;
  74     if (! (wantz || (*jobz=='N' || *jobz=='n'))) {
  75         *info = -1;
  76     } else if (! (alleig || valeig || indeig)) {
  77         *info = -2;
  78     } else if (*n < 0) {
  79         *info = -3;
  80     } else if (valeig && *n > 0 && *vu <= *vl) {
  81         *info = -7;
  82     } else if (indeig && (*il < 1 || *il > *n)) {
  83         *info = -8;
  84     } else if (indeig && (*iu < *il || *iu > *n)) {
  85         *info = -9;
  86     } else if (*ldz < 1 || (wantz && *ldz < *n)) {
  87         *info = -14;
  88     } else if (*lwork < lwmin && ! lquery) {
  89         *info = -17;
  90     } else if (*liwork < liwmin && ! lquery) {
  91         *info = -19;
  92     }
  93     if (*info == 0) {
  94         work[1] = (double) lwmin;
  95         iwork[1] = liwmin;
  96     }
  97
  98     if (*info != 0) {
  99         i__1 = -(*info);
 100         return;
 101     } else if (lquery) {
 102         return;
 103     }
 104
 105     *m = 0;
 106     if (*n == 0) {
 107         return;
 108     }
 109
 110     if (*n == 1) {
 111         if (alleig || indeig) {
 112             *m = 1;
 113             w[1] = d__[1];
 114         } else {
 115             if (*vl < d__[1] && *vu >= d__[1]) {
 116                 *m = 1;
 117                 w[1] = d__[1];
 118             }
 119         }
 120         if (wantz) {
 121             z__[z_dim1 + 1] = 1.;
 122         }
 123         return;
 124     }
 125
 126     minval = GMX_DOUBLE_MIN;
 127     safmin = minval*(1.0+GMX_DOUBLE_EPS);
 128     eps = GMX_DOUBLE_EPS;
 129     smlnum = safmin / eps;
 130     bignum = 1. / smlnum;
 131     rmin = sqrt(smlnum);
 132     d__1 = sqrt(bignum), d__2 = 1. / sqrt(sqrt(safmin));
 133     rmax = (d__1<d__2) ? d__1 : d__2;
 134     scale = 1.;
 135     tnrm = F77_FUNC(dlanst,DLANST)("M", n, &d__[1], &e[1]);
 136     if (tnrm > 0. && tnrm < rmin) {
 137         scale = rmin / tnrm;
 138     } else if (tnrm > rmax) {
 139         scale = rmax / tnrm;
 140     }
 141     if ( fabs(scale-1.0)>GMX_DOUBLE_EPS) {
 142         F77_FUNC(dscal,DSCAL)(n, &scale, &d__[1], &c__1);
 143         i__1 = *n - 1;
 144         F77_FUNC(dscal,DSCAL)(&i__1, &scale, &e[1], &c__1);
 145         tnrm *= scale;
 146     }
 147     indgrs = 1;
 148     indwrk = (*n << 1) + 1;
 149
 150     iinspl = 1;
 151     iindbl = *n + 1;
 152     iindw = (*n << 1) + 1;
 153     iindwk = *n * 3 + 1;
 154
 155     thresh = eps * tnrm;
 156     F77_FUNC(dlarrex,DLARREX)(range, n, vl, vu, il, iu, &d__[1], &e[1], &thresh, &nsplit, &
 157             iwork[iinspl], m, &w[1], &iwork[iindbl], &iwork[iindw], &work[
 158             indgrs], &work[indwrk], &iwork[iindwk], &iinfo);
 159
 160     if (iinfo != 0) {
 161         *info = 1;
 162         return;
 163     }
 164
 165     if (wantz) {
 166         d__1 = *abstol, d__2 = (double) (*n) * eps;
 167         tol = (d__1>d__2) ? d__1 : d__2;
 168         F77_FUNC(dlarrvx,DLARRVX)(n, &d__[1], &e[1], &iwork[iinspl], m, &w[1], &iwork[iindbl], &
 169                 iwork[iindw], &work[indgrs], &tol, &z__[z_offset], ldz, &
 170                 isuppz[1], &work[indwrk], &iwork[iindwk], &iinfo);
 171         if (iinfo != 0) {
 172             *info = 2;
 173             return;
 174         }
 175     }
 176
 177     i__1 = *m;
 178     for (j = 1; j <= i__1; ++j) {
 179         itmp = iwork[iindbl + j - 1];
 180         w[j] += e[iwork[iinspl + itmp - 1]];
 181     }
 182
 183     if (fabs(scale-1.0)>GMX_DOUBLE_EPS) {
 184         d__1 = 1. / scale;
 185         F77_FUNC(dscal,DSCAL)(m, &d__1, &w[1], &c__1);
 186     }
 187     if (nsplit > 1) {
 188         i__1 = *m - 1;
 189         for (j = 1; j <= i__1; ++j) {
 190             i__ = 0;
 191             tmp = w[j];
 192             i__2 = *m;
 193             for (jj = j + 1; jj <= i__2; ++jj) {
 194                 if (w[jj] < tmp) {
 195                     i__ = jj;
 196                     tmp = w[jj];
 197                 }
 198             }
 199             if (i__ != 0) {
 200                 w[i__] = w[j];
 201                 w[j] = tmp;
 202                 if (wantz) {
 203                     F77_FUNC(dswap,DSWAP)(n, &z__[i__ * z_dim1 + 1], &c__1, &z__[j * z_dim1
 204                             + 1], &c__1);
 205                     itmp = isuppz[(i__ << 1) - 1];
 206                     isuppz[(i__ << 1) - 1] = isuppz[(j << 1) - 1];
 207                     isuppz[(j << 1) - 1] = itmp;
 208                     itmp = isuppz[i__ * 2];
 209                     isuppz[i__ * 2] = isuppz[j * 2];
 210                     isuppz[j * 2] = itmp;
 211                 }
 212             }
 213         }
 214     }
 215
 216     work[1] = (double) lwmin;
 217     iwork[1] = liwmin;
 218     return;
 219
 220 }