I am a little confused with PETSc's documentation for MatSetValuesBlocked. The code below works fine for matrices when I choose small block sizes, but I get errors when using bigger block sizes for the same matrix. The two main questions:

1) What is the right way to use MatSetValuesBlocked?

2) What is the correct way to generate idxn and idxm, global indices for the row and column blocks, respectively. As an example, for a matrix with m=n=20 and block size = 2, are idxn and idxm = {0,1,2,3,4,5,6,7,8,9} ?

Not sure if there's a mistake in any other function, but I'd like to know what I am doing wrong!


    PetscInt idxm[(m/bs)] ; //row index of total number of blocks
    PetscInt idxn[(m/bs)] ; //column index of total number of blocks

    /*Values from the sparsematrix.mtx file */

    for (i=0; i<nnz; i++) 
     fscanf(file,"%d %d %le\n",&row,&col,(double*)&val);
     values[i] =  val;

/*Constructing idxm */ 

   for(m_index=0; m_index< (m/bs); m_index++) //Number of blocks rowwise
        idxm[m_index] = m_index; // For the blocks in the first row

   } //idxm construction complete

 /*Constructing idxn*/

  for(b_index = 0; b_index < (n/bs) ; b_index ++) //Each column for each block
    idxn[b_index] = b_index;

  } //idxn construction complete.

  /*Errors occur..*/
  MatSetValuesBlocked(A,m/bs,idxm,n/bs,idxn, values, INSERT_VALUES);


  PetscPrintf(PETSC_COMM_SELF,"Reading matrix completed....\n");


  MatMult(A, b, res); //Actual multiplication

  PetscPrintf(PETSC_COMM_SELF,"Matmult completed...\n");
  ierr = VecDestroy(&b);CHKERRQ(ierr);
  ierr = MatDestroy(&A);CHKERRQ(ierr);
  ierr = VecDestroy(&res);CHKERRQ(ierr);
  ierr = PetscRandomDestroy(&r);CHKERRQ(ierr);
  ierr = PetscFinalize();

MatSetValuesBlocked contributes logically dense blocks. Given the pattern shown above, you should be calling it once per block-row. You're currently calling it as though the arguments were COO format.

Also be aware that reading a matrix from a file is a workflow choice that guarantees that your software will never be scalable. Long-term, you should choose a way to avoid using the file system during matrix assembly.

  • $\begingroup$ Thanks! I revised the insertion but I am still missing something: for(i=0; i<m/bs; i++) //For each blocked row: { MatSetValuesBlocked(A,2,&idxm[i],n/bs,idxn, &values[i][0], INSERT_VALUES); } $\endgroup$
    – i8086
    Sep 13 '13 at 21:56
  • $\begingroup$ What is bs, what m, what is n, and how is values packed? $\endgroup$
    – Jed Brown
    Sep 17 '13 at 5:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.