Add matrix multiplication routine

brucehow · Sep 21, 2019 · 436a459 · 436a459
1 parent 13aa9e1
commit 436a459
Show file tree

Hide file tree

Showing 3 changed files with 99 additions and 2 deletions.
diff --git a/matrix.c b/matrix.c
@@ -370,6 +370,49 @@ int main(int argc, char *argv[]) {
             }
             break;
         case MM:
+            // Read input files
+            gettimeofday(&start, NULL);
+            type = read_mat_type(fp);
+            rows = read_mat_dim(fp);
+            cols = read_mat_dim(fp);
+            data = read_line(fp);
+            struct CSR mm1 = csr_format(rows, cols, type, data);
+            gettimeofday(&end, NULL);
+            load_time = get_time(start, end);
+
+            // Filename2 validation
+            FILE *mmfp2 = fopen(filename2, "r");
+            if (mmfp2 == NULL) {
+                fprintf(stderr, "%s: no such file\n", filename2);
+                exit(EXIT_FAILURE);
+            }
+
+            gettimeofday(&start, NULL);
+            type = read_mat_type(mmfp2);
+            rows = read_mat_dim(mmfp2);
+            cols = read_mat_dim(mmfp2);
+            data2 = read_line(mmfp2);
+            struct CSC mm2 = csc_format(rows, cols, type, data2);
+
+            gettimeofday(&end, NULL);
+            load_time += get_time(start, end);
+
+            // Check matrix constraints
+            if (mm1.cols != mm2.rows) {
+                fprintf(stderr, "matrix: the matrix multiplication routine can only be performed if the columns on matrix 1 matches the rows on matrix 2\n");
+                exit(EXIT_FAILURE);
+            } else if (mm1.type != mm2.type) {
+                fprintf(stderr, "matrix: the matrix multiplication routine should only be performed on matrices of identical variable types\n");
+                exit(EXIT_FAILURE);
+            }
+            struct COO mmresult;
+
+            gettimeofday(&start, NULL);
+            mmresult = matrix_multiply(mm1, mm2);
+            gettimeofday(&end, NULL);
+            routine_time = get_time(start, end);
+
+            write_coo_data(stdout, mmresult);
             break;
         case UNDEF: // No routines specified
             usage("no matrix algebra routine specified\n");

diff --git a/matrix.h b/matrix.h
@@ -152,6 +152,15 @@ extern struct CSR matrix_addition_f(struct CSR matrix, struct CSR matrix2);
  */
 extern struct CSR transpose(struct CSC matrix);
 
+/**
+ * Performs matrix multiplication on two given matrices
+ * 
+ * @param matrix The first matrix
+ * @param matrix2 The second matrix
+ * @return struct COO The resulting multiplied matrix
+ */
+extern struct COO matrix_multiply(struct CSR matrix, struct CSC matrix2);
+
 /**
  * Allocates memory of a given size using malloc
  * 

diff --git a/routines.c b/routines.c
@@ -231,6 +231,51 @@ struct CSR transpose(struct CSC matrix) {
     return result;
 }
 
-void matrix_multiply(struct COO matrix) {
-
+struct COO matrix_multiply(struct CSR matrix, struct CSC matrix2) {
+    struct COO result;
+    result.count = 0;
+    result.type = matrix.type;
+    size_t elements_size = MEMSIZ * sizeof(struct ELEMENT);
+    result.elements = allocate(elements_size);
+    result.rows = matrix.rows;
+    result.cols = matrix2.cols;
+
+    for (int i = 0; i < result.rows; i++) {
+        int m1count = matrix.ia[i+1] - matrix.ia[i];
+
+        for (int j = 0; j < result.cols; j++) {
+            int dp = 0; // Final dot product
+
+            int m1pos = matrix.ia[i];
+            int m1seen = 0;
+            int m2count = matrix2.ia[j+1] - matrix2.ia[j];
+            int m2pos = matrix2.ia[j];
+            int m2seen = 0;
+
+            while (m1seen != m1count && m2seen != m2count) {
+                if (matrix.ja[m1pos] == matrix2.ja[m2pos]) { // Row and col index match
+                    dp += (matrix.nnz.i[m1pos++] * matrix2.nnz.i[m2pos++]);
+                    m1seen++;
+                    m2seen++;
+                } else if (matrix.ja[m1pos] < matrix2.ja[m2pos]) {
+                    m1seen++;
+                    m1pos++;
+                } else {
+                    m2seen++;
+                    m2pos++;
+                }
+            }
+            if (dp != 0) {
+                // Dynamically allocate memory for elements struct pointer
+                if (((result.count) * sizeof(struct ELEMENT)) == elements_size) {
+                    elements_size *= 2;
+                    result.elements = reallocate(result.elements, elements_size);
+                }
+                result.elements[result.count].value.i = dp;
+                result.elements[result.count].x = i;
+                result.elements[result.count++].y = j;
+            }
+        }
+    }
+    return result;
 }