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
41 /* This module provides routines to read/write sparse or full storage
42 * matrices from/to files. It is normally used for the Hessian matrix
43 * in normal mode analysis.
46 #include "gromacs/legacyheaders/copyrite.h"
47 #include "gromacs/legacyheaders/gmx_fatal.h"
48 #include "gromacs/fileio/gmxfio.h"
49 #include "gromacs/legacyheaders/smalloc.h"
50 #include "gromacs/fileio/xdrf.h"
51 #include "gromacs/linearalgebra/sparsematrix.h"
54 /* Just a number to identify our file type */
55 #define GMX_MTXIO_MAGIC_NUMBER 0x34ce8fd2
57 #define GMX_MTXIO_FULL_MATRIX 0
58 #define GMX_MTXIO_SPARSE_MATRIX 1
62 /* Matrix file format definition:
64 * All entries are stored in XDR format.
66 * 1. Magic number integer, should be GMX_MTXIO_MAGIC_NUMBER
67 * 2. An XDR string specifying the Gromacs version used to generate the file.
68 * 3. Integer to denote precision. 1 if double, 0 if single precision.
69 * 4. Two integers specifying number of rows and columns.
70 * 5. Integer to denote storage type:
71 * GMX_MTXIO_FULL_MATRIX or GMX_MTXIO_SPARSE_MATRIX
74 * a) In case of full matrix, this is nrow*ncol floating-point values.
75 * b) In case of sparse matrix the data is:
76 * - Integer specifying compressed_symmetric format (1=yes, 0=no)
77 * - Integer specifying number of rows (again)
78 * - nrow integers specifying the number of data entries on each row ("ndata")
79 * - All the actual entries for each row, stored contiguous.
80 * Each entry consists of an integer column index and floating-point data value.
83 void gmx_mtxio_write(const char * filename,
87 gmx_sparsematrix_t * sparse_matrix)
93 gmx_bool bRead = FALSE;
96 if (full_matrix != NULL && sparse_matrix != NULL)
98 gmx_fatal(FARGS, "Both full AND sparse matrix specified to gmx_mtxio_write().\n");
101 fio = gmx_fio_open(filename, "w");
102 gmx_fio_checktype(fio);
103 xd = gmx_fio_getxdr(fio);
105 /* Write magic number */
106 i = GMX_MTXIO_MAGIC_NUMBER;
107 gmx_fio_do_int(fio, i);
109 /* Write generating Gromacs version */
110 gmx_fio_write_string(fio, GromacsVersion());
112 /* Write 1 for double, 0 for single precision */
113 if (sizeof(real) == sizeof(double))
121 gmx_fio_do_int(fio, prec);
123 gmx_fio_do_int(fio, nrow);
124 gmx_fio_do_int(fio, ncol);
126 if (full_matrix != NULL)
128 /* Full matrix storage format */
129 i = GMX_MTXIO_FULL_MATRIX;
130 gmx_fio_do_int(fio, i);
132 bDum = gmx_fio_ndo_real(fio, full_matrix, sz);
137 i = GMX_MTXIO_SPARSE_MATRIX;
138 gmx_fio_do_int(fio, i);
140 gmx_fio_do_gmx_bool(fio, sparse_matrix->compressed_symmetric);
141 gmx_fio_do_int(fio, sparse_matrix->nrow);
142 if (sparse_matrix->nrow != nrow)
144 gmx_fatal(FARGS, "Internal inconsistency in sparse matrix.\n");
146 bDum = gmx_fio_ndo_int(fio, sparse_matrix->ndata, sparse_matrix->nrow);
147 for (i = 0; i < sparse_matrix->nrow; i++)
149 for (j = 0; j < sparse_matrix->ndata[i]; j++)
151 gmx_fio_do_int(fio, sparse_matrix->data[i][j].col);
152 gmx_fio_do_real(fio, sparse_matrix->data[i][j].value);
161 gmx_mtxio_read (const char * filename,
165 gmx_sparsematrix_t ** sparse_matrix)
170 gmx_bool bDum = TRUE;
171 gmx_bool bRead = TRUE;
175 fio = gmx_fio_open(filename, "r");
176 gmx_fio_checktype(fio);
177 xd = gmx_fio_getxdr(fio);
179 /* Read and check magic number */
180 i = GMX_MTXIO_MAGIC_NUMBER;
181 gmx_fio_do_int(fio, i);
183 if (i != GMX_MTXIO_MAGIC_NUMBER)
186 "No matrix data found in file. Note that the Hessian matrix format changed\n"
187 "in Gromacs 3.3 to enable portable files and sparse matrix storage.\n");
190 /* Read generating Gromacs version */
191 gmx_fio_do_string(fio, gmxver);
193 /* Write 1 for double, 0 for single precision */
194 if (sizeof(real) == sizeof(double))
202 gmx_fio_do_int(fio, prec);
204 fprintf(stderr, "Reading %s precision matrix generated by Gromacs %s\n",
205 (prec == 1) ? "double" : "single", gmxver);
207 gmx_fio_do_int(fio, i);
209 gmx_fio_do_int(fio, i);
212 gmx_fio_do_int(fio, i);
214 if (i == GMX_MTXIO_FULL_MATRIX && NULL != full_matrix)
216 printf("Full matrix storage format, nrow=%d, ncols=%d\n", *nrow, *ncol);
218 sz = (*nrow) * (*ncol);
219 snew((*full_matrix), sz);
220 bDum = gmx_fio_ndo_real(fio, (*full_matrix), sz);
222 else if (NULL != sparse_matrix)
225 printf("Sparse matrix storage format, nrow=%d, ncols=%d\n", *nrow, *ncol);
227 snew((*sparse_matrix), 1);
228 gmx_fio_do_gmx_bool(fio, (*sparse_matrix)->compressed_symmetric);
229 gmx_fio_do_int(fio, (*sparse_matrix)->nrow);
230 if ((*sparse_matrix)->nrow != *nrow)
232 gmx_fatal(FARGS, "Internal inconsistency in sparse matrix.\n");
234 snew((*sparse_matrix)->ndata, (*sparse_matrix)->nrow);
235 snew((*sparse_matrix)->nalloc, (*sparse_matrix)->nrow);
236 snew((*sparse_matrix)->data, (*sparse_matrix)->nrow);
237 bDum = gmx_fio_ndo_int(fio, (*sparse_matrix)->ndata,
238 (*sparse_matrix)->nrow);
240 for (i = 0; i < (*sparse_matrix)->nrow; i++)
242 (*sparse_matrix)->nalloc[i] = (*sparse_matrix)->ndata[i] + 10;
243 snew(((*sparse_matrix)->data[i]), (*sparse_matrix)->nalloc[i]);
245 for (j = 0; j < (*sparse_matrix)->ndata[i]; j++)
247 gmx_fio_do_int(fio, (*sparse_matrix)->data[i][j].col);
248 gmx_fio_do_real(fio, (*sparse_matrix)->data[i][j].value);