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