# matrix vector Product/iterative solver

## matrix vector Product/iterative solver

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[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;ib[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).

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;ib[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;

}

}

2 posts / 0 new
For more complete information about compiler optimizations, see our Optimization Notice.

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[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;ib[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).

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;ib[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;

}

}