Just to give a bit more of context I am doing like this: call MatCreate(PETSC_COMM_WORLD, this%A, ierr) call MatSetSizes(this%A, lm, lm, M, M, ierr) ! lm is the local size of the matrix , M the global one call MatSetType(this%A, myType, ierr) ! myType is MATMPIBAIJ call MatSetBlockSize(this%A, 4-bdim, ierr) ! 4-bdim is equal to 3 in this case call MatMPIBAIJSetPreallocation(this%A, 4-bdim, flubioSolvers%d_nz, mesh%d_nnz, flubioSolvers%o_nz, mesh%o_nnz, ierr) ! d_nnz and o_nnz is the number of diagonal and off diagonal non zero blocks call MatSetUp(this%A, ierr) ... some non relevant code ..... call MatSetValuesBlocked(this%A, 4-bdim, mesh%cellGlobalAddr(iElement)-1, 4-bdim, mesh%cellGlobalAddr(iElement)-1, blockValues, INSERT_VALUES, ierr) mesh%cellGlobalAddr(iElement)-1 is an integer equal (not an array) to the block element number the block matrix *blockValues (3x3) *belongs to. Any evident errors?