src/matrix.c

   1 /* PSPP - computes sample statistics.
   2    Copyright (C) 1997-9, 2000 Free Software Foundation, Inc.
   3    Written by Ben Pfaff <blp@gnu.org>.
   4
   5    This program is free software; you can redistribute it and/or
   6    modify it under the terms of the GNU General Public License as
   7    published by the Free Software Foundation; either version 2 of the
   8    License, or (at your option) any later version.
   9
  10    This program is distributed in the hope that it will be useful, but
  11    WITHOUT ANY WARRANTY; without even the implied warranty of
  12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13    General Public License for more details.
  14
  15    You should have received a copy of the GNU General Public License
  16    along with this program; if not, write to the Free Software
  17    Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
  18    02111-1307, USA. */
  19
  20 #include <config.h>
  21 #include <assert.h>
  22 #include <stdlib.h>
  23 #include "alloc.h"
  24 #include "matrix.h"
  25 \f
  26 /* Kahan summation formula, Thm. 8, _What Every Computer Scientist
  27    Should Know About Floating-Point Arithmetic_, David Goldberg,
  28    orig. March 1991 issue of Computing Surveys, also at
  29    <URL:http://www.wam.umd.edu/whats_new/workshop3.0/common-tools/numerical_comp_guide/goldberg1.doc.html>.
  30    Hopefully your compiler won't try to optimize the code below too
  31    much, because that will ruin the precision. */
  32 #define KAHAN_SUMMATION_FORMULA(S)                              \
  33         do                                                      \
  34           {                                                     \
  35             double S_c;                                         \
  36             int S_j;                                            \
  37                                                                 \
  38             S = SUMMATION_ELEMENT (0);                          \
  39             S_c = 0.;                                           \
  40             for (S_j = 1; S_j < SUMMATION_COUNT; S_j++)         \
  41               {                                                 \
  42                 double S_y = SUMMATION_ELEMENT (S_j) - S_c;     \
  43                 double S_t = S + S_y;                           \
  44                 S_c = (S_t - S) - S_y;                          \
  45                 S = S_t;                                        \
  46               }                                                 \
  47           }                                                     \
  48         while (0)
  49
  50 \f
  51 /* Vectors. */
  52
  53 /* Allocate a new vector of length N. */
  54 struct vector *
  55 vec_alloc (int n)
  56 {
  57   struct vector *vec = xmalloc (sizeof *vec);
  58   vec->data = xmalloc (sizeof *vec->data * n);
  59   vec->n = vec->m = n;
  60   return vec;
  61 }
  62
  63 /* Change the length of VEC to N.  The amount of space allocated will
  64    not be lowered, but may be enlarged. */
  65 void
  66 vec_realloc (struct vector *vec, int n)
  67 {
  68   if (n < vec->m)
  69     {
  70       vec->m = n;
  71       vec->data = xrealloc (vec->data, sizeof *vec->data * n);
  72     }
  73   vec->n = n;
  74 }
  75
  76 /* Free vector VEC. */
  77 void
  78 vec_free (struct vector *vec)
  79 {
  80   free (vec->data);
  81   free (vec);
  82 }
  83
  84 /* Set the values in vector VEC to constant VALUE. */
  85 #if 0
  86 void
  87 vec_init (struct vector *vec, double value)
  88 {
  89   double *p;
  90   int i;
  91
  92   p = vec->data;
  93   for (i = 0; i < vec->n; i++)
  94     *p++ = value;
  95 }
  96 #endif
  97
  98 /* Print out vector VEC to stdout for debugging purposes. */
  99 #if GLOBAL_DEBUGGING
 100 #include <stdio.h>
 101 #include "settings.h"
 102
 103 void
 104 vec_print (const struct vector *vec)
 105 {
 106   int i;
 107
 108   for (i = 0; i < vec->n; i++)
 109     {
 110       if (i % ((set_viewwidth - 4) / 8) == 0)
 111         {
 112           if (i)
 113             putchar ('\n');
 114           printf ("%3d:", i);
 115         }
 116
 117       printf ("%8g", vec_elem (vec, i));
 118     }
 119 }
 120 #endif
 121
 122 /* Return the sum of the values in VEC. */
 123 double
 124 vec_total (const struct vector *vec)
 125 {
 126   double sum;
 127
 128 #define SUMMATION_COUNT (vec->n)
 129 #define SUMMATION_ELEMENT(INDEX) (vec_elem (vec, (INDEX)))
 130   KAHAN_SUMMATION_FORMULA (sum);
 131 #undef SUMMATION_COUNT
 132 #undef SUMMATION_ELEMENT
 133
 134   return sum;
 135 }
 136 \f
 137 /* Matrices. */
 138
 139 /* Allocate a new matrix with NR rows and NC columns. */
 140 struct matrix *
 141 mat_alloc (int nr, int nc)
 142 {
 143   struct matrix *mat = xmalloc (sizeof *mat);
 144   mat->nr = nr;
 145   mat->nc = nc;
 146   mat->m = nr * nc;
 147   mat->data = xmalloc (sizeof *mat->data * nr * nc);
 148   return mat;
 149 }
 150
 151 /* Set the size of matrix MAT to NR rows and NC columns.  The matrix
 152    data array will be enlarged if necessary but will not be shrunk. */
 153 void
 154 mat_realloc (struct matrix *mat, int nr, int nc)
 155 {
 156   if (nc * nr > mat->m)
 157     {
 158       mat->m = nc * nr;
 159       mat->data = xrealloc (mat->data, sizeof *mat->data * mat->m);
 160     }
 161   mat->nr = nr;
 162   mat->nc = nc;
 163 }
 164
 165 /* Free matrix MAT. */
 166 void
 167 mat_free (struct matrix *mat)
 168 {
 169   free (mat->data);
 170   free (mat);
 171 }
 172
 173 /* Set all matrix MAT entries to VALUE. */
 174 void
 175 mat_init (struct matrix *mat, double value)
 176 {
 177   double *p;
 178   int i;
 179
 180   p = mat->data;
 181   for (i = 0; i < mat->nr * mat->nc; i++)
 182     *p++ = value;
 183 }
 184
 185 /* Set all MAT entries in row R to VALUE. */
 186 void
 187 mat_init_row (struct matrix *mat, int r, double value)
 188 {
 189   double *p;
 190   int i;
 191
 192   p = &mat_elem (mat, r, 0);
 193   for (i = 0; i < mat->nc; i++)
 194     *p++ = value;
 195 }
 196
 197 /* Set all MAT entries in column C to VALUE. */
 198 void
 199 mat_init_col (struct matrix *mat, int c, double value)
 200 {
 201   double *p;
 202   int i;
 203
 204   p = &mat_elem (mat, 0, c);
 205   for (i = 0; i < mat->nr; i++)
 206     {
 207       *p = value;
 208       p += mat->nc;
 209     }
 210 }
 211
 212 /* Print out MAT entries to stdout, optionally with row and column
 213    labels ROW_LABELS and COL_LABELS. */
 214 #if GLOBAL_DEBUGGING
 215 void
 216 mat_print (const struct matrix *mat,
 217            const struct vector *row_labels,
 218            const struct vector *col_labels)
 219 {
 220   int r, c;
 221
 222   assert (!row_labels || row_labels->n == mat->nr);
 223   if (col_labels)
 224     {
 225       int c;
 226
 227       assert (col_labels->n == mat->nc);
 228       if (row_labels)
 229         printf ("        ");
 230       for (c = 0; c < mat->nc; c++)
 231         printf ("%8g", vec_elem (col_labels, c));
 232     }
 233
 234   for (r = 0; r < mat->nr; r++)
 235     {
 236       if (row_labels)
 237         printf ("%8g:", vec_elem (row_labels, r));
 238       for (c = 0; c < mat->nc; c++)
 239         printf ("%8g", mat_elem (mat, r, c));
 240       putchar ('\n');
 241     }
 242 }
 243 #endif /* GLOBAL_DEBUGGING */
 244
 245 /* Calculate row totals for matrix MAT into vector ROW_TOTS. */
 246 void
 247 mat_row_totals (const struct matrix *mat, struct vector *row_tots)
 248 {
 249   int r;
 250
 251   vec_realloc (row_tots, mat->nr);
 252   for (r = 0; r < mat->nr; r++)
 253     {
 254       double sum;
 255
 256 #define SUMMATION_COUNT (mat->nc)
 257 #define SUMMATION_ELEMENT(INDEX) (mat_elem (mat, r, INDEX))
 258       KAHAN_SUMMATION_FORMULA (sum);
 259 #undef SUMMATION_COUNT
 260 #undef SUMMATION_ELEMENT
 261
 262       vec_elem (row_tots, r) = sum;
 263     }
 264 }
 265
 266 /* Calculate column totals for matrix MAT into vector COL_TOTS. */
 267 void
 268 mat_col_totals (const struct matrix *mat, struct vector *col_tots)
 269 {
 270   int c;
 271
 272   vec_realloc (col_tots, mat->nc);
 273   for (c = 0; c < mat->nc; c++)
 274     {
 275       double sum;
 276
 277 #define SUMMATION_COUNT (mat->nr)
 278 #define SUMMATION_ELEMENT(INDEX) (mat_elem (mat, INDEX, c))
 279       KAHAN_SUMMATION_FORMULA (sum);
 280 #undef SUMMATION_COUNT
 281 #undef SUMMATION_ELEMENT
 282
 283       vec_elem (col_tots, c) = sum;
 284     }
 285 }
 286
 287 /* Return the grand total for matrix MAT.  Of course, if you're also
 288    calculating column or row totals, it would be faster to use
 289    vec_total on one of those sets of totals. */
 290 double
 291 mat_grand_total (const struct matrix *mat)
 292 {
 293   double sum;
 294
 295 #define SUMMATION_COUNT (mat->nr * mat->nc)
 296 #define SUMMATION_ELEMENT(INDEX) (mat->data[INDEX])
 297   KAHAN_SUMMATION_FORMULA (sum);
 298 #undef SUMMATION_COUNT
 299 #undef SUMMATION_ELEMENT
 300
 301   return sum;
 302 }