Go4Expert (http://www.go4expert.com/)
-   C (http://www.go4expert.com/forums/c/)
-   -   Problem in matrix vector multiplication code (http://www.go4expert.com/forums/matrix-vector-multiplication-code-t19520/)

 hashirkk 22Sep2009 23:52

Problem in matrix vector multiplication code

Hello all,
I am copying two different codes below for your analysis and help. The first code takes a large Sparse matrix and apply the iterative solution algorithm on that. The second code below, takes the same sparse matrix, apply a simple storage mechanism onto that and convert it into Dense matrix (containing only non-zero values ) and then apply the iterative solver on that matrix to converge.

Problem: If I set the size of the matrix size small i.e. less than 2kx2k, the results of both the code are approximatel same but for large matrices > 2k, the results of both the code donot match at all.
Note: matrix is a penta-diagonal matrix.

I would really appreciate if somebody can guide me what mistake Am I making or is there any problem in the second code Algorithm.

Code #1:

Generation of a penta-diagonal sparse matrix (A).
Code:

```for (i = 0; i < ROWS; i++) { c[i] = 0.0; for (j = 0; j < COLS; j++) { if (j==i) { if (i==0 && j==0) { matrix[i]j+1 = (float) 1.0; matrix[i][j] = (float) 2.0; matrix[i]j+GRID_SIZE = (float) 3.0; } if (i == ROWS-1) { matrix[i]j-1 = (float) 1.0; matrix[i][j] = (float) 4.0; matrix[i]j-GRID_SIZE = (float) 3.0; } else { if (i % GRID_SIZE == 0 ) matrix[i]j+1 = (float) 1.0; else { matrix[i]j-1 = (float) 1.0; if ((i+1)%GRID_SIZE !=0) matrix[i]j+1 = (float) 1.0; } matrix[i][j] = (float) 2.0; if ((j+GRID_SIZE)<= COLS) matrix[i]j+GRID_SIZE = (float) 3.0; if ((j-GRID_SIZE)>= 0) matrix[i]j-GRID_SIZE = (float) 3.0; } } // j == i } // j } //i // Generation of dummy constant vector (b) b[0] = 3; b[1]= 1; b[2] = 1;b[3] = 1; b[4]= 1; b[5] = 1;b[6] = 1; b[7]= 1; b[8] = 1; for (i=9;i<ROWS;i++) b[i] = b[8]; // Start of Iterative Solver for (it = 0; it < ROWS; it++) { ...... matrix_vector_product (matrix, x, g,ROWS); ...... matrix_vector_product (matrix, d, tmpvec,ROWS); } // end of loop //Matrix Vector product Code. void matrix_vector_product (float **a, float *b, float *c,int ROWS) { int i, j; int N = ROWS; float tmp; for (i = 0; i < N; i++) { tmp = 0.0; for (j = 0; j < N; j++) tmp += a[i][j] * b[j]; c[i] = tmp; } }```

Code #2: Converting sparse matrix to Dense matrix for effecient Data storage

Generation of a penta-diagonal sparse matrix (A).
Code:

```for (i = 0; i < ROWS; i++) { c[i] = 0.0; for (j = 0; j < COLS; j++) { if (j==i) { if (i==0 && j==0) { matrix[i]j+1 = (float) 1.0; matrix[i][j] = (float) 2.0; matrix[i]j+GRID_SIZE = (float) 3.0; } if (i == ROWS-1) { matrix[i]j-1 = (float) 1.0; matrix[i][j] = (float) 4.0; matrix[i]j-GRID_SIZE = (float) 3.0; } else { if (i % GRID_SIZE == 0 ) matrix[i]j+1 = (float) 1.0; else { matrix[i]j-1 = (float) 1.0; if ((i+1)%GRID_SIZE !=0) matrix[i]j+1 = (float) 1.0; } matrix[i][j] = (float) 2.0; if ((j+GRID_SIZE)<= COLS) matrix[i]j+GRID_SIZE = (float) 3.0; if ((j-GRID_SIZE)>= 0) matrix[i]j-GRID_SIZE = (float) 3.0; } } // j == i } // j } //i // Additional steps to conver the sparse matrix to dense matrix ROWS*DIAGONALS+1 and storing of indices in seperate matrix for (i = 0; i < ROWS; i++) { k = 0; for (j = 0; j < COLS; j++) { if (matrix[i][j] != 0.0) { matrix_min[i][k] = matrix[i][j]; // changes made for the penta_diagonal matrix.The index array // will contain the value of the vector j to which the matrix // needs to be multiplied. index[i][k] = j; k++; } } } // Generation of dummy constant vector (b) b[0] = 3; b[1]= 1; b[2] = 1;b[3] = 1; b[4]= 1; b[5] = 1;b[6] = 1; b[7]= 1; b[8] = 1; for (i=9;i<ROWS;i++) b[i] = b[8]; // Start of Iterative Solver for (it = 0; it < ROWS; it++) { ...... min_matrix_vector_product (matrix_min,index, x, g,ROWS); ...... min_matrix_vector_product (matrix_min, index,d, tmpvec,ROWS); } // end of loop //Matrix Vector product Code (modified for Dense matrix. void min_matrix_vector_product (float **a, int **d,float *b, float *c,int ROWS) { int i, j; int N = ROWS; float tmp; for (i = 0; i < N; i++) { tmp = 0.0; for (j = 0; j < DIAGONALS+1; j++){ tmp += a[i][j] * b[d[i]j]; } c[i] = tmp; } }```

 hashirkk 28Sep2009 12:30

Re: Problem in matrix vector multiplication code

Can anybody please help me in fixing the bug in the code.

Thanks

 All times are GMT +5.5. The time now is 17:48.