src/gromacs/linearalgebra/matrix.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-2008, The GROMACS development team.
   6  * Copyright (c) 2012,2013,2014, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "matrix.h"
  38
  39 #ifdef HAVE_CONFIG_H
  40 #include <config.h>
  41 #endif
  42
  43 #include <stdio.h>
  44
  45 #include "gromacs/legacyheaders/gmx_fatal.h"
  46 #include "gromacs/legacyheaders/vec.h"
  47
  48 #include "gromacs/utility/smalloc.h"
  49
  50 #include "gmx_lapack.h"
  51
  52 double **alloc_matrix(int n, int m)
  53 {
  54     double **ptr;
  55     int      i;
  56
  57     /* There's always time for more pointer arithmetic! */
  58     /* This is necessary in order to be able to work with LAPACK */
  59     snew(ptr, n);
  60     snew(ptr[0], n*m);
  61     for (i = 1; (i < n); i++)
  62     {
  63         ptr[i] = ptr[i-1]+m;
  64     }
  65     return ptr;
  66 }
  67
  68 void free_matrix(double **a)
  69 {
  70     int i;
  71
  72     sfree(a[0]);
  73     sfree(a);
  74 }
  75
  76 #define DEBUG_MATRIX
  77 void matrix_multiply(FILE *fp, int n, int m, double **x, double **y, double **z)
  78 {
  79     int i, j, k;
  80
  81 #ifdef DEBUG_MATRIX
  82     if (fp)
  83     {
  84         fprintf(fp, "Multiplying %d x %d matrix with a %d x %d matrix\n",
  85                 n, m, m, n);
  86     }
  87     if (fp)
  88     {
  89         for (i = 0; (i < n); i++)
  90         {
  91             for (j = 0; (j < m); j++)
  92             {
  93                 fprintf(fp, " %7g", x[i][j]);
  94             }
  95             fprintf(fp, "\n");
  96         }
  97     }
  98 #endif
  99     for (i = 0; (i < m); i++)
 100     {
 101         for (j = 0; (j < m); j++)
 102         {
 103             z[i][j] = 0;
 104             for (k = 0; (k < n); k++)
 105             {
 106                 z[i][j] += x[k][i]*y[j][k];
 107             }
 108         }
 109     }
 110 }
 111
 112 static void dump_matrix(FILE *fp, const char *title, int n, double **a)
 113 {
 114     double d = 1;
 115     int    i, j;
 116
 117     fprintf(fp, "%s\n", title);
 118     for (i = 0; (i < n); i++)
 119     {
 120         d = d*a[i][i];
 121         for (j = 0; (j < n); j++)
 122         {
 123             fprintf(fp, " %8.2f", a[i][j]);
 124         }
 125         fprintf(fp, "\n");
 126     }
 127     fprintf(fp, "Prod a[i][i] = %g\n", d);
 128 }
 129
 130 int matrix_invert(FILE *fp, int n, double **a)
 131 {
 132     int      i, j, m, lda, *ipiv, lwork, info;
 133     double **test = NULL, **id, *work;
 134
 135 #ifdef DEBUG_MATRIX
 136     if (fp)
 137     {
 138         fprintf(fp, "Inverting %d square matrix\n", n);
 139         test = alloc_matrix(n, n);
 140         for (i = 0; (i < n); i++)
 141         {
 142             for (j = 0; (j < n); j++)
 143             {
 144                 test[i][j] = a[i][j];
 145             }
 146         }
 147         dump_matrix(fp, "before inversion", n, a);
 148     }
 149 #endif
 150     snew(ipiv, n);
 151     lwork = n*n;
 152     snew(work, lwork);
 153     m     = lda   = n;
 154     info  = 0;
 155     F77_FUNC(dgetrf, DGETRF) (&n, &m, a[0], &lda, ipiv, &info);
 156 #ifdef DEBUG_MATRIX
 157     if (fp)
 158     {
 159         dump_matrix(fp, "after dgetrf", n, a);
 160     }
 161 #endif
 162     if (info != 0)
 163     {
 164         return info;
 165     }
 166     F77_FUNC(dgetri, DGETRI) (&n, a[0], &lda, ipiv, work, &lwork, &info);
 167 #ifdef DEBUG_MATRIX
 168     if (fp)
 169     {
 170         dump_matrix(fp, "after dgetri", n, a);
 171     }
 172 #endif
 173     if (info != 0)
 174     {
 175         return info;
 176     }
 177
 178 #ifdef DEBUG_MATRIX
 179     if (fp)
 180     {
 181         id = alloc_matrix(n, n);
 182         matrix_multiply(fp, n, n, test, a, id);
 183         dump_matrix(fp, "And here is the product of A and Ainv", n, id);
 184         free_matrix(id);
 185         free_matrix(test);
 186     }
 187 #endif
 188     sfree(ipiv);
 189     sfree(work);
 190
 191     return 0;
 192 }
 193
 194 double multi_regression(FILE *fp, int nrow, double *y, int ncol,
 195                         double **xx, double *a0)
 196 {
 197     int    row, niter, i, j;
 198     double ax, chi2, **a, **at, **ata, *atx;
 199
 200     a   = alloc_matrix(nrow, ncol);
 201     at  = alloc_matrix(ncol, nrow);
 202     ata = alloc_matrix(ncol, ncol);
 203     for (i = 0; (i < nrow); i++)
 204     {
 205         for (j = 0; (j < ncol); j++)
 206         {
 207             at[j][i] = a[i][j] = xx[j][i];
 208         }
 209     }
 210     matrix_multiply(fp, nrow, ncol, a, at, ata);
 211     if ((row = matrix_invert(fp, ncol, ata)) != 0)
 212     {
 213         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.",
 214                   row);
 215     }
 216     snew(atx, ncol);
 217
 218     for (i = 0; (i < ncol); i++)
 219     {
 220         atx[i] = 0;
 221         for (j = 0; (j < nrow); j++)
 222         {
 223             atx[i] += at[i][j]*y[j];
 224         }
 225     }
 226     for (i = 0; (i < ncol); i++)
 227     {
 228         a0[i] = 0;
 229         for (j = 0; (j < ncol); j++)
 230         {
 231             a0[i] += ata[i][j]*atx[j];
 232         }
 233     }
 234     chi2 = 0;
 235     for (j = 0; (j < nrow); j++)
 236     {
 237         ax = 0;
 238         for (i = 0; (i < ncol); i++)
 239         {
 240             ax += a0[i]*a[j][i];
 241         }
 242         chi2 += sqr(y[j]-ax);
 243     }
 244
 245     sfree(atx);
 246     free_matrix(a);
 247     free_matrix(at);
 248     free_matrix(ata);
 249
 250     return chi2;
 251 }