C++ MatGetLocalSize函数代码示例

/*@
   MatMultEqual - Compares matrix-vector products of two matrices.

   Collective on Mat

   Input Parameters:
+  A - the first matrix
-  B - the second matrix
-  n - number of random vectors to be tested

   Output Parameter:
.  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

   Level: intermediate

   Concepts: matrices^equality between
@*/
PetscErrorCode  MatMultEqual(Mat A,Mat B,PetscInt n,PetscBool  *flg)
{
  PetscErrorCode ierr;
  Vec            x,s1,s2;
  PetscRandom    rctx;
  PetscReal      r1,r2,tol=1.e-10;
  PetscInt       am,an,bm,bn,k;
  PetscScalar    none = -1.0;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(A,MAT_CLASSID,1); 
  PetscValidHeaderSpecific(B,MAT_CLASSID,2);
  ierr = MatGetLocalSize(A,&am,&an);CHKERRQ(ierr);
  ierr = MatGetLocalSize(B,&bm,&bn);CHKERRQ(ierr);
  if (am != bm || an != bn) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Mat A,Mat B: local dim %D %D %D %D",am,bm,an,bn);
  PetscCheckSameComm(A,1,B,2);

#if defined(PETSC_USE_REAL_SINGLE)
  tol  = 1.e-5;
#endif
  ierr = PetscRandomCreate(((PetscObject)A)->comm,&rctx);CHKERRQ(ierr);
  ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr);
  ierr = VecCreate(((PetscObject)A)->comm,&x);CHKERRQ(ierr);
  ierr = VecSetSizes(x,an,PETSC_DECIDE);CHKERRQ(ierr);
  ierr = VecSetFromOptions(x);CHKERRQ(ierr);
  
  ierr = VecCreate(((PetscObject)A)->comm,&s1);CHKERRQ(ierr);
  ierr = VecSetSizes(s1,am,PETSC_DECIDE);CHKERRQ(ierr);
  ierr = VecSetFromOptions(s1);CHKERRQ(ierr);
  ierr = VecDuplicate(s1,&s2);CHKERRQ(ierr);
  
  *flg = PETSC_TRUE;
  for (k=0; k<n; k++) {
    ierr = VecSetRandom(x,rctx);CHKERRQ(ierr);
    ierr = MatMult(A,x,s1);CHKERRQ(ierr);
    ierr = MatMult(B,x,s2);CHKERRQ(ierr);
    ierr = VecNorm(s2,NORM_INFINITY,&r2);CHKERRQ(ierr);
    if (r2 < tol){
      ierr = VecNorm(s1,NORM_INFINITY,&r1);CHKERRQ(ierr);
    } else {
      ierr = VecAXPY(s2,none,s1);CHKERRQ(ierr);
      ierr = VecNorm(s2,NORM_INFINITY,&r1);CHKERRQ(ierr);
      r1 /= r2;
    }
    if (r1 > tol) {
      *flg = PETSC_FALSE;
      ierr = PetscInfo2(A,"Error: %D-th MatMult() %G\n",k,r1);CHKERRQ(ierr);
      break;
    } 
  }
  ierr = PetscRandomDestroy(&rctx);CHKERRQ(ierr);
  ierr = VecDestroy(&x);CHKERRQ(ierr);
  ierr = VecDestroy(&s1);CHKERRQ(ierr);
  ierr = VecDestroy(&s2);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

