src/gromacs/linearalgebra/gmx_lapack/dsyevr.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 void
  10 F77_FUNC(dsyevr,DSYEVR)(const char *jobz, const char *range, const char *uplo, int *n,
  11         double *a, int *lda, double *vl, double *vu, int *
  12         il, int *iu, double *abstol, int *m, double *w,
  13         double *z__, int *ldz, int *isuppz, double *work,
  14         int *lwork, int *iwork, int *liwork, int *info)
  15 {
  16     /* System generated locals */
  17     int a_dim1, a_offset, z_dim1, z_offset, i__1, i__2;
  18     double d__1, d__2;
  19
  20     /* Local variables */
  21     int c__1 = 1;
  22     int i__, j, nb, jj;
  23     double eps, vll, vuu, tmp1;
  24     int indd, inde;
  25     double anrm;
  26     int imax;
  27     double rmin, rmax;
  28     int itmp1, inddd, indee;
  29     double sigma;
  30     int iinfo;
  31     int indwk;
  32     int lwmin;
  33     int lower, wantz;
  34     int alleig, indeig;
  35     int iscale, ieeeok, indibl, indifl;
  36     int valeig;
  37     double safmin,minval;
  38     double abstll, bignum;
  39     int indtau;
  40     int indwkn;
  41     int liwmin;
  42     int llwrkn, llwork;
  43     double smlnum;
  44     int lwkopt;
  45     int lquery;
  46     double posinf,neginf,negzro,newzro;
  47     double fzero = 0.0;
  48     double fone = 1.0;
  49
  50
  51     /* Parameter adjustments */
  52     a_dim1 = *lda;
  53     a_offset = 1 + a_dim1;
  54     a -= a_offset;
  55     --w;
  56     z_dim1 = *ldz;
  57     z_offset = 1 + z_dim1;
  58     z__ -= z_offset;
  59     --isuppz;
  60     --work;
  61     --iwork;
  62
  63     /* Check for IEEE-compliant FP */
  64     ieeeok = 1;
  65     posinf = fone/fzero;
  66     if(posinf<=1.0)
  67       ieeeok = 0;
  68     neginf = -fone/fzero;
  69     if(neginf>=0.0)
  70       ieeeok = 0;
  71     negzro = fone/(neginf+fone);
  72     if(negzro!=0)
  73       ieeeok = 0;
  74     neginf = fone/negzro;
  75     if(neginf>=0)
  76       ieeeok = 0;
  77     newzro = negzro + fzero;
  78     if(newzro!=fzero)
  79       ieeeok = 0;
  80     posinf = fone /newzro;
  81     if(posinf<=fone)
  82       ieeeok = 0;
  83     neginf = neginf*posinf;
  84     if(neginf>=fzero)
  85       ieeeok = 0;
  86     posinf = posinf*posinf;
  87     if(posinf<=1.0)
  88       ieeeok = 0;
  89
  90     lower = (*uplo=='L' || *uplo=='l');
  91     wantz = (*jobz=='V' || *jobz=='v');
  92     alleig = (*range=='A' || *range=='a');
  93     valeig = (*range=='V' || *range=='v');
  94     indeig = (*range=='I' || *range=='i');
  95
  96     indibl = 0;
  97     lquery = *lwork == -1 || *liwork == -1;
  98
  99     i__1 = 1;
 100     i__2 = *n * 26;
 101
 102     if(*n>0)
 103       lwmin = *n * 26;
 104     else
 105       lwmin = 1;
 106
 107     if(*n>0)
 108       liwmin = *n * 10;
 109     else
 110       liwmin = 1;
 111
 112     *info = 0;
 113     if (! (wantz || (*jobz=='N' || *jobz=='n'))) {
 114         *info = -1;
 115     } else if (! (alleig || valeig || indeig)) {
 116         *info = -2;
 117     } else if (! (lower || (*uplo=='U' || *uplo=='u'))) {
 118         *info = -3;
 119     } else if (*n < 0) {
 120         *info = -4;
 121     } else if (*lda < ((*n>1) ? *n : 1) ) {
 122         *info = -6;
 123     } else {
 124         if (valeig) {
 125             if (*n > 0 && *vu <= *vl) {
 126                 *info = -8;
 127             }
 128         } else if (indeig) {
 129           if (*il < 1 || *il > ((*n>1) ? *n : 1)) {
 130                 *info = -9;
 131             } else if (*iu < ((*n<*il) ? *n : *il) || *iu > *n) {
 132                 *info = -10;
 133             }
 134         }
 135     }
 136     if (*info == 0) {
 137       if (*ldz < 1 || (wantz && *ldz < *n)) {
 138             *info = -15;
 139         } else if (*lwork < lwmin && ! lquery) {
 140             *info = -18;
 141         } else if (*liwork < liwmin && ! lquery) {
 142             *info = -20;
 143         }
 144     }
 145
 146     if (*info == 0) {
 147       nb = 32;
 148       /* Computing MAX */
 149       i__1 = (nb + 1) * *n;
 150       lwkopt = (i__1>lwmin) ? i__1 : lwmin;
 151       work[1] = (double) lwkopt;
 152       iwork[1] = liwmin;
 153     } else
 154       return;
 155
 156     if (lquery)
 157         return;
 158
 159     *m = 0;
 160     if (*n == 0) {
 161         work[1] = 1.;
 162         return;
 163     }
 164
 165     if (*n == 1) {
 166         work[1] = 7.;
 167         if (alleig || indeig) {
 168             *m = 1;
 169             w[1] = a[a_dim1 + 1];
 170         } else {
 171             if (*vl < a[a_dim1 + 1] && *vu >= a[a_dim1 + 1]) {
 172                 *m = 1;
 173                 w[1] = a[a_dim1 + 1];
 174             }
 175         }
 176         if (wantz) {
 177             z__[z_dim1 + 1] = 1.;
 178         }
 179         return;
 180     }
 181     minval = GMX_DOUBLE_MIN;
 182     safmin = minval*(1.0+GMX_DOUBLE_EPS);
 183     eps = GMX_DOUBLE_EPS;
 184
 185     smlnum = safmin / eps;
 186     bignum = 1. / smlnum;
 187     rmin = sqrt(smlnum);
 188
 189     d__1 = sqrt(bignum), d__2 = 1. / sqrt(sqrt(safmin));
 190     rmax = (d__1<d__2) ? d__1 : d__2;
 191
 192     iscale = 0;
 193     abstll = *abstol;
 194     vll = *vl;
 195     vuu = *vu;
 196     anrm = F77_FUNC(dlansy,DLANSY)("M", uplo, n, &a[a_offset], lda, &work[1]);
 197     if (anrm > 0. && anrm < rmin) {
 198         iscale = 1;
 199         sigma = rmin / anrm;
 200     } else if (anrm > rmax) {
 201         iscale = 1;
 202         sigma = rmax / anrm;
 203     }
 204     if (iscale == 1) {
 205         if (lower) {
 206             i__1 = *n;
 207             for (j = 1; j <= i__1; ++j) {
 208                 i__2 = *n - j + 1;
 209                 F77_FUNC(dscal,DSCAL)(&i__2, &sigma, &a[j + j * a_dim1], &c__1);
 210             }
 211         } else {
 212             i__1 = *n;
 213             for (j = 1; j <= i__1; ++j) {
 214                 F77_FUNC(dscal,DSCAL)(&j, &sigma, &a[j * a_dim1 + 1], &c__1);
 215
 216             }
 217         }
 218         if (*abstol > 0.) {
 219             abstll = *abstol * sigma;
 220         }
 221         if (valeig) {
 222             vll = *vl * sigma;
 223             vuu = *vu * sigma;
 224         }
 225     }
 226
 227     indtau = 1;
 228     inde = indtau + *n;
 229     indd = inde + *n;
 230     indee = indd + *n;
 231     inddd = indee + *n;
 232     indifl = inddd + *n;
 233     indwk = indifl + *n;
 234     llwork = *lwork - indwk + 1;
 235     F77_FUNC(dsytrd,DSYTRD)(uplo, n, &a[a_offset], lda, &work[indd], &work[inde], &work[
 236             indtau], &work[indwk], &llwork, &iinfo);
 237
 238     i__1 = *n - 1;
 239     F77_FUNC(dcopy,DCOPY)(&i__1, &work[inde], &c__1, &work[indee], &c__1);
 240     F77_FUNC(dcopy,DCOPY)(n, &work[indd], &c__1, &work[inddd], &c__1);
 241
 242     F77_FUNC(dstegr,DSTEGR)(jobz, range, n, &work[inddd], &work[indee], vl, vu, il, iu,
 243             abstol, m, &w[1], &z__[z_offset], ldz, &isuppz[1],
 244             &work[indwk], lwork, &iwork[1], liwork, info);
 245     if (wantz && *info == 0) {
 246       indwkn = inde;
 247       llwrkn = *lwork - indwkn + 1;
 248       F77_FUNC(dormtr,DORMTR)("L", uplo, "N", n, m, &a[a_offset], lda, &work[indtau]
 249               , &z__[z_offset], ldz, &work[indwkn], &llwrkn, &iinfo);
 250     }
 251
 252     if (*info != 0)
 253       return;
 254
 255     if (iscale == 1) {
 256         if (*info == 0) {
 257             imax = *m;
 258         } else {
 259             imax = *info - 1;
 260         }
 261         d__1 = 1. / sigma;
 262         F77_FUNC(dscal,DSCAL)(&imax, &d__1, &w[1], &c__1);
 263     }
 264
 265     if (wantz) {
 266         i__1 = *m - 1;
 267
 268         for (j = 1; j <= i__1; ++j) {
 269             i__ = 0;
 270             tmp1 = w[j];
 271             i__2 = *m;
 272             for (jj = j + 1; jj <= i__2; ++jj) {
 273                 if (w[jj] < tmp1) {
 274                     i__ = jj;
 275                     tmp1 = w[jj];
 276                 }
 277             }
 278
 279             if (i__ != 0) {
 280                 itmp1 = iwork[indibl + i__ - 1];
 281                 w[i__] = w[j];
 282                 iwork[indibl + i__ - 1] = iwork[indibl + j - 1];
 283                 w[j] = tmp1;
 284                 iwork[indibl + j - 1] = itmp1;
 285                 F77_FUNC(dswap,DSWAP)(n, &z__[i__ * z_dim1 + 1], &c__1, &z__[j * z_dim1 + 1],
 286                          &c__1);
 287             }
 288         }
 289     }
 290
 291     work[1] = (double) lwkopt;
 292     iwork[1] = liwmin;
 293     return;
 294
 295 }