src/gromacs/linearalgebra/sparsematrix.h

   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,2018,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 #ifndef GMX_LINEARALGEBRA_SPARSEMATRIX_H
  38 #define GMX_LINEARALGEBRA_SPARSEMATRIX_H
  39
  40 #include <stdio.h>
  41
  42 #include "gromacs/utility/basedefinitions.h"
  43 #include "gromacs/utility/real.h"
  44
  45 typedef struct gmx_sparsematrix_entry
  46 {
  47     int  col;
  48     real value;
  49 } gmx_sparsematrix_entry_t;
  50
  51 /*! \brief Sparse matrix storage format
  52  *
  53  *  This structure specifies a storage format for a sparse matrix.
  54  *  The memory requirements are only proportional to the number
  55  *  of nonzero elements, and it provides a reasonably fast way to
  56  *  perform matrix-vector multiplications.
  57  *
  58  *  The data format is very similar to a neighborlist. It is optimized
  59  *  for fast access, but it is difficult to add entries. If you are
  60  *  constructing a matrix you should either do it in exactly the order
  61  *  specified here, or use some other more flexible intermediate structure.
  62  *
  63  *  The index array is of size nrow+1. All non-zero matrix elements
  64  *  on row i are stored in positions index[i] through index[i+1]-1 in
  65  *  the arrays column and value. The column array contains the column
  66  *  index for each entry, in ascending order, and the corresponding
  67  *  position in the value array contains the floating point matrix element.
  68  *
  69  *  index[nrow] should be equal to the total number of elements stored.
  70  *
  71  *  Thus, to find the value of matrix element [5,4] you should loop
  72  *  over positions index[5] to index[6]-1 in column until you either find
  73  *  the value 4, or a higher value (meaning the element was zero).
  74  *
  75  *  It is fairly easy to construct the matrix on-the-fly if you can do
  76  *  it row-by-row.
  77  *
  78  *  IMPORTANT:
  79  *  If compressed_symmetric is set to TRUE, you should only store EITHER the upper OR
  80  *  lower triangle (and the diagonal), and the other half is assumed to be
  81  *  symmetric. Otherwise, if compressed_symmetric==FALSE, no symmetry is implied and all
  82  *  elements should be stored.
  83  *
  84  *  The symmetry compression saves us a factor 2 both in storage and
  85  *  matrix multiplication CPU-time, which can be very useful for huge eigenproblems.
  86  *
  87  *  If you are unsure, just set compressed_symmetric to FALSE and list all elements. If
  88  *  you enable it but still list all elements (both upper and lower triangle) you will be sorry...
  89  *
  90  *  Internally, the sparse data is stored as a separate list for each row, where the list
  91  *  element is a structure with a column and (floating-point) data value. This makes it
  92  *  possible, although not completely transparent, to update values in random access order.
  93  *  The drawback is that the structure will allocate nrow memory regions.
  94  *  The matrix data could be stored in a single contiguous array with indices for each row,
  95  *  but then we could only insert elements at the end without copying the entire matrix.
  96  *
  97  *  After you have
  98  *
  99  *  In other words: Not perfect, but it works.
 100  */
 101 typedef struct gmx_sparsematrix
 102 {
 103     gmx_bool compressed_symmetric;   /**< Store half elements and assume symmetry. */
 104     int      nrow;                   /**< Number of rows in matrix                 */
 105     int*     ndata;                  /**< Number of entries on each row (list)     */
 106     int*     nalloc;                 /**< Allocated entry list length for each row */
 107     gmx_sparsematrix_entry_t** data; /**< data[i] is a list with entries on row i  */
 108 } gmx_sparsematrix_t;
 109
 110
 111 /*! \brief Allocate a new sparse matrix structure
 112  *
 113  *  The number of rows is used to allocate the index array entry. Obviously you
 114  *  can reallocate these later yourself if necessary - this is a
 115  *  convenience routine.
 116  *
 117  *  By default, the compressed_symmetric flag in the structure will
 118  *  be FALSE. Set it to TRUE manually if you are only storing either the
 119  *  upper or lower half of the matrix.
 120  */
 121 gmx_sparsematrix_t* gmx_sparsematrix_init(int nrow);
 122
 123
 124 /*! \brief Release all resources used by a sparse matrix structure
 125  *
 126  *  All arrays in the structure will be freed, and the structure itself.
 127  */
 128 void gmx_sparsematrix_destroy(gmx_sparsematrix_t* A);
 129
 130
 131 /*! \brief Print sparse matrix to a stream.
 132  *
 133  *  Mainly used for debugging. Be warned that the real sparse matrices used
 134  *  in Gromacs runs can be HUGE (think 100,000 rows).
 135  */
 136 void gmx_sparsematrix_print(FILE* stream, gmx_sparsematrix_t* A);
 137
 138 /* Adds value at row,col. If the value did not exist
 139  * previously it is added, otherwise it is incremented with difference.
 140  *
 141  * The column sort order might change, so you need to run fix_sparsematrix
 142  * once you are done changing the matrix.
 143  */
 144 real gmx_sparsematrix_value(gmx_sparsematrix_t* A, int row, int col);
 145
 146
 147 /* Adds value at row,col. If the value did not exist
 148  * previously it is added, otherwise it is incremented with difference.
 149  *
 150  * The column sort order might change, so you need to run fix_sparsematrix
 151  * once you are done changing the matrix.
 152  */
 153 void gmx_sparsematrix_increment_value(gmx_sparsematrix_t* A, int row, int col, real difference);
 154
 155
 156 /*! \brief Sort elements in each column and remove zeros.
 157  *
 158  *  Sparse matrix access is faster when the elements are stored in
 159  *  increasing column order in each row. In some cases previously non-zero
 160  *  elements will be zero after adding more data, and this routine also removes
 161  *  those entries to reduce the storage requirements.
 162  *
 163  *  It never hurts to run this routine if you have been updating the matrix...
 164  */
 165 void gmx_sparsematrix_compress(gmx_sparsematrix_t* A);
 166
 167
 168 /*! \brief Sparse matrix vector multiplication
 169  *
 170  * Calculate y = A * x for a sparse matrix A.
 171  */
 172 void gmx_sparsematrix_vector_multiply(gmx_sparsematrix_t* A, real* x, real* y);
 173
 174
 175 #endif