3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
35 #ifndef GMX_LINEARALGEBRA_SPARSEMATRIX_H
36 #define GMX_LINEARALGEBRA_SPARSEMATRIX_H
40 #include "../legacyheaders/types/simple.h"
47 gmx_sparsematrix_entry
51 } gmx_sparsematrix_entry_t;
53 /*! \brief Sparse matrix storage format
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.
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.
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.
71 * index[nrow] should be equal to the total number of elements stored.
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).
77 * It is fairly easy to construct the matrix on-the-fly if you can do
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.
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.
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...
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.
101 * In other words: Not perfect, but it works.
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 */
115 /*! \brief Allocate a new sparse matrix structure
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.
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.
126 gmx_sparsematrix_init (int nrow);
129 /*! \brief Release all resources used by a sparse matrix structure
131 * All arrays in the structure will be freed, and the structure itself.
134 gmx_sparsematrix_destroy (gmx_sparsematrix_t * A);
137 /*! \brief Print sparse matrix to a stream.
139 * Mainly used for debugging. Be warned that the real sparse matrices used
140 * in Gromacs runs can be HUGE (think 100,000 rows).
143 gmx_sparsematrix_print (FILE * stream,
144 gmx_sparsematrix_t * A);
146 /* Adds value at row,col. If the value did not exist
147 * previously it is added, otherwise it is incremented with difference.
149 * The column sort order might change, so you need to run fix_sparsematrix
150 * once you are done changing the matrix.
153 gmx_sparsematrix_value (gmx_sparsematrix_t * A,
158 /* Adds value at row,col. If the value did not exist
159 * previously it is added, otherwise it is incremented with difference.
161 * The column sort order might change, so you need to run fix_sparsematrix
162 * once you are done changing the matrix.
165 gmx_sparsematrix_increment_value(gmx_sparsematrix_t * A,
172 /*! \brief Sort elements in each column and remove zeros.
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.
179 * It never hurts to run this routine if you have been updating the matrix...
182 gmx_sparsematrix_compress (gmx_sparsematrix_t * A);
186 /*! \brief Sparse matrix vector multiplication
188 * Calculate y = A * x for a sparse matrix A.
191 gmx_sparsematrix_vector_multiply(gmx_sparsematrix_t * A,