Dynamic Memory Allocation for Matrix Multiplication

Question

Currently I am working on a code that computes the following equation with two matrices, X and Y, to return the value of matrix W.

W = (XT * X)^-1 * XT * Y

Input Matrix train:

4
10
3.000000,1.000000,1180.000000,1955.000000,221900.000000
3.000000,2.250000,2570.000000,1951.000000,538000.000000
2.000000,1.000000,770.000000,1933.000000,180000.000000
4.000000,3.000000,1960.000000,1965.000000,604000.000000
3.000000,2.000000,1680.000000,1987.000000,510000.000000
4.000000,4.500000,5420.000000,2001.000000,1230000.000000
3.000000,2.250000,1715.000000,1995.000000,257500.000000
3.000000,1.500000,1060.000000,1963.000000,291850.000000
3.000000,1.000000,1780.000000,1960.000000,229500.000000
3.000000,2.500000,1890.000000,2003.000000,323000.000000

Input Matrix test:

3
3.000000,2.500000,3560.000000,1965.000000
2.000000,1.000000,1160.000000,1942.000000
3.000000,1.000000,1430.000000,1927.000000

Result Matrix:

716559
194430
323391

My code returns the proper values for the testcases with the exception of matrices over the size of 1000. I know this is because the size is not dynamically allocated, but I am not sure what the best approach to doing this in my code would be:

#include <stdlib.h>
#include <stdio.h>

int main(int argc, char* argv[]){
        if(argc < 3){
                printf("error.");
                return 0;
        }
        FILE *fptrain = fopen(argv[1], "r");
        if(fptrain == NULL)
        {
                printf("error.");
                return 0;
        }
        int row, col, i, j;
        fscanf(fptrain, "%d", &col);
        col = col+1;
        fscanf(fptrain, "%d", &row);
        char ch;
        //creates the original X and Y matrix

        double trainX[row][col];
        double trainY[row][1];
        for(i=0; i<row; i++)
        {
                trainX[i][0] = 1.000000;
                for(j=1; j<col; j++)
                {
                        fscanf(fptrain, "%lf%c", &trainX[i][j], &ch);
                }
                        fscanf(fptrain, "%lf%c", &trainY[i][0], &ch);
        }
        //creates the X transposed matrix
        double trainXtrans[col][row];
        for(i=0; i<row; i++)
        {
                for(j=0; j<col; j++)
                {
                        trainXtrans[j][i] = trainX[i][j];
                }
        }
        //multiplies X and X transposed
        double trainXtemp[row][row];
        int s;
        double num=0;
        for(i=0; i<col; i++)
        {
                for(j=0; j<col; j++)
                {
                        for(s=0; s<row; s++)
                        {
                                num = num +  trainX[s][j]*trainXtrans[i][s];
                        }
                        trainXtemp[i][j] = num;
                        num = 0;
                }
        }
        //finds the identity matrix of X times X transposed
        double trainXinden[col][col*2];
        for(i=0; i<col; i++)
        {
                for(j=0; j<col; j++)
                {
                        trainXinden[i][j] = trainXtemp[i][j];
                }
                for(j=col; j<col*2; j++)
                {
                        if(j==i+col)
                        {
                                trainXinden[i][j] = 1.000000;
                        }
                        else{
                                trainXinden[i][j] = 0.000000;
                        }
                }
        }
        //finds the inverse of X times X transposed through Gauss Jordan Elimination
        int k;
        double divscalar;
        for(i=0; i<col; i++)
        {
                divscalar = trainXinden[i][i];
                for(j=0; j<col*2; j++)
                {
                        if(trainXinden[i][j] != 0)
                        {
                                trainXinden[i][j] = trainXinden[i][j]/divscalar;
                        }
                }
                for(k=0; k<col; k++)
                {
                        if(i!=k)
                        {
                                double subscalar = trainXinden[k][i];
                                for(j=0; j<col*2; j++)
                                {
                                        trainXinden[k][j] = trainXinden[k][j] - subscalar*trainXinden[i][j];
                                }
                        }
                }
        }
        double trainXinverse[row][row];       
        for(i=0; i<row; i++)          
        {        
                for(j=0; j<col; j++)              
                {        
                        trainXinverse[i][j] = trainXinden[i][j+col];            
                }        
        }        
        double trainXinvXt[col][row];
        for(i=0; i<col; i++)
        {
                for(j=0; j<row; j++)                  
                {                
                        for(s=0; s<col; s++)
                        {        
                                num = num + trainXinverse[i][s]*trainXtrans[s][j];
                        }
                                trainXinvXt[i][j] = num;
                                num = 0;
                }
        }
        //multiples (trainXinvXt) by Y
        double weight[row][1];
        for(i=0; i<col; i++)
        {
                for(s=0; s<row; s++)
                {
                        weight[i][0] +=  trainXinvXt[i][s]*trainY[s][0];
                }
        }

        FILE *fptest = fopen(argv[2], "r");
        if(fptest == NULL)
        {
                printf("error.");
                return 0;
        }
        int testrows;
        fscanf(fptest, "%d", &testrows);
        //creates the test file matrix

        double testM[testrows][col];
        for(i=0; i<testrows; i++)
        {
                testM[i][0] = 1.000000;
                for(j=1; j<col; j++)
                {
                        fscanf(fptest, "%lf%c", &testM[i][j], &ch);
                }
        }


        double prices[testrows][1];
        for(i=0; i<testrows; i++)
        {
                for(s=0; s<col; s++)
                {
                        num = num + testM[i][s]*weight[s][0];
                }
                prices[i][0] = num;
                num = 0;
        }

        for(i=0; i<testrows; i++)
        {
                printf("%0.0lf", prices[i][0]);
                printf("\n");
        }
return 0;
}

