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,2017, 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           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 != nullptr) ? "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 #if 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  = static_cast<int>(w0);
 114     liwork = iw0;
 115
 116     snew(work, lwork);
 117     snew(iwork, liwork);
 118
 119     abstol = 0;
 120
 121 #if 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 #if 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         fflush(stderr);
 226     }
 227     while (info == 0 && (ido == -1 || ido == 1));
 228
 229     fprintf(stderr, "\n");
 230     if (info == 1)
 231     {
 232         gmx_fatal(FARGS,
 233                   "Maximum number of iterations (%d) reached in Arnoldi\n"
 234                   "diagonalization, but only %d of %d eigenvectors converged.\n",
 235                   maxiter, iparam[4], neig);
 236     }
 237     else if (info != 0)
 238     {
 239         gmx_fatal(FARGS, "Unspecified error from Arnoldi diagonalization:%d\n", info);
 240     }
 241
 242     info = 0;
 243     /* Extract eigenvalues and vectors from data */
 244     fprintf(stderr, "Calculating eigenvalues and eigenvectors...\n");
 245
 246 #if GMX_DOUBLE
 247     F77_FUNC(pdseupd, PDSEUPD) (&dovec, "A", select, eigenvalues, eigenvectors,
 248                                 &n, NULL, "I", &n, "SA", &neig, &abstol,
 249                                 resid, &ncv, v, &n, iparam, ipntr,
 250                                 workd, workl, &lworkl, &info);
 251 #else
 252     F77_FUNC(psseupd, PSSEUPD) (&dovec, "A", select, eigenvalues, eigenvectors,
 253                                 &n, NULL, "I", &n, "SA", &neig, &abstol,
 254                                 resid, &ncv, v, &n, iparam, ipntr,
 255                                 workd, workl, &lworkl, &info);
 256 #endif
 257
 258     sfree(v);
 259     sfree(resid);
 260     sfree(workd);
 261     sfree(workl);
 262     sfree(select);
 263 }
 264 #endif
 265
 266
 267 void
 268 sparse_eigensolver(gmx_sparsematrix_t *    A,
 269                    int                     neig,
 270                    real *                  eigenvalues,
 271                    real *                  eigenvectors,
 272                    int                     maxiter)
 273 {
 274     int      iwork[80];
 275     int      iparam[11];
 276     int      ipntr[11];
 277     real *   resid;
 278     real *   workd;
 279     real *   workl;
 280     real *   v;
 281     int      n;
 282     int      ido, info, lworkl, i, ncv, dovec;
 283     real     abstol;
 284     int *    select;
 285     int      iter;
 286
 287 #ifdef GMX_MPI_NOT
 288     MPI_Comm_size( MPI_COMM_WORLD, &n );
 289     if (n > 1)
 290     {
 291         sparse_parallel_eigensolver(A, neig, eigenvalues, eigenvectors, maxiter);
 292         return;
 293     }
 294 #endif
 295
 296     if (eigenvectors != nullptr)
 297     {
 298         dovec = 1;
 299     }
 300     else
 301     {
 302         dovec = 0;
 303     }
 304
 305     n   = A->nrow;
 306     ncv = 2*neig;
 307
 308     if (ncv > n)
 309     {
 310         ncv = n;
 311     }
 312
 313     for (i = 0; i < 11; i++)
 314     {
 315         iparam[i] = ipntr[i] = 0;
 316     }
 317
 318     iparam[0] = 1;       /* Don't use explicit shifts */
 319     iparam[2] = maxiter; /* Max number of iterations */
 320     iparam[6] = 1;       /* Standard symmetric eigenproblem */
 321
 322     lworkl = ncv*(8+ncv);
 323     snew(resid, n);
 324     snew(workd, (3*n+4));
 325     snew(workl, lworkl);
 326     snew(select, ncv);
 327     snew(v, n*ncv);
 328
 329     /* Use machine tolerance - roughly 1e-16 in double precision */
 330     abstol = 0;
 331
 332     ido = info = 0;
 333     fprintf(stderr, "Calculation Ritz values and Lanczos vectors, max %d iterations...\n", maxiter);
 334
 335     iter = 1;
 336     do
 337     {
 338 #if GMX_DOUBLE
 339         F77_FUNC(dsaupd, DSAUPD) (&ido, "I", &n, "SA", &neig, &abstol,
 340                                   resid, &ncv, v, &n, iparam, ipntr,
 341                                   workd, iwork, workl, &lworkl, &info);
 342 #else
 343         F77_FUNC(ssaupd, SSAUPD) (&ido, "I", &n, "SA", &neig, &abstol,
 344                                   resid, &ncv, v, &n, iparam, ipntr,
 345                                   workd, iwork, workl, &lworkl, &info);
 346 #endif
 347         if (ido == -1 || ido == 1)
 348         {
 349             gmx_sparsematrix_vector_multiply(A, workd+ipntr[0]-1, workd+ipntr[1]-1);
 350         }
 351
 352         fprintf(stderr, "\rIteration %4d: %3d out of %3d Ritz values converged.", iter++, iparam[4], neig);
 353         fflush(stderr);
 354     }
 355     while (info == 0 && (ido == -1 || ido == 1));
 356
 357     fprintf(stderr, "\n");
 358     if (info == 1)
 359     {
 360         gmx_fatal(FARGS,
 361                   "Maximum number of iterations (%d) reached in Arnoldi\n"
 362                   "diagonalization, but only %d of %d eigenvectors converged.\n",
 363                   maxiter, iparam[4], neig);
 364     }
 365     else if (info != 0)
 366     {
 367         gmx_fatal(FARGS, "Unspecified error from Arnoldi diagonalization:%d\n", info);
 368     }
 369
 370     info = 0;
 371     /* Extract eigenvalues and vectors from data */
 372     fprintf(stderr, "Calculating eigenvalues and eigenvectors...\n");
 373
 374 #if GMX_DOUBLE
 375     F77_FUNC(dseupd, DSEUPD) (&dovec, "A", select, eigenvalues, eigenvectors,
 376                               &n, nullptr, "I", &n, "SA", &neig, &abstol,
 377                               resid, &ncv, v, &n, iparam, ipntr,
 378                               workd, workl, &lworkl, &info);
 379 #else
 380     F77_FUNC(sseupd, SSEUPD) (&dovec, "A", select, eigenvalues, eigenvectors,
 381                               &n, nullptr, "I", &n, "SA", &neig, &abstol,
 382                               resid, &ncv, v, &n, iparam, ipntr,
 383                               workd, workl, &lworkl, &info);
 384 #endif
 385
 386     sfree(v);
 387     sfree(resid);
 388     sfree(workd);
 389     sfree(workl);
 390     sfree(select);
 391 }