2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2012,2014, 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.
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.
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.
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.
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.
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.
39 #include "sparsematrix.h"
45 #include "gromacs/utility/smalloc.h"
48 gmx_sparsematrix_init(int nrow
)
51 gmx_sparsematrix_t
*A
;
57 snew(A
->nalloc
, nrow
);
60 for (i
= 0; i
< nrow
; i
++)
72 gmx_sparsematrix_destroy(gmx_sparsematrix_t
* A
)
76 /* Release each row */
77 for (i
= 0; i
< A
->nrow
; i
++)
79 if (A
->data
[i
] != NULL
)
84 /* Release the rowdata arrays */
88 /* Release matrix structure itself */
95 gmx_sparsematrix_print(FILE * stream
,
96 gmx_sparsematrix_t
* A
)
100 for (i
= 0; i
< A
->nrow
; i
++)
102 if (A
->ndata
[i
] == 0)
104 for (j
= 0; j
< A
->nrow
; j
++)
106 fprintf(stream
, " %6.3f", 0.0);
113 for (j
= 0; j
< A
->ndata
[i
]; j
++)
115 while (k
++ < A
->data
[i
][j
].col
)
117 fprintf(stream
, " %6.3f", 0.0);
119 fprintf(stream
, " %6.3f", A
->data
[i
][j
].value
);
121 while (k
++ < A
->nrow
)
123 fprintf(stream
, " %6.3f", 0.0);
126 fprintf(stream
, "\n");
133 gmx_sparsematrix_value(gmx_sparsematrix_t
* A
,
137 gmx_bool found
= FALSE
;
141 assert(row
< A
->nrow
);
145 /* Find previous value */
146 for (i
= 0; i
< A
->ndata
[row
] && (found
== FALSE
); i
++)
148 if (A
->data
[row
][i
].col
== col
)
151 value
= A
->data
[row
][i
].value
;
155 /* value=0 if we didn't find any match */
162 gmx_sparsematrix_increment_value(gmx_sparsematrix_t
* A
,
167 gmx_bool found
= FALSE
;
170 assert(row
< A
->nrow
);
172 /* Try to find a previous entry with this row/col */
173 for (i
= 0; i
< A
->ndata
[row
] && !found
; i
++)
175 if (A
->data
[row
][i
].col
== col
)
178 A
->data
[row
][i
].value
+= difference
;
182 /* Add a new entry if nothing was found */
185 /* add the value at the end of the row */
186 if (A
->ndata
[row
] == A
->nalloc
[row
])
188 A
->nalloc
[row
] += 100;
189 if (A
->data
[row
] == NULL
)
191 snew(A
->data
[row
], A
->nalloc
[row
]);
195 srenew(A
->data
[row
], A
->nalloc
[row
]);
198 A
->data
[row
][A
->ndata
[row
]].col
= col
;
199 /* Previous value was 0.0 */
200 A
->data
[row
][A
->ndata
[row
]].value
= difference
;
206 /* Routine to compare column values of two entries, used for quicksort of each row.
208 * The data entries to compare are of the type gmx_sparsematrix_entry_t, but quicksort
209 * uses void pointers as arguments, so we cast them back internally.
212 compare_columns(const void *v1
, const void *v2
)
214 int c1
= ((gmx_sparsematrix_entry_t
*)v1
)->col
;
215 int c2
= ((gmx_sparsematrix_entry_t
*)v2
)->col
;
233 gmx_sparsematrix_compress(gmx_sparsematrix_t
* A
)
237 for (i
= 0; i
< A
->nrow
; i
++)
239 /* Remove last value on this row while it is zero */
240 while (A
->ndata
[i
] > 0 && A
->data
[i
][A
->ndata
[i
]-1].value
== 0)
245 /* Go through values on this row and look for more zero elements */
246 for (j
= 0; j
< A
->ndata
[i
]; j
++)
248 /* If this element was zero, exchange it with the last non-zero
249 * element on the row (yes, this will invalidate the sort order)
251 if (A
->data
[i
][j
].value
== 0)
253 A
->data
[i
][j
].value
= A
->data
[i
][A
->ndata
[i
]-1].value
;
254 A
->data
[i
][j
].col
= A
->data
[i
][A
->ndata
[i
]-1].col
;
258 /* Only non-zero elements remaining on this row. Sort them after column index */
259 qsort((void *)(A
->data
[i
]),
261 sizeof(gmx_sparsematrix_entry_t
),
268 gmx_sparsematrix_vector_multiply(gmx_sparsematrix_t
* A
,
274 gmx_sparsematrix_entry_t
* data
; /* pointer to simplify data access */
276 for (i
= 0; i
< A
->nrow
; i
++)
281 if (A
->compressed_symmetric
)
283 for (i
= 0; i
< A
->nrow
; i
++)
289 for (k
= 0; k
< A
->ndata
[i
]; k
++)
304 /* not compressed symmetric storage */
305 for (i
= 0; i
< A
->nrow
; i
++)
311 for (k
= 0; k
< A
->ndata
[i
]; k
++)