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