htswanalysis/MACS/lib/gsl/gsl-1.11/linalg/multiply.c

   1 /* linalg/multiply.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 /* Author:  G. Jungman */
  21
  22 #include <config.h>
  23 #include <stdlib.h>
  24 #include <gsl/gsl_math.h>
  25 #include <gsl/gsl_linalg.h>
  26
  27 #define SWAP_SIZE_T(a, b)  do { size_t swap_tmp = a; a = b; b = swap_tmp; } while(0)
  28
  29 int
  30 gsl_linalg_matmult (const gsl_matrix * A, const gsl_matrix * B, gsl_matrix * C)
  31 {
  32   if (A->size2 != B->size1 || A->size1 != C->size1 || B->size2 != C->size2)
  33     {
  34       GSL_ERROR ("matrix sizes are not conformant", GSL_EBADLEN);
  35     }
  36   else
  37     {
  38       double a, b;
  39       double temp;
  40       size_t i, j, k;
  41
  42       for (i = 0; i < C->size1; i++)
  43         {
  44           for (j = 0; j < C->size2; j++)
  45             {
  46               a = gsl_matrix_get (A, i, 0);
  47               b = gsl_matrix_get (B, 0, j);
  48               temp = a * b;
  49               for (k = 1; k < A->size2; k++)
  50                 {
  51                   a = gsl_matrix_get (A, i, k);
  52                   b = gsl_matrix_get (B, k, j);
  53                   temp += a * b;
  54                 }
  55               gsl_matrix_set (C, i, j, temp);
  56             }
  57         }
  58
  59       return GSL_SUCCESS;
  60     }
  61 }
  62
  63
  64 int
  65 gsl_linalg_matmult_mod (const gsl_matrix * A, gsl_linalg_matrix_mod_t modA,
  66                     const gsl_matrix * B, gsl_linalg_matrix_mod_t modB,
  67                     gsl_matrix * C)
  68 {
  69   if (modA == GSL_LINALG_MOD_NONE && modB == GSL_LINALG_MOD_NONE)
  70     {
  71       return gsl_linalg_matmult (A, B, C);
  72     }
  73   else
  74     {
  75       size_t dim1_A = A->size1;
  76       size_t dim2_A = A->size2;
  77       size_t dim1_B = B->size1;
  78       size_t dim2_B = B->size2;
  79       size_t dim1_C = C->size1;
  80       size_t dim2_C = C->size2;
  81
  82       if (modA & GSL_LINALG_MOD_TRANSPOSE)
  83         SWAP_SIZE_T (dim1_A, dim2_A);
  84       if (modB & GSL_LINALG_MOD_TRANSPOSE)
  85         SWAP_SIZE_T (dim1_B, dim2_B);
  86
  87       if (dim2_A != dim1_B || dim1_A != dim1_C || dim2_B != dim2_C)
  88         {
  89           GSL_ERROR ("matrix sizes are not conformant", GSL_EBADLEN);
  90         }
  91       else
  92         {
  93           double a, b;
  94           double temp;
  95           size_t i, j, k;
  96           size_t a1, a2, b1, b2;
  97
  98           for (i = 0; i < dim1_C; i++)
  99             {
 100               for (j = 0; j < dim2_C; j++)
 101                 {
 102                   a1 = i;
 103                   a2 = 0;
 104                   b1 = 0;
 105                   b2 = j;
 106                   if (modA & GSL_LINALG_MOD_TRANSPOSE)
 107                     SWAP_SIZE_T (a1, a2);
 108                   if (modB & GSL_LINALG_MOD_TRANSPOSE)
 109                     SWAP_SIZE_T (b1, b2);
 110
 111                   a = gsl_matrix_get (A, a1, a2);
 112                   b = gsl_matrix_get (B, b1, b2);
 113                   temp = a * b;
 114
 115                   for (k = 1; k < dim2_A; k++)
 116                     {
 117                       a1 = i;
 118                       a2 = k;
 119                       b1 = k;
 120                       b2 = j;
 121                       if (modA & GSL_LINALG_MOD_TRANSPOSE)
 122                         SWAP_SIZE_T (a1, a2);
 123                       if (modB & GSL_LINALG_MOD_TRANSPOSE)
 124                         SWAP_SIZE_T (b1, b2);
 125                       a = gsl_matrix_get (A, a1, a2);
 126                       b = gsl_matrix_get (B, b1, b2);
 127                       temp += a * b;
 128                     }
 129
 130                   gsl_matrix_set (C, i, j, temp);
 131                 }
 132             }
 133
 134           return GSL_SUCCESS;
 135         }
 136     }
 137 }