Changed the pV term to the correct one, with the reference pressure. Also made it
[gromacs/rigid-bodies.git] / src / tools / eigensolver.c
blob0708a507240e9dd1326f8ac4ecada11fd011466f
1 /*
2 *
3 * This source code is part of
4 *
5 * G R O M A C S
6 *
7 * GROningen MAchine for Chemical Simulations
8 *
9 * VERSION 3.2.0
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
32 * And Hey:
33 * Green Red Orange Magenta Azure Cyan Skyblue
36 #ifdef HAVE_CONFIG_H
37 #include <config.h>
38 #endif
40 #include "types/simple.h"
41 #include "smalloc.h"
42 #include "gmx_fatal.h"
44 #include "sparsematrix.h"
45 #include "eigensolver.h"
47 #ifndef F77_FUNC
48 #define F77_FUNC(name,NAME) name ## _
49 #endif
51 #include "gmx_lapack.h"
52 #include "gmx_arpack.h"
54 void
55 eigensolver(real * a,
56 int n,
57 int index_lower,
58 int index_upper,
59 real * eigenvalues,
60 real * eigenvectors)
62 int * isuppz;
63 int lwork,liwork;
64 int il,iu,m,iw0,info;
65 real w0,abstol;
66 int * iwork;
67 real * work;
68 real vl,vu;
69 const char * jobz;
71 if(index_lower<0)
72 index_lower = 0;
74 if(index_upper>=n)
75 index_upper = n-1;
77 /* Make jobz point to the character "V" if eigenvectors
78 * should be calculated, otherwise "N" (only eigenvalues).
79 */
80 jobz = (eigenvectors != NULL) ? "V" : "N";
82 /* allocate lapack stuff */
83 snew(isuppz,2*n);
84 vl = vu = 0;
86 /* First time we ask the routine how much workspace it needs */
87 lwork = -1;
88 liwork = -1;
89 abstol = 0;
91 /* Convert indices to fortran standard */
92 index_lower++;
93 index_upper++;
95 /* Call LAPACK routine using fortran interface. Note that we use upper storage,
96 * but this corresponds to lower storage ("L") in Fortran.
97 */
98 #ifdef GMX_DOUBLE
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);
102 #else
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);
106 #endif
108 if(info != 0)
110 sfree(isuppz);
111 gmx_fatal(FARGS,"Internal errror in LAPACK diagonalization.");
114 lwork = w0;
115 liwork = iw0;
117 snew(work,lwork);
118 snew(iwork,liwork);
120 abstol = 0;
122 #ifdef GMX_DOUBLE
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);
126 #else
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);
130 #endif
132 sfree(isuppz);
133 sfree(work);
134 sfree(iwork);
136 if(info != 0)
138 gmx_fatal(FARGS,"Internal errror in LAPACK diagonalization.");
144 #ifdef GMX_MPI_NOT
145 void
146 sparse_parallel_eigensolver(gmx_sparsematrix_t * A,
147 int neig,
148 real * eigenvalues,
149 real * eigenvectors,
150 int maxiter)
152 int iwork[80];
153 int iparam[11];
154 int ipntr[11];
155 real * resid;
156 real * workd;
157 real * workl;
158 real * v;
159 int n;
160 int ido,info,lworkl,i,ncv,dovec;
161 real abstol;
162 int * select;
163 int iter;
164 int nnodes,rank;
166 MPI_Comm_size( MPI_COMM_WORLD, &nnodes );
167 MPI_Comm_rank( MPI_COMM_WORLD, &rank );
169 if(eigenvectors != NULL)
170 dovec = 1;
171 else
172 dovec = 0;
174 n = A->nrow;
175 ncv = 2*neig;
177 if(ncv>n)
178 ncv=n;
180 for(i=0;i<11;i++)
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);
188 snew(resid,n);
189 snew(workd,(3*n+4));
190 snew(workl,lworkl);
191 snew(select,ncv);
192 snew(v,n*ncv);
194 /* Use machine tolerance - roughly 1e-16 in double precision */
195 abstol = 0;
197 ido = info = 0;
198 fprintf(stderr,"Calculation Ritz values and Lanczos vectors, max %d iterations...\n",maxiter);
200 iter = 1;
201 do {
202 #ifdef GMX_DOUBLE
203 F77_FUNC(pdsaupd,PDSAUPD)(&ido, "I", &n, "SA", &neig, &abstol,
204 resid, &ncv, v, &n, iparam, ipntr,
205 workd, iwork, workl, &lworkl, &info);
206 #else
207 F77_FUNC(pssaupd,PSSAUPD)(&ido, "I", &n, "SA", &neig, &abstol,
208 resid, &ncv, v, &n, iparam, ipntr,
209 workd, iwork, workl, &lworkl, &info);
210 #endif
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");
218 if(info==1)
220 gmx_fatal(FARGS,
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);
225 else if(info!=0)
227 gmx_fatal(FARGS,"Unspecified error from Arnoldi diagonalization:%d\n",info);
230 info = 0;
231 /* Extract eigenvalues and vectors from data */
232 fprintf(stderr,"Calculating eigenvalues and eigenvectors...\n");
234 #ifdef GMX_DOUBLE
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);
239 #else
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);
244 #endif
246 sfree(v);
247 sfree(resid);
248 sfree(workd);
249 sfree(workl);
250 sfree(select);
252 #endif
255 void
256 sparse_eigensolver(gmx_sparsematrix_t * A,
257 int neig,
258 real * eigenvalues,
259 real * eigenvectors,
260 int maxiter)
262 int iwork[80];
263 int iparam[11];
264 int ipntr[11];
265 real * resid;
266 real * workd;
267 real * workl;
268 real * v;
269 int n;
270 int ido,info,lworkl,i,ncv,dovec;
271 real abstol;
272 int * select;
273 int iter;
275 #ifdef GMX_MPI_NOT
276 MPI_Comm_size( MPI_COMM_WORLD, &n );
277 if(n > 1)
279 sparse_parallel_eigensolver(A,neig,eigenvalues,eigenvectors,maxiter);
280 return;
282 #endif
284 if(eigenvectors != NULL)
285 dovec = 1;
286 else
287 dovec = 0;
289 n = A->nrow;
290 ncv = 2*neig;
292 if(ncv>n)
293 ncv=n;
295 for(i=0;i<11;i++)
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);
303 snew(resid,n);
304 snew(workd,(3*n+4));
305 snew(workl,lworkl);
306 snew(select,ncv);
307 snew(v,n*ncv);
309 /* Use machine tolerance - roughly 1e-16 in double precision */
310 abstol = 0;
312 ido = info = 0;
313 fprintf(stderr,"Calculation Ritz values and Lanczos vectors, max %d iterations...\n",maxiter);
315 iter = 1;
316 do {
317 #ifdef GMX_DOUBLE
318 F77_FUNC(dsaupd,DSAUPD)(&ido, "I", &n, "SA", &neig, &abstol,
319 resid, &ncv, v, &n, iparam, ipntr,
320 workd, iwork, workl, &lworkl, &info);
321 #else
322 F77_FUNC(ssaupd,SSAUPD)(&ido, "I", &n, "SA", &neig, &abstol,
323 resid, &ncv, v, &n, iparam, ipntr,
324 workd, iwork, workl, &lworkl, &info);
325 #endif
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");
333 if(info==1)
335 gmx_fatal(FARGS,
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);
340 else if(info!=0)
342 gmx_fatal(FARGS,"Unspecified error from Arnoldi diagonalization:%d\n",info);
345 info = 0;
346 /* Extract eigenvalues and vectors from data */
347 fprintf(stderr,"Calculating eigenvalues and eigenvectors...\n");
349 #ifdef GMX_DOUBLE
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);
354 #else
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);
359 #endif
361 sfree(v);
362 sfree(resid);
363 sfree(workd);
364 sfree(workl);
365 sfree(select);