计算科学 - PETSc：使用 MatCreateSeqBAIJ 和 MatSetValuesBlocked 阻塞矩阵 - 吾爱随笔录

我对 PETSc 的MatSetValuesBlocked. 当我选择小块大小时，下面的代码适用于矩阵，但在对同一矩阵使用更大块大小时会出错。两个主要问题：

1）什么是正确的使用方法MatSetValuesBlocked？

2)分别为行块和列块生成idxn和全局索引的正确方法是什么。idxm例如，对于 m=n=20 且块大小 = 2 的矩阵，idxn 和 idxm = {0,1,2,3,4,5,6,7,8,9} 是什么？

不确定其他功能是否有错误，但我想知道我做错了什么！

谢谢！

    MatCreateSeqBAIJ(PETSC_COMM_SELF,bs,m,n,nnz*2/m,0,&A);
    MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
    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;
    }
    fclose(file);



/*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);

  MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
  MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

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

  PetscOptionsGetString(NULL,"-fout",fileout,PETSC_MAX_PATH_LEN,NULL);
  PetscViewerBinaryOpen(PETSC_COMM_WORLD,fileout,FILE_MODE_WRITE,&view);
  MatView(A,view);
  VecView(b,view);

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

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