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