matrix vector Product/iterative solver

[hashirkk]

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).

for (i = 0; i < ROWS; i++) {
c = 0.0;
for (j = 0; j < COLS; j++)
{
if (j==i)
{
if (i==0 && j==0)
{
matrixj+1 = (float) 1.0;
matrix = (float) 2.0;
matrixj+GRID_SIZE = (float) 3.0;
}
if (i == ROWS-1)
{
matrixj-1 = (float) 1.0;
matrix = (float) 4.0;
matrixj-GRID_SIZE = (float) 3.0;
}
else
{
if (i % GRID_SIZE == 0 )
matrixj+1 = (float) 1.0;
else
{
matrixj-1 = (float) 1.0;
if ((i+1)%GRID_SIZE !=0)
matrixj+1 = (float) 1.0;
}

matrix = (float) 2.0;
if ((j+GRID_SIZE)<= COLS)
matrixj+GRID_SIZE = (float) 3.0;
if ((j-GRID_SIZE)>= 0)
matrixj-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;ib = 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 * b;
c = tmp;

}

}


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

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

for (i = 0; i < ROWS; i++) {
c = 0.0;
for (j = 0; j < COLS; j++)
{
if (j==i)
{
if (i==0 && j==0)
{
matrixj+1 = (float) 1.0;
matrix = (float) 2.0;
matrixj+GRID_SIZE = (float) 3.0;
}
if (i == ROWS-1)
{
matrixj-1 = (float) 1.0;
matrix = (float) 4.0;
matrixj-GRID_SIZE = (float) 3.0;
}
else
{
if (i % GRID_SIZE == 0 )
matrixj+1 = (float) 1.0;
else
{
matrixj-1 = (float) 1.0;
if ((i+1)%GRID_SIZE !=0)
matrixj+1 = (float) 1.0;
}

matrix = (float) 2.0;
if ((j+GRID_SIZE)<= COLS)
matrixj+GRID_SIZE = (float) 3.0;
if ((j-GRID_SIZE)>= 0)
matrixj-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 != 0.0)
{
matrix_min = matrix;

// 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 = 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;ib = 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 * b[dj];
}
c = tmp;

}

}

[hashirkk]

Can please somebody help me with the above code problem? Let me know if the question or the pasted code is not clear

Thanks

Quoting - hashirkk

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).

for (i = 0; i < ROWS; i++) {
c = 0.0;
for (j = 0; j < COLS; j++)
{
if (j==i)
{
if (i==0 && j==0)
{
matrixj+1 = (float) 1.0;
matrix = (float) 2.0;
matrixj+GRID_SIZE = (float) 3.0;
}
if (i == ROWS-1)
{
matrixj-1 = (float) 1.0;
matrix = (float) 4.0;
matrixj-GRID_SIZE = (float) 3.0;
}
else
{
if (i % GRID_SIZE == 0 )
matrixj+1 = (float) 1.0;
else
{
matrixj-1 = (float) 1.0;
if ((i+1)%GRID_SIZE !=0)
matrixj+1 = (float) 1.0;
}

matrix = (float) 2.0;
if ((j+GRID_SIZE)<= COLS)
matrixj+GRID_SIZE = (float) 3.0;
if ((j-GRID_SIZE)>= 0)
matrixj-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;ib = 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 * b;
c = tmp;

}

}


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

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

for (i = 0; i < ROWS; i++) {
c = 0.0;
for (j = 0; j < COLS; j++)
{
if (j==i)
{
if (i==0 && j==0)
{
matrixj+1 = (float) 1.0;
matrix = (float) 2.0;
matrixj+GRID_SIZE = (float) 3.0;
}
if (i == ROWS-1)
{
matrixj-1 = (float) 1.0;
matrix = (float) 4.0;
matrixj-GRID_SIZE = (float) 3.0;
}
else
{
if (i % GRID_SIZE == 0 )
matrixj+1 = (float) 1.0;
else
{
matrixj-1 = (float) 1.0;
if ((i+1)%GRID_SIZE !=0)
matrixj+1 = (float) 1.0;
}

matrix = (float) 2.0;
if ((j+GRID_SIZE)<= COLS)
matrixj+GRID_SIZE = (float) 3.0;
if ((j-GRID_SIZE)>= 0)
matrixj-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 != 0.0)
{
matrix_min = matrix;

// 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 = 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;ib = 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 * b[dj];
}
c = tmp;

}

}