Changeset 15785 for trunk/psLib/src/math/psMatrix.c

Timestamp:

Dec 11, 2007, 12:10:13 PM (18 years ago)

Author:

Paul Price

Message:

Fixing matrix multiplication. The check of dimensions was wrong, so I
went ahead and fixed the lot. Matrix multiplication is sufficiently
simple to implement that we shouldn't have to wrap GSL.

File:

• trunk/psLib/src/math/psMatrix.c (modified) (5 diffs)

trunk/psLib/src/math/psMatrix.c

@@ -22 +22 @@
  *  @author Andy Becker, University of Washington (SVD).
  *
- *  @version $Revision: 1.51 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-12-08 01:48:34 $
+ *  @version $Revision: 1.52 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-12-11 22:10:13 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -555 +555 @@
     gsl_linalg_LU_invert(lu, perm, inv);
     if (determinant) {
-        *determinant = (float)gsl_linalg_LU_det(lu, signum);
+        *determinant = (float)gsl_linalg_LU_lndet(lu);
     }
 
@@ -609 +609 @@
                           const psImage* in2)
 {
-    gsl_matrix *m1 = NULL;
-    gsl_matrix *m2 = NULL;
-    gsl_matrix *m3 = NULL;
-
     #define MULTIPLY_CLEANUP { psFree(out); return NULL; }
 
@@ -625 +621 @@
     PS_CHECK_POINTERS(in1, in2, MULTIPLY_CLEANUP);
 
-    if (in1->numRows != in2->numCols) {
-        psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Invalid operation: number of rows of in1 != number of cols of in2.");
+    int rows1 = in1->numRows, cols1 = in1->numCols; // Size of input 1
+    int rows2 = in2->numRows, cols2 = in2->numCols; // Size of input 2
+    if (cols1 != rows2) {
+        psError(PS_ERR_BAD_PARAMETER_SIZE, true,
+                "Incompatible dimensions for matrix multiplication: %dx%d * %dx%d (row x col)",
+                rows1, cols1, rows2, cols2);
         MULTIPLY_CLEANUP;
     }
+    int common = cols1;                 // Common dimension
+
     if (in1->type.type != in2->type.type) {
         psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Invalid operation: data types of in1 and in2 must match.");
@@ -634 +636 @@
     }
 
-    out = psImageRecycle(out, in2->numCols, in1->numRows, in2->type.type);
-
-    // Initialize GSL data
-    m1 = gsl_matrix_alloc(in1->numRows, in1->numCols);
-    psImageToGslMatrix(m1, in1);
-    m2 = gsl_matrix_alloc(in2->numRows, in2->numCols);
-    psImageToGslMatrix(m2, in2);
-    m3 = gsl_matrix_alloc(out->numRows, out->numCols);
-
-    // Perform multiplication
-    // gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m1, m2, 0.0, m3);
-    gsl_linalg_matmult(m1, m2, m3);
-
-    // Copy GSL matrix data to psImage data
-    gslMatrixToPsImage(out, m3);
-
-    // Free GSL structs
-    gsl_matrix_free(m1);
-    gsl_matrix_free(m2);
-    gsl_matrix_free(m3);
+    psElemType type = in1->type.type;   // Data type
+    int outRows = rows1, outCols = cols2; // Size of output
+    out = psImageRecycle(out, outCols, outRows, type);
+
+#define MATRIX_MULTIPLY_CASE(TYPE) \
+  case PS_TYPE_##TYPE: \
+    for (int i = 0; i < outRows; i++) { \
+        for (int j = 0; j < outCols; j++) { \
+            ps##TYPE value = 0.0; \
+            for (int k = 0; k < common; k++) { \
+                value += in1->data.TYPE[i][k] * in2->data.TYPE[k][j]; \
+            } \
+            out->data.TYPE[i][j] = value; \
+        } \
+    } \
+    break;
+
+    switch (type) {
+        MATRIX_MULTIPLY_CASE(F32);
+        MATRIX_MULTIPLY_CASE(F64);
+      default:
+        psAbort("Should never get here.  Unsupported type: %x", type);
+    }
 
     return out;