src/gmxlib/gmx_matrix.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 4.0.99
  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-2008, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  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.
  19  *
  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.
  26  *
  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.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Groningen Machine for Chemical Simulation
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <stdlib.h>
  40 #include "typedefs.h"
  41 #include "smalloc.h"
  42 #include "vec.h"
  43 #include "gmx_fatal.h"
  44 #include "gmx_matrix.h"
  45 #include "gmx_lapack.h"
  46
  47 double **alloc_matrix(int n,int m)
  48 {
  49     double **ptr;
  50     int i;
  51
  52     /* There's always time for more pointer arithmetic! */
  53     /* This is necessary in order to be able to work with LAPACK */
  54     snew(ptr,n);
  55     snew(ptr[0],n*m);
  56     for(i=1; (i<n); i++)
  57         ptr[i] = ptr[i-1]+m;
  58     return ptr;
  59 }
  60
  61 void free_matrix(double **a,int n)
  62 {
  63     int i;
  64
  65     sfree(a[0]);
  66     sfree(a);
  67 }
  68
  69 #define DEBUG_MATRIX
  70 void matrix_multiply(FILE *fp,int n,int m,double **x,double **y,double **z)
  71 {
  72     int i,j,k;
  73
  74 #ifdef DEBUG_MATRIX
  75     if (fp)
  76         fprintf(fp,"Multiplying %d x %d matrix with a %d x %d matrix\n",
  77                 n,m,m,n);
  78     if (fp)
  79         for(i=0; (i<n); i++)
  80         {
  81             for(j=0; (j<m); j++)
  82                 fprintf(fp," %7g",x[i][j]);
  83             fprintf(fp,"\n");
  84         }
  85 #endif
  86     for(i=0; (i<m); i++)
  87     {
  88         for(j=0; (j<m); j++)
  89         {
  90             z[i][j] = 0;
  91             for(k=0; (k<n); k++)
  92                 z[i][j] += x[k][i]*y[j][k];
  93         }
  94     }
  95 }
  96
  97 static void dump_matrix(FILE *fp,const char *title,int n,double **a)
  98 {
  99     double d=1;
 100     int i,j;
 101
 102     fprintf(fp,"%s\n",title);
 103     for(i=0; (i<n); i++)
 104     {
 105         d = d*a[i][i];
 106         for(j=0; (j<n); j++)
 107         {
 108             fprintf(fp," %8.2f",a[i][j]);
 109         }
 110         fprintf(fp,"\n");
 111     }
 112     fprintf(fp,"Prod a[i][i] = %g\n",d);
 113 }
 114
 115 int matrix_invert(FILE *fp,int n,double **a)
 116 {
 117     int i,j,m,lda,*ipiv,lwork,info;
 118     double **test=NULL,**id,*work;
 119
 120 #ifdef DEBUG_MATRIX
 121     if (fp)
 122     {
 123         fprintf(fp,"Inverting %d square matrix\n",n);
 124         test = alloc_matrix(n,n);
 125         for(i=0; (i<n); i++)
 126         {
 127             for(j=0; (j<n); j++)
 128             {
 129                 test[i][j] = a[i][j];
 130             }
 131         }
 132         dump_matrix(fp,"before inversion",n,a);
 133     }
 134 #endif
 135     snew(ipiv,n);
 136     lwork = n*n;
 137     snew(work,lwork);
 138     m = lda   = n;
 139     info  = 0;
 140     F77_FUNC(dgetrf,DGETRF)(&n,&m,a[0],&lda,ipiv,&info);
 141 #ifdef DEBUG_MATRIX
 142     if (fp)
 143         dump_matrix(fp,"after dgetrf",n,a);
 144 #endif
 145     if (info != 0)
 146         return info;
 147     F77_FUNC(dgetri,DGETRI)(&n,a[0],&lda,ipiv,work,&lwork,&info);
 148 #ifdef DEBUG_MATRIX
 149     if (fp)
 150         dump_matrix(fp,"after dgetri",n,a);
 151 #endif
 152     if (info != 0)
 153         return info;
 154
 155 #ifdef DEBUG_MATRIX
 156     if (fp)
 157     {
 158         id = alloc_matrix(n,n);
 159         matrix_multiply(fp,n,n,test,a,id);
 160         dump_matrix(fp,"And here is the product of A and Ainv",n,id);
 161         free_matrix(id,n);
 162         free_matrix(test,n);
 163     }
 164 #endif
 165     sfree(ipiv);
 166     sfree(work);
 167
 168     return 0;
 169 }
 170
 171 double multi_regression(FILE *fp,int nrow,double *y,int ncol,
 172                         double **xx,double *a0)
 173 {
 174   int    row,niter,i,j;
 175   double ax,chi2,**a,**at,**ata,*atx;
 176
 177   a   = alloc_matrix(nrow,ncol);
 178   at  = alloc_matrix(ncol,nrow);
 179   ata = alloc_matrix(ncol,ncol);
 180   for(i=0; (i<nrow); i++)
 181       for(j=0; (j<ncol); j++)
 182           at[j][i] = a[i][j] = xx[j][i];
 183   matrix_multiply(fp,nrow,ncol,a,at,ata);
 184   if ((row = matrix_invert(fp,ncol,ata)) != 0) {
 185       gmx_fatal(FARGS,"Matrix inversion failed. Incorrect row = %d.\nThis probably indicates that you do not have sufficient data points, or that some parameters are linearly dependent.",
 186                 row);
 187   }
 188   snew(atx,ncol);
 189
 190   for(i=0; (i<ncol); i++)
 191   {
 192       atx[i] = 0;
 193       for(j=0; (j<nrow); j++)
 194           atx[i] += at[i][j]*y[j];
 195   }
 196   for(i=0; (i<ncol); i++)
 197   {
 198       a0[i] = 0;
 199       for(j=0; (j<ncol); j++)
 200           a0[i] += ata[i][j]*atx[j];
 201   }
 202   chi2 = 0;
 203   for(j=0; (j<nrow); j++)
 204   {
 205       ax = 0;
 206       for(i=0; (i<ncol); i++)
 207           ax += a0[i]*a[j][i];
 208       chi2 += sqr(y[j]-ax);
 209   }
 210
 211   sfree(atx);
 212   free_matrix(a,nrow);
 213   free_matrix(at,ncol);
 214   free_matrix(ata,ncol);
 215
 216   return chi2;
 217 }
 218