src/gmxlib/gmx_lapack/sstegr.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 #include "gmx_blas.h"
  37 #include "gmx_lapack.h"
  38 #include "lapack_limits.h"
  39
  40 #include <types/simple.h>
  41
  42 void
  43 F77_FUNC(sstegr,SSTEGR)(const char *jobz,
  44         const char *range,
  45         int *n,
  46         float *d__,
  47         float *e,
  48         float *vl,
  49         float *vu,
  50         int *il,
  51         int *iu,
  52         float *abstol,
  53         int *m,
  54         float *w,
  55         float *z__,
  56         int *ldz,
  57         int *isuppz,
  58         float *work,
  59         int *lwork,
  60         int *iwork,
  61         int *liwork,
  62         int *info)
  63 {
  64     int z_dim1, z_offset, i__1, i__2;
  65     float d__1, d__2;
  66     int c__1 = 1;
  67
  68     int i__, j;
  69     int jj;
  70     float eps, tol, tmp, rmin, rmax;
  71     int itmp;
  72     float tnrm;
  73     float scale;
  74     int iinfo, iindw;
  75     int lwmin;
  76     int wantz;
  77     int iindbl;
  78     int valeig,alleig,indeig;
  79     float safmin,minval;
  80     float bignum;
  81     int iindwk, indgrs;
  82     float thresh;
  83     int iinspl, indwrk, liwmin, nsplit;
  84     float smlnum;
  85     int lquery;
  86
  87
  88     --d__;
  89     --e;
  90     --w;
  91     z_dim1 = *ldz;
  92     z_offset = 1 + z_dim1;
  93     z__ -= z_offset;
  94     --isuppz;
  95     --work;
  96     --iwork;
  97
  98     wantz = (*jobz=='V' || *jobz=='v');
  99     alleig = (*range=='A' || *range=='a');
 100     valeig = (*range=='V' || *range=='v');
 101     indeig = (*range=='I' || *range=='i');
 102
 103     lquery = *lwork == -1 || *liwork == -1;
 104     lwmin = *n * 17;
 105     liwmin = *n * 10;
 106
 107     *info = 0;
 108     if (! (wantz || (*jobz=='N' || *jobz=='n'))) {
 109         *info = -1;
 110     } else if (! (alleig || valeig || indeig)) {
 111         *info = -2;
 112     } else if (*n < 0) {
 113         *info = -3;
 114     } else if (valeig && *n > 0 && *vu <= *vl) {
 115         *info = -7;
 116     } else if (indeig && (*il < 1 || *il > *n)) {
 117         *info = -8;
 118     } else if (indeig && (*iu < *il || *iu > *n)) {
 119         *info = -9;
 120     } else if (*ldz < 1 || (wantz && *ldz < *n)) {
 121         *info = -14;
 122     } else if (*lwork < lwmin && ! lquery) {
 123         *info = -17;
 124     } else if (*liwork < liwmin && ! lquery) {
 125         *info = -19;
 126     }
 127     if (*info == 0) {
 128         work[1] = (float) lwmin;
 129         iwork[1] = liwmin;
 130     }
 131
 132     if (*info != 0) {
 133         i__1 = -(*info);
 134         return;
 135     } else if (lquery) {
 136         return;
 137     }
 138
 139     *m = 0;
 140     if (*n == 0) {
 141         return;
 142     }
 143
 144     if (*n == 1) {
 145         if (alleig || indeig) {
 146             *m = 1;
 147             w[1] = d__[1];
 148         } else {
 149             if (*vl < d__[1] && *vu >= d__[1]) {
 150                 *m = 1;
 151                 w[1] = d__[1];
 152             }
 153         }
 154         if (wantz) {
 155             z__[z_dim1 + 1] = 1.;
 156         }
 157         return;
 158     }
 159
 160     minval = GMX_FLOAT_MIN;
 161     safmin = minval*(1.0+GMX_FLOAT_EPS);
 162     eps = GMX_FLOAT_EPS;
 163     smlnum = safmin / eps;
 164     bignum = 1. / smlnum;
 165     rmin = sqrt(smlnum);
 166     d__1 = sqrt(bignum), d__2 = 1. / sqrt(sqrt(safmin));
 167     rmax = (d__1<d__2) ? d__1 : d__2;
 168     scale = 1.;
 169     tnrm = F77_FUNC(slanst,SLANST)("M", n, &d__[1], &e[1]);
 170     if (tnrm > 0. && tnrm < rmin) {
 171         scale = rmin / tnrm;
 172     } else if (tnrm > rmax) {
 173         scale = rmax / tnrm;
 174     }
 175     if ( fabs(scale-1.0)>GMX_FLOAT_EPS) {
 176         F77_FUNC(sscal,SSCAL)(n, &scale, &d__[1], &c__1);
 177         i__1 = *n - 1;
 178         F77_FUNC(sscal,SSCAL)(&i__1, &scale, &e[1], &c__1);
 179         tnrm *= scale;
 180     }
 181     indgrs = 1;
 182     indwrk = (*n << 1) + 1;
 183
 184     iinspl = 1;
 185     iindbl = *n + 1;
 186     iindw = (*n << 1) + 1;
 187     iindwk = *n * 3 + 1;
 188
 189     thresh = eps * tnrm;
 190     F77_FUNC(slarrex,SLARREX)(range, n, vl, vu, il, iu, &d__[1], &e[1], &thresh, &nsplit, &
 191             iwork[iinspl], m, &w[1], &iwork[iindbl], &iwork[iindw], &work[
 192             indgrs], &work[indwrk], &iwork[iindwk], &iinfo);
 193
 194     if (iinfo != 0) {
 195         *info = 1;
 196         return;
 197     }
 198
 199     if (wantz) {
 200         d__1 = *abstol, d__2 = (float) (*n) * eps;
 201         tol = (d__1>d__2) ? d__1 : d__2;
 202         F77_FUNC(slarrvx,SLARRVX)(n, &d__[1], &e[1], &iwork[iinspl], m, &w[1], &iwork[iindbl], &
 203                 iwork[iindw], &work[indgrs], &tol, &z__[z_offset], ldz, &
 204                 isuppz[1], &work[indwrk], &iwork[iindwk], &iinfo);
 205         if (iinfo != 0) {
 206             *info = 2;
 207             return;
 208         }
 209     }
 210
 211     i__1 = *m;
 212     for (j = 1; j <= i__1; ++j) {
 213         itmp = iwork[iindbl + j - 1];
 214         w[j] += e[iwork[iinspl + itmp - 1]];
 215     }
 216
 217     if (fabs(scale-1.0)>GMX_FLOAT_EPS) {
 218         d__1 = 1. / scale;
 219         F77_FUNC(sscal,SSCAL)(m, &d__1, &w[1], &c__1);
 220     }
 221     if (nsplit > 1) {
 222         i__1 = *m - 1;
 223         for (j = 1; j <= i__1; ++j) {
 224             i__ = 0;
 225             tmp = w[j];
 226             i__2 = *m;
 227             for (jj = j + 1; jj <= i__2; ++jj) {
 228                 if (w[jj] < tmp) {
 229                     i__ = jj;
 230                     tmp = w[jj];
 231                 }
 232             }
 233             if (i__ != 0) {
 234                 w[i__] = w[j];
 235                 w[j] = tmp;
 236                 if (wantz) {
 237                     F77_FUNC(sswap,SSWAP)(n, &z__[i__ * z_dim1 + 1], &c__1, &z__[j * z_dim1
 238                             + 1], &c__1);
 239                     itmp = isuppz[(i__ << 1) - 1];
 240                     isuppz[(i__ << 1) - 1] = isuppz[(j << 1) - 1];
 241                     isuppz[(j << 1) - 1] = itmp;
 242                     itmp = isuppz[i__ * 2];
 243                     isuppz[i__ * 2] = isuppz[j * 2];
 244                     isuppz[j * 2] = itmp;
 245                 }
 246             }
 247         }
 248     }
 249
 250     work[1] = (float) lwmin;
 251     iwork[1] = liwmin;
 252     return;
 253
 254 }