src/gromacs/linearalgebra/sparsematrix.h

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Green Red Orange Magenta Azure Cyan Skyblue
  34  */
  35 #ifndef GMX_LINEARALGEBRA_SPARSEMATRIX_H
  36 #define GMX_LINEARALGEBRA_SPARSEMATRIX_H
  37
  38 #include <stdio.h>
  39
  40 #include "../legacyheaders/types/simple.h"
  41
  42 #ifdef __cplusplus
  43 extern "C" {
  44 #endif
  45
  46 typedef struct
  47     gmx_sparsematrix_entry
  48 {
  49     int      col;
  50     real     value;
  51 } gmx_sparsematrix_entry_t;
  52
  53 /*! \brief Sparse matrix storage format
  54  *
  55  *  This structure specifies a storage format for a sparse matrix.
  56  *  The memory requirements are only proportional to the number
  57  *  of nonzero elements, and it provides a reasonably fast way to
  58  *  perform matrix-vector multiplications.
  59  *
  60  *  The data format is very similar to a neighborlist. It is optimized
  61  *  for fast access, but it is difficult to add entries. If you are
  62  *  constructing a matrix you should either do it in exactly the order
  63  *  specified here, or use some other more flexible intermediate structure.
  64  *
  65  *  The index array is of size nrow+1. All non-zero matrix elements
  66  *  on row i are stored in positions index[i] through index[i+1]-1 in
  67  *  the arrays column and value. The column array contains the column
  68  *  index for each entry, in ascending order, and the corresponding
  69  *  position in the value array contains the floating point matrix element.
  70  *
  71  *  index[nrow] should be equal to the total number of elements stored.
  72  *
  73  *  Thus, to find the value of matrix element [5,4] you should loop
  74  *  over positions index[5] to index[6]-1 in column until you either find
  75  *  the value 4, or a higher value (meaning the element was zero).
  76  *
  77  *  It is fairly easy to construct the matrix on-the-fly if you can do
  78  *  it row-by-row.
  79  *
  80  *  IMPORTANT:
  81  *  If compressed_symmetric is set to TRUE, you should only store EITHER the upper OR
  82  *  lower triangle (and the diagonal), and the other half is assumed to be
  83  *  symmetric. Otherwise, if compressed_symmetric==FALSE, no symmetry is implied and all
  84  *  elements should be stored.
  85  *
  86  *  The symmetry compression saves us a factor 2 both in storage and
  87  *  matrix multiplication CPU-time, which can be very useful for huge eigenproblems.
  88  *
  89  *  If you are unsure, just set compressed_symmetric to FALSE and list all elements. If
  90  *  you enable it but still list all elements (both upper and lower triangle) you will be sorry...
  91  *
  92  *  Internally, the sparse data is stored as a separate list for each row, where the list
  93  *  element is a structure with a column and (floating-point) data value. This makes it
  94  *  possible, although not completely transparent, to update values in random access order.
  95  *  The drawback is that the structure will allocate nrow memory regions.
  96  *  The matrix data could be stored in a single contiguous array with indices for each row,
  97  *  but then we could only insert elements at the end without copying the entire matrix.
  98  *
  99  *  After you have
 100  *
 101  *  In other words: Not perfect, but it works.
 102  */
 103 typedef struct
 104     gmx_sparsematrix
 105 {
 106     gmx_bool                     compressed_symmetric; /**< Store half elements and assume symmetry. */
 107     int                          nrow;                 /**< Number of rows in matrix                 */
 108     int *                        ndata;                /**< Number of entries on each row (list)     */
 109     int *                        nalloc;               /**< Allocated entry list length for each row */
 110     gmx_sparsematrix_entry_t **  data;                 /**< data[i] is a list with entries on row i  */
 111 }
 112 gmx_sparsematrix_t;
 113
 114
 115 /*! \brief Allocate a new sparse matrix structure
 116  *
 117  *  The number of rows is used to allocate the index array entry. Obviously you
 118  *  can reallocate these later yourself if necessary - this is a
 119  *  convenience routine.
 120  *
 121  *  By default, the compressed_symmetric flag in the structure will
 122  *  be FALSE. Set it to TRUE manually if you are only storing either the
 123  *  upper or lower half of the matrix.
 124  */
 125 gmx_sparsematrix_t *
 126 gmx_sparsematrix_init            (int                    nrow);
 127
 128
 129 /*! \brief Release all resources used by a sparse matrix structure
 130  *
 131  *  All arrays in the structure will be freed, and the structure itself.
 132  */
 133 void
 134 gmx_sparsematrix_destroy         (gmx_sparsematrix_t *   A);
 135
 136
 137 /*! \brief Print sparse matrix to a stream.
 138  *
 139  *  Mainly used for debugging. Be warned that the real sparse matrices used
 140  *  in Gromacs runs can be HUGE (think 100,000 rows).
 141  */
 142 void
 143 gmx_sparsematrix_print           (FILE *                 stream,
 144                                   gmx_sparsematrix_t *   A);
 145
 146 /* Adds value at row,col. If the value did not exist
 147  * previously it is added, otherwise it is incremented with difference.
 148  *
 149  * The column sort order might change, so you need to run fix_sparsematrix
 150  * once you are done changing the matrix.
 151  */
 152 real
 153 gmx_sparsematrix_value          (gmx_sparsematrix_t *    A,
 154                                  int                     row,
 155                                  int                     col);
 156
 157
 158 /* Adds value at row,col. If the value did not exist
 159  * previously it is added, otherwise it is incremented with difference.
 160  *
 161  * The column sort order might change, so you need to run fix_sparsematrix
 162  * once you are done changing the matrix.
 163  */
 164 void
 165 gmx_sparsematrix_increment_value(gmx_sparsematrix_t *    A,
 166                                  int                     row,
 167                                  int                     col,
 168                                  real                    difference);
 169
 170
 171
 172 /*! \brief Sort elements in each column and remove zeros.
 173  *
 174  *  Sparse matrix access is faster when the elements are stored in
 175  *  increasing column order in each row. In some cases previously non-zero
 176  *  elements will be zero after adding more data, and this routine also removes
 177  *  those entries to reduce the storage requirements.
 178  *
 179  *  It never hurts to run this routine if you have been updating the matrix...
 180  */
 181 void
 182 gmx_sparsematrix_compress       (gmx_sparsematrix_t *    A);
 183
 184
 185
 186 /*! \brief Sparse matrix vector multiplication
 187  *
 188  * Calculate y = A * x for a sparse matrix A.
 189  */
 190 void
 191 gmx_sparsematrix_vector_multiply(gmx_sparsematrix_t *    A,
 192                                  real *                  x,
 193                                  real *                  y);
 194
 195 #ifdef __cplusplus
 196 }
 197 #endif
 198
 199
 200 #endif