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"
42 #include "gmx_lapack.h"
43 #include "gmx_arpack.h"
55 int il, iu, m, iw0, info;
72 /* Make jobz point to the character "V" if eigenvectors
73 * should be calculated, otherwise "N" (only eigenvalues).
75 jobz = (eigenvectors != NULL) ? "V" : "N";
77 /* allocate lapack stuff */
81 /* First time we ask the routine how much workspace it needs */
86 /* Convert indices to fortran standard */
90 /* Call LAPACK routine using fortran interface. Note that we use upper storage,
91 * but this corresponds to lower storage ("L") in Fortran.
94 F77_FUNC(dsyevr, DSYEVR) (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);
98 F77_FUNC(ssyevr, SSYEVR) (jobz, "I", "L", &n, a, &n, &vl, &vu, &index_lower, &index_upper,
99 &abstol, &m, eigenvalues, eigenvectors, &n,
100 isuppz, &w0, &lwork, &iw0, &liwork, &info);
106 gmx_fatal(FARGS, "Internal errror in LAPACK diagonalization.");
118 F77_FUNC(dsyevr, DSYEVR) (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);
122 F77_FUNC(ssyevr, SSYEVR) (jobz, "I", "L", &n, a, &n, &vl, &vu, &index_lower, &index_upper,
123 &abstol, &m, eigenvalues, eigenvectors, &n,
124 isuppz, work, &lwork, iwork, &liwork, &info);
133 gmx_fatal(FARGS, "Internal errror in LAPACK diagonalization.");
141 sparse_parallel_eigensolver(gmx_sparsematrix_t * A,
155 int ido, info, lworkl, i, ncv, dovec;
161 MPI_Comm_size( MPI_COMM_WORLD, &nnodes );
162 MPI_Comm_rank( MPI_COMM_WORLD, &rank );
164 if (eigenvectors != NULL)
181 for (i = 0; i < 11; i++)
183 iparam[i] = ipntr[i] = 0;
186 iparam[0] = 1; /* Don't use explicit shifts */
187 iparam[2] = maxiter; /* Max number of iterations */
188 iparam[6] = 1; /* Standard symmetric eigenproblem */
190 lworkl = ncv*(8+ncv);
192 snew(workd, (3*n+4));
197 /* Use machine tolerance - roughly 1e-16 in double precision */
201 fprintf(stderr, "Calculation Ritz values and Lanczos vectors, max %d iterations...\n", maxiter);
207 F77_FUNC(pdsaupd, PDSAUPD) (&ido, "I", &n, "SA", &neig, &abstol,
208 resid, &ncv, v, &n, iparam, ipntr,
209 workd, iwork, workl, &lworkl, &info);
211 F77_FUNC(pssaupd, PSSAUPD) (&ido, "I", &n, "SA", &neig, &abstol,
212 resid, &ncv, v, &n, iparam, ipntr,
213 workd, iwork, workl, &lworkl, &info);
215 if (ido == -1 || ido == 1)
217 gmx_sparsematrix_vector_multiply(A, workd+ipntr[0]-1, workd+ipntr[1]-1);
220 fprintf(stderr, "\rIteration %4d: %3d out of %3d Ritz values converged.", iter++, iparam[4], neig);
222 while (info == 0 && (ido == -1 || ido == 1));
224 fprintf(stderr, "\n");
228 "Maximum number of iterations (%d) reached in Arnoldi\n"
229 "diagonalization, but only %d of %d eigenvectors converged.\n",
230 maxiter, iparam[4], neig);
234 gmx_fatal(FARGS, "Unspecified error from Arnoldi diagonalization:%d\n", info);
238 /* Extract eigenvalues and vectors from data */
239 fprintf(stderr, "Calculating eigenvalues and eigenvectors...\n");
242 F77_FUNC(pdseupd, PDSEUPD) (&dovec, "A", select, eigenvalues, eigenvectors,
243 &n, NULL, "I", &n, "SA", &neig, &abstol,
244 resid, &ncv, v, &n, iparam, ipntr,
245 workd, workl, &lworkl, &info);
247 F77_FUNC(psseupd, PSSEUPD) (&dovec, "A", select, eigenvalues, eigenvectors,
248 &n, NULL, "I", &n, "SA", &neig, &abstol,
249 resid, &ncv, v, &n, iparam, ipntr,
250 workd, workl, &lworkl, &info);
263 sparse_eigensolver(gmx_sparsematrix_t * A,
277 int ido, info, lworkl, i, ncv, dovec;
283 MPI_Comm_size( MPI_COMM_WORLD, &n );
286 sparse_parallel_eigensolver(A, neig, eigenvalues, eigenvectors, maxiter);
291 if (eigenvectors != NULL)
308 for (i = 0; i < 11; i++)
310 iparam[i] = ipntr[i] = 0;
313 iparam[0] = 1; /* Don't use explicit shifts */
314 iparam[2] = maxiter; /* Max number of iterations */
315 iparam[6] = 1; /* Standard symmetric eigenproblem */
317 lworkl = ncv*(8+ncv);
319 snew(workd, (3*n+4));
324 /* Use machine tolerance - roughly 1e-16 in double precision */
328 fprintf(stderr, "Calculation Ritz values and Lanczos vectors, max %d iterations...\n", maxiter);
334 F77_FUNC(dsaupd, DSAUPD) (&ido, "I", &n, "SA", &neig, &abstol,
335 resid, &ncv, v, &n, iparam, ipntr,
336 workd, iwork, workl, &lworkl, &info);
338 F77_FUNC(ssaupd, SSAUPD) (&ido, "I", &n, "SA", &neig, &abstol,
339 resid, &ncv, v, &n, iparam, ipntr,
340 workd, iwork, workl, &lworkl, &info);
342 if (ido == -1 || ido == 1)
344 gmx_sparsematrix_vector_multiply(A, workd+ipntr[0]-1, workd+ipntr[1]-1);
347 fprintf(stderr, "\rIteration %4d: %3d out of %3d Ritz values converged.", iter++, iparam[4], neig);
349 while (info == 0 && (ido == -1 || ido == 1));
351 fprintf(stderr, "\n");
355 "Maximum number of iterations (%d) reached in Arnoldi\n"
356 "diagonalization, but only %d of %d eigenvectors converged.\n",
357 maxiter, iparam[4], neig);
361 gmx_fatal(FARGS, "Unspecified error from Arnoldi diagonalization:%d\n", info);
365 /* Extract eigenvalues and vectors from data */
366 fprintf(stderr, "Calculating eigenvalues and eigenvectors...\n");
369 F77_FUNC(dseupd, DSEUPD) (&dovec, "A", select, eigenvalues, eigenvectors,
370 &n, NULL, "I", &n, "SA", &neig, &abstol,
371 resid, &ncv, v, &n, iparam, ipntr,
372 workd, workl, &lworkl, &info);
374 F77_FUNC(sseupd, SSEUPD) (&dovec, "A", select, eigenvalues, eigenvectors,
375 &n, NULL, "I", &n, "SA", &neig, &abstol,
376 resid, &ncv, v, &n, iparam, ipntr,
377 workd, workl, &lworkl, &info);