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

   1 /* matrix/oper_complex_source.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 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 int
  21 FUNCTION (gsl_matrix, add) (TYPE (gsl_matrix) * a,
  22                             const TYPE (gsl_matrix) * b)
  23 {
  24   const size_t M = a->size1;
  25   const size_t N = a->size2;
  26
  27   if (b->size1 != M || b->size2 != N)
  28     {
  29       GSL_ERROR ("matrices must have same dimensions", GSL_EBADLEN);
  30     }
  31   else
  32     {
  33       const size_t tda_a = a->tda;
  34       const size_t tda_b = b->tda;
  35
  36       size_t i, j;
  37
  38       for (i = 0; i < M; i++)
  39         {
  40           for (j = 0; j < N; j++)
  41             {
  42               const size_t aij = 2 * (i * tda_a + j);
  43               const size_t bij = 2 * (i * tda_b + j);
  44
  45               a->data[aij] += b->data[bij];
  46               a->data[aij + 1] += b->data[bij + 1];
  47             }
  48         }
  49
  50       return GSL_SUCCESS;
  51     }
  52 }
  53
  54 int
  55 FUNCTION (gsl_matrix, sub) (TYPE (gsl_matrix) * a,
  56                             const TYPE (gsl_matrix) * b)
  57 {
  58   const size_t M = a->size1;
  59   const size_t N = a->size2;
  60
  61   if (b->size1 != M || b->size2 != N)
  62     {
  63       GSL_ERROR ("matrices must have same dimensions", GSL_EBADLEN);
  64     }
  65   else
  66     {
  67       const size_t tda_a = a->tda;
  68       const size_t tda_b = b->tda;
  69
  70       size_t i, j;
  71
  72       for (i = 0; i < M; i++)
  73         {
  74           for (j = 0; j < N; j++)
  75             {
  76               const size_t aij = 2 * (i * tda_a + j);
  77               const size_t bij = 2 * (i * tda_b + j);
  78
  79               a->data[aij] -= b->data[bij];
  80               a->data[aij + 1] -= b->data[bij + 1];
  81             }
  82         }
  83
  84       return GSL_SUCCESS;
  85     }
  86 }
  87
  88 int
  89 FUNCTION (gsl_matrix, mul_elements) (TYPE (gsl_matrix) * a,
  90                                      const TYPE (gsl_matrix) * b)
  91 {
  92   const size_t M = a->size1;
  93   const size_t N = a->size2;
  94
  95   if (b->size1 != M || b->size2 != N)
  96     {
  97       GSL_ERROR ("matrices must have same dimensions", GSL_EBADLEN);
  98     }
  99   else
 100     {
 101       const size_t tda_a = a->tda;
 102       const size_t tda_b = b->tda;
 103
 104       size_t i, j;
 105
 106       for (i = 0; i < M; i++)
 107         {
 108           for (j = 0; j < N; j++)
 109             {
 110               const size_t aij = 2 * (i * tda_a + j);
 111               const size_t bij = 2 * (i * tda_b + j);
 112
 113               ATOMIC ar = a->data[aij];
 114               ATOMIC ai = a->data[aij + 1];
 115
 116               ATOMIC br = b->data[bij];
 117               ATOMIC bi = b->data[bij + 1];
 118
 119               a->data[aij] = ar * br - ai * bi;
 120               a->data[aij + 1] = ar * bi + ai * br;
 121             }
 122         }
 123
 124       return GSL_SUCCESS;
 125     }
 126 }
 127
 128 int
 129 FUNCTION (gsl_matrix, div_elements) (TYPE (gsl_matrix) * a,
 130                                      const TYPE (gsl_matrix) * b)
 131 {
 132   const size_t M = a->size1;
 133   const size_t N = a->size2;
 134
 135   if (b->size1 != M || b->size2 != N)
 136     {
 137       GSL_ERROR ("matrices must have same dimensions", GSL_EBADLEN);
 138     }
 139   else
 140     {
 141       const size_t tda_a = a->tda;
 142       const size_t tda_b = b->tda;
 143
 144       size_t i, j;
 145
 146       for (i = 0; i < M; i++)
 147         {
 148           for (j = 0; j < N; j++)
 149             {
 150               const size_t aij = 2 * (i * tda_a + j);
 151               const size_t bij = 2 * (i * tda_b + j);
 152
 153               ATOMIC ar = a->data[aij];
 154               ATOMIC ai = a->data[aij + 1];
 155
 156               ATOMIC br = b->data[bij];
 157               ATOMIC bi = b->data[bij + 1];
 158
 159               ATOMIC s = 1.0 / hypot(br, bi);
 160
 161               ATOMIC sbr = s * br;
 162               ATOMIC sbi = s * bi;
 163
 164               a->data[aij] = (ar * sbr + ai * sbi) * s;
 165               a->data[aij + 1] = (ai * sbr - ar * sbi) * s;
 166             }
 167         }
 168
 169       return GSL_SUCCESS;
 170     }
 171 }
 172
 173 int FUNCTION (gsl_matrix, scale) (TYPE (gsl_matrix) * a, const BASE x)
 174 {
 175   const size_t M = a->size1;
 176   const size_t N = a->size2;
 177   const size_t tda = a->tda;
 178
 179   size_t i, j;
 180
 181   ATOMIC xr = GSL_REAL(x);
 182   ATOMIC xi = GSL_IMAG(x);
 183
 184   for (i = 0; i < M; i++)
 185     {
 186       for (j = 0; j < N; j++)
 187         {
 188           const size_t aij = 2 * (i * tda + j);
 189
 190           ATOMIC ar = a->data[aij];
 191           ATOMIC ai = a->data[aij + 1];
 192
 193           a->data[aij] = ar * xr - ai * xi;
 194           a->data[aij + 1] = ar * xi + ai * xr;
 195         }
 196     }
 197
 198   return GSL_SUCCESS;
 199 }
 200
 201 int FUNCTION (gsl_matrix, add_constant) (TYPE (gsl_matrix) * a, const BASE x)
 202 {
 203   const size_t M = a->size1;
 204   const size_t N = a->size2;
 205   const size_t tda = a->tda;
 206
 207   size_t i, j;
 208
 209   for (i = 0; i < M; i++)
 210     {
 211       for (j = 0; j < N; j++)
 212         {
 213           a->data[2 * (i * tda + j)] += GSL_REAL (x);
 214           a->data[2 * (i * tda + j) + 1] += GSL_IMAG (x);
 215         }
 216     }
 217
 218   return GSL_SUCCESS;
 219 }
 220
 221
 222 int FUNCTION (gsl_matrix, add_diagonal) (TYPE (gsl_matrix) * a, const BASE x)
 223 {
 224   const size_t M = a->size1;
 225   const size_t N = a->size2;
 226   const size_t tda = a->tda;
 227   const size_t loop_lim = (M < N ? M : N);
 228   size_t i;
 229   for (i = 0; i < loop_lim; i++)
 230     {
 231       a->data[2 * (i * tda + i)] += GSL_REAL (x);
 232       a->data[2 * (i * tda + i) + 1] += GSL_IMAG (x);
 233     }
 234
 235   return GSL_SUCCESS;
 236 }