src/gromacs/linearalgebra/eigensolver.c

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