2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2012,2013, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
37 #include "eigensolver.h"
39 #include "gromacs/legacyheaders/types/simple.h"
40 #include "gromacs/legacyheaders/gmx_fatal.h"
41 #include "gromacs/legacyheaders/smalloc.h"
43 #include "gromacs/linearalgebra/sparsematrix.h"
44 #include "gmx_lapack.h"
45 #include "gmx_arpack.h"
57 int il, iu, m, iw0, info;
74 /* Make jobz point to the character "V" if eigenvectors
75 * should be calculated, otherwise "N" (only eigenvalues).
77 jobz = (eigenvectors != NULL) ? "V" : "N";
79 /* allocate lapack stuff */
83 /* First time we ask the routine how much workspace it needs */
88 /* Convert indices to fortran standard */
92 /* Call LAPACK routine using fortran interface. Note that we use upper storage,
93 * but this corresponds to lower storage ("L") in Fortran.
96 F77_FUNC(dsyevr, DSYEVR) (jobz, "I", "L", &n, a, &n, &vl, &vu, &index_lower, &index_upper,
97 &abstol, &m, eigenvalues, eigenvectors, &n,
98 isuppz, &w0, &lwork, &iw0, &liwork, &info);
100 F77_FUNC(ssyevr, SSYEVR) (jobz, "I", "L", &n, a, &n, &vl, &vu, &index_lower, &index_upper,
101 &abstol, &m, eigenvalues, eigenvectors, &n,
102 isuppz, &w0, &lwork, &iw0, &liwork, &info);
108 gmx_fatal(FARGS, "Internal errror in LAPACK diagonalization.");
120 F77_FUNC(dsyevr, DSYEVR) (jobz, "I", "L", &n, a, &n, &vl, &vu, &index_lower, &index_upper,
121 &abstol, &m, eigenvalues, eigenvectors, &n,
122 isuppz, work, &lwork, iwork, &liwork, &info);
124 F77_FUNC(ssyevr, SSYEVR) (jobz, "I", "L", &n, a, &n, &vl, &vu, &index_lower, &index_upper,
125 &abstol, &m, eigenvalues, eigenvectors, &n,
126 isuppz, work, &lwork, iwork, &liwork, &info);
135 gmx_fatal(FARGS, "Internal errror in LAPACK diagonalization.");
143 sparse_parallel_eigensolver(gmx_sparsematrix_t * A,
157 int ido, info, lworkl, i, ncv, dovec;
163 MPI_Comm_size( MPI_COMM_WORLD, &nnodes );
164 MPI_Comm_rank( MPI_COMM_WORLD, &rank );
166 if (eigenvectors != NULL)
183 for (i = 0; i < 11; i++)
185 iparam[i] = ipntr[i] = 0;
188 iparam[0] = 1; /* Don't use explicit shifts */
189 iparam[2] = maxiter; /* Max number of iterations */
190 iparam[6] = 1; /* Standard symmetric eigenproblem */
192 lworkl = ncv*(8+ncv);
194 snew(workd, (3*n+4));
199 /* Use machine tolerance - roughly 1e-16 in double precision */
203 fprintf(stderr, "Calculation Ritz values and Lanczos vectors, max %d iterations...\n", maxiter);
209 F77_FUNC(pdsaupd, PDSAUPD) (&ido, "I", &n, "SA", &neig, &abstol,
210 resid, &ncv, v, &n, iparam, ipntr,
211 workd, iwork, workl, &lworkl, &info);
213 F77_FUNC(pssaupd, PSSAUPD) (&ido, "I", &n, "SA", &neig, &abstol,
214 resid, &ncv, v, &n, iparam, ipntr,
215 workd, iwork, workl, &lworkl, &info);
217 if (ido == -1 || ido == 1)
219 gmx_sparsematrix_vector_multiply(A, workd+ipntr[0]-1, workd+ipntr[1]-1);
222 fprintf(stderr, "\rIteration %4d: %3d out of %3d Ritz values converged.", iter++, iparam[4], neig);
224 while (info == 0 && (ido == -1 || ido == 1));
226 fprintf(stderr, "\n");
230 "Maximum number of iterations (%d) reached in Arnoldi\n"
231 "diagonalization, but only %d of %d eigenvectors converged.\n",
232 maxiter, iparam[4], neig);
236 gmx_fatal(FARGS, "Unspecified error from Arnoldi diagonalization:%d\n", info);
240 /* Extract eigenvalues and vectors from data */
241 fprintf(stderr, "Calculating eigenvalues and eigenvectors...\n");
244 F77_FUNC(pdseupd, PDSEUPD) (&dovec, "A", select, eigenvalues, eigenvectors,
245 &n, NULL, "I", &n, "SA", &neig, &abstol,
246 resid, &ncv, v, &n, iparam, ipntr,
247 workd, workl, &lworkl, &info);
249 F77_FUNC(psseupd, PSSEUPD) (&dovec, "A", select, eigenvalues, eigenvectors,
250 &n, NULL, "I", &n, "SA", &neig, &abstol,
251 resid, &ncv, v, &n, iparam, ipntr,
252 workd, workl, &lworkl, &info);
265 sparse_eigensolver(gmx_sparsematrix_t * A,
279 int ido, info, lworkl, i, ncv, dovec;
285 MPI_Comm_size( MPI_COMM_WORLD, &n );
288 sparse_parallel_eigensolver(A, neig, eigenvalues, eigenvectors, maxiter);
293 if (eigenvectors != NULL)
310 for (i = 0; i < 11; i++)
312 iparam[i] = ipntr[i] = 0;
315 iparam[0] = 1; /* Don't use explicit shifts */
316 iparam[2] = maxiter; /* Max number of iterations */
317 iparam[6] = 1; /* Standard symmetric eigenproblem */
319 lworkl = ncv*(8+ncv);
321 snew(workd, (3*n+4));
326 /* Use machine tolerance - roughly 1e-16 in double precision */
330 fprintf(stderr, "Calculation Ritz values and Lanczos vectors, max %d iterations...\n", maxiter);
336 F77_FUNC(dsaupd, DSAUPD) (&ido, "I", &n, "SA", &neig, &abstol,
337 resid, &ncv, v, &n, iparam, ipntr,
338 workd, iwork, workl, &lworkl, &info);
340 F77_FUNC(ssaupd, SSAUPD) (&ido, "I", &n, "SA", &neig, &abstol,
341 resid, &ncv, v, &n, iparam, ipntr,
342 workd, iwork, workl, &lworkl, &info);
344 if (ido == -1 || ido == 1)
346 gmx_sparsematrix_vector_multiply(A, workd+ipntr[0]-1, workd+ipntr[1]-1);
349 fprintf(stderr, "\rIteration %4d: %3d out of %3d Ritz values converged.", iter++, iparam[4], neig);
351 while (info == 0 && (ido == -1 || ido == 1));
353 fprintf(stderr, "\n");
357 "Maximum number of iterations (%d) reached in Arnoldi\n"
358 "diagonalization, but only %d of %d eigenvectors converged.\n",
359 maxiter, iparam[4], neig);
363 gmx_fatal(FARGS, "Unspecified error from Arnoldi diagonalization:%d\n", info);
367 /* Extract eigenvalues and vectors from data */
368 fprintf(stderr, "Calculating eigenvalues and eigenvectors...\n");
371 F77_FUNC(dseupd, DSEUPD) (&dovec, "A", select, eigenvalues, eigenvectors,
372 &n, NULL, "I", &n, "SA", &neig, &abstol,
373 resid, &ncv, v, &n, iparam, ipntr,
374 workd, workl, &lworkl, &info);
376 F77_FUNC(sseupd, SSEUPD) (&dovec, "A", select, eigenvalues, eigenvectors,
377 &n, NULL, "I", &n, "SA", &neig, &abstol,
378 resid, &ncv, v, &n, iparam, ipntr,
379 workd, workl, &lworkl, &info);