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