src/gromacs/linearalgebra/matrix.cpp

   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,2015,2017,2019, 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 "gmxpre.h"
  38
  39 #include "matrix.h"
  40
  41 #include "config.h"
  42
  43 #include <stdio.h>
  44
  45 #include "gromacs/utility/fatalerror.h"
  46 #include "gromacs/utility/smalloc.h"
  47
  48 #include "gmx_lapack.h"
  49
  50 double** alloc_matrix(int n, int m)
  51 {
  52     double** ptr;
  53     int      i;
  54
  55     /* There's always time for more pointer arithmetic! */
  56     /* This is necessary in order to be able to work with LAPACK */
  57     snew(ptr, n);
  58     snew(ptr[0], n * m);
  59     for (i = 1; (i < n); i++)
  60     {
  61         ptr[i] = ptr[i - 1] + m;
  62     }
  63     return ptr;
  64 }
  65
  66 void free_matrix(double** a)
  67 {
  68     sfree(a[0]);
  69     sfree(a);
  70 }
  71
  72 #define DEBUG_MATRIX
  73 void matrix_multiply(FILE* fp, int n, int m, double** x, double** y, double** z)
  74 {
  75     int i, j, k;
  76
  77 #ifdef DEBUG_MATRIX
  78     if (fp)
  79     {
  80         fprintf(fp, "Multiplying %d x %d matrix with a %d x %d matrix\n", n, m, m, n);
  81     }
  82     if (fp)
  83     {
  84         for (i = 0; (i < n); i++)
  85         {
  86             for (j = 0; (j < m); j++)
  87             {
  88                 fprintf(fp, " %7g", x[i][j]);
  89             }
  90             fprintf(fp, "\n");
  91         }
  92     }
  93 #endif
  94     for (i = 0; (i < m); i++)
  95     {
  96         for (j = 0; (j < m); j++)
  97         {
  98             z[i][j] = 0;
  99             for (k = 0; (k < n); k++)
 100             {
 101                 z[i][j] += x[k][i] * y[j][k];
 102             }
 103         }
 104     }
 105 }
 106
 107 static void dump_matrix(FILE* fp, const char* title, int n, double** a)
 108 {
 109     double d = 1;
 110     int    i, j;
 111
 112     fprintf(fp, "%s\n", title);
 113     for (i = 0; (i < n); i++)
 114     {
 115         d = d * a[i][i];
 116         for (j = 0; (j < n); j++)
 117         {
 118             fprintf(fp, " %8.2f", a[i][j]);
 119         }
 120         fprintf(fp, "\n");
 121     }
 122     fprintf(fp, "Prod a[i][i] = %g\n", d);
 123 }
 124
 125 int matrix_invert(FILE* fp, int n, double** a)
 126 {
 127     int      i, j, m, lda, *ipiv, lwork, info;
 128     double **test = nullptr, **id, *work;
 129
 130 #ifdef DEBUG_MATRIX
 131     if (fp)
 132     {
 133         fprintf(fp, "Inverting %d square matrix\n", n);
 134         test = alloc_matrix(n, n);
 135         for (i = 0; (i < n); i++)
 136         {
 137             for (j = 0; (j < n); j++)
 138             {
 139                 test[i][j] = a[i][j];
 140             }
 141         }
 142         dump_matrix(fp, "before inversion", n, a);
 143     }
 144 #endif
 145     snew(ipiv, n);
 146     lwork = n * n;
 147     snew(work, lwork);
 148     m = lda = n;
 149     info    = 0;
 150     F77_FUNC(dgetrf, DGETRF)(&n, &m, a[0], &lda, ipiv, &info);
 151 #ifdef DEBUG_MATRIX
 152     if (fp)
 153     {
 154         dump_matrix(fp, "after dgetrf", n, a);
 155     }
 156 #endif
 157     if (info != 0)
 158     {
 159         return info;
 160     }
 161     F77_FUNC(dgetri, DGETRI)(&n, a[0], &lda, ipiv, work, &lwork, &info);
 162 #ifdef DEBUG_MATRIX
 163     if (fp)
 164     {
 165         dump_matrix(fp, "after dgetri", n, a);
 166     }
 167 #endif
 168     if (info != 0)
 169     {
 170         return info;
 171     }
 172
 173 #ifdef DEBUG_MATRIX
 174     if (fp)
 175     {
 176         id = alloc_matrix(n, n);
 177         matrix_multiply(fp, n, n, test, a, id);
 178         dump_matrix(fp, "And here is the product of A and Ainv", n, id);
 179         free_matrix(id);
 180         free_matrix(test);
 181     }
 182 #endif
 183     sfree(ipiv);
 184     sfree(work);
 185
 186     return 0;
 187 }
 188
 189 double multi_regression(FILE* fp, int nrow, double* y, int ncol, double** xx, double* a0)
 190 {
 191     int    row, i, j;
 192     double ax, chi2, **a, **at, **ata, *atx;
 193
 194     a   = alloc_matrix(nrow, ncol);
 195     at  = alloc_matrix(ncol, nrow);
 196     ata = alloc_matrix(ncol, ncol);
 197     for (i = 0; (i < nrow); i++)
 198     {
 199         for (j = 0; (j < ncol); j++)
 200         {
 201             at[j][i] = a[i][j] = xx[j][i];
 202         }
 203     }
 204     matrix_multiply(fp, nrow, ncol, a, at, ata);
 205     if ((row = matrix_invert(fp, ncol, ata)) != 0)
 206     {
 207         gmx_fatal(
 208                 FARGS,
 209                 "Matrix inversion failed. Incorrect row = %d.\nThis probably indicates that you do "
 210                 "not have sufficient data points, or that some parameters are linearly dependent.",
 211                 row);
 212     }
 213     snew(atx, ncol);
 214
 215     for (i = 0; (i < ncol); i++)
 216     {
 217         atx[i] = 0;
 218         for (j = 0; (j < nrow); j++)
 219         {
 220             atx[i] += at[i][j] * y[j];
 221         }
 222     }
 223     for (i = 0; (i < ncol); i++)
 224     {
 225         a0[i] = 0;
 226         for (j = 0; (j < ncol); j++)
 227         {
 228             a0[i] += ata[i][j] * atx[j];
 229         }
 230     }
 231     chi2 = 0;
 232     for (j = 0; (j < nrow); j++)
 233     {
 234         ax = 0;
 235         for (i = 0; (i < ncol); i++)
 236         {
 237             ax += a0[i] * a[j][i];
 238         }
 239         chi2 += (y[j] - ax) * (y[j] - ax);
 240     }
 241
 242     sfree(atx);
 243     free_matrix(a);
 244     free_matrix(at);
 245     free_matrix(ata);
 246
 247     return chi2;
 248 }