htswanalysis/MACS/lib/gsl/gsl-1.11/matrix/init_source.c

   1 /* matrix/init_source.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Gerard Jungman, Brian Gough
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 TYPE (gsl_matrix) *
  21 FUNCTION (gsl_matrix, alloc) (const size_t n1, const size_t n2)
  22 {
  23   TYPE (gsl_block) * block;
  24   TYPE (gsl_matrix) * m;
  25
  26   if (n1 == 0)
  27     {
  28       GSL_ERROR_VAL ("matrix dimension n1 must be positive integer",
  29                         GSL_EINVAL, 0);
  30     }
  31   else if (n2 == 0)
  32     {
  33       GSL_ERROR_VAL ("matrix dimension n2 must be positive integer",
  34                         GSL_EINVAL, 0);
  35     }
  36
  37   m = (TYPE (gsl_matrix) *) malloc (sizeof (TYPE (gsl_matrix)));
  38
  39   if (m == 0)
  40     {
  41       GSL_ERROR_VAL ("failed to allocate space for matrix struct",
  42                         GSL_ENOMEM, 0);
  43     }
  44
  45   /* FIXME: n1*n2 could overflow for large dimensions */
  46
  47   block = FUNCTION(gsl_block, alloc) (n1 * n2) ;
  48
  49   if (block == 0)
  50     {
  51       GSL_ERROR_VAL ("failed to allocate space for block",
  52                         GSL_ENOMEM, 0);
  53     }
  54
  55   m->data = block->data;
  56   m->size1 = n1;
  57   m->size2 = n2;
  58   m->tda = n2;
  59   m->block = block;
  60   m->owner = 1;
  61
  62   return m;
  63 }
  64
  65 TYPE (gsl_matrix) *
  66 FUNCTION (gsl_matrix, calloc) (const size_t n1, const size_t n2)
  67 {
  68   size_t i;
  69
  70   TYPE (gsl_matrix) * m = FUNCTION (gsl_matrix, alloc) (n1, n2);
  71
  72   if (m == 0)
  73     return 0;
  74
  75   /* initialize matrix to zero */
  76
  77   for (i = 0; i < MULTIPLICITY * n1 * n2; i++)
  78     {
  79       m->data[i] = 0;
  80     }
  81
  82   return m;
  83 }
  84
  85 TYPE (gsl_matrix) *
  86 FUNCTION (gsl_matrix, alloc_from_block) (TYPE(gsl_block) * block,
  87                                          const size_t offset,
  88                                          const size_t n1,
  89                                          const size_t n2,
  90                                          const size_t d2)
  91 {
  92   TYPE (gsl_matrix) * m;
  93
  94   if (n1 == 0)
  95     {
  96       GSL_ERROR_VAL ("matrix dimension n1 must be positive integer",
  97                         GSL_EINVAL, 0);
  98     }
  99   else if (n2 == 0)
 100     {
 101       GSL_ERROR_VAL ("matrix dimension n2 must be positive integer",
 102                         GSL_EINVAL, 0);
 103     }
 104   else if (d2 < n2)
 105     {
 106       GSL_ERROR_VAL ("matrix dimension d2 must be greater than n2",
 107                         GSL_EINVAL, 0);
 108     }
 109   else if (block->size < offset + n1 * d2)
 110     {
 111       GSL_ERROR_VAL ("matrix size exceeds available block size",
 112                         GSL_EINVAL, 0);
 113     }
 114
 115   m = (TYPE (gsl_matrix) *) malloc (sizeof (TYPE (gsl_matrix)));
 116
 117   if (m == 0)
 118     {
 119       GSL_ERROR_VAL ("failed to allocate space for matrix struct",
 120                         GSL_ENOMEM, 0);
 121     }
 122
 123   m->data = block->data + MULTIPLICITY * offset;
 124   m->size1 = n1;
 125   m->size2 = n2;
 126   m->tda = d2;
 127   m->block = block;
 128   m->owner = 0;
 129
 130   return m;
 131 }
 132
 133
 134 TYPE (gsl_matrix) *
 135 FUNCTION (gsl_matrix, alloc_from_matrix) (TYPE(gsl_matrix) * mm,
 136                                           const size_t k1,
 137                                           const size_t k2,
 138                                           const size_t n1,
 139                                           const size_t n2)
 140 {
 141   TYPE (gsl_matrix) * m;
 142
 143   if (n1 == 0)
 144     {
 145       GSL_ERROR_VAL ("matrix dimension n1 must be positive integer",
 146                         GSL_EINVAL, 0);
 147     }
 148   else if (n2 == 0)
 149     {
 150       GSL_ERROR_VAL ("matrix dimension n2 must be positive integer",
 151                         GSL_EINVAL, 0);
 152     }
 153   else if (k1 + n1 > mm->size1)
 154     {
 155       GSL_ERROR_VAL ("submatrix dimension 1 exceeds size of original",
 156                         GSL_EINVAL, 0);
 157     }
 158   else if (k2 + n2 > mm->size2)
 159     {
 160       GSL_ERROR_VAL ("submatrix dimension 2 exceeds size of original",
 161                         GSL_EINVAL, 0);
 162     }
 163
 164   m = (TYPE (gsl_matrix) *) malloc (sizeof (TYPE (gsl_matrix)));
 165
 166   if (m == 0)
 167     {
 168       GSL_ERROR_VAL ("failed to allocate space for matrix struct",
 169                         GSL_ENOMEM, 0);
 170     }
 171
 172   m->data = mm->data + k1 * mm->tda + k2 ;
 173   m->size1 = n1;
 174   m->size2 = n2;
 175   m->tda = mm->tda;
 176   m->block = mm->block;
 177   m->owner = 0;
 178
 179   return m;
 180 }
 181
 182 void
 183 FUNCTION (gsl_matrix, free) (TYPE (gsl_matrix) * m)
 184 {
 185   if (m->owner)
 186     {
 187       FUNCTION(gsl_block, free) (m->block);
 188     }
 189
 190   free (m);
 191 }
 192 void
 193 FUNCTION (gsl_matrix, set_identity) (TYPE (gsl_matrix) * m)
 194 {
 195   size_t i, j;
 196   ATOMIC * const data = m->data;
 197   const size_t p = m->size1 ;
 198   const size_t q = m->size2 ;
 199   const size_t tda = m->tda ;
 200
 201   const BASE zero = ZERO;
 202   const BASE one = ONE;
 203
 204   for (i = 0; i < p; i++)
 205     {
 206       for (j = 0; j < q; j++)
 207         {
 208           *(BASE *) (data + MULTIPLICITY * (i * tda + j)) = ((i == j) ? one : zero);
 209         }
 210     }
 211 }
 212
 213 void
 214 FUNCTION (gsl_matrix, set_zero) (TYPE (gsl_matrix) * m)
 215 {
 216   size_t i, j;
 217   ATOMIC * const data = m->data;
 218   const size_t p = m->size1 ;
 219   const size_t q = m->size2 ;
 220   const size_t tda = m->tda ;
 221
 222   const BASE zero = ZERO;
 223
 224   for (i = 0; i < p; i++)
 225     {
 226       for (j = 0; j < q; j++)
 227         {
 228           *(BASE *) (data + MULTIPLICITY * (i * tda + j)) = zero;
 229         }
 230     }
 231 }
 232
 233 void
 234 FUNCTION (gsl_matrix, set_all) (TYPE (gsl_matrix) * m, BASE x)
 235 {
 236   size_t i, j;
 237   ATOMIC * const data = m->data;
 238   const size_t p = m->size1 ;
 239   const size_t q = m->size2 ;
 240   const size_t tda = m->tda ;
 241
 242   for (i = 0; i < p; i++)
 243     {
 244       for (j = 0; j < q; j++)
 245         {
 246           *(BASE *) (data + MULTIPLICITY * (i * tda + j)) = x;
 247         }
 248     }
 249 }
 250