src/gmxlib/gmx_lapack/dsyevr.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012, by the GROMACS development team, led by
   5  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   6  * others, as listed in the AUTHORS file in the top-level source
   7  * directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 #include <math.h>
  36
  37 #include "gmx_blas.h"
  38 #include "gmx_lapack.h"
  39 #include "lapack_limits.h"
  40
  41 void
  42 F77_FUNC(dsyevr,DSYEVR)(const char *jobz, const char *range, const char *uplo, int *n,
  43         double *a, int *lda, double *vl, double *vu, int *
  44         il, int *iu, double *abstol, int *m, double *w,
  45         double *z__, int *ldz, int *isuppz, double *work,
  46         int *lwork, int *iwork, int *liwork, int *info)
  47 {
  48     /* System generated locals */
  49     int a_dim1, a_offset, z_dim1, z_offset, i__1, i__2;
  50     double d__1, d__2;
  51
  52     /* Local variables */
  53     int c__1 = 1;
  54     int i__, j, nb, jj;
  55     double eps, vll, vuu, tmp1;
  56     int indd, inde;
  57     double anrm;
  58     int imax;
  59     double rmin, rmax;
  60     int itmp1, inddd, indee;
  61     double sigma;
  62     int iinfo;
  63     int indwk;
  64     int lwmin;
  65     int lower, wantz;
  66     int alleig, indeig;
  67     int iscale, ieeeok, indibl, indifl;
  68     int valeig;
  69     double safmin,minval;
  70     double abstll, bignum;
  71     int indtau;
  72     int indwkn;
  73     int liwmin;
  74     int llwrkn, llwork;
  75     double smlnum;
  76     int lwkopt;
  77     int lquery;
  78     double posinf,neginf,negzro,newzro;
  79     double fzero = 0.0;
  80     double fone = 1.0;
  81
  82
  83     /* Parameter adjustments */
  84     a_dim1 = *lda;
  85     a_offset = 1 + a_dim1;
  86     a -= a_offset;
  87     --w;
  88     z_dim1 = *ldz;
  89     z_offset = 1 + z_dim1;
  90     z__ -= z_offset;
  91     --isuppz;
  92     --work;
  93     --iwork;
  94
  95     /* Check for IEEE-compliant FP */
  96     ieeeok = 1;
  97     posinf = fone/fzero;
  98     if(posinf<=1.0)
  99       ieeeok = 0;
 100     neginf = -fone/fzero;
 101     if(neginf>=0.0)
 102       ieeeok = 0;
 103     negzro = fone/(neginf+fone);
 104     if(negzro!=0)
 105       ieeeok = 0;
 106     neginf = fone/negzro;
 107     if(neginf>=0)
 108       ieeeok = 0;
 109     newzro = negzro + fzero;
 110     if(newzro!=fzero)
 111       ieeeok = 0;
 112     posinf = fone /newzro;
 113     if(posinf<=fone)
 114       ieeeok = 0;
 115     neginf = neginf*posinf;
 116     if(neginf>=fzero)
 117       ieeeok = 0;
 118     posinf = posinf*posinf;
 119     if(posinf<=1.0)
 120       ieeeok = 0;
 121
 122     lower = (*uplo=='L' || *uplo=='l');
 123     wantz = (*jobz=='V' || *jobz=='v');
 124     alleig = (*range=='A' || *range=='a');
 125     valeig = (*range=='V' || *range=='v');
 126     indeig = (*range=='I' || *range=='i');
 127
 128     indibl = 0;
 129     lquery = *lwork == -1 || *liwork == -1;
 130
 131     i__1 = 1;
 132     i__2 = *n * 26;
 133
 134     if(*n>0)
 135       lwmin = *n * 26;
 136     else
 137       lwmin = 1;
 138
 139     if(*n>0)
 140       liwmin = *n * 10;
 141     else
 142       liwmin = 1;
 143
 144     *info = 0;
 145     if (! (wantz || (*jobz=='N' || *jobz=='n'))) {
 146         *info = -1;
 147     } else if (! (alleig || valeig || indeig)) {
 148         *info = -2;
 149     } else if (! (lower || (*uplo=='U' || *uplo=='u'))) {
 150         *info = -3;
 151     } else if (*n < 0) {
 152         *info = -4;
 153     } else if (*lda < ((*n>1) ? *n : 1) ) {
 154         *info = -6;
 155     } else {
 156         if (valeig) {
 157             if (*n > 0 && *vu <= *vl) {
 158                 *info = -8;
 159             }
 160         } else if (indeig) {
 161           if (*il < 1 || *il > ((*n>1) ? *n : 1)) {
 162                 *info = -9;
 163             } else if (*iu < ((*n<*il) ? *n : *il) || *iu > *n) {
 164                 *info = -10;
 165             }
 166         }
 167     }
 168     if (*info == 0) {
 169       if (*ldz < 1 || (wantz && *ldz < *n)) {
 170             *info = -15;
 171         } else if (*lwork < lwmin && ! lquery) {
 172             *info = -18;
 173         } else if (*liwork < liwmin && ! lquery) {
 174             *info = -20;
 175         }
 176     }
 177
 178     if (*info == 0) {
 179       nb = 32;
 180       /* Computing MAX */
 181       i__1 = (nb + 1) * *n;
 182       lwkopt = (i__1>lwmin) ? i__1 : lwmin;
 183       work[1] = (double) lwkopt;
 184       iwork[1] = liwmin;
 185     } else
 186       return;
 187
 188     if (lquery)
 189         return;
 190
 191     *m = 0;
 192     if (*n == 0) {
 193         work[1] = 1.;
 194         return;
 195     }
 196
 197     if (*n == 1) {
 198         work[1] = 7.;
 199         if (alleig || indeig) {
 200             *m = 1;
 201             w[1] = a[a_dim1 + 1];
 202         } else {
 203             if (*vl < a[a_dim1 + 1] && *vu >= a[a_dim1 + 1]) {
 204                 *m = 1;
 205                 w[1] = a[a_dim1 + 1];
 206             }
 207         }
 208         if (wantz) {
 209             z__[z_dim1 + 1] = 1.;
 210         }
 211         return;
 212     }
 213     minval = GMX_DOUBLE_MIN;
 214     safmin = minval*(1.0+GMX_DOUBLE_EPS);
 215     eps = GMX_DOUBLE_EPS;
 216
 217     smlnum = safmin / eps;
 218     bignum = 1. / smlnum;
 219     rmin = sqrt(smlnum);
 220
 221     d__1 = sqrt(bignum), d__2 = 1. / sqrt(sqrt(safmin));
 222     rmax = (d__1<d__2) ? d__1 : d__2;
 223
 224     iscale = 0;
 225     abstll = *abstol;
 226     vll = *vl;
 227     vuu = *vu;
 228     anrm = F77_FUNC(dlansy,DLANSY)("M", uplo, n, &a[a_offset], lda, &work[1]);
 229     if (anrm > 0. && anrm < rmin) {
 230         iscale = 1;
 231         sigma = rmin / anrm;
 232     } else if (anrm > rmax) {
 233         iscale = 1;
 234         sigma = rmax / anrm;
 235     }
 236     if (iscale == 1) {
 237         if (lower) {
 238             i__1 = *n;
 239             for (j = 1; j <= i__1; ++j) {
 240                 i__2 = *n - j + 1;
 241                 F77_FUNC(dscal,DSCAL)(&i__2, &sigma, &a[j + j * a_dim1], &c__1);
 242             }
 243         } else {
 244             i__1 = *n;
 245             for (j = 1; j <= i__1; ++j) {
 246                 F77_FUNC(dscal,DSCAL)(&j, &sigma, &a[j * a_dim1 + 1], &c__1);
 247
 248             }
 249         }
 250         if (*abstol > 0.) {
 251             abstll = *abstol * sigma;
 252         }
 253         if (valeig) {
 254             vll = *vl * sigma;
 255             vuu = *vu * sigma;
 256         }
 257     }
 258
 259     indtau = 1;
 260     inde = indtau + *n;
 261     indd = inde + *n;
 262     indee = indd + *n;
 263     inddd = indee + *n;
 264     indifl = inddd + *n;
 265     indwk = indifl + *n;
 266     llwork = *lwork - indwk + 1;
 267     F77_FUNC(dsytrd,DSYTRD)(uplo, n, &a[a_offset], lda, &work[indd], &work[inde], &work[
 268             indtau], &work[indwk], &llwork, &iinfo);
 269
 270     i__1 = *n - 1;
 271     F77_FUNC(dcopy,DCOPY)(&i__1, &work[inde], &c__1, &work[indee], &c__1);
 272     F77_FUNC(dcopy,DCOPY)(n, &work[indd], &c__1, &work[inddd], &c__1);
 273
 274     F77_FUNC(dstegr,DSTEGR)(jobz, range, n, &work[inddd], &work[indee], vl, vu, il, iu,
 275             abstol, m, &w[1], &z__[z_offset], ldz, &isuppz[1],
 276             &work[indwk], lwork, &iwork[1], liwork, info);
 277     if (wantz && *info == 0) {
 278       indwkn = inde;
 279       llwrkn = *lwork - indwkn + 1;
 280       F77_FUNC(dormtr,DORMTR)("L", uplo, "N", n, m, &a[a_offset], lda, &work[indtau]
 281               , &z__[z_offset], ldz, &work[indwkn], &llwrkn, &iinfo);
 282     }
 283
 284     if (*info != 0)
 285       return;
 286
 287     if (iscale == 1) {
 288         if (*info == 0) {
 289             imax = *m;
 290         } else {
 291             imax = *info - 1;
 292         }
 293         d__1 = 1. / sigma;
 294         F77_FUNC(dscal,DSCAL)(&imax, &d__1, &w[1], &c__1);
 295     }
 296
 297     if (wantz) {
 298         i__1 = *m - 1;
 299
 300         for (j = 1; j <= i__1; ++j) {
 301             i__ = 0;
 302             tmp1 = w[j];
 303             i__2 = *m;
 304             for (jj = j + 1; jj <= i__2; ++jj) {
 305                 if (w[jj] < tmp1) {
 306                     i__ = jj;
 307                     tmp1 = w[jj];
 308                 }
 309             }
 310
 311             if (i__ != 0) {
 312                 itmp1 = iwork[indibl + i__ - 1];
 313                 w[i__] = w[j];
 314                 iwork[indibl + i__ - 1] = iwork[indibl + j - 1];
 315                 w[j] = tmp1;
 316                 iwork[indibl + j - 1] = itmp1;
 317                 F77_FUNC(dswap,DSWAP)(n, &z__[i__ * z_dim1 + 1], &c__1, &z__[j * z_dim1 + 1],
 318                          &c__1);
 319             }
 320         }
 321     }
 322
 323     work[1] = (double) lwkopt;
 324     iwork[1] = liwmin;
 325     return;
 326
 327 }