src/tools/ql77.c

   1 /*
   2  * $Id$
   3  *
   4  *       This source code is part of
   5  *
   6  *        G   R   O   M   A   C   S
   7  *
   8  * GROningen MAchine for Chemical Simulations
   9  *
  10  *               VERSION 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * Great Red Oystrich Makes All Chemists Sane
  28  */
  29 static char *SRCID_ql77_c = "$Id$";
  30
  31 #include <math.h>
  32 #include "typedefs.h"
  33 #include "vec.h"
  34 #include "smalloc.h"
  35
  36 void ql77 (int n,real *x,real *d)
  37      /*c
  38        c     matrix diagonalization routine
  39        c
  40        c     code derived from eispack
  41        c     update : fmp 16.2.86 (conversion to fortran 77)
  42        c     update : fmp 09.3.89 (get rid of double precision)
  43        c
  44        c     update : bh  16.5.99 (now works, real C code iso f2c)
  45        c
  46        c
  47        c     input
  48        c
  49        c     n           order of the matrix, not larger than nmax
  50        c     x(n*n)      the real symmetric matrix to be diagonalized as a
  51        c                 two-indexed square array of which only the lower left
  52        c                 triangle is used and need be supplied
  53        c
  54        c     output
  55        c
  56        c     d(n)        computed eigenvalues in descending order
  57        c     x(n*n)      eigenvectors (in columns) in the same order as the
  58        c                 eigenvalues
  59        c
  60        */
  61 {
  62   int i,j,k,l,ni;
  63   real *e,h,g,f,b,s,p,r,c,absp;
  64   real totwork,work;
  65
  66 #define pr_pr(a,b,c) fprintf(stderr,"\rreduction: %g%%  accumulation:  %g%%  accumulation: %g%%",a,b,c)
  67
  68   const real eps=7.e-14,tol=1.e-30;
  69
  70   snew(e,n);
  71
  72   /*
  73    *  householder reduction to tridiagonal form
  74    */
  75
  76   totwork = 0;
  77   for(ni=1; ni<n; ni++)
  78     totwork += pow(n-ni,3);
  79
  80   work=0;
  81   pr_pr(0.,0.,0.);
  82   for(ni=1; (ni < n); ni++) {
  83     i=n-ni;
  84     l=i-1;
  85     h=0.0;
  86     g=x[i+n*(i-1)];
  87     if (l > 0) {
  88       for(k=0; (k<l); k++)
  89         h=h+sqr(x[i+n*k]);
  90       s=h+g*g;
  91       if (s < tol)
  92         h=0.0;
  93       else if (h > 0.0) {
  94         l=l+1;
  95         f=g;
  96         g=sqrt(s);
  97         g=-g;
  98         h=s-f*g;
  99         x[i+n*(i-1)]=f-g;
 100         f=0.0;
 101         for(j=0; (j<l); j++) {
 102           x[j+n*i]=x[i+n*j]/h;
 103           s=0.0;
 104           for(k=0; (k<=j); k++)
 105             s=s+x[j+n*k]*x[i+n*k];
 106           for(k=j+1; (k<l); k++)
 107             s=s+x[k+n*j]*x[i+n*k];
 108           e[j]=s/h;
 109           f=f+s*x[j+n*i];
 110         }
 111         f=f/(h+h);
 112         for(j=0; (j<l); j++)
 113           e[j]=e[j]-f*x[i+n*j];
 114         for(j=0; (j<l); j++) {
 115           f=x[i+n*j];
 116           s=e[j];
 117           for(k=0; (k<=j); k++)
 118             x[j+n*k]=x[j+n*k]-f*e[k]-x[i+n*k]*s;
 119         }
 120       }
 121     }
 122     d[i]=h;
 123     e[i-1]=g;
 124
 125     work += pow(n-ni,3);
 126     if (ni % 5 == 0)
 127       pr_pr(floor(100*work/totwork+0.5),0.,0.);
 128   }
 129
 130   /*
 131    *  accumulation of transformation matrix and intermediate d vector
 132    */
 133
 134   d[0]=x[0];
 135   x[0]=1.0;
 136
 137   work=0;
 138   pr_pr(100.,0.,0.);
 139   for(i=1; (i<n); i++) {
 140     if (d[i] > 0.0) {
 141       for(j=0; (j<i); j++) {
 142         s=0.0;
 143         for(k=0; (k<i); k++)
 144           s=s+x[i+n*k]*x[k+n*j];
 145         for(k=0; (k<i); k++)
 146           x[k+n*j]=x[k+n*j]-s*x[k+n*i];
 147       }
 148     }
 149     d[i]=x[i+n*i];
 150     x[i+n*i]=1.0;
 151     for(j=0; (j<i); j++) {
 152       x[i+n*j]=0.0;
 153       x[j+n*i]=0.0;
 154     }
 155     work += pow(i,3);
 156     if (i % 5 == 0)
 157       pr_pr(100.,floor(100*work/totwork+0.5),0.);
 158   }
 159
 160   /*
 161    *  ql iterates
 162    */
 163
 164   b=0.0;
 165   f=0.0;
 166   e[n-1]=0.0;
 167   totwork += pow(n,3);
 168   work=0;
 169   pr_pr(100.,100.,0.);
 170   for(l=0; (l<n); l++) {
 171     h=eps*(fabs(d[l])+fabs(e[l]));
 172     if (h > b)
 173       b=h;
 174     for(j=l; (j<n); j++) {
 175       if(fabs(e[j]) <= b)
 176         break;
 177     }
 178     if (j != l) {
 179       do {
 180         g=d[l];
 181         p=(d[l+1]-g)*0.5/e[l];
 182         r=sqrt(p*p+1.0);
 183         if(p < 0.0)
 184           p=p-r;
 185         else
 186           p=p+r;
 187         d[l]=e[l]/p;
 188         h=g-d[l];
 189         for(i=l+1; (i<n); i++)
 190           d[i]=d[i]-h;
 191         f=f+h;
 192         p=d[j];
 193         c=1.0;
 194         s=0.0;
 195         for(ni=l; (ni<j); ni++) {
 196           i=l+j-1-ni;
 197           g=c*e[i];
 198           h=c*p;
 199           if(fabs(p) >= fabs(e[i])) {
 200             c=e[i]/p;
 201             r=sqrt(c*c+1.0);
 202             e[i+1]=s*p*r;
 203             s=c/r;
 204             c=1.0/r;
 205           } else {
 206             c=p/e[i];
 207             r=sqrt(c*c+1.0);
 208             e[i+1]=s*e[i]*r;
 209             s=1.0/r;
 210             c=c/r;
 211           }
 212           p=c*d[i]-s*g;
 213           d[i+1]=h+s*(c*g+s*d[i]);
 214           for(k=0; (k<n); k++) {
 215             h=x[k+n*(i+1)];
 216             x[k+n*(i+1)]=x[k+n*i]*s+h*c;
 217             x[k+n*i]=x[k+n*i]*c-h*s;
 218           }
 219         }
 220         e[l]=s*p;
 221         d[l]=c*p;
 222       } while (fabs(e[l]) > b);
 223     }
 224     d[l]=d[l]+f;
 225
 226     work += pow(n-l,3);
 227     if (l % 5 == 0)
 228       pr_pr(100.,100.,floor(100*work/totwork+0.5));
 229   }
 230   fprintf(stderr,"\n");
 231
 232   /*
 233    *  put eigenvalues and eigenvectors in
 234    *  desired ascending order
 235    */
 236
 237   for(i=0; (i<n-1); i++) {
 238     k    = i;
 239     p    = d[i];
 240     absp = fabs(d[i]);
 241     for(j=i+1; (j<n); j++) {
 242       if(fabs(d[j]) < absp) {
 243         k    = j;
 244         p    = d[j];
 245         absp = fabs(d[j]);
 246       }
 247     }
 248     if (k != i) {
 249       d[k] = d[i];
 250       d[i] = p;
 251       for(j=0; (j<n); j++) {
 252         p        = x[j+n*i];
 253         x[j+n*i] = x[j+n*k];
 254         x[j+n*k] = p;
 255       }
 256     }
 257   }
 258
 259   sfree(e);
 260
 261   /*c
 262     c     fmp
 263     c     last but not least i have to confess that there is an original
 264     c     g. binsch remark on this routine: 'ql is purported to be one of
 265     c     the fastest and most compact routines of its kind presently known
 266     c     to mankind.'
 267     c
 268
 269     SIC! DvdS
 270
 271     */
 272
 273 }