4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
31 * For more info, check our website at http://www.gromacs.org
34 * GROningen Mixture of Alchemy and Childrens' Stories
36 /* This file is completely threadsafe - keep it that way! */
45 #define ROTATE(a,i,j,k,l) g=a[i][j];h=a[k][l];a[i][j]=g-s*(h+g*tau);\
46 a[k][l]=h+s*(g-h*tau);
48 void jacobi(double **a
,int n
,double d
[],double **v
,int *nrot
)
52 double tresh
,theta
,tau
,t
,sm
,s
,h
,g
,c
,*b
,*z
;
56 for (ip
=0; ip
<n
; ip
++) {
57 for (iq
=0; iq
<n
; iq
++) v
[ip
][iq
]=0.0;
60 for (ip
=0; ip
<n
;ip
++) {
61 b
[ip
]=d
[ip
]=a
[ip
][ip
];
65 for (i
=1; i
<=50; i
++) {
67 for (ip
=0; ip
<n
-1; ip
++) {
68 for (iq
=ip
+1; iq
<n
; iq
++)
69 sm
+= fabs(a
[ip
][iq
]);
80 for (ip
=0; ip
<n
-1; ip
++) {
81 for (iq
=ip
+1; iq
<n
; iq
++) {
82 g
=100.0*fabs(a
[ip
][iq
]);
83 if (i
> 4 && fabs(d
[ip
])+g
== fabs(d
[ip
])
84 && fabs(d
[iq
])+g
== fabs(d
[iq
]))
86 else if (fabs(a
[ip
][iq
]) > tresh
) {
88 if (fabs(h
)+g
== fabs(h
))
91 theta
=0.5*h
/(a
[ip
][iq
]);
92 t
=1.0/(fabs(theta
)+sqrt(1.0+theta
*theta
));
93 if (theta
< 0.0) t
= -t
;
104 for (j
=0; j
<ip
; j
++) {
107 for (j
=ip
+1; j
<iq
; j
++) {
110 for (j
=iq
+1; j
<n
; j
++) {
113 for (j
=0; j
<n
; j
++) {
120 for (ip
=0; ip
<n
; ip
++) {
126 fatal_error(0,"Error: Too many iterations in routine JACOBI\n");
129 int m_inv_gen(real
**m
,int n
,real
**minv
)
131 double **md
,**v
,*eig
,tol
,s
;
132 int nzero
,i
,j
,k
,nrot
;
147 tol
+= fabs(md
[i
][i
]);
150 jacobi(md
,n
,eig
,v
,&nrot
);
154 if (fabs(eig
[i
]) < tol
) {
164 s
+= eig
[k
]*v
[i
][k
]*v
[j
][k
];