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