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