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