When I use malloc on each matrix, for some reason it appears that it is not allowing me to create the augmented matrix or perform my gauss-jordan reduction, which is ruining my final answer.

Ugh! Haven't you heard of functions? You should be using them. — Jonathan Leffler, Oct 25 '17 at 03:51
Anyway, you allocate the memory with `malloc()`. What's the trouble? What have you tried? — Jonathan Leffler, Oct 25 '17 at 03:55
@JonathanLeffler I have tried using a for loop with matrix[i] = (double*)malloc(sizeof(double)*row) for each matrix, but this broke my code and caused more segmentation faults — Joe Cur, Oct 25 '17 at 03:56
I can't even compile this code on VS2013. What are you using to compile it? — jwdonahue, Oct 25 '17 at 04:09
@jwdonahue: The code is using C99 features (notably variable length arrays, or VLAs) that may not be supported by MS Visual Studio, which is still (depressingly) primarily a C90 compiler with some extra features from C99 and/or C11. Note that in C11, VLA support is optional; it is mandatory in C99. — Jonathan Leffler, Oct 25 '17 at 04:39
@JoeCur: It would be helpful to have sample data files in a small size (say 10x10 or thereabouts) so that we don't have to guess or deduce the format of the data. Please read about creating an MCVE ([MCVE]) — the input data and expected output is an important part of an MCVE. — Jonathan Leffler, Oct 25 '17 at 04:41
Is it just me, or when I see `4` read as `col` and `10` read as `row` I would expect a `10x4` *Input Matrix train* data set, but I see `10x5`, and in your code, you create a `10x4` 2D VLA `trainX`, but then fill the first col with `1.000` leaving you reading 3-cols of data out of the `5` from your input data into `trainX` -- WTF? `fscanf` — David C. Rankin, Oct 25 '17 at 06:50
@DavidC.Rankin yes, 10x4 initially, and then we are meant to add a column of 1s to make it a 10x5, but that is besides the point — Joe Cur, Oct 25 '17 at 06:52
But do you know your `fscanf(fptrain, "%lf%c", &trainX[i][j], &ch);` would be wrapping from row 1 to row 2 to read `i = 1`, etc.? I guess it's just an example data mismatch, but is sure looks wonky. — David C. Rankin, Oct 25 '17 at 06:54

