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