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,2013,2014,2015,2017,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.
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.
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/fileio/gmxfio.h"
47 #include "gromacs/fileio/gmxfio_xdr.h"
48 #include "gromacs/linearalgebra/sparsematrix.h"
49 #include "gromacs/utility/baseversion.h"
50 #include "gromacs/utility/fatalerror.h"
51 #include "gromacs/utility/smalloc.h"
53 /* Just a number to identify our file type */
54 #define GMX_MTXIO_MAGIC_NUMBER 0x34ce8fd2
56 #define GMX_MTXIO_FULL_MATRIX 0
57 #define GMX_MTXIO_SPARSE_MATRIX 1
60 /* Matrix file format definition:
62 * All entries are stored in XDR format.
64 * 1. Magic number integer, should be GMX_MTXIO_MAGIC_NUMBER
65 * 2. An XDR string specifying the Gromacs version used to generate the file.
66 * 3. Integer to denote precision. 1 if double, 0 if single precision.
67 * 4. Two integers specifying number of rows and columns.
68 * 5. Integer to denote storage type:
69 * GMX_MTXIO_FULL_MATRIX or GMX_MTXIO_SPARSE_MATRIX
72 * a) In case of full matrix, this is nrow*ncol floating-point values.
73 * b) In case of sparse matrix the data is:
74 * - Integer specifying compressed_symmetric format (1=yes, 0=no)
75 * - Integer specifying number of rows (again)
76 * - nrow integers specifying the number of data entries on each row ("ndata")
77 * - All the actual entries for each row, stored contiguous.
78 * Each entry consists of an integer column index and floating-point data value.
81 void gmx_mtxio_write(const char* filename
, int nrow
, int ncol
, real
* full_matrix
, gmx_sparsematrix_t
* sparse_matrix
)
87 if (full_matrix
!= nullptr && sparse_matrix
!= nullptr)
89 gmx_fatal(FARGS
, "Both full AND sparse matrix specified to gmx_mtxio_write().\n");
92 fio
= gmx_fio_open(filename
, "w");
94 /* Write magic number */
95 i
= GMX_MTXIO_MAGIC_NUMBER
;
96 gmx_fio_do_int(fio
, i
);
98 /* Write generating Gromacs version */
99 gmx_fio_write_string(fio
, gmx_version());
101 /* Write 1 for double, 0 for single precision */
102 if (sizeof(real
) == sizeof(double))
110 gmx_fio_do_int(fio
, prec
);
112 gmx_fio_do_int(fio
, nrow
);
113 gmx_fio_do_int(fio
, ncol
);
115 if (full_matrix
!= nullptr)
117 /* Full matrix storage format */
118 i
= GMX_MTXIO_FULL_MATRIX
;
119 gmx_fio_do_int(fio
, i
);
121 gmx_fio_ndo_real(fio
, full_matrix
, sz
);
126 i
= GMX_MTXIO_SPARSE_MATRIX
;
127 gmx_fio_do_int(fio
, i
);
129 gmx_fio_do_gmx_bool(fio
, sparse_matrix
->compressed_symmetric
);
130 gmx_fio_do_int(fio
, sparse_matrix
->nrow
);
131 if (sparse_matrix
->nrow
!= nrow
)
133 gmx_fatal(FARGS
, "Internal inconsistency in sparse matrix.\n");
135 gmx_fio_ndo_int(fio
, sparse_matrix
->ndata
, sparse_matrix
->nrow
);
136 for (i
= 0; i
< sparse_matrix
->nrow
; i
++)
138 for (j
= 0; j
< sparse_matrix
->ndata
[i
]; j
++)
140 gmx_fio_do_int(fio
, sparse_matrix
->data
[i
][j
].col
);
141 gmx_fio_do_real(fio
, sparse_matrix
->data
[i
][j
].value
);
149 void gmx_mtxio_read(const char* filename
, int* nrow
, int* ncol
, real
** full_matrix
, gmx_sparsematrix_t
** sparse_matrix
)
156 fio
= gmx_fio_open(filename
, "r");
158 /* Read and check magic number */
159 i
= GMX_MTXIO_MAGIC_NUMBER
;
160 gmx_fio_do_int(fio
, i
);
162 if (i
!= GMX_MTXIO_MAGIC_NUMBER
)
165 "No matrix data found in file. Note that the Hessian matrix format changed\n"
166 "in GROMACS 3.3 to enable portable files and sparse matrix storage.\n");
169 /* Read generating Gromacs version */
170 gmx_fio_do_string(fio
, gmxver
);
172 /* Write 1 for double, 0 for single precision */
173 if (sizeof(real
) == sizeof(double))
181 gmx_fio_do_int(fio
, prec
);
183 fprintf(stderr
, "Reading %s precision matrix generated by GROMACS %s\n",
184 (prec
== 1) ? "double" : "single", gmxver
);
186 gmx_fio_do_int(fio
, i
);
188 gmx_fio_do_int(fio
, i
);
191 gmx_fio_do_int(fio
, i
);
193 if (i
== GMX_MTXIO_FULL_MATRIX
&& nullptr != full_matrix
)
195 printf("Full matrix storage format, nrow=%d, ncols=%d\n", *nrow
, *ncol
);
197 sz
= (*nrow
) * (*ncol
);
198 snew((*full_matrix
), sz
);
199 gmx_fio_ndo_real(fio
, (*full_matrix
), sz
);
201 else if (nullptr != sparse_matrix
)
204 printf("Sparse matrix storage format, nrow=%d, ncols=%d\n", *nrow
, *ncol
);
206 snew((*sparse_matrix
), 1);
207 gmx_fio_do_gmx_bool(fio
, (*sparse_matrix
)->compressed_symmetric
);
208 gmx_fio_do_int(fio
, (*sparse_matrix
)->nrow
);
209 if ((*sparse_matrix
)->nrow
!= *nrow
)
211 gmx_fatal(FARGS
, "Internal inconsistency in sparse matrix.\n");
213 snew((*sparse_matrix
)->ndata
, (*sparse_matrix
)->nrow
);
214 snew((*sparse_matrix
)->nalloc
, (*sparse_matrix
)->nrow
);
215 snew((*sparse_matrix
)->data
, (*sparse_matrix
)->nrow
);
216 gmx_fio_ndo_int(fio
, (*sparse_matrix
)->ndata
, (*sparse_matrix
)->nrow
);
218 for (i
= 0; i
< (*sparse_matrix
)->nrow
; i
++)
220 (*sparse_matrix
)->nalloc
[i
] = (*sparse_matrix
)->ndata
[i
] + 10;
221 snew(((*sparse_matrix
)->data
[i
]), (*sparse_matrix
)->nalloc
[i
]);
223 for (j
= 0; j
< (*sparse_matrix
)->ndata
[i
]; j
++)
225 gmx_fio_do_int(fio
, (*sparse_matrix
)->data
[i
][j
].col
);
226 gmx_fio_do_real(fio
, (*sparse_matrix
)->data
[i
][j
].value
);