0

I am trying to write a program that will do matrix multiplication using Strassen's method using a 2D array represented as a 1D array in column major order. Here is what I have right now.

This is the main method for the calculation.

vector<double> strassen_mult(int n, vector<double> A, vector<double> B) {
    vector<double> C(n*n);
    if (n == 1) {
        C.push_back(A[0] * B[0]);
    } else {
        vector<double> S1 = matrix_partitioner(n/2, 1, B) - matrix_partitioner(n/2, 3, B);
        vector<double> S2 = matrix_partitioner(n/2, 0, A) + matrix_partitioner(n/2, 1, A);
        vector<double> S3 = matrix_partitioner(n/2, 2, A) + matrix_partitioner(n/2, 3, A);
        vector<double> S4 = matrix_partitioner(n/2, 2, B) - matrix_partitioner(n/2, 0, B);
        vector<double> S5 = matrix_partitioner(n/2, 0, A) + matrix_partitioner(n/2, 3, A);
        vector<double> S6 = matrix_partitioner(n/2, 0, B) + matrix_partitioner(n/2, 3, B);
        vector<double> S7 = matrix_partitioner(n/2, 1, A) - matrix_partitioner(n/2, 3, A);
        vector<double> S8 = matrix_partitioner(n/2, 2, B) + matrix_partitioner(n/2, 3, B);
        vector<double> S9 = matrix_partitioner(n/2, 0, A) - matrix_partitioner(n/2, 2, A);
        vector<double> S10 = matrix_partitioner(n/2, 0, B) + matrix_partitioner(n/2, 1, B);
        vector<double> P1 = strassen_mult(n/2, matrix_partitioner(n/2, 0 ,A), S1);
        vector<double> P2 = strassen_mult(n/2, S2, matrix_partitioner(n/2, 3, B));
        vector<double> P3 = strassen_mult(n/2, S3, matrix_partitioner(n/2, 0, B));
        vector<double> P4 = strassen_mult(n/2, matrix_partitioner(n/2, 3 ,A), S4);
        vector<double> P5 = strassen_mult(n/2, S5, S6);
        vector<double> P6 = strassen_mult(n/2, S7, S8);
        vector<double> P7 = strassen_mult(n/2, S9, S10);
        C = equals(C, P5 + P4 - P2 + P6);
        C = equals(C, P1 + P2);
        C = equals(C, P3 + P4);
        C = equals(C, P5 + P1 - P3 - P7);
    }
    return C;
}

Matrix partitioner gets a specific quadrant from the given matrix.

vector<double> matrix_partitioner(int n, int section, vector<double> A) {
    vector<double> C(n*n);
    int start_i, start_j, tmp_j;
    if (section == 0) {
        start_i = 0;
        tmp_j = 0;
    }
    else if (section == 1) {
        start_i = 0;
        tmp_j = n;
    }
    else if (section == 2) {
        start_i = n;
        tmp_j = 0;
    }
    else if (section == 3) {
        start_i = n;
        tmp_j = n;
    }
    for (int i = 0; i < n; i++) {
        start_j = tmp_j;
        for (int j = 0; j < n; j++) {
            C[i+(j*n)] = A[start_i+(start_j*(n*2))];
            start_j++;
        }
        start_i++;
    }
    return C;
}

Equals puts the result from adding the two matrices together into C

vector<double> equals(vector<double> A, vector<double> B) {
    for (int i = 0; i < B.size(); i++) {
        A.push_back(B[i]);
    }
    return A;
}

I have also overloaded the '+' and '-' operator to make adding and subtracting the matrices easier.

These are the results I am getting (I have the iterative method results to compare to, they both use the same print method):

Iterative | 250 260 270 280 | | 618 644 670 696 | | 986 1028 1070 1112 | | 1354 1412 1470 1528 | Strassen | 0 0 0 0 | | 0 0 0 0 | | 0 0 0 0 | | 0 0 0 0 |

Clearly my results are not correct. Can someone help me fix this code?

I have tried moving things around but have not found a good way to do it

Banthrall
  • 27
  • 1
  • Your `equals` is NOT recombining matrix quadrants. You'll end up with an (N/2,N*2) array. You would need to interleave the first row of one array, then the first row of the other array, then the second row of the first array, etc. – Tim Roberts Mar 26 '23 at 05:00
  • 1
    You have too many dynamic allocations. Their overhead is likely to far outweigh any potential speedup you can get from Strassen. – n. m. could be an AI Mar 26 '23 at 05:01
  • Yes, and you are making 24 calls to `matrix_partitioner`, when there are only 4 quadrants in each matrix. – Tim Roberts Mar 26 '23 at 05:03
  • Are you really asking the same question as in your [earlier post](https://stackoverflow.com/q/75845695/10871073)? There seems to be rather a lot of code duplication between the two questions. – Adrian Mole Mar 26 '23 at 05:09
  • Your final additions are not correct -- there should only be one subtraction in that last summation. What is your input data? – Tim Roberts Mar 26 '23 at 05:10
  • `vector equals(vector A, vector B) ` -- Why are you passing the vectors by value? This looks like a Python or Java programmer's attempt at C++ without learning C++ properly. You are making unnecessary copies of vectors. You are *not* passing references to the vectors -- C++ does not work this way. This should be `vector equals(vector& A, const vector& B);`. The same thing here: `vector matrix_partitioner(int n, int section, vector A)` --> `vector matrix_partitioner(int n, int section, const vector& A)` – PaulMcKenzie Mar 26 '23 at 09:04

0 Answers0