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 #include "sparsematrix.h"
41 #include "gromacs/legacyheaders/smalloc.h"
47 gmx_sparsematrix_init(int nrow)
50 gmx_sparsematrix_t *A;
71 gmx_sparsematrix_destroy(gmx_sparsematrix_t * A)
75 /* Release each row */
76 for(i=0;i<A->nrow;i++)
81 /* Release the rowdata arrays */
85 /* Release matrix structure itself */
92 gmx_sparsematrix_print(FILE * stream,
93 gmx_sparsematrix_t * A)
97 for(i=0;i<A->nrow;i++)
101 for(j=0;j<A->nrow;j++)
102 fprintf(stream," %6.3f",0.0);
108 for(j=0;j<A->ndata[i];j++)
110 while(k++<A->data[i][j].col)
111 fprintf(stream," %6.3f",0.0);
112 fprintf(stream," %6.3f",A->data[i][j].value);
115 fprintf(stream," %6.3f",0.0);
117 fprintf(stream,"\n");
124 gmx_sparsematrix_value(gmx_sparsematrix_t * A,
128 gmx_bool found = FALSE;
136 /* Find previous value */
137 for(i=0;i<A->ndata[row] && (found==FALSE);i++)
139 if(A->data[row][i].col==col)
142 value = A->data[row][i].value;
146 /* value=0 if we didn't find any match */
153 gmx_sparsematrix_increment_value(gmx_sparsematrix_t * A,
158 gmx_bool found = FALSE;
163 /* Try to find a previous entry with this row/col */
164 for(i=0;i<A->ndata[row] && !found;i++)
166 if(A->data[row][i].col==col)
169 A->data[row][i].value += difference;
173 /* Add a new entry if nothing was found */
176 /* add the value at the end of the row */
177 if(A->ndata[row] == A->nalloc[row])
179 A->nalloc[row] += 100;
180 if(A->data[row]==NULL)
182 snew(A->data[row],A->nalloc[row]);
186 srenew(A->data[row],A->nalloc[row]);
189 A->data[row][A->ndata[row]].col = col;
190 /* Previous value was 0.0 */
191 A->data[row][A->ndata[row]].value = difference;
197 /* Routine to compare column values of two entries, used for quicksort of each row.
199 * The data entries to compare are of the type gmx_sparsematrix_entry_t, but quicksort
200 * uses void pointers as arguments, so we cast them back internally.
203 compare_columns(const void *v1, const void *v2)
205 int c1 = ((gmx_sparsematrix_entry_t *)v1)->col;
206 int c2 = ((gmx_sparsematrix_entry_t *)v2)->col;
218 gmx_sparsematrix_compress(gmx_sparsematrix_t * A)
222 for (i=0;i<A->nrow;i++)
224 /* Remove last value on this row while it is zero */
225 while(A->ndata[i]>0 && A->data[i][A->ndata[i]-1].value==0)
228 /* Go through values on this row and look for more zero elements */
229 for(j=0;j<A->ndata[i];j++)
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)
234 if(A->data[i][j].value==0)
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;
241 /* Only non-zero elements remaining on this row. Sort them after column index */
242 qsort((void *)(A->data[i]),
244 sizeof(gmx_sparsematrix_entry_t),
251 gmx_sparsematrix_vector_multiply(gmx_sparsematrix_t * A,
257 gmx_sparsematrix_entry_t * data; /* pointer to simplify data access */
259 for (i = 0; i < A->nrow; i ++)
262 if(A->compressed_symmetric)
264 for (i = 0; i < A->nrow; i ++)
270 for (k=0;k<A->ndata[i];k++)
283 /* not compressed symmetric storage */
284 for (i = 0; i < A->nrow; i ++)
290 for (k=0;k<A->ndata[i];k++)