src/gromacs/linearalgebra/mtxio.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  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.
  14
  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.
  19  *
  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.
  26  *
  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.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Green Red Orange Magenta Azure Cyan Skyblue
  34  */
  35 #include "mtxio.h"
  36
  37 #ifdef HAVE_CONFIG_H
  38 #include <config.h>
  39 #endif
  40
  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.
  44  */
  45
  46 #include "gromacs/legacyheaders/copyrite.h"
  47 #include "gromacs/legacyheaders/gmx_fatal.h"
  48 #include "gromacs/legacyheaders/gmxfio.h"
  49 #include "gromacs/legacyheaders/smalloc.h"
  50 #include "gromacs/legacyheaders/xdrf.h"
  51 #include "gromacs/linearalgebra/sparsematrix.h"
  52
  53
  54 /* Just a number to identify our file type */
  55 #define GMX_MTXIO_MAGIC_NUMBER  0x34ce8fd2
  56
  57 #define GMX_MTXIO_FULL_MATRIX     0
  58 #define GMX_MTXIO_SPARSE_MATRIX   1
  59
  60
  61
  62 /* Matrix file format definition:
  63  *
  64  * All entries are stored in XDR format.
  65  *
  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
  72  *
  73  * 6. Matrix data.
  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.
  81  */
  82
  83 void gmx_mtxio_write(const char *             filename,
  84                      int                      nrow,
  85                      int                      ncol,
  86                      real *                   full_matrix,
  87                      gmx_sparsematrix_t *     sparse_matrix)
  88 {
  89     t_fileio   *fio;
  90     XDR     *   xd;
  91     int         i, j, prec;
  92     gmx_bool    bDum  = TRUE;
  93     gmx_bool    bRead = FALSE;
  94     size_t      sz;
  95
  96     if (full_matrix != NULL && sparse_matrix != NULL)
  97     {
  98         gmx_fatal(FARGS, "Both full AND sparse matrix specified to gmx_mtxio_write().\n");
  99     }
 100
 101     fio = gmx_fio_open(filename, "w");
 102     gmx_fio_checktype(fio);
 103     xd = gmx_fio_getxdr(fio);
 104
 105     /* Write magic number */
 106     i = GMX_MTXIO_MAGIC_NUMBER;
 107     gmx_fio_do_int(fio, i);
 108
 109     /* Write generating Gromacs version */
 110     gmx_fio_write_string(fio, GromacsVersion());
 111
 112     /* Write 1 for double, 0 for single precision */
 113     if (sizeof(real) == sizeof(double))
 114     {
 115         prec = 1;
 116     }
 117     else
 118     {
 119         prec = 0;
 120     }
 121     gmx_fio_do_int(fio, prec);
 122
 123     gmx_fio_do_int(fio, nrow);
 124     gmx_fio_do_int(fio, ncol);
 125
 126     if (full_matrix != NULL)
 127     {
 128         /* Full matrix storage format */
 129         i = GMX_MTXIO_FULL_MATRIX;
 130         gmx_fio_do_int(fio, i);
 131         sz   = nrow*ncol;
 132         bDum = gmx_fio_ndo_real(fio, full_matrix, sz);
 133     }
 134     else
 135     {
 136         /* Sparse storage */
 137         i = GMX_MTXIO_SPARSE_MATRIX;
 138         gmx_fio_do_int(fio, i);
 139
 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)
 143         {
 144             gmx_fatal(FARGS, "Internal inconsistency in sparse matrix.\n");
 145         }
 146         bDum = gmx_fio_ndo_int(fio, sparse_matrix->ndata, sparse_matrix->nrow);
 147         for (i = 0; i < sparse_matrix->nrow; i++)
 148         {
 149             for (j = 0; j < sparse_matrix->ndata[i]; j++)
 150             {
 151                 gmx_fio_do_int(fio, sparse_matrix->data[i][j].col);
 152                 gmx_fio_do_real(fio, sparse_matrix->data[i][j].value);
 153             }
 154         }
 155     }
 156     gmx_fio_close(fio);
 157 }
 158
 159
 160 void
 161 gmx_mtxio_read (const char *            filename,
 162                 int *                   nrow,
 163                 int *                   ncol,
 164                 real **                 full_matrix,
 165                 gmx_sparsematrix_t **   sparse_matrix)
 166 {
 167     t_fileio   *fio;
 168     XDR     *   xd;
 169     int         i, j, prec;
 170     gmx_bool    bDum  = TRUE;
 171     gmx_bool    bRead = TRUE;
 172     char        gmxver[256];
 173     size_t      sz;
 174
 175     fio = gmx_fio_open(filename, "r");
 176     gmx_fio_checktype(fio);
 177     xd = gmx_fio_getxdr(fio);
 178
 179     /* Read and check magic number */
 180     i = GMX_MTXIO_MAGIC_NUMBER;
 181     gmx_fio_do_int(fio, i);
 182
 183     if (i != GMX_MTXIO_MAGIC_NUMBER)
 184     {
 185         gmx_fatal(FARGS,
 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");
 188     }
 189
 190     /* Read generating Gromacs version */
 191     gmx_fio_do_string(fio, gmxver);
 192
 193     /* Write 1 for double, 0 for single precision */
 194     if (sizeof(real) == sizeof(double))
 195     {
 196         prec = 1;
 197     }
 198     else
 199     {
 200         prec = 0;
 201     }
 202     gmx_fio_do_int(fio, prec);
 203
 204     fprintf(stderr, "Reading %s precision matrix generated by Gromacs %s\n",
 205             (prec == 1) ? "double" : "single", gmxver);
 206
 207     gmx_fio_do_int(fio, i);
 208     *nrow = i;
 209     gmx_fio_do_int(fio, i);
 210     *ncol = i;
 211
 212     gmx_fio_do_int(fio, i);
 213
 214     if (i == GMX_MTXIO_FULL_MATRIX && NULL != full_matrix)
 215     {
 216         printf("Full matrix storage format, nrow=%d, ncols=%d\n", *nrow, *ncol);
 217
 218         sz = (*nrow) * (*ncol);
 219         snew((*full_matrix), sz);
 220         bDum = gmx_fio_ndo_real(fio, (*full_matrix), sz);
 221     }
 222     else if (NULL != sparse_matrix)
 223     {
 224         /* Sparse storage */
 225         printf("Sparse matrix storage format, nrow=%d, ncols=%d\n", *nrow, *ncol);
 226
 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)
 231         {
 232             gmx_fatal(FARGS, "Internal inconsistency in sparse matrix.\n");
 233         }
 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);
 239
 240         for (i = 0; i < (*sparse_matrix)->nrow; i++)
 241         {
 242             (*sparse_matrix)->nalloc[i] = (*sparse_matrix)->ndata[i] + 10;
 243             snew(((*sparse_matrix)->data[i]), (*sparse_matrix)->nalloc[i]);
 244
 245             for (j = 0; j < (*sparse_matrix)->ndata[i]; j++)
 246             {
 247                 gmx_fio_do_int(fio, (*sparse_matrix)->data[i][j].col);
 248                 gmx_fio_do_real(fio, (*sparse_matrix)->data[i][j].value);
 249             }
 250         }
 251     }
 252     gmx_fio_close(fio);
 253 }