3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
35 #include "eigensolver.h"
37 #include "gromacs/legacyheaders/types/simple.h"
38 #include "gromacs/legacyheaders/gmx_fatal.h"
39 #include "gromacs/legacyheaders/smalloc.h"
41 #include "gromacs/linearalgebra/sparsematrix.h"
43 #include "gmx_lapack.h"
44 #include "gmx_arpack.h"
69 /* Make jobz point to the character "V" if eigenvectors
70 * should be calculated, otherwise "N" (only eigenvalues).
72 jobz = (eigenvectors != NULL) ? "V" : "N";
74 /* allocate lapack stuff */
78 /* First time we ask the routine how much workspace it needs */
83 /* Convert indices to fortran standard */
87 /* Call LAPACK routine using fortran interface. Note that we use upper storage,
88 * but this corresponds to lower storage ("L") in Fortran.
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);
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);
103 gmx_fatal(FARGS,"Internal errror in LAPACK diagonalization.");
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);
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);
130 gmx_fatal(FARGS,"Internal errror in LAPACK diagonalization.");
138 sparse_parallel_eigensolver(gmx_sparsematrix_t * A,
152 int ido,info,lworkl,i,ncv,dovec;
158 MPI_Comm_size( MPI_COMM_WORLD, &nnodes );
159 MPI_Comm_rank( MPI_COMM_WORLD, &rank );
161 if(eigenvectors != NULL)
173 iparam[i]=ipntr[i]=0;
175 iparam[0] = 1; /* Don't use explicit shifts */
176 iparam[2] = maxiter; /* Max number of iterations */
177 iparam[6] = 1; /* Standard symmetric eigenproblem */
179 lworkl = ncv*(8+ncv);
186 /* Use machine tolerance - roughly 1e-16 in double precision */
190 fprintf(stderr,"Calculation Ritz values and Lanczos vectors, max %d iterations...\n",maxiter);
195 F77_FUNC(pdsaupd,PDSAUPD)(&ido, "I", &n, "SA", &neig, &abstol,
196 resid, &ncv, v, &n, iparam, ipntr,
197 workd, iwork, workl, &lworkl, &info);
199 F77_FUNC(pssaupd,PSSAUPD)(&ido, "I", &n, "SA", &neig, &abstol,
200 resid, &ncv, v, &n, iparam, ipntr,
201 workd, iwork, workl, &lworkl, &info);
203 if(ido==-1 || ido==1)
204 gmx_sparsematrix_vector_multiply(A,workd+ipntr[0]-1, workd+ipntr[1]-1);
206 fprintf(stderr,"\rIteration %4d: %3d out of %3d Ritz values converged.",iter++,iparam[4],neig);
207 } while(info==0 && (ido==-1 || ido==1));
209 fprintf(stderr,"\n");
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);
219 gmx_fatal(FARGS,"Unspecified error from Arnoldi diagonalization:%d\n",info);
223 /* Extract eigenvalues and vectors from data */
224 fprintf(stderr,"Calculating eigenvalues and eigenvectors...\n");
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);
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);
248 sparse_eigensolver(gmx_sparsematrix_t * A,
262 int ido,info,lworkl,i,ncv,dovec;
268 MPI_Comm_size( MPI_COMM_WORLD, &n );
271 sparse_parallel_eigensolver(A,neig,eigenvalues,eigenvectors,maxiter);
276 if(eigenvectors != NULL)
288 iparam[i]=ipntr[i]=0;
290 iparam[0] = 1; /* Don't use explicit shifts */
291 iparam[2] = maxiter; /* Max number of iterations */
292 iparam[6] = 1; /* Standard symmetric eigenproblem */
294 lworkl = ncv*(8+ncv);
301 /* Use machine tolerance - roughly 1e-16 in double precision */
305 fprintf(stderr,"Calculation Ritz values and Lanczos vectors, max %d iterations...\n",maxiter);
310 F77_FUNC(dsaupd,DSAUPD)(&ido, "I", &n, "SA", &neig, &abstol,
311 resid, &ncv, v, &n, iparam, ipntr,
312 workd, iwork, workl, &lworkl, &info);
314 F77_FUNC(ssaupd,SSAUPD)(&ido, "I", &n, "SA", &neig, &abstol,
315 resid, &ncv, v, &n, iparam, ipntr,
316 workd, iwork, workl, &lworkl, &info);
318 if(ido==-1 || ido==1)
319 gmx_sparsematrix_vector_multiply(A,workd+ipntr[0]-1, workd+ipntr[1]-1);
321 fprintf(stderr,"\rIteration %4d: %3d out of %3d Ritz values converged.",iter++,iparam[4],neig);
322 } while(info==0 && (ido==-1 || ido==1));
324 fprintf(stderr,"\n");
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);
334 gmx_fatal(FARGS,"Unspecified error from Arnoldi diagonalization:%d\n",info);
338 /* Extract eigenvalues and vectors from data */
339 fprintf(stderr,"Calculating eigenvalues and eigenvectors...\n");
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);
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);