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
40 #include "types/simple.h"
42 #include "gmx_fatal.h"
44 #include "sparsematrix.h"
45 #include "eigensolver.h"
48 #define F77_FUNC(name,NAME) name ## _
51 #include "gmx_lapack.h"
52 #include "gmx_arpack.h"
77 /* Make jobz point to the character "V" if eigenvectors
78 * should be calculated, otherwise "N" (only eigenvalues).
80 jobz = (eigenvectors != NULL) ? "V" : "N";
82 /* allocate lapack stuff */
86 /* First time we ask the routine how much workspace it needs */
91 /* Convert indices to fortran standard */
95 /* Call LAPACK routine using fortran interface. Note that we use upper storage,
96 * but this corresponds to lower storage ("L") in Fortran.
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);
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);
111 gmx_fatal(FARGS,"Internal errror in LAPACK diagonalization.");
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);
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);
138 gmx_fatal(FARGS,"Internal errror in LAPACK diagonalization.");
146 sparse_parallel_eigensolver(gmx_sparsematrix_t * A,
160 int ido,info,lworkl,i,ncv,dovec;
166 MPI_Comm_size( MPI_COMM_WORLD, &nnodes );
167 MPI_Comm_rank( MPI_COMM_WORLD, &rank );
169 if(eigenvectors != NULL)
181 iparam[i]=ipntr[i]=0;
183 iparam[0] = 1; /* Don't use explicit shifts */
184 iparam[2] = maxiter; /* Max number of iterations */
185 iparam[6] = 1; /* Standard symmetric eigenproblem */
187 lworkl = ncv*(8+ncv);
194 /* Use machine tolerance - roughly 1e-16 in double precision */
198 fprintf(stderr,"Calculation Ritz values and Lanczos vectors, max %d iterations...\n",maxiter);
203 F77_FUNC(pdsaupd,PDSAUPD)(&ido, "I", &n, "SA", &neig, &abstol,
204 resid, &ncv, v, &n, iparam, ipntr,
205 workd, iwork, workl, &lworkl, &info);
207 F77_FUNC(pssaupd,PSSAUPD)(&ido, "I", &n, "SA", &neig, &abstol,
208 resid, &ncv, v, &n, iparam, ipntr,
209 workd, iwork, workl, &lworkl, &info);
211 if(ido==-1 || ido==1)
212 gmx_sparsematrix_vector_multiply(A,workd+ipntr[0]-1, workd+ipntr[1]-1);
214 fprintf(stderr,"\rIteration %4d: %3d out of %3d Ritz values converged.",iter++,iparam[4],neig);
215 } while(info==0 && (ido==-1 || ido==1));
217 fprintf(stderr,"\n");
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);
227 gmx_fatal(FARGS,"Unspecified error from Arnoldi diagonalization:%d\n",info);
231 /* Extract eigenvalues and vectors from data */
232 fprintf(stderr,"Calculating eigenvalues and eigenvectors...\n");
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);
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);
256 sparse_eigensolver(gmx_sparsematrix_t * A,
270 int ido,info,lworkl,i,ncv,dovec;
276 MPI_Comm_size( MPI_COMM_WORLD, &n );
279 sparse_parallel_eigensolver(A,neig,eigenvalues,eigenvectors,maxiter);
284 if(eigenvectors != NULL)
296 iparam[i]=ipntr[i]=0;
298 iparam[0] = 1; /* Don't use explicit shifts */
299 iparam[2] = maxiter; /* Max number of iterations */
300 iparam[6] = 1; /* Standard symmetric eigenproblem */
302 lworkl = ncv*(8+ncv);
309 /* Use machine tolerance - roughly 1e-16 in double precision */
313 fprintf(stderr,"Calculation Ritz values and Lanczos vectors, max %d iterations...\n",maxiter);
318 F77_FUNC(dsaupd,DSAUPD)(&ido, "I", &n, "SA", &neig, &abstol,
319 resid, &ncv, v, &n, iparam, ipntr,
320 workd, iwork, workl, &lworkl, &info);
322 F77_FUNC(ssaupd,SSAUPD)(&ido, "I", &n, "SA", &neig, &abstol,
323 resid, &ncv, v, &n, iparam, ipntr,
324 workd, iwork, workl, &lworkl, &info);
326 if(ido==-1 || ido==1)
327 gmx_sparsematrix_vector_multiply(A,workd+ipntr[0]-1, workd+ipntr[1]-1);
329 fprintf(stderr,"\rIteration %4d: %3d out of %3d Ritz values converged.",iter++,iparam[4],neig);
330 } while(info==0 && (ido==-1 || ido==1));
332 fprintf(stderr,"\n");
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);
342 gmx_fatal(FARGS,"Unspecified error from Arnoldi diagonalization:%d\n",info);
346 /* Extract eigenvalues and vectors from data */
347 fprintf(stderr,"Calculating eigenvalues and eigenvectors...\n");
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);
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);