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 "eigensolver.h"
  38
  39 #include "gromacs/legacyheaders/types/simple.h"
  40
  41 #include "gromacs/linearalgebra/sparsematrix.h"
  42 #include "gromacs/utility/fatalerror.h"
  43 #include "gromacs/utility/smalloc.h"
  44
  45 #include "gmx_lapack.h"
  46 #include "gmx_arpack.h"
  47
  48 void
  49 eigensolver(real *   a,
  50             int      n,
  51             int      index_lower,
  52             int      index_upper,
  53             real *   eigenvalues,
  54             real *   eigenvectors)
  55 {
  56     int       *   isuppz;
  57     int           lwork, liwork;
  58     int           il, iu, m, iw0, info;
  59     real          w0, abstol;
  60     int       *   iwork;
  61     real       *  work;
  62     real          vl, vu;
  63     const char *  jobz;
  64
  65     if (index_lower < 0)
  66     {
  67         index_lower = 0;
  68     }
  69
  70     if (index_upper >= n)
  71     {
  72         index_upper = n-1;
  73     }
  74
  75     /* Make jobz point to the character "V" if eigenvectors
  76      * should be calculated, otherwise "N" (only eigenvalues).
  77      */
  78     jobz = (eigenvectors != NULL) ? "V" : "N";
  79
  80     /* allocate lapack stuff */
  81     snew(isuppz, 2*n);
  82     vl = vu = 0;
  83
  84     /* First time we ask the routine how much workspace it needs */
  85     lwork  = -1;
  86     liwork = -1;
  87     abstol = 0;
  88
  89     /* Convert indices to fortran standard */
  90     index_lower++;
  91     index_upper++;
  92
  93     /* Call LAPACK routine using fortran interface. Note that we use upper storage,
  94      * but this corresponds to lower storage ("L") in Fortran.
  95      */
  96 #ifdef GMX_DOUBLE
  97     F77_FUNC(dsyevr, DSYEVR) (jobz, "I", "L", &n, a, &n, &vl, &vu, &index_lower, &index_upper,
  98                               &abstol, &m, eigenvalues, eigenvectors, &n,
  99                               isuppz, &w0, &lwork, &iw0, &liwork, &info);
 100 #else
 101     F77_FUNC(ssyevr, SSYEVR) (jobz, "I", "L", &n, a, &n, &vl, &vu, &index_lower, &index_upper,
 102                               &abstol, &m, eigenvalues, eigenvectors, &n,
 103                               isuppz, &w0, &lwork, &iw0, &liwork, &info);
 104 #endif
 105
 106     if (info != 0)
 107     {
 108         sfree(isuppz);
 109         gmx_fatal(FARGS, "Internal errror in LAPACK diagonalization.");
 110     }
 111
 112     lwork  = w0;
 113     liwork = iw0;
 114
 115     snew(work, lwork);
 116     snew(iwork, liwork);
 117
 118     abstol = 0;
 119
 120 #ifdef GMX_DOUBLE
 121     F77_FUNC(dsyevr, DSYEVR) (jobz, "I", "L", &n, a, &n, &vl, &vu, &index_lower, &index_upper,
 122                               &abstol, &m, eigenvalues, eigenvectors, &n,
 123                               isuppz, work, &lwork, iwork, &liwork, &info);
 124 #else
 125     F77_FUNC(ssyevr, SSYEVR) (jobz, "I", "L", &n, a, &n, &vl, &vu, &index_lower, &index_upper,
 126                               &abstol, &m, eigenvalues, eigenvectors, &n,
 127                               isuppz, work, &lwork, iwork, &liwork, &info);
 128 #endif
 129
 130     sfree(isuppz);
 131     sfree(work);
 132     sfree(iwork);
 133
 134     if (info != 0)
 135     {
 136         gmx_fatal(FARGS, "Internal errror in LAPACK diagonalization.");
 137     }
 138
 139 }
 140
 141
 142 #ifdef GMX_MPI_NOT
 143 void
 144 sparse_parallel_eigensolver(gmx_sparsematrix_t *    A,
 145                             int                     neig,
 146                             real *                  eigenvalues,
 147                             real *                  eigenvectors,
 148                             int                     maxiter)
 149 {
 150     int      iwork[80];
 151     int      iparam[11];
 152     int      ipntr[11];
 153     real *   resid;
 154     real *   workd;
 155     real *   workl;
 156     real *   v;
 157     int      n;
 158     int      ido, info, lworkl, i, ncv, dovec;
 159     real     abstol;
 160     int *    select;
 161     int      iter;
 162     int      nnodes, rank;
 163
 164     MPI_Comm_size( MPI_COMM_WORLD, &nnodes );
 165     MPI_Comm_rank( MPI_COMM_WORLD, &rank );
 166
 167     if (eigenvectors != NULL)
 168     {
 169         dovec = 1;
 170     }
 171     else
 172     {
 173         dovec = 0;
 174     }
 175
 176     n   = A->nrow;
 177     ncv = 2*neig;
 178
 179     if (ncv > n)
 180     {
 181         ncv = n;
 182     }
 183
 184     for (i = 0; i < 11; i++)
 185     {
 186         iparam[i] = ipntr[i] = 0;
 187     }
 188
 189     iparam[0] = 1;       /* Don't use explicit shifts */
 190     iparam[2] = maxiter; /* Max number of iterations */
 191     iparam[6] = 1;       /* Standard symmetric eigenproblem */
 192
 193     lworkl = ncv*(8+ncv);
 194     snew(resid, n);
 195     snew(workd, (3*n+4));
 196     snew(workl, lworkl);
 197     snew(select, ncv);
 198     snew(v, n*ncv);
 199
 200     /* Use machine tolerance - roughly 1e-16 in double precision */
 201     abstol = 0;
 202
 203     ido = info = 0;
 204     fprintf(stderr, "Calculation Ritz values and Lanczos vectors, max %d iterations...\n", maxiter);
 205
 206     iter = 1;
 207     do
 208     {
 209 #ifdef GMX_DOUBLE
 210         F77_FUNC(pdsaupd, PDSAUPD) (&ido, "I", &n, "SA", &neig, &abstol,
 211                                     resid, &ncv, v, &n, iparam, ipntr,
 212                                     workd, iwork, workl, &lworkl, &info);
 213 #else
 214         F77_FUNC(pssaupd, PSSAUPD) (&ido, "I", &n, "SA", &neig, &abstol,
 215                                     resid, &ncv, v, &n, iparam, ipntr,
 216                                     workd, iwork, workl, &lworkl, &info);
 217 #endif
 218         if (ido == -1 || ido == 1)
 219         {
 220             gmx_sparsematrix_vector_multiply(A, workd+ipntr[0]-1, workd+ipntr[1]-1);
 221         }
 222
 223         fprintf(stderr, "\rIteration %4d: %3d out of %3d Ritz values converged.", iter++, iparam[4], neig);
 224     }
 225     while (info == 0 && (ido == -1 || ido == 1));
 226
 227     fprintf(stderr, "\n");
 228     if (info == 1)
 229     {
 230         gmx_fatal(FARGS,
 231                   "Maximum number of iterations (%d) reached in Arnoldi\n"
 232                   "diagonalization, but only %d of %d eigenvectors converged.\n",
 233                   maxiter, iparam[4], neig);
 234     }
 235     else if (info != 0)
 236     {
 237         gmx_fatal(FARGS, "Unspecified error from Arnoldi diagonalization:%d\n", info);
 238     }
 239
 240     info = 0;
 241     /* Extract eigenvalues and vectors from data */
 242     fprintf(stderr, "Calculating eigenvalues and eigenvectors...\n");
 243
 244 #ifdef GMX_DOUBLE
 245     F77_FUNC(pdseupd, PDSEUPD) (&dovec, "A", select, eigenvalues, eigenvectors,
 246                                 &n, NULL, "I", &n, "SA", &neig, &abstol,
 247                                 resid, &ncv, v, &n, iparam, ipntr,
 248                                 workd, workl, &lworkl, &info);
 249 #else
 250     F77_FUNC(psseupd, PSSEUPD) (&dovec, "A", select, eigenvalues, eigenvectors,
 251                                 &n, NULL, "I", &n, "SA", &neig, &abstol,
 252                                 resid, &ncv, v, &n, iparam, ipntr,
 253                                 workd, workl, &lworkl, &info);
 254 #endif
 255
 256     sfree(v);
 257     sfree(resid);
 258     sfree(workd);
 259     sfree(workl);
 260     sfree(select);
 261 }
 262 #endif
 263
 264
 265 void
 266 sparse_eigensolver(gmx_sparsematrix_t *    A,
 267                    int                     neig,
 268                    real *                  eigenvalues,
 269                    real *                  eigenvectors,
 270                    int                     maxiter)
 271 {
 272     int      iwork[80];
 273     int      iparam[11];
 274     int      ipntr[11];
 275     real *   resid;
 276     real *   workd;
 277     real *   workl;
 278     real *   v;
 279     int      n;
 280     int      ido, info, lworkl, i, ncv, dovec;
 281     real     abstol;
 282     int *    select;
 283     int      iter;
 284
 285 #ifdef GMX_MPI_NOT
 286     MPI_Comm_size( MPI_COMM_WORLD, &n );
 287     if (n > 1)
 288     {
 289         sparse_parallel_eigensolver(A, neig, eigenvalues, eigenvectors, maxiter);
 290         return;
 291     }
 292 #endif
 293
 294     if (eigenvectors != NULL)
 295     {
 296         dovec = 1;
 297     }
 298     else
 299     {
 300         dovec = 0;
 301     }
 302
 303     n   = A->nrow;
 304     ncv = 2*neig;
 305
 306     if (ncv > n)
 307     {
 308         ncv = n;
 309     }
 310
 311     for (i = 0; i < 11; i++)
 312     {
 313         iparam[i] = ipntr[i] = 0;
 314     }
 315
 316     iparam[0] = 1;       /* Don't use explicit shifts */
 317     iparam[2] = maxiter; /* Max number of iterations */
 318     iparam[6] = 1;       /* Standard symmetric eigenproblem */
 319
 320     lworkl = ncv*(8+ncv);
 321     snew(resid, n);
 322     snew(workd, (3*n+4));
 323     snew(workl, lworkl);
 324     snew(select, ncv);
 325     snew(v, n*ncv);
 326
 327     /* Use machine tolerance - roughly 1e-16 in double precision */
 328     abstol = 0;
 329
 330     ido = info = 0;
 331     fprintf(stderr, "Calculation Ritz values and Lanczos vectors, max %d iterations...\n", maxiter);
 332
 333     iter = 1;
 334     do
 335     {
 336 #ifdef GMX_DOUBLE
 337         F77_FUNC(dsaupd, DSAUPD) (&ido, "I", &n, "SA", &neig, &abstol,
 338                                   resid, &ncv, v, &n, iparam, ipntr,
 339                                   workd, iwork, workl, &lworkl, &info);
 340 #else
 341         F77_FUNC(ssaupd, SSAUPD) (&ido, "I", &n, "SA", &neig, &abstol,
 342                                   resid, &ncv, v, &n, iparam, ipntr,
 343                                   workd, iwork, workl, &lworkl, &info);
 344 #endif
 345         if (ido == -1 || ido == 1)
 346         {
 347             gmx_sparsematrix_vector_multiply(A, workd+ipntr[0]-1, workd+ipntr[1]-1);
 348         }
 349
 350         fprintf(stderr, "\rIteration %4d: %3d out of %3d Ritz values converged.", iter++, iparam[4], neig);
 351     }
 352     while (info == 0 && (ido == -1 || ido == 1));
 353
 354     fprintf(stderr, "\n");
 355     if (info == 1)
 356     {
 357         gmx_fatal(FARGS,
 358                   "Maximum number of iterations (%d) reached in Arnoldi\n"
 359                   "diagonalization, but only %d of %d eigenvectors converged.\n",
 360                   maxiter, iparam[4], neig);
 361     }
 362     else if (info != 0)
 363     {
 364         gmx_fatal(FARGS, "Unspecified error from Arnoldi diagonalization:%d\n", info);
 365     }
 366
 367     info = 0;
 368     /* Extract eigenvalues and vectors from data */
 369     fprintf(stderr, "Calculating eigenvalues and eigenvectors...\n");
 370
 371 #ifdef GMX_DOUBLE
 372     F77_FUNC(dseupd, DSEUPD) (&dovec, "A", select, eigenvalues, eigenvectors,
 373                               &n, NULL, "I", &n, "SA", &neig, &abstol,
 374                               resid, &ncv, v, &n, iparam, ipntr,
 375                               workd, workl, &lworkl, &info);
 376 #else
 377     F77_FUNC(sseupd, SSEUPD) (&dovec, "A", select, eigenvalues, eigenvectors,
 378                               &n, NULL, "I", &n, "SA", &neig, &abstol,
 379                               resid, &ncv, v, &n, iparam, ipntr,
 380                               workd, workl, &lworkl, &info);
 381 #endif
 382
 383     sfree(v);
 384     sfree(resid);
 385     sfree(workd);
 386     sfree(workl);
 387     sfree(select);
 388 }