1. This site uses cookies. By continuing to use this site, you are agreeing to our use of cookies. Learn More.

Problem in matrix vector multiplication code

Discussion in 'C' started by hashirkk, Sep 22, 2009.

  1. hashirkk

    hashirkk New Member

    Joined:
    Sep 22, 2009
    Messages:
    2
    Likes Received:
    0
    Trophy Points:
    0
    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;
    
    }
    
    }
     
    Last edited by a moderator: Sep 23, 2009
  2. hashirkk

    hashirkk New Member

    Joined:
    Sep 22, 2009
    Messages:
    2
    Likes Received:
    0
    Trophy Points:
    0
    Can anybody please help me in fixing the bug in the code.

    Thanks
     

Share This Page