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
);