src/gromacs/linearalgebra/sparsematrix.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 "sparsematrix.h"
  36
  37 #include <assert.h>
  38 #include <stdlib.h>
  39 #include <stdio.h>
  40
  41 #include "gromacs/legacyheaders/smalloc.h"
  42
  43
  44
  45
  46 gmx_sparsematrix_t *
  47 gmx_sparsematrix_init(int                    nrow)
  48 {
  49     int                 i;
  50     gmx_sparsematrix_t *A;
  51
  52     snew(A, 1);
  53
  54     A->nrow = nrow;
  55     snew(A->ndata, nrow);
  56     snew(A->nalloc, nrow);
  57     snew(A->data, nrow);
  58
  59     for (i = 0; i < nrow; i++)
  60     {
  61         A->ndata[i]    = 0;
  62         A->nalloc[i]   = 0;
  63         A->data[i]     = NULL;
  64     }
  65     return A;
  66 }
  67
  68
  69
  70 void
  71 gmx_sparsematrix_destroy(gmx_sparsematrix_t *   A)
  72 {
  73     int i;
  74
  75     /* Release each row */
  76     for (i = 0; i < A->nrow; i++)
  77     {
  78         if (A->data[i] != NULL)
  79         {
  80             sfree(A->data[i]);
  81         }
  82     }
  83     /* Release the rowdata arrays */
  84     sfree(A->ndata);
  85     sfree(A->nalloc);
  86     sfree(A->data);
  87     /* Release matrix structure itself */
  88     sfree(A);
  89 }
  90
  91
  92
  93 void
  94 gmx_sparsematrix_print(FILE *                 stream,
  95                        gmx_sparsematrix_t *   A)
  96 {
  97     int i, j, k;
  98
  99     for (i = 0; i < A->nrow; i++)
 100     {
 101         if (A->ndata[i] == 0)
 102         {
 103             for (j = 0; j < A->nrow; j++)
 104             {
 105                 fprintf(stream, " %6.3f", 0.0);
 106             }
 107         }
 108         else
 109         {
 110             k = 0;
 111             j = 0;
 112             for (j = 0; j < A->ndata[i]; j++)
 113             {
 114                 while (k++ < A->data[i][j].col)
 115                 {
 116                     fprintf(stream, " %6.3f", 0.0);
 117                 }
 118                 fprintf(stream, " %6.3f", A->data[i][j].value);
 119             }
 120             while (k++ < A->nrow)
 121             {
 122                 fprintf(stream, " %6.3f", 0.0);
 123             }
 124         }
 125         fprintf(stream, "\n");
 126     }
 127
 128 }
 129
 130
 131 real
 132 gmx_sparsematrix_value(gmx_sparsematrix_t *    A,
 133                        int                     row,
 134                        int                     col)
 135 {
 136     gmx_bool found  = FALSE;
 137     int      i;
 138     real     value;
 139
 140     assert(row < A->nrow);
 141
 142     value = 0;
 143
 144     /* Find previous value */
 145     for (i = 0; i < A->ndata[row] && (found == FALSE); i++)
 146     {
 147         if (A->data[row][i].col == col)
 148         {
 149             found = TRUE;
 150             value = A->data[row][i].value;
 151         }
 152     }
 153
 154     /* value=0 if we didn't find any match */
 155     return value;
 156 }
 157
 158
 159
 160 void
 161 gmx_sparsematrix_increment_value(gmx_sparsematrix_t *    A,
 162                                  int                     row,
 163                                  int                     col,
 164                                  real                    difference)
 165 {
 166     gmx_bool found  = FALSE;
 167     int      i;
 168
 169     assert(row < A->nrow);
 170
 171     /* Try to find a previous entry with this row/col */
 172     for (i = 0; i < A->ndata[row] && !found; i++)
 173     {
 174         if (A->data[row][i].col == col)
 175         {
 176             found                  = TRUE;
 177             A->data[row][i].value += difference;
 178         }
 179     }
 180
 181     /* Add a new entry if nothing was found */
 182     if (!found)
 183     {
 184         /* add the value at the end of the row */
 185         if (A->ndata[row] == A->nalloc[row])
 186         {
 187             A->nalloc[row] += 100;
 188             if (A->data[row] == NULL)
 189             {
 190                 snew(A->data[row], A->nalloc[row]);
 191             }
 192             else
 193             {
 194                 srenew(A->data[row], A->nalloc[row]);
 195             }
 196         }
 197         A->data[row][A->ndata[row]].col   = col;
 198         /* Previous value was 0.0 */
 199         A->data[row][A->ndata[row]].value = difference;
 200         A->ndata[row]++;
 201     }
 202 }
 203
 204
 205 /* Routine to compare column values of two entries, used for quicksort of each row.
 206  *
 207  * The data entries to compare are of the type gmx_sparsematrix_entry_t, but quicksort
 208  * uses void pointers as arguments, so we cast them back internally.
 209  */
 210 static int
 211 compare_columns(const void *v1, const void *v2)
 212 {
 213     int c1 = ((gmx_sparsematrix_entry_t *)v1)->col;
 214     int c2 = ((gmx_sparsematrix_entry_t *)v2)->col;
 215
 216     if (c1 < c2)
 217     {
 218         return -1;
 219     }
 220     else if (c1 > c2)
 221     {
 222         return 1;
 223     }
 224     else
 225     {
 226         return 0;
 227     }
 228 }
 229
 230
 231 void
 232 gmx_sparsematrix_compress(gmx_sparsematrix_t *    A)
 233 {
 234     int i, j;
 235
 236     for (i = 0; i < A->nrow; i++)
 237     {
 238         /* Remove last value on this row while it is zero */
 239         while (A->ndata[i] > 0 && A->data[i][A->ndata[i]-1].value == 0)
 240         {
 241             A->ndata[i]--;
 242         }
 243
 244         /* Go through values on this row and look for more zero elements */
 245         for (j = 0; j < A->ndata[i]; j++)
 246         {
 247             /* If this element was zero, exchange it with the last non-zero
 248              * element on the row (yes, this will invalidate the sort order)
 249              */
 250             if (A->data[i][j].value == 0)
 251             {
 252                 A->data[i][j].value = A->data[i][A->ndata[i]-1].value;
 253                 A->data[i][j].col   = A->data[i][A->ndata[i]-1].col;
 254                 A->ndata[i]--;
 255             }
 256         }
 257         /* Only non-zero elements remaining on this row. Sort them after column index */
 258         qsort((void *)(A->data[i]),
 259               A->ndata[i],
 260               sizeof(gmx_sparsematrix_entry_t),
 261               compare_columns);
 262     }
 263 }
 264
 265
 266 void
 267 gmx_sparsematrix_vector_multiply(gmx_sparsematrix_t *    A,
 268                                  real *                  x,
 269                                  real *                  y)
 270 {
 271     real                        s, v, xi;
 272     int                         i, j, k;
 273     gmx_sparsematrix_entry_t *  data; /* pointer to simplify data access */
 274
 275     for (i = 0; i < A->nrow; i++)
 276     {
 277         y[i] = 0;
 278     }
 279
 280     if (A->compressed_symmetric)
 281     {
 282         for (i = 0; i < A->nrow; i++)
 283         {
 284             xi   = x[i];
 285             s    = 0.0;
 286             data = A->data[i];
 287
 288             for (k = 0; k < A->ndata[i]; k++)
 289             {
 290                 j  = data[k].col;
 291                 v  = data[k].value;
 292                 s += v * x[j];
 293                 if (i != j)
 294                 {
 295                     y[j] += v * xi;
 296                 }
 297             }
 298             y[i] += s;
 299         }
 300     }
 301     else
 302     {
 303         /* not compressed symmetric storage */
 304         for (i = 0; i < A->nrow; i++)
 305         {
 306             xi   = x[i];
 307             s    = 0.0;
 308             data = A->data[i];
 309
 310             for (k = 0; k < A->ndata[i]; k++)
 311             {
 312                 j  = data[k].col;
 313                 v  = data[k].value;
 314                 s += v * x[j];
 315             }
 316             y[i] += s;
 317         }
 318     }
 319 }