H.S. · Answer 1 · 2017-10-25T10:40:01.960

1

The best approach to allocate memory to a multi-dimensional array (taking e.g. of 2D array, since you are using matrix in your program):

int(*matrix)[col] = malloc (sizeof(int[row][col]));

If you want to write a function for it then:

void* allocMatrix (int row, int col)
{
    return malloc (sizeof(int[row][col]));
}

If you are not familiar with this approach of allocating memory dynamically to a multi-dimensional array then read this (very nicely explained by Lundin).

In your program, you are having the matrix of type double, so the allocMatrix() will be -

void* allocMatrix (int row, int col)
{
    return malloc (sizeof(double[row][col]));
}

and in main(), you need to make following changes to create matrix dynamically -

double (*trainX)[col] = allocMatrix (row, col);

Do same changes for other matrices also and make sure to free() the dynamically allocated memory to matrices at appropriate places in your program.

Just for the knowledge purpose, this practice of allocating memory to the multi-dimensional array has been followed commonly, though it's not the best.

edited Oct 25 '17 at 10:40

answered Oct 25 '17 at 04:27

H.S.

11,654
2
15
32

I am being told that the initialization makes a pointer from an integer without a cast, that there are conflicting types for allocMatrix, and previous implicit declaration, when commenting out trainX and trainY and placing int(*trainX)[col] = allocMatrix(row, col); along with its trainY counterpart beneath the commented out matrices – Joe Cur Oct 25 '17 at 04:46
@JoeCur I tried making changes in your program for `trainX` and I am not getting any warning message. Please let me know which compiler you are using. Also, in your program, `trainX` is of type `double` but in your comment you have mentioned that you are placing `int(*trainX)[col]...`. Change `int` to `double`. – H.S. Oct 25 '17 at 06:26
@JoeCur: this should work. But it is hard for anyone to test your code because we don't have a minimal data set, and the layout is curious, at best. – Jonathan Leffler Oct 25 '17 at 06:27
@JonathanLeffler Ill add it just for <3 – Joe Cur Oct 25 '17 at 06:37
Thanks! That'll help. And size 3 is fine; anything small enough to be easy is fine. – Jonathan Leffler Oct 25 '17 at 06:37
@JonathanLeffler edited it with input data and included the problems I am having now with malloc – Joe Cur Oct 25 '17 at 06:41
@JoeCur: Thanks. I'm getting some odd results from my code, with it sometimes working and sometimes producing spurious numbers. I suspect it is 'address space layout randomization' kicking in and doing its stuff, but I've had the same binary run consecutively giving me correct and incorrect results. That's infuriating — I'm only grateful that I spotted it before sending the code to you. – Jonathan Leffler Oct 25 '17 at 07:36
@JoeCur: Please make sure that the `allocMatrix()` should allocate memory to the matrix of type `double`. So, the statement will be - `return malloc (sizeof(double[row][col]));`. With this, I am getting the desired output for the given input that you have specified in your question when I am allocating memory dynamically to matrices in your program. – H.S. Oct 25 '17 at 09:21
@H.S. by adding malloc I am able to handle the oversized arrays, but now I am having issues with creating the augmented matrix as it is not adding the 1.000000s as it is supposed to – Joe Cur Oct 25 '17 at 13:45
@H.S. it is also no longer multiplying properly, why is this? – Joe Cur Oct 25 '17 at 13:56
@JoeCur: I am getting the exact same "Result Matrix" for the given input after adding code for dynamic memory allocation. The `malloc` has nothing to do with the multiplying logic. – H.S. Oct 25 '17 at 16:22

Dynamic Memory Allocation for Matrix Multiplication

1 Answers1