src/tools/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  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012,2013, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * 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
  39 #ifdef HAVE_CONFIG_H
  40 #include <config.h>
  41 #endif
  42
  43 #include "types/simple.h"
  44 #include "smalloc.h"
  45 #include "gmx_fatal.h"
  46
  47 #include "sparsematrix.h"
  48 #include "eigensolver.h"
  49
  50 #include "gmx_lapack.h"
  51 #include "gmx_arpack.h"
  52
  53 void
  54 eigensolver(real *   a,
  55             int      n,
  56             int      index_lower,
  57             int      index_upper,
  58             real *   eigenvalues,
  59             real *   eigenvectors)
  60 {
  61     int *   isuppz;
  62     int     lwork,liwork;
  63     int     il,iu,m,iw0,info;
  64     real    w0,abstol;
  65     int *   iwork;
  66     real *  work;
  67     real    vl,vu;
  68     const char *  jobz;
  69
  70     if(index_lower<0)
  71         index_lower = 0;
  72
  73     if(index_upper>=n)
  74         index_upper = n-1;
  75
  76     /* Make jobz point to the character "V" if eigenvectors
  77      * should be calculated, otherwise "N" (only eigenvalues).
  78      */
  79     jobz = (eigenvectors != NULL) ? "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 #ifdef 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 = w0;
 114     liwork = iw0;
 115
 116     snew(work,lwork);
 117     snew(iwork,liwork);
 118
 119     abstol = 0;
 120
 121 #ifdef 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         dovec = 1;
 170     else
 171         dovec = 0;
 172
 173     n   = A->nrow;
 174     ncv = 2*neig;
 175
 176     if(ncv>n)
 177         ncv=n;
 178
 179     for(i=0;i<11;i++)
 180         iparam[i]=ipntr[i]=0;
 181
 182         iparam[0] = 1;       /* Don't use explicit shifts */
 183         iparam[2] = maxiter; /* Max number of iterations */
 184         iparam[6] = 1;       /* Standard symmetric eigenproblem */
 185
 186         lworkl = ncv*(8+ncv);
 187     snew(resid,n);
 188     snew(workd,(3*n+4));
 189     snew(workl,lworkl);
 190     snew(select,ncv);
 191     snew(v,n*ncv);
 192
 193     /* Use machine tolerance - roughly 1e-16 in double precision */
 194     abstol = 0;
 195
 196         ido = info = 0;
 197     fprintf(stderr,"Calculation Ritz values and Lanczos vectors, max %d iterations...\n",maxiter);
 198
 199     iter = 1;
 200         do {
 201 #ifdef GMX_DOUBLE
 202                 F77_FUNC(pdsaupd,PDSAUPD)(&ido, "I", &n, "SA", &neig, &abstol,
 203                                                                   resid, &ncv, v, &n, iparam, ipntr,
 204                                                                   workd, iwork, workl, &lworkl, &info);
 205 #else
 206                 F77_FUNC(pssaupd,PSSAUPD)(&ido, "I", &n, "SA", &neig, &abstol,
 207                                                                   resid, &ncv, v, &n, iparam, ipntr,
 208                                                                   workd, iwork, workl, &lworkl, &info);
 209 #endif
 210         if(ido==-1 || ido==1)
 211             gmx_sparsematrix_vector_multiply(A,workd+ipntr[0]-1, workd+ipntr[1]-1);
 212
 213         fprintf(stderr,"\rIteration %4d: %3d out of %3d Ritz values converged.",iter++,iparam[4],neig);
 214         } while(info==0 && (ido==-1 || ido==1));
 215
 216     fprintf(stderr,"\n");
 217         if(info==1)
 218     {
 219             gmx_fatal(FARGS,
 220                   "Maximum number of iterations (%d) reached in Arnoldi\n"
 221                   "diagonalization, but only %d of %d eigenvectors converged.\n",
 222                   maxiter,iparam[4],neig);
 223     }
 224         else if(info!=0)
 225     {
 226         gmx_fatal(FARGS,"Unspecified error from Arnoldi diagonalization:%d\n",info);
 227     }
 228
 229         info = 0;
 230         /* Extract eigenvalues and vectors from data */
 231     fprintf(stderr,"Calculating eigenvalues and eigenvectors...\n");
 232
 233 #ifdef GMX_DOUBLE
 234     F77_FUNC(pdseupd,PDSEUPD)(&dovec, "A", select, eigenvalues, eigenvectors,
 235                                                           &n, NULL, "I", &n, "SA", &neig, &abstol,
 236                                                           resid, &ncv, v, &n, iparam, ipntr,
 237                                                           workd, workl, &lworkl, &info);
 238 #else
 239     F77_FUNC(psseupd,PSSEUPD)(&dovec, "A", select, eigenvalues, eigenvectors,
 240                                                           &n, NULL, "I", &n, "SA", &neig, &abstol,
 241                                                           resid, &ncv, v, &n, iparam, ipntr,
 242                                                           workd, workl, &lworkl, &info);
 243 #endif
 244
 245     sfree(v);
 246     sfree(resid);
 247     sfree(workd);
 248     sfree(workl);
 249     sfree(select);
 250 }
 251 #endif
 252
 253
 254 void
 255 sparse_eigensolver(gmx_sparsematrix_t *    A,
 256                    int                     neig,
 257                    real *                  eigenvalues,
 258                    real *                  eigenvectors,
 259                    int                     maxiter)
 260 {
 261     int      iwork[80];
 262     int      iparam[11];
 263     int      ipntr[11];
 264     real *   resid;
 265     real *   workd;
 266     real *   workl;
 267     real *   v;
 268     int      n;
 269     int      ido,info,lworkl,i,ncv,dovec;
 270     real     abstol;
 271     int *    select;
 272     int      iter;
 273
 274 #ifdef GMX_MPI_NOT
 275         MPI_Comm_size( MPI_COMM_WORLD, &n );
 276         if(n > 1)
 277         {
 278                 sparse_parallel_eigensolver(A,neig,eigenvalues,eigenvectors,maxiter);
 279                 return;
 280         }
 281 #endif
 282
 283     if(eigenvectors != NULL)
 284         dovec = 1;
 285     else
 286         dovec = 0;
 287
 288     n   = A->nrow;
 289     ncv = 2*neig;
 290
 291     if(ncv>n)
 292         ncv=n;
 293
 294     for(i=0;i<11;i++)
 295         iparam[i]=ipntr[i]=0;
 296
 297         iparam[0] = 1;       /* Don't use explicit shifts */
 298         iparam[2] = maxiter; /* Max number of iterations */
 299         iparam[6] = 1;       /* Standard symmetric eigenproblem */
 300
 301         lworkl = ncv*(8+ncv);
 302     snew(resid,n);
 303     snew(workd,(3*n+4));
 304     snew(workl,lworkl);
 305     snew(select,ncv);
 306     snew(v,n*ncv);
 307
 308     /* Use machine tolerance - roughly 1e-16 in double precision */
 309     abstol = 0;
 310
 311         ido = info = 0;
 312     fprintf(stderr,"Calculation Ritz values and Lanczos vectors, max %d iterations...\n",maxiter);
 313
 314     iter = 1;
 315         do {
 316 #ifdef GMX_DOUBLE
 317             F77_FUNC(dsaupd,DSAUPD)(&ido, "I", &n, "SA", &neig, &abstol,
 318                                     resid, &ncv, v, &n, iparam, ipntr,
 319                                     workd, iwork, workl, &lworkl, &info);
 320 #else
 321             F77_FUNC(ssaupd,SSAUPD)(&ido, "I", &n, "SA", &neig, &abstol,
 322                                     resid, &ncv, v, &n, iparam, ipntr,
 323                                     workd, iwork, workl, &lworkl, &info);
 324 #endif
 325         if(ido==-1 || ido==1)
 326             gmx_sparsematrix_vector_multiply(A,workd+ipntr[0]-1, workd+ipntr[1]-1);
 327
 328         fprintf(stderr,"\rIteration %4d: %3d out of %3d Ritz values converged.",iter++,iparam[4],neig);
 329         } while(info==0 && (ido==-1 || ido==1));
 330
 331     fprintf(stderr,"\n");
 332         if(info==1)
 333     {
 334             gmx_fatal(FARGS,
 335                   "Maximum number of iterations (%d) reached in Arnoldi\n"
 336                   "diagonalization, but only %d of %d eigenvectors converged.\n",
 337                   maxiter,iparam[4],neig);
 338     }
 339         else if(info!=0)
 340     {
 341         gmx_fatal(FARGS,"Unspecified error from Arnoldi diagonalization:%d\n",info);
 342     }
 343
 344         info = 0;
 345         /* Extract eigenvalues and vectors from data */
 346     fprintf(stderr,"Calculating eigenvalues and eigenvectors...\n");
 347
 348 #ifdef GMX_DOUBLE
 349     F77_FUNC(dseupd,DSEUPD)(&dovec, "A", select, eigenvalues, eigenvectors,
 350                             &n, NULL, "I", &n, "SA", &neig, &abstol,
 351                             resid, &ncv, v, &n, iparam, ipntr,
 352                             workd, workl, &lworkl, &info);
 353 #else
 354     F77_FUNC(sseupd,SSEUPD)(&dovec, "A", select, eigenvalues, eigenvectors,
 355                             &n, NULL, "I", &n, "SA", &neig, &abstol,
 356                             resid, &ncv, v, &n, iparam, ipntr,
 357                             workd, workl, &lworkl, &info);
 358 #endif
 359
 360     sfree(v);
 361     sfree(resid);
 362     sfree(workd);
 363     sfree(workl);
 364     sfree(select);
 365 }
 366
 367