src/gmxlib/nrjac.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 /* This file is completely threadsafe - keep it that way! */
  39 #ifdef HAVE_CONFIG_H
  40 #include <config.h>
  41 #endif
  42
  43 #include <stdio.h>
  44 #include <stdlib.h>
  45 #include <math.h>
  46 #include "gstat.h"
  47 #include "smalloc.h"
  48 #include "gmx_fatal.h"
  49 #include "nrjac.h"
  50
  51
  52 static inline
  53 void do_rotate(double **a, int i, int j, int k, int l, double tau, double s)
  54 {
  55     double g, h;
  56     g = a[i][j];
  57     h = a[k][l];
  58     a[i][j] = g - s * (h + g * tau);
  59     a[k][l] = h + s * (g - h * tau);
  60 }
  61
  62 void jacobi(double **a,int n,double d[],double **v,int *nrot)
  63 {
  64   int j,i;
  65   int iq,ip;
  66   double tresh,theta,tau,t,sm,s,h,g,c,*b,*z;
  67
  68   snew(b,n);
  69   snew(z,n);
  70   for (ip=0; ip<n; ip++) {
  71     for (iq=0; iq<n; iq++) v[ip][iq]=0.0;
  72     v[ip][ip]=1.0;
  73  }
  74   for (ip=0; ip<n;ip++) {
  75     b[ip]=d[ip]=a[ip][ip];
  76     z[ip]=0.0;
  77   }
  78   *nrot=0;
  79   for (i=1; i<=50; i++) {
  80     sm=0.0;
  81     for (ip=0; ip<n-1; ip++) {
  82       for (iq=ip+1; iq<n; iq++)
  83         sm += fabs(a[ip][iq]);
  84     }
  85     if (sm == 0.0) {
  86       sfree(z);
  87       sfree(b);
  88       return;
  89     }
  90     if (i < 4)
  91       tresh=0.2*sm/(n*n);
  92     else
  93       tresh=0.0;
  94     for (ip=0; ip<n-1; ip++) {
  95       for (iq=ip+1; iq<n; iq++) {
  96         g=100.0*fabs(a[ip][iq]);
  97         if (i > 4 && fabs(d[ip])+g == fabs(d[ip])
  98             && fabs(d[iq])+g == fabs(d[iq]))
  99           a[ip][iq]=0.0;
 100         else if (fabs(a[ip][iq]) > tresh) {
 101           h=d[iq]-d[ip];
 102           if (fabs(h)+g == fabs(h))
 103             t=(a[ip][iq])/h;
 104           else {
 105             theta=0.5*h/(a[ip][iq]);
 106             t=1.0/(fabs(theta)+sqrt(1.0+theta*theta));
 107             if (theta < 0.0) t = -t;
 108           }
 109           c=1.0/sqrt(1+t*t);
 110           s=t*c;
 111           tau=s/(1.0+c);
 112           h=t*a[ip][iq];
 113           z[ip] -= h;
 114           z[iq] += h;
 115           d[ip] -= h;
 116           d[iq] += h;
 117           a[ip][iq]=0.0;
 118           for (j=0; j<ip; j++)
 119           {
 120               do_rotate(a,j,ip,j,iq,tau,s);
 121           }
 122           for (j=ip+1; j<iq; j++)
 123           {
 124               do_rotate(a,ip,j,j,iq,tau,s);
 125           }
 126           for (j=iq+1; j<n; j++)
 127           {
 128               do_rotate(a,ip,j,iq,j,tau,s);
 129           }
 130           for (j=0; j<n; j++)
 131           {
 132               do_rotate(v,j,ip,j,iq,tau,s);
 133           }
 134           ++(*nrot);
 135         }
 136       }
 137     }
 138     for (ip=0; ip<n; ip++) {
 139       b[ip] +=  z[ip];
 140       d[ip]  =  b[ip];
 141       z[ip]  =  0.0;
 142     }
 143   }
 144   gmx_fatal(FARGS,"Error: Too many iterations in routine JACOBI\n");
 145 }
 146
 147 int m_inv_gen(real **m,int n,real **minv)
 148 {
 149   double **md,**v,*eig,tol,s;
 150   int    nzero,i,j,k,nrot;
 151
 152   snew(md,n);
 153   for(i=0; i<n; i++)
 154     snew(md[i],n);
 155   snew(v,n);
 156   for(i=0; i<n; i++)
 157     snew(v[i],n);
 158   snew(eig,n);
 159   for(i=0; i<n; i++)
 160     for(j=0; j<n; j++)
 161       md[i][j] = m[i][j];
 162
 163   tol = 0;
 164   for(i=0; i<n; i++)
 165     tol += fabs(md[i][i]);
 166   tol = 1e-6*tol/n;
 167
 168   jacobi(md,n,eig,v,&nrot);
 169
 170   nzero = 0;
 171   for(i=0; i<n; i++)
 172     if (fabs(eig[i]) < tol) {
 173       eig[i] = 0;
 174       nzero++;
 175     } else
 176       eig[i] = 1.0/eig[i];
 177
 178   for(i=0; i<n; i++)
 179     for(j=0; j<n; j++) {
 180       s = 0;
 181       for(k=0; k<n; k++)
 182         s += eig[k]*v[i][k]*v[j][k];
 183       minv[i][j] = s;
 184     }
 185
 186   sfree(eig);
 187   for(i=0; i<n; i++)
 188     sfree(v[i]);
 189   sfree(v);
 190   for(i=0; i<n; i++)
 191     sfree(md[i]);
 192   sfree(md);
 193
 194   return nzero;
 195 }