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! */
45 #include "gmx_fatal.h"
50 void do_rotate(double **a, int i, int j, int k, int l, double tau, double s)
55 a[i][j] = g - s * (h + g * tau);
56 a[k][l] = h + s * (g - h * tau);
59 void jacobi(double **a, int n, double d[], double **v, int *nrot)
63 double tresh, theta, tau, t, sm, s, h, g, c, *b, *z;
67 for (ip = 0; ip < n; ip++)
69 for (iq = 0; iq < n; iq++)
75 for (ip = 0; ip < n; ip++)
77 b[ip] = d[ip] = a[ip][ip];
81 for (i = 1; i <= 50; i++)
84 for (ip = 0; ip < n-1; ip++)
86 for (iq = ip+1; iq < n; iq++)
88 sm += fabs(a[ip][iq]);
105 for (ip = 0; ip < n-1; ip++)
107 for (iq = ip+1; iq < n; iq++)
109 g = 100.0*fabs(a[ip][iq]);
110 if (i > 4 && fabs(d[ip])+g == fabs(d[ip])
111 && fabs(d[iq])+g == fabs(d[iq]))
115 else if (fabs(a[ip][iq]) > tresh)
118 if (fabs(h)+g == fabs(h))
124 theta = 0.5*h/(a[ip][iq]);
125 t = 1.0/(fabs(theta)+sqrt(1.0+theta*theta));
140 for (j = 0; j < ip; j++)
142 do_rotate(a, j, ip, j, iq, tau, s);
144 for (j = ip+1; j < iq; j++)
146 do_rotate(a, ip, j, j, iq, tau, s);
148 for (j = iq+1; j < n; j++)
150 do_rotate(a, ip, j, iq, j, tau, s);
152 for (j = 0; j < n; j++)
154 do_rotate(v, j, ip, j, iq, tau, s);
160 for (ip = 0; ip < n; ip++)
167 gmx_fatal(FARGS, "Error: Too many iterations in routine JACOBI\n");
170 int m_inv_gen(real **m, int n, real **minv)
172 double **md, **v, *eig, tol, s;
173 int nzero, i, j, k, nrot;
176 for (i = 0; i < n; i++)
181 for (i = 0; i < n; i++)
186 for (i = 0; i < n; i++)
188 for (j = 0; j < n; j++)
195 for (i = 0; i < n; i++)
197 tol += fabs(md[i][i]);
201 jacobi(md, n, eig, v, &nrot);
204 for (i = 0; i < n; i++)
206 if (fabs(eig[i]) < tol)
217 for (i = 0; i < n; i++)
219 for (j = 0; j < n; j++)
222 for (k = 0; k < n; k++)
224 s += eig[k]*v[i][k]*v[j][k];
231 for (i = 0; i < n; i++)
236 for (i = 0; i < n; i++)