src/gromacs/linearalgebra/sparsematrix.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-2004, The GROMACS development team.
   6  * Copyright (c) 2012,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 "sparsematrix.h"
  40
  41 #include <stdio.h>
  42 #include <stdlib.h>
  43
  44 #include <cassert>
  45
  46 #include "gromacs/utility/smalloc.h"
  47
  48 gmx_sparsematrix_t* gmx_sparsematrix_init(int nrow)
  49 {
  50     int                 i;
  51     gmx_sparsematrix_t* A;
  52
  53     snew(A, 1);
  54
  55     A->nrow = nrow;
  56     snew(A->ndata, nrow);
  57     snew(A->nalloc, nrow);
  58     snew(A->data, nrow);
  59
  60     for (i = 0; i < nrow; i++)
  61     {
  62         A->ndata[i]  = 0;
  63         A->nalloc[i] = 0;
  64         A->data[i]   = nullptr;
  65     }
  66     return A;
  67 }
  68
  69
  70 void gmx_sparsematrix_destroy(gmx_sparsematrix_t* A)
  71 {
  72     int i;
  73
  74     /* Release each row */
  75     for (i = 0; i < A->nrow; i++)
  76     {
  77         if (A->data[i] != nullptr)
  78         {
  79             sfree(A->data[i]);
  80         }
  81     }
  82     /* Release the rowdata arrays */
  83     sfree(A->ndata);
  84     sfree(A->nalloc);
  85     sfree(A->data);
  86     /* Release matrix structure itself */
  87     sfree(A);
  88 }
  89
  90
  91 void gmx_sparsematrix_print(FILE* stream, gmx_sparsematrix_t* A)
  92 {
  93     int i, j, k;
  94
  95     for (i = 0; i < A->nrow; i++)
  96     {
  97         if (A->ndata[i] == 0)
  98         {
  99             for (j = 0; j < A->nrow; j++)
 100             {
 101                 fprintf(stream, " %6.3f", 0.0);
 102             }
 103         }
 104         else
 105         {
 106             k = 0;
 107             for (j = 0; j < A->ndata[i]; j++)
 108             {
 109                 while (k++ < A->data[i][j].col)
 110                 {
 111                     fprintf(stream, " %6.3f", 0.0);
 112                 }
 113                 fprintf(stream, " %6.3f", A->data[i][j].value);
 114             }
 115             while (k++ < A->nrow)
 116             {
 117                 fprintf(stream, " %6.3f", 0.0);
 118             }
 119         }
 120         fprintf(stream, "\n");
 121     }
 122 }
 123
 124
 125 real gmx_sparsematrix_value(gmx_sparsematrix_t* A, int row, int col)
 126 {
 127     gmx_bool found = FALSE;
 128     int      i;
 129     real     value;
 130
 131     assert(row < A->nrow);
 132
 133     value = 0;
 134
 135     /* Find previous value */
 136     for (i = 0; i < A->ndata[row] && (found == FALSE); i++)
 137     {
 138         if (A->data[row][i].col == col)
 139         {
 140             found = TRUE;
 141             value = A->data[row][i].value;
 142         }
 143     }
 144
 145     /* value=0 if we didn't find any match */
 146     return value;
 147 }
 148
 149
 150 void gmx_sparsematrix_increment_value(gmx_sparsematrix_t* A, int row, int col, real difference)
 151 {
 152     gmx_bool found = FALSE;
 153     int      i;
 154
 155     assert(row < A->nrow);
 156
 157     /* Try to find a previous entry with this row/col */
 158     for (i = 0; i < A->ndata[row] && !found; i++)
 159     {
 160         if (A->data[row][i].col == col)
 161         {
 162             found = TRUE;
 163             A->data[row][i].value += difference;
 164         }
 165     }
 166
 167     /* Add a new entry if nothing was found */
 168     if (!found)
 169     {
 170         /* add the value at the end of the row */
 171         if (A->ndata[row] == A->nalloc[row])
 172         {
 173             A->nalloc[row] += 100;
 174             if (A->data[row] == nullptr)
 175             {
 176                 snew(A->data[row], A->nalloc[row]);
 177             }
 178             else
 179             {
 180                 srenew(A->data[row], A->nalloc[row]);
 181             }
 182         }
 183         A->data[row][A->ndata[row]].col = col;
 184         /* Previous value was 0.0 */
 185         A->data[row][A->ndata[row]].value = difference;
 186         A->ndata[row]++;
 187     }
 188 }
 189
 190
 191 /* Routine to compare column values of two entries, used for quicksort of each row.
 192  *
 193  * The data entries to compare are of the type gmx_sparsematrix_entry_t, but quicksort
 194  * uses void pointers as arguments, so we cast them back internally.
 195  */
 196 static int compare_columns(const void* v1, const void* v2)
 197 {
 198     int c1 = ((gmx_sparsematrix_entry_t*)v1)->col;
 199     int c2 = ((gmx_sparsematrix_entry_t*)v2)->col;
 200
 201     if (c1 < c2)
 202     {
 203         return -1;
 204     }
 205     else if (c1 > c2)
 206     {
 207         return 1;
 208     }
 209     else
 210     {
 211         return 0;
 212     }
 213 }
 214
 215
 216 void gmx_sparsematrix_compress(gmx_sparsematrix_t* A)
 217 {
 218     int i, j;
 219
 220     for (i = 0; i < A->nrow; i++)
 221     {
 222         /* Remove last value on this row while it is zero */
 223         while (A->ndata[i] > 0 && A->data[i][A->ndata[i] - 1].value == 0)
 224         {
 225             A->ndata[i]--;
 226         }
 227
 228         /* Go through values on this row and look for more zero elements */
 229         for (j = 0; j < A->ndata[i]; j++)
 230         {
 231             /* If this element was zero, exchange it with the last non-zero
 232              * element on the row (yes, this will invalidate the sort order)
 233              */
 234             if (A->data[i][j].value == 0)
 235             {
 236                 A->data[i][j].value = A->data[i][A->ndata[i] - 1].value;
 237                 A->data[i][j].col   = A->data[i][A->ndata[i] - 1].col;
 238                 A->ndata[i]--;
 239             }
 240         }
 241         /* Only non-zero elements remaining on this row. Sort them after column index */
 242         qsort((void*)(A->data[i]), A->ndata[i], sizeof(gmx_sparsematrix_entry_t), compare_columns);
 243     }
 244 }
 245
 246
 247 void gmx_sparsematrix_vector_multiply(gmx_sparsematrix_t* A, real* x, real* y)
 248 {
 249     real                      s, v, xi;
 250     int                       i, j, k;
 251     gmx_sparsematrix_entry_t* data; /* pointer to simplify data access */
 252
 253     for (i = 0; i < A->nrow; i++)
 254     {
 255         y[i] = 0;
 256     }
 257
 258     if (A->compressed_symmetric)
 259     {
 260         for (i = 0; i < A->nrow; i++)
 261         {
 262             xi   = x[i];
 263             s    = 0.0;
 264             data = A->data[i];
 265
 266             for (k = 0; k < A->ndata[i]; k++)
 267             {
 268                 j = data[k].col;
 269                 v = data[k].value;
 270                 s += v * x[j];
 271                 if (i != j)
 272                 {
 273                     y[j] += v * xi;
 274                 }
 275             }
 276             y[i] += s;
 277         }
 278     }
 279     else
 280     {
 281         /* not compressed symmetric storage */
 282         for (i = 0; i < A->nrow; i++)
 283         {
 284             s    = 0.0;
 285             data = A->data[i];
 286
 287             for (k = 0; k < A->ndata[i]; k++)
 288             {
 289                 j = data[k].col;
 290                 v = data[k].value;
 291                 s += v * x[j];
 292             }
 293             y[i] += s;
 294         }
 295     }
 296 }