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

   1 /* matrix/swap_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 int
  21 FUNCTION (gsl_matrix, swap_rows) (TYPE (gsl_matrix) * m,
  22                                  const size_t i, const size_t j)
  23 {
  24   const size_t size1 = m->size1;
  25   const size_t size2 = m->size2;
  26
  27   if (i >= size1)
  28     {
  29       GSL_ERROR ("first row index is out of range", GSL_EINVAL);
  30     }
  31
  32   if (j >= size1)
  33     {
  34       GSL_ERROR ("second row index is out of range", GSL_EINVAL);
  35     }
  36
  37   if (i != j)
  38     {
  39       ATOMIC *row1 = m->data + MULTIPLICITY * i * m->tda;
  40       ATOMIC *row2 = m->data + MULTIPLICITY * j * m->tda;
  41
  42       size_t k;
  43
  44       for (k = 0; k < MULTIPLICITY * size2; k++)
  45         {
  46           ATOMIC tmp = row1[k] ;
  47           row1[k] = row2[k] ;
  48           row2[k] = tmp ;
  49         }
  50     }
  51
  52   return GSL_SUCCESS;
  53 }
  54
  55 int
  56 FUNCTION (gsl_matrix, swap_columns) (TYPE (gsl_matrix) * m,
  57                                      const size_t i, const size_t j)
  58 {
  59   const size_t size1 = m->size1;
  60   const size_t size2 = m->size2;
  61
  62   if (i >= size2)
  63     {
  64       GSL_ERROR ("first column index is out of range", GSL_EINVAL);
  65     }
  66
  67   if (j >= size2)
  68     {
  69       GSL_ERROR ("second column index is out of range", GSL_EINVAL);
  70     }
  71
  72   if (i != j)
  73     {
  74       ATOMIC *col1 = m->data + MULTIPLICITY * i;
  75       ATOMIC *col2 = m->data + MULTIPLICITY * j;
  76
  77       size_t p;
  78
  79       for (p = 0; p < size1; p++)
  80         {
  81           size_t k;
  82           size_t n = p * MULTIPLICITY * m->tda;
  83
  84           for (k = 0; k < MULTIPLICITY; k++)
  85             {
  86               ATOMIC tmp = col1[n+k] ;
  87               col1[n+k] = col2[n+k] ;
  88               col2[n+k] = tmp ;
  89             }
  90         }
  91     }
  92
  93   return GSL_SUCCESS;
  94 }
  95
  96
  97 int
  98 FUNCTION (gsl_matrix, swap_rowcol) (TYPE (gsl_matrix) * m,
  99                                     const size_t i, const size_t j)
 100 {
 101   const size_t size1 = m->size1;
 102   const size_t size2 = m->size2;
 103
 104   if (size1 != size2)
 105     {
 106       GSL_ERROR ("matrix must be square to swap row and column", GSL_ENOTSQR);
 107     }
 108
 109   if (i >= size1)
 110     {
 111       GSL_ERROR ("row index is out of range", GSL_EINVAL);
 112     }
 113
 114   if (j >= size2)
 115     {
 116       GSL_ERROR ("column index is out of range", GSL_EINVAL);
 117     }
 118
 119   {
 120     ATOMIC *row = m->data + MULTIPLICITY * i * m->tda;
 121     ATOMIC *col = m->data + MULTIPLICITY * j;
 122
 123     size_t p;
 124
 125     for (p = 0; p < size1; p++)
 126       {
 127         size_t k;
 128
 129         size_t r = p * MULTIPLICITY;
 130         size_t c = p * MULTIPLICITY * m->tda;
 131
 132           for (k = 0; k < MULTIPLICITY; k++)
 133             {
 134               ATOMIC tmp = col[c+k] ;
 135               col[c+k] = row[r+k] ;
 136               row[r+k] = tmp ;
 137             }
 138         }
 139     }
 140
 141   return GSL_SUCCESS;
 142 }
 143
 144
 145 int
 146 FUNCTION (gsl_matrix, transpose) (TYPE (gsl_matrix) * m)
 147 {
 148   const size_t size1 = m->size1;
 149   const size_t size2 = m->size2;
 150   size_t i, j, k;
 151
 152   if (size1 != size2)
 153     {
 154       GSL_ERROR ("matrix must be square to take transpose", GSL_ENOTSQR);
 155     }
 156
 157   for (i = 0; i < size1; i++)
 158     {
 159       for (j = i + 1 ; j < size2 ; j++)
 160         {
 161           for (k = 0; k < MULTIPLICITY; k++)
 162             {
 163               size_t e1 = (i *  m->tda + j) * MULTIPLICITY + k ;
 164               size_t e2 = (j *  m->tda + i) * MULTIPLICITY + k ;
 165               {
 166                 ATOMIC tmp = m->data[e1] ;
 167                 m->data[e1] = m->data[e2] ;
 168                 m->data[e2] = tmp ;
 169               }
 170             }
 171         }
 172     }
 173
 174   return GSL_SUCCESS;
 175 }
 176
 177 int
 178 FUNCTION (gsl_matrix, transpose_memcpy) (TYPE (gsl_matrix) * dest,
 179                                          const TYPE (gsl_matrix) * src)
 180 {
 181   const size_t src_size1 = src->size1;
 182   const size_t src_size2 = src->size2;
 183
 184   const size_t dest_size1 = dest->size1;
 185   const size_t dest_size2 = dest->size2;
 186
 187   size_t i, j, k;
 188
 189   if (dest_size2 != src_size1 || dest_size1 != src_size2)
 190     {
 191       GSL_ERROR ("dimensions of dest matrix must be transpose of src matrix",
 192                  GSL_EBADLEN);
 193     }
 194
 195   for (i = 0; i < dest_size1; i++)
 196     {
 197       for (j = 0 ; j < dest_size2; j++)
 198         {
 199           for (k = 0; k < MULTIPLICITY; k++)
 200             {
 201               size_t e1 = (i *  dest->tda + j) * MULTIPLICITY + k ;
 202               size_t e2 = (j *  src->tda + i) * MULTIPLICITY + k ;
 203
 204               dest->data[e1] = src->data[e2] ;
 205             }
 206         }
 207     }
 208
 209   return GSL_SUCCESS;
 210 }