src/gromacs/linearalgebra/eigensolver.cpp

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