src/gmxlib/gmx_matrix.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  * $Id: gmx_matrix.c,v 1.4 2008/12/02 18:27:57 spoel Exp $
   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 4.0.99
  11  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  12  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  13  * Copyright (c) 2001-2008, The GROMACS development team,
  14  * check out http://www.gromacs.org for more information.
  15
  16  * This program is free software; you can redistribute it and/or
  17  * modify it under the terms of the GNU General Public License
  18  * as published by the Free Software Foundation; either version 2
  19  * of the License, or (at your option) any later version.
  20  *
  21  * If you want to redistribute modifications, please consider that
  22  * scientific software is very special. Version control is crucial -
  23  * bugs must be traceable. We will be happy to consider code for
  24  * inclusion in the official distribution, but derived work must not
  25  * be called official GROMACS. Details are found in the README & COPYING
  26  * files - if they are missing, get the official version at www.gromacs.org.
  27  *
  28  * To help us fund GROMACS development, we humbly ask that you cite
  29  * the papers on the package - you can find them in the top README file.
  30  *
  31  * For more info, check our website at http://www.gromacs.org
  32  *
  33  * And Hey:
  34  * Groningen Machine for Chemical Simulation
  35  */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include <stdlib.h>
  41 #include "typedefs.h"
  42 #include "smalloc.h"
  43 #include "vec.h"
  44 #include "gmx_fatal.h"
  45 #include "gmx_matrix.h"
  46 #include "gmx_lapack.h"
  47
  48 double **alloc_matrix(int n,int m)
  49 {
  50     double **ptr;
  51     int i;
  52
  53     /* There's always time for more pointer arithmetic! */
  54     /* This is necessary in order to be able to work with LAPACK */
  55     snew(ptr,n);
  56     snew(ptr[0],n*m);
  57     for(i=1; (i<n); i++)
  58         ptr[i] = ptr[i-1]+m;
  59     return ptr;
  60 }
  61
  62 void free_matrix(double **a,int n)
  63 {
  64     int i;
  65
  66     sfree(a[0]);
  67     sfree(a);
  68 }
  69
  70 #define DEBUG_MATRIX
  71 void matrix_multiply(FILE *fp,int n,int m,double **x,double **y,double **z)
  72 {
  73     int i,j,k;
  74
  75 #ifdef DEBUG_MATRIX
  76     if (fp)
  77         fprintf(fp,"Multiplying %d x %d matrix with a %d x %d matrix\n",
  78                 n,m,m,n);
  79     if (fp)
  80         for(i=0; (i<n); i++)
  81         {
  82             for(j=0; (j<m); j++)
  83                 fprintf(fp," %7g",x[i][j]);
  84             fprintf(fp,"\n");
  85         }
  86 #endif
  87     for(i=0; (i<m); i++)
  88     {
  89         for(j=0; (j<m); j++)
  90         {
  91             z[i][j] = 0;
  92             for(k=0; (k<n); k++)
  93                 z[i][j] += x[k][i]*y[j][k];
  94         }
  95     }
  96 }
  97
  98 static void dump_matrix(FILE *fp,const char *title,int n,double **a)
  99 {
 100     double d=1;
 101     int i,j;
 102
 103     fprintf(fp,"%s\n",title);
 104     for(i=0; (i<n); i++)
 105     {
 106         d = d*a[i][i];
 107         for(j=0; (j<n); j++)
 108         {
 109             fprintf(fp," %8.2f",a[i][j]);
 110         }
 111         fprintf(fp,"\n");
 112     }
 113     fprintf(fp,"Prod a[i][i] = %g\n",d);
 114 }
 115
 116 int matrix_invert(FILE *fp,int n,double **a)
 117 {
 118     int i,j,m,lda,*ipiv,lwork,info;
 119     double **test=NULL,**id,*work;
 120
 121 #ifdef DEBUG_MATRIX
 122     if (fp)
 123     {
 124         fprintf(fp,"Inverting %d square matrix\n",n);
 125         test = alloc_matrix(n,n);
 126         for(i=0; (i<n); i++)
 127         {
 128             for(j=0; (j<n); j++)
 129             {
 130                 test[i][j] = a[i][j];
 131             }
 132         }
 133         dump_matrix(fp,"before inversion",n,a);
 134     }
 135 #endif
 136     snew(ipiv,n);
 137     lwork = n*n;
 138     snew(work,lwork);
 139     m = lda   = n;
 140     info  = 0;
 141     F77_FUNC(dgetrf,DGETRF)(&n,&m,a[0],&lda,ipiv,&info);
 142 #ifdef DEBUG_MATRIX
 143     if (fp)
 144         dump_matrix(fp,"after dgetrf",n,a);
 145 #endif
 146     if (info != 0)
 147         return info;
 148     F77_FUNC(dgetri,DGETRI)(&n,a[0],&lda,ipiv,work,&lwork,&info);
 149 #ifdef DEBUG_MATRIX
 150     if (fp)
 151         dump_matrix(fp,"after dgetri",n,a);
 152 #endif
 153     if (info != 0)
 154         return info;
 155
 156 #ifdef DEBUG_MATRIX
 157     if (fp)
 158     {
 159         id = alloc_matrix(n,n);
 160         matrix_multiply(fp,n,n,test,a,id);
 161         dump_matrix(fp,"And here is the product of A and Ainv",n,id);
 162         free_matrix(id,n);
 163         free_matrix(test,n);
 164     }
 165 #endif
 166     sfree(ipiv);
 167     sfree(work);
 168
 169     return 0;
 170 }
 171
 172 double multi_regression(FILE *fp,int nrow,double *y,int ncol,
 173                         double **xx,double *a0)
 174 {
 175   int    row,niter,i,j;
 176   double ax,chi2,**a,**at,**ata,*atx;
 177
 178   a   = alloc_matrix(nrow,ncol);
 179   at  = alloc_matrix(ncol,nrow);
 180   ata = alloc_matrix(ncol,ncol);
 181   for(i=0; (i<nrow); i++)
 182       for(j=0; (j<ncol); j++)
 183           at[j][i] = a[i][j] = xx[j][i];
 184   matrix_multiply(fp,nrow,ncol,a,at,ata);
 185   if ((row = matrix_invert(fp,ncol,ata)) != 0) {
 186       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.",
 187                 row);
 188   }
 189   snew(atx,ncol);
 190
 191   for(i=0; (i<ncol); i++)
 192   {
 193       atx[i] = 0;
 194       for(j=0; (j<nrow); j++)
 195           atx[i] += at[i][j]*y[j];
 196   }
 197   for(i=0; (i<ncol); i++)
 198   {
 199       a0[i] = 0;
 200       for(j=0; (j<ncol); j++)
 201           a0[i] += ata[i][j]*atx[j];
 202   }
 203   chi2 = 0;
 204   for(j=0; (j<nrow); j++)
 205   {
 206       ax = 0;
 207       for(i=0; (i<ncol); i++)
 208           ax += a0[i]*a[j][i];
 209       chi2 += sqr(y[j]-ax);
 210   }
 211
 212   sfree(atx);
 213   free_matrix(a,nrow);
 214   free_matrix(at,ncol);
 215   free_matrix(ata,ncol);
 216
 217   return chi2;
 218 }
 219