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 * GROningen Mixture of Alchemy and Childrens' Stories
35 /* This file is completely threadsafe - keep it that way! */
47 #define ROTATE(a,i,j,k,l) g=a[i][j];h=a[k][l];a[i][j]=g-s*(h+g*tau);\
48 a[k][l]=h+s*(g-h*tau);
50 void jacobi(double **a
,int n
,double d
[],double **v
,int *nrot
)
54 double tresh
,theta
,tau
,t
,sm
,s
,h
,g
,c
,*b
,*z
;
58 for (ip
=0; ip
<n
; ip
++) {
59 for (iq
=0; iq
<n
; iq
++) v
[ip
][iq
]=0.0;
62 for (ip
=0; ip
<n
;ip
++) {
63 b
[ip
]=d
[ip
]=a
[ip
][ip
];
67 for (i
=1; i
<=50; i
++) {
69 for (ip
=0; ip
<n
-1; ip
++) {
70 for (iq
=ip
+1; iq
<n
; iq
++)
71 sm
+= fabs(a
[ip
][iq
]);
82 for (ip
=0; ip
<n
-1; ip
++) {
83 for (iq
=ip
+1; iq
<n
; iq
++) {
84 g
=100.0*fabs(a
[ip
][iq
]);
85 if (i
> 4 && fabs(d
[ip
])+g
== fabs(d
[ip
])
86 && fabs(d
[iq
])+g
== fabs(d
[iq
]))
88 else if (fabs(a
[ip
][iq
]) > tresh
) {
90 if (fabs(h
)+g
== fabs(h
))
93 theta
=0.5*h
/(a
[ip
][iq
]);
94 t
=1.0/(fabs(theta
)+sqrt(1.0+theta
*theta
));
95 if (theta
< 0.0) t
= -t
;
106 for (j
=0; j
<ip
; j
++) {
109 for (j
=ip
+1; j
<iq
; j
++) {
112 for (j
=iq
+1; j
<n
; j
++) {
115 for (j
=0; j
<n
; j
++) {
122 for (ip
=0; ip
<n
; ip
++) {
128 gmx_fatal(FARGS
,"Error: Too many iterations in routine JACOBI\n");
131 int m_inv_gen(real
**m
,int n
,real
**minv
)
133 double **md
,**v
,*eig
,tol
,s
;
134 int nzero
,i
,j
,k
,nrot
;
149 tol
+= fabs(md
[i
][i
]);
152 jacobi(md
,n
,eig
,v
,&nrot
);
156 if (fabs(eig
[i
]) < tol
) {
166 s
+= eig
[k
]*v
[i
][k
]*v
[j
][k
];