| /* linalg/multiply.c |
| * |
| * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman, Brian Gough |
| * |
| * This program is free software; you can redistribute it and/or modify |
| * it under the terms of the GNU General Public License as published by |
| * the Free Software Foundation; either version 2 of the License, or (at |
| * your option) any later version. |
| * |
| * This program is distributed in the hope that it will be useful, but |
| * WITHOUT ANY WARRANTY; without even the implied warranty of |
| * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
| * General Public License for more details. |
| * |
| * You should have received a copy of the GNU General Public License |
| * along with this program; if not, write to the Free Software |
| * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. |
| */ |
| |
| /* Author: G. Jungman */ |
| |
| #include <config.h> |
| #include <stdlib.h> |
| #include <gsl/gsl_math.h> |
| #include <gsl/gsl_linalg.h> |
| |
| #define SWAP_SIZE_T(a, b) do { size_t swap_tmp = a; a = b; b = swap_tmp; } while(0) |
| |
| int |
| gsl_linalg_matmult (const gsl_matrix * A, const gsl_matrix * B, gsl_matrix * C) |
| { |
| if (A->size2 != B->size1 || A->size1 != C->size1 || B->size2 != C->size2) |
| { |
| GSL_ERROR ("matrix sizes are not conformant", GSL_EBADLEN); |
| } |
| else |
| { |
| double a, b; |
| double temp; |
| size_t i, j, k; |
| |
| for (i = 0; i < C->size1; i++) |
| { |
| for (j = 0; j < C->size2; j++) |
| { |
| a = gsl_matrix_get (A, i, 0); |
| b = gsl_matrix_get (B, 0, j); |
| temp = a * b; |
| for (k = 1; k < A->size2; k++) |
| { |
| a = gsl_matrix_get (A, i, k); |
| b = gsl_matrix_get (B, k, j); |
| temp += a * b; |
| } |
| gsl_matrix_set (C, i, j, temp); |
| } |
| } |
| |
| return GSL_SUCCESS; |
| } |
| } |
| |
| |
| int |
| gsl_linalg_matmult_mod (const gsl_matrix * A, gsl_linalg_matrix_mod_t modA, |
| const gsl_matrix * B, gsl_linalg_matrix_mod_t modB, |
| gsl_matrix * C) |
| { |
| if (modA == GSL_LINALG_MOD_NONE && modB == GSL_LINALG_MOD_NONE) |
| { |
| return gsl_linalg_matmult (A, B, C); |
| } |
| else |
| { |
| size_t dim1_A = A->size1; |
| size_t dim2_A = A->size2; |
| size_t dim1_B = B->size1; |
| size_t dim2_B = B->size2; |
| size_t dim1_C = C->size1; |
| size_t dim2_C = C->size2; |
| |
| if (modA & GSL_LINALG_MOD_TRANSPOSE) |
| SWAP_SIZE_T (dim1_A, dim2_A); |
| if (modB & GSL_LINALG_MOD_TRANSPOSE) |
| SWAP_SIZE_T (dim1_B, dim2_B); |
| |
| if (dim2_A != dim1_B || dim1_A != dim1_C || dim2_B != dim2_C) |
| { |
| GSL_ERROR ("matrix sizes are not conformant", GSL_EBADLEN); |
| } |
| else |
| { |
| double a, b; |
| double temp; |
| size_t i, j, k; |
| size_t a1, a2, b1, b2; |
| |
| for (i = 0; i < dim1_C; i++) |
| { |
| for (j = 0; j < dim2_C; j++) |
| { |
| a1 = i; |
| a2 = 0; |
| b1 = 0; |
| b2 = j; |
| if (modA & GSL_LINALG_MOD_TRANSPOSE) |
| SWAP_SIZE_T (a1, a2); |
| if (modB & GSL_LINALG_MOD_TRANSPOSE) |
| SWAP_SIZE_T (b1, b2); |
| |
| a = gsl_matrix_get (A, a1, a2); |
| b = gsl_matrix_get (B, b1, b2); |
| temp = a * b; |
| |
| for (k = 1; k < dim2_A; k++) |
| { |
| a1 = i; |
| a2 = k; |
| b1 = k; |
| b2 = j; |
| if (modA & GSL_LINALG_MOD_TRANSPOSE) |
| SWAP_SIZE_T (a1, a2); |
| if (modB & GSL_LINALG_MOD_TRANSPOSE) |
| SWAP_SIZE_T (b1, b2); |
| a = gsl_matrix_get (A, a1, a2); |
| b = gsl_matrix_get (B, b1, b2); |
| temp += a * b; |
| } |
| |
| gsl_matrix_set (C, i, j, temp); |
| } |
| } |
| |
| return GSL_SUCCESS; |
| } |
| } |
| } |