This is the mail archive of the
gsl-discuss@sourceware.cygnus.com
mailing list for the GSL project.
matrix multiplication broken
- To: gsl-discuss at sourceware dot cygnus dot com
- Subject: matrix multiplication broken
- From: Paul Walmsley <pjw42 at eng dot cam dot ac dot uk>
- Date: Tue, 15 Feb 2000 15:31:46 +0000 (GMT)
- Reply-To: Paul Walmsley <pwalmsley at iee dot org>
I've recently downloaded and started using the gsl library (primarily with
C++ wrappers for matrix classes) but I've found that the generalised
multiply operation gsl_la_matmult_mod_impl() (in multiply.c) is broken.
It works fine for NxN matrices but not MxN.
I've attached a patch for this. Many apologies: it isn't very elegant, and
has a separate subroutine for each combination of A*B, A'*B, A*B' and
A'*B', and only deals with the case for matrices of doubles. I've also
included a patch to the test_la.c program to test MxN multiplication,
which will fail on the original (release 0.5) source.
The GSL is a great resource -- I'm very glad to have come across it. Keep
up the good work.
Paul Walmsley
--------------------------------------------------------------------------
Signal Processing Group
Cambridge University Engineering Department
pwalmsley@iee.org pjw42@cam.ac.uk http://www-sigproc.eng.cam.ac.uk/~pjw42
*** test_la.c 2000/02/15 14:52:49 1.1
--- test_la.c 2000/02/15 14:59:34
***************
*** 1,5 ****
/* Author: G. Jungman
! * RCS: $Id: test_la.c,v 1.1 2000/02/15 14:52:49 pjw42 Exp $
*/
#include <gsl_test.h>
#include <gsl_math.h>
--- 1,5 ----
/* Author: G. Jungman
! * RCS: $Id: test_la.c,v 1.2 2000/02/15 14:59:23 pjw42 Exp pjw42 $
*/
#include <gsl_test.h>
#include <gsl_math.h>
***************
*** 104,109 ****
--- 104,112 ----
gsl_matrix * A = gsl_matrix_calloc(3, 3);
gsl_matrix * B = gsl_matrix_calloc(3, 3);
gsl_matrix * C = gsl_matrix_calloc(3, 3);
+ gsl_matrix * D = gsl_matrix_calloc(2, 3);
+ gsl_matrix * E = gsl_matrix_calloc(2, 3);
+ gsl_matrix * F = gsl_matrix_calloc(2, 2);
gsl_matrix_set(A, 0, 0, 10.0);
gsl_matrix_set(A, 0, 1, 5.0);
***************
*** 125,130 ****
--- 128,147 ----
gsl_matrix_set(B, 2, 1, 3.0);
gsl_matrix_set(B, 2, 2, 2.0);
+ gsl_matrix_set(D, 0, 0, 10.0);
+ gsl_matrix_set(D, 0, 1, 5.0);
+ gsl_matrix_set(D, 0, 2, 1.0);
+ gsl_matrix_set(D, 1, 0, 1.0);
+ gsl_matrix_set(D, 1, 1, 20.0);
+ gsl_matrix_set(D, 1, 2, 5.0);
+
+ gsl_matrix_set(E, 0, 0, 10.0);
+ gsl_matrix_set(E, 0, 1, 5.0);
+ gsl_matrix_set(E, 0, 2, 2.0);
+ gsl_matrix_set(E, 1, 0, 1.0);
+ gsl_matrix_set(E, 1, 1, 3.0);
+ gsl_matrix_set(E, 1, 2, 2.0);
+
gsl_la_matmult_mod_impl(A, GSL_LA_MOD_NONE, B, GSL_LA_MOD_NONE, C);
s += ( fabs(gsl_matrix_get(C, 0, 0) - 106.0) > GSL_DBL_EPSILON );
s += ( fabs(gsl_matrix_get(C, 0, 1) - 68.0) > GSL_DBL_EPSILON );
***************
*** 169,177 ****
--- 186,217 ----
s += ( fabs(gsl_matrix_get(C, 2, 1) - 30.0) > GSL_DBL_EPSILON );
s += ( fabs(gsl_matrix_get(C, 2, 2) - 30.0) > GSL_DBL_EPSILON );
+ /* now try for non-symmetric matrices */
+ gsl_la_matmult_mod_impl(D, GSL_LA_MOD_TRANSPOSE, E, GSL_LA_MOD_NONE, C);
+ s += ( fabs(gsl_matrix_get(C, 0, 0) - 101.0) > GSL_DBL_EPSILON );
+ s += ( fabs(gsl_matrix_get(C, 0, 1) - 53.0) > GSL_DBL_EPSILON );
+ s += ( fabs(gsl_matrix_get(C, 0, 2) - 22.0) > GSL_DBL_EPSILON );
+ s += ( fabs(gsl_matrix_get(C, 1, 0) - 70.0) > GSL_DBL_EPSILON );
+ s += ( fabs(gsl_matrix_get(C, 1, 1) - 85.0) > GSL_DBL_EPSILON );
+ s += ( fabs(gsl_matrix_get(C, 1, 2) - 50.0) > GSL_DBL_EPSILON );
+ s += ( fabs(gsl_matrix_get(C, 2, 0) - 15.0) > GSL_DBL_EPSILON );
+ s += ( fabs(gsl_matrix_get(C, 2, 1) - 20.0) > GSL_DBL_EPSILON );
+ s += ( fabs(gsl_matrix_get(C, 2, 2) - 12.0) > GSL_DBL_EPSILON );
+
+
+ gsl_la_matmult_mod_impl(D, GSL_LA_MOD_NONE, E, GSL_LA_MOD_TRANSPOSE, F);
+ s += ( fabs(gsl_matrix_get(F, 0, 0) - 127.0) > GSL_DBL_EPSILON );
+ s += ( fabs(gsl_matrix_get(F, 0, 1) - 27.0) > GSL_DBL_EPSILON );
+ s += ( fabs(gsl_matrix_get(F, 1, 0) - 120.0) > GSL_DBL_EPSILON );
+ s += ( fabs(gsl_matrix_get(F, 1, 1) - 71.0) > GSL_DBL_EPSILON );
+
+
gsl_matrix_free(A);
gsl_matrix_free(B);
gsl_matrix_free(C);
+ gsl_matrix_free(D);
+ gsl_matrix_free(E);
+ gsl_matrix_free(F);
return s;
}
*** multiply.c 2000/02/15 14:11:16 1.1
--- multiply.c 2000/02/15 14:59:58
***************
*** 1,5 ****
/* Author: G. Jungman
! * RCS: $Id: multiply.c,v 1.1 2000/02/15 14:11:16 pjw42 Exp $
*/
#include <config.h>
#include <stdlib.h>
--- 1,5 ----
/* Author: G. Jungman
! * RCS: $Id: multiply.c,v 1.2 2000/02/15 14:59:40 pjw42 Exp pjw42 $
*/
#include <config.h>
#include <stdlib.h>
***************
*** 42,48 ****
--- 42,170 ----
}
}
+ // multiply: A'*B
+ int
+ gsl_la_matmult_mod1_impl2(const gsl_matrix * A, const gsl_matrix * B, gsl_matrix * C)
+ {
+ if(A == 0 || B == 0 || C == 0) {
+ return GSL_EFAULT;
+ }
+ else if(A->size1 != B->size1 || A->size2 != C->size1 || B->size2 != C->size2) {
+ return 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 = A->data[0*A->size2 + i];
+ b = B->data[0 + j];
+ temp = a * b;
+ for(k=1; k<A->size1; k++) {
+ a = A->data[k*A->size2 + i];
+ b = B->data[k*B->size2 + j];
+ temp += a * b;
+ }
+ C->data[i*C->size2 + j] = temp;
+ }
+ }
+ return GSL_SUCCESS;
+ }
+ }
+
+
+ // multiply: A*B'
+ int
+ gsl_la_matmult_mod2_impl2(const gsl_matrix * A, const gsl_matrix * B, gsl_matrix * C)
+ {
+ if(A == 0 || B == 0 || C == 0) {
+ return GSL_EFAULT;
+ }
+ else if(A->size2 != B->size2 || A->size1 != C->size1 || B->size1 != C->size2) {
+ return 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 = A->data[i*A->size2 + 0];
+ b = B->data[j*B->size2 + 0];
+ temp = a * b;
+ for(k=1; k<A->size2; k++) {
+ a = A->data[i*A->size2 + k];
+ b = B->data[j*B->size2 + k];
+ temp += a * b;
+ }
+ C->data[i*C->size2 + j] = temp;
+ }
+ }
+ return GSL_SUCCESS;
+ }
+ }
+
+ // multiply: A'*B'
+ int
+ gsl_la_matmult_mod3_impl2(const gsl_matrix * A, const gsl_matrix * B, gsl_matrix * C)
+ {
+ if(A == 0 || B == 0 || C == 0) {
+ return GSL_EFAULT;
+ }
+ else if(A->size1 != B->size2 || A->size2 != C->size1 || B->size1 != C->size2) {
+ return 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 = A->data[0*A->size2 + i];
+ b = B->data[j*B->size2 + 0];
+ temp = a * b;
+ for(k=1; k<A->size1; k++) {
+ a = A->data[k*A->size2 + i];
+ b = B->data[j*B->size2 + k];
+ temp += a * b;
+ }
+ C->data[i*C->size2 + j] = temp;
+ }
+ }
+ return GSL_SUCCESS;
+ }
+ }
+
+
+
+ int
+ gsl_la_matmult_mod_impl(const gsl_matrix * A, gsl_la_matrix_mod_t modA,
+ const gsl_matrix * B, gsl_la_matrix_mod_t modB,
+ gsl_matrix * C)
+ {
+ if(A == 0 || B == 0 || C == 0) {
+ return GSL_EFAULT;
+ }
+ else if(modA == GSL_LA_MOD_NONE && modB == GSL_LA_MOD_NONE) {
+ return gsl_la_matmult_impl(A, B, C);
+ }
+ else if (modA == GSL_LA_MOD_TRANSPOSE && modB == GSL_LA_MOD_NONE) {
+ return gsl_la_matmult_mod1_impl2(A,B,C);
+ }
+ else if (modA == GSL_LA_MOD_NONE && modB == GSL_LA_MOD_TRANSPOSE) {
+ return gsl_la_matmult_mod2_impl2(A,B,C);
+ } else { /* transpose A and B */
+ return gsl_la_matmult_mod3_impl2(A,B,C);
+ }
+ }
+
+
+ /*
int
gsl_la_matmult_mod_impl(const gsl_matrix * A, gsl_la_matrix_mod_t modA,
const gsl_matrix * B, gsl_la_matrix_mod_t modB,
***************
*** 107,109 ****
--- 229,232 ----
}
}
}
+ */