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 176