src/gmxlib/sparsematrix.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012,2013, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include <stdlib.h>
  43 #include <stdio.h>
  44 #include <assert.h>
  45
  46 #include "smalloc.h"
  47
  48 #include "sparsematrix.h"
  49
  50
  51
  52
  53 gmx_sparsematrix_t *
  54 gmx_sparsematrix_init(int                    nrow)
  55 {
  56     int i;
  57     gmx_sparsematrix_t *A;
  58
  59     snew(A,1);
  60
  61     A->nrow=nrow;
  62     snew(A->ndata,nrow);
  63     snew(A->nalloc,nrow);
  64     snew(A->data,nrow);
  65
  66     for(i=0;i<nrow;i++)
  67     {
  68         A->ndata[i]    = 0;
  69         A->nalloc[i]   = 0;
  70         A->data[i]     = NULL;
  71     }
  72     return A;
  73 }
  74
  75
  76
  77 void
  78 gmx_sparsematrix_destroy(gmx_sparsematrix_t *   A)
  79 {
  80     int i;
  81
  82     /* Release each row */
  83     for(i=0;i<A->nrow;i++)
  84     {
  85         if(A->data[i]!=NULL)
  86             sfree(A->data[i]);
  87     }
  88     /* Release the rowdata arrays */
  89     sfree(A->ndata);
  90     sfree(A->nalloc);
  91     sfree(A->data);
  92     /* Release matrix structure itself */
  93     sfree(A);
  94 }
  95
  96
  97
  98 void
  99 gmx_sparsematrix_print(FILE *                 stream,
 100                        gmx_sparsematrix_t *   A)
 101 {
 102     int i,j,k;
 103
 104     for(i=0;i<A->nrow;i++)
 105     {
 106         if(A->ndata[i]==0)
 107         {
 108             for(j=0;j<A->nrow;j++)
 109                 fprintf(stream," %6.3f",0.0);
 110         }
 111         else
 112         {
 113             k=0;
 114             j=0;
 115             for(j=0;j<A->ndata[i];j++)
 116             {
 117                 while(k++<A->data[i][j].col)
 118                     fprintf(stream," %6.3f",0.0);
 119                 fprintf(stream," %6.3f",A->data[i][j].value);
 120             }
 121             while(k++<A->nrow)
 122                 fprintf(stream," %6.3f",0.0);
 123         }
 124         fprintf(stream,"\n");
 125     }
 126
 127 }
 128
 129
 130 real
 131 gmx_sparsematrix_value(gmx_sparsematrix_t *    A,
 132                        int                     row,
 133                        int                     col)
 134 {
 135     gmx_bool found  = FALSE;
 136     int  i;
 137     real value;
 138
 139     assert(row<A->nrow);
 140
 141     value = 0;
 142
 143     /* Find previous value */
 144     for(i=0;i<A->ndata[row] && (found==FALSE);i++)
 145     {
 146         if(A->data[row][i].col==col)
 147         {
 148             found = TRUE;
 149             value = A->data[row][i].value;
 150         }
 151     }
 152
 153     /* value=0 if we didn't find any match */
 154     return value;
 155 }
 156
 157
 158
 159 void
 160 gmx_sparsematrix_increment_value(gmx_sparsematrix_t *    A,
 161                                  int                     row,
 162                                  int                     col,
 163                                  real                    difference)
 164 {
 165     gmx_bool found  = FALSE;
 166     int i;
 167
 168     assert(row<A->nrow);
 169
 170     /* Try to find a previous entry with this row/col */
 171     for(i=0;i<A->ndata[row] && !found;i++)
 172     {
 173         if(A->data[row][i].col==col)
 174         {
 175             found = TRUE;
 176             A->data[row][i].value += difference;
 177         }
 178     }
 179
 180     /* Add a new entry if nothing was found */
 181     if(!found)
 182     {
 183         /* add the value at the end of the row */
 184         if(A->ndata[row] == A->nalloc[row])
 185         {
 186             A->nalloc[row] += 100;
 187             if(A->data[row]==NULL)
 188             {
 189                 snew(A->data[row],A->nalloc[row]);
 190             }
 191             else
 192             {
 193                 srenew(A->data[row],A->nalloc[row]);
 194             }
 195         }
 196         A->data[row][A->ndata[row]].col   = col;
 197         /* Previous value was 0.0 */
 198         A->data[row][A->ndata[row]].value = difference;
 199         A->ndata[row]++;
 200     }
 201 }
 202
 203
 204 /* Routine to compare column values of two entries, used for quicksort of each row.
 205  *
 206  * The data entries to compare are of the type gmx_sparsematrix_entry_t, but quicksort
 207  * uses void pointers as arguments, so we cast them back internally.
 208  */
 209 static int
 210 compare_columns(const void *v1, const void *v2)
 211 {
 212     int c1 = ((gmx_sparsematrix_entry_t *)v1)->col;
 213     int c2 = ((gmx_sparsematrix_entry_t *)v2)->col;
 214
 215     if(c1<c2)
 216         return -1;
 217     else if(c1>c2)
 218         return 1;
 219     else
 220         return 0;
 221 }
 222
 223
 224 void
 225 gmx_sparsematrix_compress(gmx_sparsematrix_t *    A)
 226 {
 227     int i,j;
 228
 229     for (i=0;i<A->nrow;i++)
 230     {
 231         /* Remove last value on this row while it is zero */
 232         while(A->ndata[i]>0 && A->data[i][A->ndata[i]-1].value==0)
 233             A->ndata[i]--;
 234
 235         /* Go through values on this row and look for more zero elements */
 236         for(j=0;j<A->ndata[i];j++)
 237         {
 238             /* If this element was zero, exchange it with the last non-zero
 239              * element on the row (yes, this will invalidate the sort order)
 240              */
 241             if(A->data[i][j].value==0)
 242             {
 243                 A->data[i][j].value = A->data[i][A->ndata[i]-1].value;
 244                 A->data[i][j].col   = A->data[i][A->ndata[i]-1].col;
 245                 A->ndata[i]--;
 246             }
 247         }
 248         /* Only non-zero elements remaining on this row. Sort them after column index */
 249         qsort((void *)(A->data[i]),
 250               A->ndata[i],
 251               sizeof(gmx_sparsematrix_entry_t),
 252               compare_columns);
 253     }
 254 }
 255
 256
 257 void
 258 gmx_sparsematrix_vector_multiply(gmx_sparsematrix_t *    A,
 259                                  real *                  x,
 260                                  real *                  y)
 261 {
 262     real                        s,v,xi;
 263     int                         i,j,k;
 264     gmx_sparsematrix_entry_t *  data; /* pointer to simplify data access */
 265
 266     for (i = 0; i < A->nrow; i ++)
 267         y[i] = 0;
 268
 269     if(A->compressed_symmetric)
 270     {
 271         for (i = 0; i < A->nrow; i ++)
 272         {
 273             xi = x[i];
 274             s = 0.0;
 275             data = A->data[i];
 276
 277             for (k=0;k<A->ndata[i];k++)
 278             {
 279                 j = data[k].col;
 280                 v = data[k].value;
 281                 s += v * x[j];
 282                 if(i!=j)
 283                     y[j] += v * xi;
 284             }
 285             y[i] += s;
 286         }
 287     }
 288     else
 289     {
 290         /* not compressed symmetric storage */
 291         for (i = 0; i < A->nrow; i ++)
 292         {
 293             xi = x[i];
 294             s = 0.0;
 295             data = A->data[i];
 296
 297             for (k=0;k<A->ndata[i];k++)
 298             {
 299                 j = data[k].col;
 300                 v = data[k].value;
 301                 s += v * x[j];
 302             }
 303             y[i] += s;
 304         }
 305     }
 306 }
 307
 308
 309