src/gromacs/linearalgebra/eigensolver.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2012,2013,2014, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include "eigensolver.h"
  40
  41 #include "gromacs/linearalgebra/sparsematrix.h"
  42 #include "gromacs/utility/fatalerror.h"
  43 #include "gromacs/utility/real.h"
  44 #include "gromacs/utility/smalloc.h"
  45
  46 #include "gmx_arpack.h"
  47 #include "gmx_lapack.h"
  48
  49 void
  50 eigensolver(real *   a,
  51             int      n,
  52             int      index_lower,
  53             int      index_upper,
  54             real *   eigenvalues,
  55             real *   eigenvectors)
  56 {
  57     int       *   isuppz;
  58     int           lwork, liwork;
  59     int           il, iu, m, iw0, info;
  60     real          w0, abstol;
  61     int       *   iwork;
  62     real       *  work;
  63     real          vl, vu;
  64     const char *  jobz;
  65
  66     if (index_lower < 0)
  67     {
  68         index_lower = 0;
  69     }
  70
  71     if (index_upper >= n)
  72     {
  73         index_upper = n-1;
  74     }
  75
  76     /* Make jobz point to the character "V" if eigenvectors
  77      * should be calculated, otherwise "N" (only eigenvalues).
  78      */
  79     jobz = (eigenvectors != NULL) ? "V" : "N";
  80
  81     /* allocate lapack stuff */
  82     snew(isuppz, 2*n);
  83     vl = vu = 0;
  84
  85     /* First time we ask the routine how much workspace it needs */
  86     lwork  = -1;
  87     liwork = -1;
  88     abstol = 0;
  89
  90     /* Convert indices to fortran standard */
  91     index_lower++;
  92     index_upper++;
  93
  94     /* Call LAPACK routine using fortran interface. Note that we use upper storage,
  95      * but this corresponds to lower storage ("L") in Fortran.
  96      */
  97 #ifdef GMX_DOUBLE
  98     F77_FUNC(dsyevr, DSYEVR) (jobz, "I", "L", &n, a, &n, &vl, &vu, &index_lower, &index_upper,
  99                               &abstol, &m, eigenvalues, eigenvectors, &n,
 100                               isuppz, &w0, &lwork, &iw0, &liwork, &info);
 101 #else
 102     F77_FUNC(ssyevr, SSYEVR) (jobz, "I", "L", &n, a, &n, &vl, &vu, &index_lower, &index_upper,
 103                               &abstol, &m, eigenvalues, eigenvectors, &n,
 104                               isuppz, &w0, &lwork, &iw0, &liwork, &info);
 105 #endif
 106
 107     if (info != 0)
 108     {
 109         sfree(isuppz);
 110         gmx_fatal(FARGS, "Internal errror in LAPACK diagonalization.");
 111     }
 112
 113     lwork  = w0;
 114     liwork = iw0;
 115
 116     snew(work, lwork);
 117     snew(iwork, liwork);
 118
 119     abstol = 0;
 120
 121 #ifdef GMX_DOUBLE
 122     F77_FUNC(dsyevr, DSYEVR) (jobz, "I", "L", &n, a, &n, &vl, &vu, &index_lower, &index_upper,
 123                               &abstol, &m, eigenvalues, eigenvectors, &n,
 124                               isuppz, work, &lwork, iwork, &liwork, &info);
 125 #else
 126     F77_FUNC(ssyevr, SSYEVR) (jobz, "I", "L", &n, a, &n, &vl, &vu, &index_lower, &index_upper,
 127                               &abstol, &m, eigenvalues, eigenvectors, &n,
 128                               isuppz, work, &lwork, iwork, &liwork, &info);
 129 #endif
 130
 131     sfree(isuppz);
 132     sfree(work);
 133     sfree(iwork);
 134
 135     if (info != 0)
 136     {
 137         gmx_fatal(FARGS, "Internal errror in LAPACK diagonalization.");
 138     }
 139
 140 }
 141
 142
 143 #ifdef GMX_MPI_NOT
 144 void
 145 sparse_parallel_eigensolver(gmx_sparsematrix_t *    A,
 146                             int                     neig,
 147                             real *                  eigenvalues,
 148                             real *                  eigenvectors,
 149                             int                     maxiter)
 150 {
 151     int      iwork[80];
 152     int      iparam[11];
 153     int      ipntr[11];
 154     real *   resid;
 155     real *   workd;
 156     real *   workl;
 157     real *   v;
 158     int      n;
 159     int      ido, info, lworkl, i, ncv, dovec;
 160     real     abstol;
 161     int *    select;
 162     int      iter;
 163     int      nnodes, rank;
 164
 165     MPI_Comm_size( MPI_COMM_WORLD, &nnodes );
 166     MPI_Comm_rank( MPI_COMM_WORLD, &rank );
 167
 168     if (eigenvectors != NULL)
 169     {
 170         dovec = 1;
 171     }
 172     else
 173     {
 174         dovec = 0;
 175     }
 176
 177     n   = A->nrow;
 178     ncv = 2*neig;
 179
 180     if (ncv > n)
 181     {
 182         ncv = n;
 183     }
 184
 185     for (i = 0; i < 11; i++)
 186     {
 187         iparam[i] = ipntr[i] = 0;
 188     }
 189
 190     iparam[0] = 1;       /* Don't use explicit shifts */
 191     iparam[2] = maxiter; /* Max number of iterations */
 192     iparam[6] = 1;       /* Standard symmetric eigenproblem */
 193
 194     lworkl = ncv*(8+ncv);
 195     snew(resid, n);
 196     snew(workd, (3*n+4));
 197     snew(workl, lworkl);
 198     snew(select, ncv);
 199     snew(v, n*ncv);
 200
 201     /* Use machine tolerance - roughly 1e-16 in double precision */
 202     abstol = 0;
 203
 204     ido = info = 0;
 205     fprintf(stderr, "Calculation Ritz values and Lanczos vectors, max %d iterations...\n", maxiter);
 206
 207     iter = 1;
 208     do
 209     {
 210 #ifdef GMX_DOUBLE
 211         F77_FUNC(pdsaupd, PDSAUPD) (&ido, "I", &n, "SA", &neig, &abstol,
 212                                     resid, &ncv, v, &n, iparam, ipntr,
 213                                     workd, iwork, workl, &lworkl, &info);
 214 #else
 215         F77_FUNC(pssaupd, PSSAUPD) (&ido, "I", &n, "SA", &neig, &abstol,
 216                                     resid, &ncv, v, &n, iparam, ipntr,
 217                                     workd, iwork, workl, &lworkl, &info);
 218 #endif
 219         if (ido == -1 || ido == 1)
 220         {
 221             gmx_sparsematrix_vector_multiply(A, workd+ipntr[0]-1, workd+ipntr[1]-1);
 222         }
 223
 224         fprintf(stderr, "\rIteration %4d: %3d out of %3d Ritz values converged.", iter++, iparam[4], neig);
 225     }
 226     while (info == 0 && (ido == -1 || ido == 1));
 227
 228     fprintf(stderr, "\n");
 229     if (info == 1)
 230     {
 231         gmx_fatal(FARGS,
 232                   "Maximum number of iterations (%d) reached in Arnoldi\n"
 233                   "diagonalization, but only %d of %d eigenvectors converged.\n",
 234                   maxiter, iparam[4], neig);
 235     }
 236     else if (info != 0)
 237     {
 238         gmx_fatal(FARGS, "Unspecified error from Arnoldi diagonalization:%d\n", info);
 239     }
 240
 241     info = 0;
 242     /* Extract eigenvalues and vectors from data */
 243     fprintf(stderr, "Calculating eigenvalues and eigenvectors...\n");
 244
 245 #ifdef GMX_DOUBLE
 246     F77_FUNC(pdseupd, PDSEUPD) (&dovec, "A", select, eigenvalues, eigenvectors,
 247                                 &n, NULL, "I", &n, "SA", &neig, &abstol,
 248                                 resid, &ncv, v, &n, iparam, ipntr,
 249                                 workd, workl, &lworkl, &info);
 250 #else
 251     F77_FUNC(psseupd, PSSEUPD) (&dovec, "A", select, eigenvalues, eigenvectors,
 252                                 &n, NULL, "I", &n, "SA", &neig, &abstol,
 253                                 resid, &ncv, v, &n, iparam, ipntr,
 254                                 workd, workl, &lworkl, &info);
 255 #endif
 256
 257     sfree(v);
 258     sfree(resid);
 259     sfree(workd);
 260     sfree(workl);
 261     sfree(select);
 262 }
 263 #endif
 264
 265
 266 void
 267 sparse_eigensolver(gmx_sparsematrix_t *    A,
 268                    int                     neig,
 269                    real *                  eigenvalues,
 270                    real *                  eigenvectors,
 271                    int                     maxiter)
 272 {
 273     int      iwork[80];
 274     int      iparam[11];
 275     int      ipntr[11];
 276     real *   resid;
 277     real *   workd;
 278     real *   workl;
 279     real *   v;
 280     int      n;
 281     int      ido, info, lworkl, i, ncv, dovec;
 282     real     abstol;
 283     int *    select;
 284     int      iter;
 285
 286 #ifdef GMX_MPI_NOT
 287     MPI_Comm_size( MPI_COMM_WORLD, &n );
 288     if (n > 1)
 289     {
 290         sparse_parallel_eigensolver(A, neig, eigenvalues, eigenvectors, maxiter);
 291         return;
 292     }
 293 #endif
 294
 295     if (eigenvectors != NULL)
 296     {
 297         dovec = 1;
 298     }
 299     else
 300     {
 301         dovec = 0;
 302     }
 303
 304     n   = A->nrow;
 305     ncv = 2*neig;
 306
 307     if (ncv > n)
 308     {
 309         ncv = n;
 310     }
 311
 312     for (i = 0; i < 11; i++)
 313     {
 314         iparam[i] = ipntr[i] = 0;
 315     }
 316
 317     iparam[0] = 1;       /* Don't use explicit shifts */
 318     iparam[2] = maxiter; /* Max number of iterations */
 319     iparam[6] = 1;       /* Standard symmetric eigenproblem */
 320
 321     lworkl = ncv*(8+ncv);
 322     snew(resid, n);
 323     snew(workd, (3*n+4));
 324     snew(workl, lworkl);
 325     snew(select, ncv);
 326     snew(v, n*ncv);
 327
 328     /* Use machine tolerance - roughly 1e-16 in double precision */
 329     abstol = 0;
 330
 331     ido = info = 0;
 332     fprintf(stderr, "Calculation Ritz values and Lanczos vectors, max %d iterations...\n", maxiter);
 333
 334     iter = 1;
 335     do
 336     {
 337 #ifdef GMX_DOUBLE
 338         F77_FUNC(dsaupd, DSAUPD) (&ido, "I", &n, "SA", &neig, &abstol,
 339                                   resid, &ncv, v, &n, iparam, ipntr,
 340                                   workd, iwork, workl, &lworkl, &info);
 341 #else
 342         F77_FUNC(ssaupd, SSAUPD) (&ido, "I", &n, "SA", &neig, &abstol,
 343                                   resid, &ncv, v, &n, iparam, ipntr,
 344                                   workd, iwork, workl, &lworkl, &info);
 345 #endif
 346         if (ido == -1 || ido == 1)
 347         {
 348             gmx_sparsematrix_vector_multiply(A, workd+ipntr[0]-1, workd+ipntr[1]-1);
 349         }
 350
 351         fprintf(stderr, "\rIteration %4d: %3d out of %3d Ritz values converged.", iter++, iparam[4], neig);
 352     }
 353     while (info == 0 && (ido == -1 || ido == 1));
 354
 355     fprintf(stderr, "\n");
 356     if (info == 1)
 357     {
 358         gmx_fatal(FARGS,
 359                   "Maximum number of iterations (%d) reached in Arnoldi\n"
 360                   "diagonalization, but only %d of %d eigenvectors converged.\n",
 361                   maxiter, iparam[4], neig);
 362     }
 363     else if (info != 0)
 364     {
 365         gmx_fatal(FARGS, "Unspecified error from Arnoldi diagonalization:%d\n", info);
 366     }
 367
 368     info = 0;
 369     /* Extract eigenvalues and vectors from data */
 370     fprintf(stderr, "Calculating eigenvalues and eigenvectors...\n");
 371
 372 #ifdef GMX_DOUBLE
 373     F77_FUNC(dseupd, DSEUPD) (&dovec, "A", select, eigenvalues, eigenvectors,
 374                               &n, NULL, "I", &n, "SA", &neig, &abstol,
 375                               resid, &ncv, v, &n, iparam, ipntr,
 376                               workd, workl, &lworkl, &info);
 377 #else
 378     F77_FUNC(sseupd, SSEUPD) (&dovec, "A", select, eigenvalues, eigenvectors,
 379                               &n, NULL, "I", &n, "SA", &neig, &abstol,
 380                               resid, &ncv, v, &n, iparam, ipntr,
 381                               workd, workl, &lworkl, &info);
 382 #endif
 383
 384     sfree(v);
 385     sfree(resid);
 386     sfree(workd);
 387     sfree(workl);
 388     sfree(select);
 389 }