Go4Expert

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 18:41.