/*@
      MatSchurComplementSet - Sets the matrices that define the Schur complement

   Collective on Mat

   Input Parameter:
+   N - matrix obtained with MatCreate() and MatSetType(MATSCHURCOMPLEMENT);
-   A00,A01,A10,A11  - the four parts of the original matrix (A00 is optional)

   Level: intermediate

   Notes: The Schur complement is NOT actually formed! Rather this
          object performs the matrix-vector product by using the the formula for
          the Schur complement and a KSP solver to approximate the action of inv(A)

          All four matrices must have the same MPI communicator

          A00 and  A11 must be square matrices

.seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatSchurComplementUpdate(), MatCreateTranspose(), MatGetSchurComplement()

@*/
PetscErrorCode  MatSchurComplementSet(Mat N,Mat A00,Mat Ap00,Mat A01,Mat A10,Mat A11)
{
  PetscErrorCode       ierr;
  PetscInt             m,n;
  Mat_SchurComplement  *Na = (Mat_SchurComplement*)N->data;

  PetscFunctionBegin;
  if (N->assembled) SETERRQ(((PetscObject)N)->comm,PETSC_ERR_ARG_WRONGSTATE,"Use MatSchurComplementUpdate() for already used matrix");
  PetscValidHeaderSpecific(A00,MAT_CLASSID,1);
  PetscValidHeaderSpecific(Ap00,MAT_CLASSID,2);
  PetscValidHeaderSpecific(A01,MAT_CLASSID,3);
  PetscValidHeaderSpecific(A10,MAT_CLASSID,4);
  PetscCheckSameComm(A00,1,Ap00,2);
  PetscCheckSameComm(A00,1,A01,3);
  PetscCheckSameComm(A00,1,A10,4);
  if (A00->rmap->n != A00->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local rows of A00 %D do not equal local columns %D",A00->rmap->n,A00->cmap->n);
  if (A00->rmap->n != Ap00->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local rows of A00 %D do not equal local rows of Ap00 %D",A00->rmap->n,Ap00->rmap->n);
  if (Ap00->rmap->n != Ap00->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local rows of Ap00 %D do not equal local columns %D",Ap00->rmap->n,Ap00->cmap->n);
  if (A00->cmap->n != A01->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local columns of A00 %D do not equal local rows of A01 %D",A00->cmap->n,A01->rmap->n);
  if (A10->cmap->n != A00->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local columns of A10 %D do not equal local rows of A00 %D",A10->cmap->n,A00->rmap->n);
  if (A11) {
    PetscValidHeaderSpecific(A11,MAT_CLASSID,5);
    PetscCheckSameComm(A00,1,A11,5);
    if (A11->rmap->n != A11->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local rows of A11 %D do not equal local columns %D",A11->rmap->n,A11->cmap->n);
    if (A10->rmap->n != A11->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local rows of A10 %D do not equal local rows A11 %D",A10->rmap->n,A11->rmap->n);
  }

  ierr = MatGetLocalSize(A01,PETSC_NULL,&n);CHKERRQ(ierr);
  ierr = MatGetLocalSize(A10,&m,PETSC_NULL);CHKERRQ(ierr);
  ierr = MatSetSizes(N,m,n,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
  ierr      = PetscObjectReference((PetscObject)A00);CHKERRQ(ierr);
  ierr      = PetscObjectReference((PetscObject)Ap00);CHKERRQ(ierr);
  ierr      = PetscObjectReference((PetscObject)A01);CHKERRQ(ierr);
  ierr      = PetscObjectReference((PetscObject)A10);CHKERRQ(ierr);
  Na->A     = A00;
  Na->Ap    = Ap00;
  Na->B     = A01;
  Na->C     = A10;
  Na->D     = A11;
  if (A11) {
    ierr = PetscObjectReference((PetscObject)A11);CHKERRQ(ierr);
  }
  N->assembled           = PETSC_TRUE;
  N->preallocated        = PETSC_TRUE;

  ierr = PetscLayoutSetUp((N)->rmap);CHKERRQ(ierr);
  ierr = PetscLayoutSetUp((N)->cmap);CHKERRQ(ierr);
  ierr = KSPSetOperators(Na->ksp,A00,Ap00,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

static void ApplyPCPrecPETSCGen(void *x, PRIMME_INT *ldx, void *y,
      PRIMME_INT *ldy, int *blockSize, int trans, PC *pc, MPI_Comm comm) {
   int i;
   Vec xvec, yvec;
   Mat matrix;
   PetscErrorCode ierr;
   PetscInt mLocal, nLocal;
   
   ierr = PCGetOperators(pc[0],&matrix,NULL); CHKERRABORT(comm, ierr);

   assert(sizeof(PetscScalar) == sizeof(SCALAR));   
   ierr = MatGetLocalSize(matrix, &mLocal, &nLocal); CHKERRABORT(comm, ierr);
   assert(mLocal == nLocal && nLocal <= *ldx && mLocal <= *ldy);

   #if PETSC_VERSION_LT(3,6,0)
      ierr = MatGetVecs(matrix, &xvec, &yvec); CHKERRABORT(comm, ierr);
   #else
      ierr = MatCreateVecs(matrix, &xvec, &yvec); CHKERRABORT(comm, ierr);
   #endif
   for (i=0; i<*blockSize; i++) {
      ierr = VecPlaceArray(xvec, ((SCALAR*)x) + (*ldx)*i); CHKERRABORT(comm, ierr);
      ierr = VecPlaceArray(yvec, ((SCALAR*)y) + (*ldy)*i); CHKERRABORT(comm, ierr);
      if (trans == 0) {
         ierr = PCApply(*pc, xvec, yvec); CHKERRABORT(comm, ierr);
      } else if (pc[1]) {
         ierr = PCApply(pc[1], xvec, yvec); CHKERRABORT(comm, ierr);
      } else {
         ierr = PCApplyTranspose(pc[0], xvec, yvec);
      }
      ierr = VecResetArray(xvec); CHKERRABORT(comm, ierr);
      ierr = VecResetArray(yvec); CHKERRABORT(comm, ierr);
   }
   ierr = VecDestroy(&xvec); CHKERRABORT(comm, ierr);
   ierr = VecDestroy(&yvec); CHKERRABORT(comm, ierr);
}

static void PETScMatvecGenNoBlock(void *x, PRIMME_INT ldx, void *y, PRIMME_INT ldy,
      int blockSize, int trans, Mat matrix, MPI_Comm comm) {
   int i;
   Vec xvec, yvec;
   PetscInt m, n, mLocal, nLocal;
   PetscErrorCode ierr;
   
   assert(sizeof(PetscScalar) == sizeof(SCALAR));   
   ierr = MatGetSize(matrix, &m, &n); CHKERRABORT(comm, ierr);
   ierr = MatGetLocalSize(matrix, &mLocal, &nLocal); CHKERRABORT(comm, ierr);

   #if PETSC_VERSION_LT(3,6,0)
      ierr = MatGetVecs(matrix, &xvec, &yvec); CHKERRABORT(comm, ierr);
   #else
      ierr = MatCreateVecs(matrix, &xvec, &yvec); CHKERRABORT(comm, ierr);
   #endif
   if (trans == 1) {
      Vec aux = xvec; xvec = yvec; yvec = aux;
   }
   for (i=0; i<blockSize; i++) {
      ierr = VecPlaceArray(xvec, ((SCALAR*)x) + ldx*i); CHKERRABORT(comm, ierr);
      ierr = VecPlaceArray(yvec, ((SCALAR*)y) + ldy*i); CHKERRABORT(comm, ierr);
      if (trans == 0) {
         ierr = MatMult(matrix, xvec, yvec); CHKERRABORT(comm, ierr);
      } else {
         ierr = MatMultHermitianTranspose(matrix, xvec, yvec); CHKERRABORT(comm, ierr);
      }
      ierr = VecResetArray(xvec); CHKERRABORT(comm, ierr);
      ierr = VecResetArray(yvec); CHKERRABORT(comm, ierr);
   }
   ierr = VecDestroy(&xvec); CHKERRABORT(comm, ierr);
   ierr = VecDestroy(&yvec); CHKERRABORT(comm, ierr);
}

static PetscErrorCode MatCoarsenApply_MIS(MatCoarsen coarse)
{
    /* MatCoarsen_MIS *MIS = (MatCoarsen_MIS*)coarse->; */
    PetscErrorCode ierr;
    Mat            mat = coarse->graph;

    PetscFunctionBegin;
    PetscValidHeaderSpecific(coarse,MAT_COARSEN_CLASSID,1);
    if (!coarse->perm) {
        IS       perm;
        PetscInt n,m;
        MPI_Comm comm;
        ierr = PetscObjectGetComm((PetscObject)mat,&comm);
        CHKERRQ(ierr);
        ierr = MatGetLocalSize(mat, &m, &n);
        CHKERRQ(ierr);
        ierr = ISCreateStride(comm, m, 0, 1, &perm);
        CHKERRQ(ierr);
        ierr = maxIndSetAgg(perm, mat, coarse->strict_aggs, &coarse->agg_lists);
        CHKERRQ(ierr);
        ierr = ISDestroy(&perm);
        CHKERRQ(ierr);
    } else {
        ierr = maxIndSetAgg(coarse->perm, mat, coarse->strict_aggs,  &coarse->agg_lists);
        CHKERRQ(ierr);
    }
    PetscFunctionReturn(0);
}

/*@
      MatCreateNormalHermitian - Creates a new matrix object that behaves like (A*)'*A.

   Collective on Mat

   Input Parameter:
.   A  - the (possibly rectangular complex) matrix

   Output Parameter:
.   N - the matrix that represents (A*)'*A

   Level: intermediate

   Notes: The product (A*)'*A is NOT actually formed! Rather the new matrix
          object performs the matrix-vector product by first multiplying by
          A and then (A*)'
@*/
PetscErrorCode  MatCreateNormalHermitian(Mat A,Mat *N)
{
  PetscErrorCode ierr;
  PetscInt       m,n;
  Mat_Normal     *Na;

  PetscFunctionBegin;
  ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr);
  ierr = MatCreate(PetscObjectComm((PetscObject)A),N);CHKERRQ(ierr);
  ierr = MatSetSizes(*N,n,n,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
  ierr = PetscObjectChangeTypeName((PetscObject)*N,MATNORMALHERMITIAN);CHKERRQ(ierr);

  ierr       = PetscNewLog(*N,&Na);CHKERRQ(ierr);
  (*N)->data = (void*) Na;
  ierr       = PetscObjectReference((PetscObject)A);CHKERRQ(ierr);
  Na->A      = A;
  Na->scale  = 1.0;

  ierr = VecCreateMPI(PetscObjectComm((PetscObject)A),m,PETSC_DECIDE,&Na->w);CHKERRQ(ierr);

  (*N)->ops->destroy          = MatDestroyHermitian_Normal;
  (*N)->ops->mult             = MatMultHermitian_Normal;
  (*N)->ops->multtranspose    = MatMultHermitianTranspose_Normal;
  (*N)->ops->multtransposeadd = MatMultHermitianTransposeAdd_Normal;
  (*N)->ops->multadd          = MatMultHermitianAdd_Normal;
  (*N)->ops->getdiagonal      = MatGetDiagonalHermitian_Normal;
  (*N)->ops->scale            = MatScaleHermitian_Normal;
  (*N)->ops->diagonalscale    = MatDiagonalScaleHermitian_Normal;
  (*N)->assembled             = PETSC_TRUE;
  (*N)->cmap->N               = A->cmap->N;
  (*N)->rmap->N               = A->cmap->N;
  (*N)->cmap->n               = A->cmap->n;
  (*N)->rmap->n               = A->cmap->n;
  PetscFunctionReturn(0);
}

int main(int argc,char **args)
{
  Mat            A,B,MA;
  PetscViewer    fd;
  char           file[PETSC_MAX_PATH_LEN];
  PetscInt       m,n,M,N,dof=1;
  PetscMPIInt    rank,size;
  PetscErrorCode ierr;
  PetscBool      flg;

  PetscInitialize(&argc,&args,(char*)0,help);
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);

#if defined(PETSC_USE_COMPLEX)
  SETERRQ(PETSC_COMM_WORLD,1,"This example does not work with complex numbers");
#else

  /* Load aij matrix A */
  ierr = PetscOptionsGetString(NULL,"-f",file,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr);
  if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Must indicate binary file with the -f option");
  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);CHKERRQ(ierr);
  ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
  ierr = MatLoad(A,fd);CHKERRQ(ierr);
  ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr);

  /* Get dof, then create maij matrix MA */
  ierr = PetscOptionsGetInt(NULL,"-dof",&dof,NULL);CHKERRQ(ierr);
  ierr = MatCreateMAIJ(A,dof,&MA);CHKERRQ(ierr);
  ierr = MatGetLocalSize(MA,&m,&n);CHKERRQ(ierr);
  ierr = MatGetSize(MA,&M,&N);CHKERRQ(ierr);

  if (size == 1) {
    ierr = MatConvert(MA,MATSEQAIJ,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr);
  } else {
    ierr = MatConvert(MA,MATMPIAIJ,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr);
  }

  /* Test MatMult() */
  ierr = MatMultEqual(MA,B,10,&flg);CHKERRQ(ierr);
  if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_CONV_FAILED,"Error: MatMul() for MAIJ matrix");
  /* Test MatMultAdd() */
  ierr = MatMultAddEqual(MA,B,10,&flg);CHKERRQ(ierr);
  if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_CONV_FAILED,"Error: MatMulAdd() for MAIJ matrix");

  /* Test MatMultTranspose() */
  ierr = MatMultTransposeEqual(MA,B,10,&flg);CHKERRQ(ierr);
  if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_CONV_FAILED,"Error: MatMulAdd() for MAIJ matrix");

  /* Test MatMultTransposeAdd() */
  ierr = MatMultTransposeAddEqual(MA,B,10,&flg);CHKERRQ(ierr);
  if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_CONV_FAILED,"Error: MatMulTransposeAdd() for MAIJ matrix");

  ierr = MatDestroy(&MA);CHKERRQ(ierr);
  ierr = MatDestroy(&A);CHKERRQ(ierr);
  ierr = MatDestroy(&B);CHKERRQ(ierr);
  ierr = PetscFinalize();
#endif
  return 0;
}

PetscErrorCode cHamiltonianMatrix::hamiltonianConstruction(){
	ierr = MatCreate(PETSC_COMM_WORLD,&Hpolaron);CHKERRQ(ierr);
	  ierr = MatSetType(Hpolaron,MATMPIAIJ);CHKERRQ(ierr);
	  ierr = MatSetSizes(Hpolaron,PETSC_DECIDE,PETSC_DECIDE,DIM,DIM);CHKERRQ(ierr);
	  // TODO: should be able to set the symmetric/hermitian option and
	  // only do upper-right triangle part of matrix construction .
	  // and perform corresponding operations thereon.
	  // ierr = MatSetOption(Hpolaron,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr);
	  // ierr = MatSetOption(Hpolaron,MAT_HERMITIAN,PETSC_TRUE);CHKERRQ(ierr);

	  // TODO: what is the estimate of the pre-allocation?
	  // -- number of nonzeros per row in DIAGONAL portion of local submatrix
	  // (same value is used for all local rows) ? I put dim temporarily here.
	  // number of nonzeros per row in the OFF-DIAGONAL portion of local submatrix
	  // (same value is used for all local rows) ?  I put dim temporarily here..
	  // More details at http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatMPIAIJSetPreallocation.html
	  ierr = MatMPIAIJSetPreallocation(Hpolaron,DIM,NULL,DIM,NULL);CHKERRQ(ierr);
	  ierr = MatSeqAIJSetPreallocation(Hpolaron,DIM,NULL);CHKERRQ(ierr);

	  ierr = MatGetOwnershipRange(Hpolaron,&rstart,&rend);CHKERRQ(ierr);
	  ierr = MatGetLocalSize(Hpolaron,&nlocal, NULL);CHKERRQ(ierr);

	  ierr = assemblance();CHKERRQ(ierr);

	  ierr = MatAssemblyBegin(Hpolaron,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
      ierr = MatAssemblyEnd(Hpolaron,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

//	  ierr = PetscViewerSetFormat(PETSC_VIEWER_STDOUT_WORLD,	PETSC_VIEWER_ASCII_DENSE  );CHKERRQ(ierr);
//      ierr = PetscViewerSetFormat(PETSC_VIEWER_STDOUT_WORLD,	PETSC_VIEWER_ASCII_MATLAB  );CHKERRQ(ierr);
//	  ierr = MatView(Hpolaron,	PETSC_VIEWER_STDOUT_WORLD );CHKERRQ(ierr);

	  return ierr;
}

static PetscErrorCode PCSetUp_Eisenstat(PC pc)
{
  PetscErrorCode ierr;
  PetscInt       M,N,m,n;
  PC_Eisenstat   *eis = (PC_Eisenstat*)pc->data;

  PetscFunctionBegin;
  if (!pc->setupcalled) {
    ierr = MatGetSize(pc->mat,&M,&N);CHKERRQ(ierr);
    ierr = MatGetLocalSize(pc->mat,&m,&n);CHKERRQ(ierr);
    ierr = MatCreate(((PetscObject)pc)->comm,&eis->shell);CHKERRQ(ierr);
    ierr = MatSetSizes(eis->shell,m,n,M,N);CHKERRQ(ierr);
    ierr = MatSetType(eis->shell,MATSHELL);CHKERRQ(ierr);
    ierr = MatSetUp(eis->shell);CHKERRQ(ierr);
    ierr = MatShellSetContext(eis->shell,(void*)pc);CHKERRQ(ierr);
    ierr = PetscLogObjectParent(pc,eis->shell);CHKERRQ(ierr);
    ierr = MatShellSetOperation(eis->shell,MATOP_MULT,(void(*)(void))PCMult_Eisenstat);CHKERRQ(ierr);
  }
  if (!eis->usediag) PetscFunctionReturn(0);
  if (!pc->setupcalled) {
    ierr = MatGetVecs(pc->pmat,&eis->diag,0);CHKERRQ(ierr);
    ierr = PetscLogObjectParent(pc,eis->diag);CHKERRQ(ierr);
  }
  ierr = MatGetDiagonal(pc->pmat,eis->diag);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PetscErrorCode MatSetUp_LMVM(Mat B)
{
  Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
  PetscErrorCode    ierr;
  PetscInt          m, n, M, N;
  PetscMPIInt       size;
  MPI_Comm          comm = PetscObjectComm((PetscObject)B);
  
  PetscFunctionBegin;
  ierr = MatGetSize(B, &M, &N);CHKERRQ(ierr);
  if (M == 0 && N == 0) SETERRQ(comm, PETSC_ERR_ORDER, "MatSetSizes() must be called before MatSetUp()");
  if (!lmvm->allocated) {
    ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr);
    if (size == 1) {
      ierr = VecCreateSeq(comm, N, &lmvm->Xprev);CHKERRQ(ierr);
      ierr = VecCreateSeq(comm, M, &lmvm->Fprev);CHKERRQ(ierr);
    } else {
      ierr = MatGetLocalSize(B, &m, &n);CHKERRQ(ierr);
      ierr = VecCreateMPI(comm, n, N, &lmvm->Xprev);CHKERRQ(ierr);
      ierr = VecCreateMPI(comm, m, M, &lmvm->Fprev);CHKERRQ(ierr);
    }
    if (lmvm->m > 0) {
      ierr = VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lmvm->S);CHKERRQ(ierr);
      ierr = VecDuplicateVecs(lmvm->Fprev, lmvm->m, &lmvm->Y);CHKERRQ(ierr);
    }
    lmvm->allocated = PETSC_TRUE;
    B->preallocated = PETSC_TRUE;
    B->assembled = PETSC_TRUE;
  }
  PetscFunctionReturn(0);
}

/*@
      MatCreateTranspose - Creates a new matrix object that behaves like A'

   Collective on Mat

   Input Parameter:
.   A  - the (possibly rectangular) matrix

   Output Parameter:
.   N - the matrix that represents A'

   Level: intermediate

   Notes:
    The transpose A' is NOT actually formed! Rather the new matrix
          object performs the matrix-vector product by using the MatMultTranspose() on
          the original matrix

.seealso: MatCreateNormal(), MatMult(), MatMultTranspose(), MatCreate()

@*/
PetscErrorCode  MatCreateTranspose(Mat A,Mat *N)
{
  PetscErrorCode ierr;
  PetscInt       m,n;
  Mat_Transpose  *Na;

  PetscFunctionBegin;
  ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr);
  ierr = MatCreate(PetscObjectComm((PetscObject)A),N);CHKERRQ(ierr);
  ierr = MatSetSizes(*N,n,m,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
  ierr = PetscLayoutSetUp((*N)->rmap);CHKERRQ(ierr);
  ierr = PetscLayoutSetUp((*N)->cmap);CHKERRQ(ierr);
  ierr = PetscObjectChangeTypeName((PetscObject)*N,MATTRANSPOSEMAT);CHKERRQ(ierr);

  ierr       = PetscNewLog(*N,&Na);CHKERRQ(ierr);
  (*N)->data = (void*) Na;
  ierr       = PetscObjectReference((PetscObject)A);CHKERRQ(ierr);
  Na->A      = A;

  (*N)->ops->destroy          = MatDestroy_Transpose;
  (*N)->ops->mult             = MatMult_Transpose;
  (*N)->ops->multadd          = MatMultAdd_Transpose;
  (*N)->ops->multtranspose    = MatMultTranspose_Transpose;
  (*N)->ops->multtransposeadd = MatMultTransposeAdd_Transpose;
  (*N)->ops->duplicate        = MatDuplicate_Transpose;
  (*N)->ops->getvecs          = MatCreateVecs_Transpose;
  (*N)->ops->axpy             = MatAXPY_Transpose;
  (*N)->assembled             = PETSC_TRUE;

  ierr = PetscObjectComposeFunction((PetscObject)(*N),"MatTransposeGetMat_C",MatTransposeGetMat_Transpose);CHKERRQ(ierr);
  ierr = MatSetBlockSizes(*N,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));CHKERRQ(ierr);
  ierr = MatSetUp(*N);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PetscErrorCode MatGetDiagonal_Brussel(Mat A,Vec diag)
{
  Vec            d1,d2;
  PetscInt       n;
  PetscScalar    *pd;
  MPI_Comm       comm;
  CTX_BRUSSEL    *ctx;
  PetscErrorCode ierr;

  PetscFunctionBeginUser;
  ierr = MatShellGetContext(A,(void**)&ctx);CHKERRQ(ierr);
  ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr);
  ierr = MatGetLocalSize(ctx->T,&n,NULL);CHKERRQ(ierr);
  ierr = VecGetArray(diag,&pd);CHKERRQ(ierr);
  ierr = VecCreateMPIWithArray(comm,1,n,PETSC_DECIDE,pd,&d1);CHKERRQ(ierr);
  ierr = VecCreateMPIWithArray(comm,1,n,PETSC_DECIDE,pd+n,&d2);CHKERRQ(ierr);

  ierr = VecSet(d1,-2.0*ctx->tau1 + ctx->beta - 1.0 + ctx->sigma);CHKERRQ(ierr);
  ierr = VecSet(d2,-2.0*ctx->tau2 - ctx->alpha*ctx->alpha + ctx->sigma);CHKERRQ(ierr);

  ierr = VecDestroy(&d1);CHKERRQ(ierr);
  ierr = VecDestroy(&d2);CHKERRQ(ierr);
  ierr = VecRestoreArray(diag,&pd);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

// -------------------------------------------------------------
// MatTranspose_DenseGA
// -------------------------------------------------------------
static
PetscErrorCode
MatTranspose_DenseGA(Mat mat, MatReuse reuse, Mat *B)
{
  PetscErrorCode ierr = 0;
  
  MPI_Comm comm;
  ierr = PetscObjectGetComm((PetscObject)mat, &comm); CHKERRQ(ierr);

  struct MatGACtx *ctx, *newctx;
  ierr = MatShellGetContext(mat, &ctx); CHKERRQ(ierr);

  PetscInt lrows, grows, lcols, gcols;
  ierr = MatGetSize(mat, &grows, &gcols); CHKERRQ(ierr);
  ierr = MatGetLocalSize(mat, &lrows, &lcols); CHKERRQ(ierr);

  ierr = PetscMalloc(sizeof(struct MatGACtx), &newctx); CHKERRQ(ierr);
  newctx->gaGroup = ctx->gaGroup;
  ierr = CreateMatGA(newctx->gaGroup, lcols, lrows, gcols, grows, &(newctx->ga)); CHKERRQ(ierr);
  GA_Transpose(ctx->ga, newctx->ga);

  ierr = MatCreateShell(comm, lcols, lrows, gcols, grows, newctx, B); CHKERRQ(ierr);
  ierr = MatSetOperations_DenseGA(*B);

  return ierr;
}

/*@
   MatXAIJSetPreallocation - set preallocation for serial and parallel AIJ, BAIJ, and SBAIJ matrices

   Collective on Mat

   Input Arguments:
+  A - matrix being preallocated
.  bs - block size
.  dnnz - number of nonzero blocks per block row of diagonal part of parallel matrix
.  onnz - number of nonzero blocks per block row of off-diagonal part of parallel matrix
.  dnnzu - number of nonzero blocks per block row of upper-triangular part of diagonal part of parallel matrix
-  onnzu - number of nonzero blocks per block row of upper-triangular part of off-diagonal part of parallel matrix

   Level: beginner

.seealso: MatSeqAIJSetPreallocation(), MatMPIAIJSetPreallocation(), MatSeqBAIJSetPreallocation(), MatMPIBAIJSetPreallocation(), MatSeqSBAIJSetPreallocation(), MatMPISBAIJSetPreallocation(),
          PetscSplitOwnership()
@*/
PetscErrorCode MatXAIJSetPreallocation(Mat A,PetscInt bs,const PetscInt dnnz[],const PetscInt onnz[],const PetscInt dnnzu[],const PetscInt onnzu[])
{
    PetscErrorCode ierr;
    void           (*aij)(void);

    PetscFunctionBegin;
    ierr = MatSetBlockSize(A,bs);
    CHKERRQ(ierr);
    ierr = PetscLayoutSetUp(A->rmap);
    CHKERRQ(ierr);
    ierr = PetscLayoutSetUp(A->cmap);
    CHKERRQ(ierr);
    ierr = MatSeqBAIJSetPreallocation(A,bs,0,dnnz);
    CHKERRQ(ierr);
    ierr = MatMPIBAIJSetPreallocation(A,bs,0,dnnz,0,onnz);
    CHKERRQ(ierr);
    ierr = MatSeqSBAIJSetPreallocation(A,bs,0,dnnzu);
    CHKERRQ(ierr);
    ierr = MatMPISBAIJSetPreallocation(A,bs,0,dnnzu,0,onnzu);
    CHKERRQ(ierr);
    /*
      In general, we have to do extra work to preallocate for scalar (AIJ) matrices so we check whether it will do any
      good before going on with it.
    */
    ierr = PetscObjectQueryFunction((PetscObject)A,"MatMPIAIJSetPreallocation_C",&aij);
    CHKERRQ(ierr);
    if (!aij) {
        ierr = PetscObjectQueryFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",&aij);
        CHKERRQ(ierr);
    }
    if (aij) {
        if (bs == 1) {
            ierr = MatSeqAIJSetPreallocation(A,0,dnnz);
            CHKERRQ(ierr);
            ierr = MatMPIAIJSetPreallocation(A,0,dnnz,0,onnz);
            CHKERRQ(ierr);
        } else {                    /* Convert block-row precallocation to scalar-row */
            PetscInt i,m,*sdnnz,*sonnz;
            ierr = MatGetLocalSize(A,&m,NULL);
            CHKERRQ(ierr);
            ierr = PetscMalloc2((!!dnnz)*m,PetscInt,&sdnnz,(!!onnz)*m,PetscInt,&sonnz);
            CHKERRQ(ierr);
            for (i=0; i<m; i++) {
                if (dnnz) sdnnz[i] = dnnz[i/bs] * bs;
                if (onnz) sonnz[i] = onnz[i/bs] * bs;
            }
            ierr = MatSeqAIJSetPreallocation(A,0,dnnz ? sdnnz : NULL);
            CHKERRQ(ierr);
            ierr = MatMPIAIJSetPreallocation(A,0,dnnz ? sdnnz : NULL,0,onnz ? sonnz : NULL);
            CHKERRQ(ierr);
            ierr = PetscFree2(sdnnz,sonnz);
            CHKERRQ(ierr);
        }
    }
    PetscFunctionReturn(0);
}

HYPRE_Int EuclidReadLocalNz(void *Ain)
{
  START_FUNC_DH
  Mat A = (Mat)Ain;
  HYPRE_Int m, n, ierr;

  ierr = MatGetLocalSize(Ain, &m, &n); 
  if (ierr) SET_ERROR(-1, "PETSc::MatGetLocalSize failed!\n");
  END_FUNC_VAL(m)
}