C++ MatGetRow函数代码示例

PetscErrorCode MatAXPY_Basic(Mat Y,PetscScalar a,Mat X,MatStructure str)
{
  PetscInt          i,start,end,j,ncols,m,n;
  PetscErrorCode    ierr;
  const PetscInt    *row;
  PetscScalar       *val;
  const PetscScalar *vals;

  PetscFunctionBegin;
  ierr = MatGetSize(X,&m,&n);CHKERRQ(ierr);
  ierr = MatGetOwnershipRange(X,&start,&end);CHKERRQ(ierr);
  if (a == 1.0) {
    for (i = start; i < end; i++) {
      ierr = MatGetRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr);
      ierr = MatSetValues(Y,1,&i,ncols,row,vals,ADD_VALUES);CHKERRQ(ierr);
      ierr = MatRestoreRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr);
    }
  } else {
    ierr = PetscMalloc((n+1)*sizeof(PetscScalar),&val);CHKERRQ(ierr);
    for (i=start; i<end; i++) {
      ierr = MatGetRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr);
      for (j=0; j<ncols; j++) {
        val[j] = a*vals[j];
      }
      ierr = MatSetValues(Y,1,&i,ncols,row,val,ADD_VALUES);CHKERRQ(ierr);
      ierr = MatRestoreRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr);
    }
    ierr = PetscFree(val);CHKERRQ(ierr);
  }
  ierr = MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

int testBSplineSetD2R1Mat() {
  PrintTimeStamp(PETSC_COMM_SELF, "D2", NULL);

  MPI_Comm comm = PETSC_COMM_SELF;
  BPS bps; BPSCreate(comm, &bps); BPSSetLine(bps, 5.0, 6);
  int order = 3;
  BSS bss; BSSCreate(comm, &bss); BSSSetKnots(bss, order, bps);
  BSSSetUp(bss);

  // compute matrix
  Mat M; BSSCreateR1Mat(bss, &M); BSSD2R1Mat(bss, M);

  // value check
  const PetscScalar *row;
  MatGetRow(M, 1, NULL, NULL, &row);
  ASSERT_DOUBLE_NEAR(0.1666666666, row[0], pow(10.0, -10.0));
  ASSERT_DOUBLE_NEAR(-1.0, row[1], pow(10.0, -10.0));
  ASSERT_DOUBLE_NEAR(0.3333333333, row[2], pow(10.0, -10.0));
  ASSERT_DOUBLE_NEAR(0.16666666666, row[3], pow(10.0, -10.0));
  MatRestoreRow(M, 1, NULL, NULL, &row);

  MatGetRow(M, 2, NULL, NULL, &row);
  ASSERT_DOUBLE_NEAR(0.1666666666666, row[0], pow(10.0, -10.0));
  ASSERT_DOUBLE_NEAR(0.3333333333333, row[1], pow(10.0, -10.0));
  ASSERT_DOUBLE_NEAR(-1.0, row[2] , pow(10.0, -10.0));
  ASSERT_DOUBLE_NEAR(0.3333333333333, row[3], pow(10.0, -10.0));
  ASSERT_DOUBLE_NEAR(0.1666666666666, row[4], pow(10.0, -10.0));
  MatRestoreRow(M, 2, NULL, NULL, &row);

  // Finalize
  MatDestroy(&M);
  BSSDestroy(&bss);

  return 0;
}

TEST(schur, complement) {
    Profiler::initialize();
   

	IS set;
	// vytvorit rozdeleni bloku na procesory ve tvaru "part" (tj. indexy prvnich radku na procesorech)
    int np, rank;

    MPI_Comm_size(PETSC_COMM_WORLD, &np);
    MPI_Comm_rank(PETSC_COMM_WORLD, &rank);

    Distribution ds(block_size, MPI_COMM_WORLD);
    Distribution block_ds(block_count, MPI_COMM_WORLD);
    Distribution all_ds(block_size + block_count, MPI_COMM_WORLD);
    /*if (rank == 0) {
        cout << all_ds;
        cout << ds;
        cout << block_ds;
    }*/

	ISCreateStride(PETSC_COMM_WORLD, ds.lsize(), all_ds.begin(), 1, &set);
	ISView(set, PETSC_VIEWER_STDOUT_WORLD);

    // volat s lokalni velkosti = pocet radku na lokalnim proc.
	SchurComplementTest * schurComplement = new SchurComplementTest(set, &all_ds);
	schurComplement->set_solution(NULL);
	schurComplement->set_positive_definite();
	schurComplement->start_allocation();
	schurComplement->fill_matrix( rank, ds, block_ds); // preallocate matrix
	schurComplement->start_add_assembly();
	schurComplement->fill_matrix( rank, ds, block_ds); // fill matrix
	schurComplement->finish_assembly();
	MatView(*(schurComplement->get_matrix()),PETSC_VIEWER_STDOUT_WORLD);

	LinSys * lin_sys = new LinSysPetscTest( schurComplement->make_complement_distribution() );
	schurComplement->set_complement( (LinSys_PETSC *)lin_sys );
	schurComplement->solve();

	// test of computed values
	{
		PetscInt ncols;
		const PetscInt *cols;
		const PetscScalar *vals;
		for (unsigned int i=0; i<block_size; i++) {
			MatGetRow(schurComplement->get_a_inv(), i + rank*block_size, &ncols, &cols, &vals);
			EXPECT_FLOAT_EQ( (1.0 / (double)(rank + 2)), vals[i] );
			MatRestoreRow(schurComplement->get_a_inv(), i + rank*block_size, &ncols, &cols, &vals);
		}
		MatGetRow(*(schurComplement->get_system()->get_matrix()), rank, &ncols, &cols, &vals);
		EXPECT_FLOAT_EQ( ((double)block_size / (double)(rank + 2)), vals[0] );
		MatRestoreRow(*(schurComplement->get_system()->get_matrix()), rank, &ncols, &cols, &vals);
	}
}

int testBSplineSetSR1Mat() {
  PrintTimeStamp(PETSC_COMM_SELF, "S", NULL);

  MPI_Comm comm = PETSC_COMM_SELF;
  BPS bps; BPSCreate(comm, &bps); BPSSetLine(bps, 5.0, 6);
  int order = 3;
  BSS bss; BSSCreate(comm, &bss); BSSSetKnots(bss, order, bps);
  BSSSetUp(bss);

  // compute S matrix
  Mat S;
  PetscErrorCode ierr;
  ierr = BSSCreateR1Mat(bss, &S);
  ierr = BSSSR1Mat(bss, S);  CHKERRQ(ierr);

  // Get structure and check
  int m, n;
  ierr = MatGetSize(S, &m, &n); CHKERRQ(ierr);
  ASSERT_EQ(5, m);
  ASSERT_EQ(5, n);

  // value check
  const PetscScalar *row;
  PetscInt ncols;
  const PetscInt *cols;
  ierr =   MatGetRow(S, 0, &ncols, &cols, &row); CHKERRQ(ierr);
  ASSERT_EQ(3, ncols);
  ASSERT_EQ(0, cols[0]);
  ASSERT_DOUBLE_EQ(1.0/3.0, row[0]);
  ASSERT_EQ(1, cols[1]);
  ASSERT_DOUBLE_NEAR(0.2083333333, row[1], pow(10.0, -10.0));
  ASSERT_EQ(2, cols[2]);
  ASSERT_DOUBLE_NEAR(0.0083333333, row[2], pow(10.0, -10.0));
  ierr = MatRestoreRow(S, 0, &ncols, &cols, &row); CHKERRQ(ierr);

  ierr = MatGetRow(S, 2, &ncols, &cols, &row); CHKERRQ(ierr);
  ASSERT_EQ(5, ncols);
  for(int i = 0; i < 5; i++)
    ASSERT_EQ(i, cols[i]);
  ASSERT_DOUBLE_NEAR(0.0083333333333, row[0], pow(10.0, -10.0));
  ASSERT_DOUBLE_NEAR(0.2166666666666, row[1], pow(10.0, -10.0));
  ASSERT_DOUBLE_NEAR(0.55, row[2], pow(10.0, -10.0));  
  ASSERT_DOUBLE_NEAR(0.2166666666666, row[3], pow(10.0, -10.0));
  ASSERT_DOUBLE_NEAR(0.0083333333, row[4], pow(10.0, -10.0));
  ierr = MatRestoreRow(S, 3, &ncols, &cols, &row); CHKERRQ(ierr);

  // Finalize
  MatDestroy(&S);
  BSSDestroy(&bss);

  return 0;
}

/*@
  MatChop - Set all values in the matrix less than the tolerance to zero

  Input Parameters:
+ A   - The matrix
- tol - The zero tolerance

  Output Parameters:
. A - The chopped matrix

  Level: intermediate

.seealso: MatCreate(), MatZeroEntries()
 @*/
PetscErrorCode MatChop(Mat A, PetscReal tol)
{
    PetscScalar    *newVals;
    PetscInt       *newCols;
    PetscInt       rStart, rEnd, numRows, maxRows, r, colMax = 0;
    PetscErrorCode ierr;

    PetscFunctionBegin;
    ierr = MatGetOwnershipRange(A, &rStart, &rEnd);
    CHKERRQ(ierr);
    for (r = rStart; r < rEnd; ++r) {
        PetscInt ncols;

        ierr   = MatGetRow(A, r, &ncols, NULL, NULL);
        CHKERRQ(ierr);
        colMax = PetscMax(colMax, ncols);
        CHKERRQ(ierr);
        ierr   = MatRestoreRow(A, r, &ncols, NULL, NULL);
        CHKERRQ(ierr);
    }
    numRows = rEnd - rStart;
    ierr    = MPI_Allreduce(&numRows, &maxRows, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)A));
    CHKERRQ(ierr);
    ierr    = PetscMalloc2(colMax,&newCols,colMax,&newVals);
    CHKERRQ(ierr);
    for (r = rStart; r < rStart+maxRows; ++r) {
        const PetscScalar *vals;
        const PetscInt    *cols;
        PetscInt           ncols, newcols, c;

        if (r < rEnd) {
            ierr = MatGetRow(A, r, &ncols, &cols, &vals);
            CHKERRQ(ierr);
            for (c = 0; c < ncols; ++c) {
                newCols[c] = cols[c];
                newVals[c] = PetscAbsScalar(vals[c]) < tol ? 0.0 : vals[c];
            }
            newcols = ncols;
            ierr = MatRestoreRow(A, r, &ncols, &cols, &vals);
            CHKERRQ(ierr);
            ierr = MatSetValues(A, 1, &r, newcols, newCols, newVals, INSERT_VALUES);
            CHKERRQ(ierr);
        }
        ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
        CHKERRQ(ierr);
        ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
        CHKERRQ(ierr);
    }
    ierr = PetscFree2(newCols,newVals);
    CHKERRQ(ierr);
    PetscFunctionReturn(0);
}

int Epetra_PETScAIJMatrix::NumMyRowEntries(int Row, int & NumEntries) const 
{
  int globalRow = PetscRowStart_ + Row;
  MatGetRow(Amat_, globalRow, &NumEntries, PETSC_NULL, PETSC_NULL);
  MatRestoreRow(Amat_,globalRow,&NumEntries, PETSC_NULL, PETSC_NULL);
  return(0);
}

void PETSc::Get_row_of_A(PetscInt i, PetscInt *ncol, PetscInt **cols, double **row)
{	
	ierr = MatGetRow(A,i,ncol,(const PetscInt**)cols,
			(const PetscScalar**)row);
	ierr = MatRestoreRow(A,i,ncol,(const PetscInt**)cols,
			(const PetscScalar**)row);
}

void PetscSparseStorage::fromGetSpRow( int row, int col,
				       double A[], int /*lenA */,
				       int jcolA[], int& nnz,
				       int colExtent, int& info )
{
  const double * B;
  const int    * jcolB;
  int k, ierr; 

  int m, n;
  this->getSize(m,n);
  ierr = MatGetRow(M, row, &nnz, &jcolB, &B ); assert( ierr == 0 );
  if( col == 0 && colExtent == n ) {
    // We want the whole row, just do a simple copy
    for( k = 0; k < nnz; k++ ) {
      A[k]     = B[k];
      jcolA[k] = jcolB[k];
    }
  } else { // Copy only those in range
    int i = 0;
    for( k = 0; k < nnz; k++ ) {
      if( jcolB[k] >= col && jcolB[k] < col + colExtent ) {
	A[i]     = B[k];
	jcolA[i] = jcolB[k];
	i++;
      }
    }
  }
  ierr = MatRestoreRow( M, row, &nnz, &jcolB, &B ); assert( ierr == 0 );
  info = 0;
}

int testBSplineSetENMatR1Mat() {
  PrintTimeStamp(PETSC_COMM_SELF, "EN", NULL);
 
  MPI_Comm comm = PETSC_COMM_SELF;
  BPS bps; BPSCreate(comm, &bps); BPSSetLine(bps, 5.0, 6);
  int order = 3;
  BSS bss; BSSCreate(comm, &bss); BSSSetKnots(bss, order, bps);
  BSSSetUp(bss);

  // compute matrix
  Mat M; BSSCreateR1Mat(bss, &M); BSSENR1Mat(bss, 2, 0.7, M);

  // value check
  const PetscScalar *row;
  MatGetRow(M, 4, NULL, NULL, &row);
  ASSERT_DOUBLE_NEAR(0.0000964408077, row[0], pow(10.0, -8));
  ASSERT_DOUBLE_NEAR(0.00168615366, row[1], pow(10.0, -8));
  ASSERT_DOUBLE_NEAR(0.00211136230, row[2], pow(10.0, -8));
  MatRestoreRow(M, 4, NULL, NULL, &row);

  // Finalize
  MatDestroy(&M);
  BSSDestroy(&bss);

  return 0;
}

/*
     MatDumpSPAI - Dumps a PETSc matrix to a file in an ASCII format
  suitable for the SPAI code of Stephen Barnard to solve. This routine
  is simply here to allow testing of matrices directly with the SPAI
  code, rather then through the PETSc interface.

*/
PetscErrorCode  MatDumpSPAI(Mat A,FILE *file)
{
  const PetscScalar *vals;
  PetscErrorCode ierr;
  int               i,j,n,size,nz;
  const int         *cols;
  MPI_Comm          comm;

  PetscObjectGetComm((PetscObject)A,&comm);

  MPI_Comm_size(comm,&size);
  if (size > 1) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_SUP,"Only single processor dumps");

  ierr = MatGetSize(A,&n,&n);CHKERRQ(ierr);

  /* print the matrix */
  fprintf(file,"%d\n",n);
  for (i=0; i<n; i++) {
    ierr     = MatGetRow(A,i,&nz,&cols,&vals);CHKERRQ(ierr);
    for (j=0; j<nz; j++) {
      fprintf(file,"%d %d %16.14e\n",i+1,cols[j]+1,vals[j]);
    }
    ierr     = MatRestoreRow(A,i,&nz,&cols,&vals);CHKERRQ(ierr);
  }

  PetscFunctionReturn(0);
}

/*@
  MatComputeBandwidth - Calculate the full bandwidth of the matrix, meaning the width 2k+1 where k diagonals on either side are sufficient to contain all the matrix nonzeros.

  Collective on Mat

  Input Parameters:
+ A - The Mat
- fraction - An optional percentage of the Frobenius norm of the matrix that the bandwidth should enclose

  Output Parameter:
. bw - The matrix bandwidth

  Level: beginner

.seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexSetChart()
@*/
PetscErrorCode MatComputeBandwidth(Mat A, PetscReal fraction, PetscInt *bw)
{
  PetscInt       lbw[2] = {0, 0}, gbw[2];
  PetscInt       rStart, rEnd, r;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(A, MAT_CLASSID, 1);
  PetscValidLogicalCollectiveReal(A,fraction,2);
  PetscValidPointer(bw, 3);
  if ((fraction > 0.0) && (fraction < 1.0)) SETERRQ(PetscObjectComm((PetscObject) A), PETSC_ERR_SUP, "We do not yet support a fractional bandwidth");
  ierr = MatGetOwnershipRange(A, &rStart, &rEnd);CHKERRQ(ierr);
  for (r = rStart; r < rEnd; ++r) {
    const PetscInt *cols;
    PetscInt        ncols;

    ierr = MatGetRow(A, r, &ncols, &cols, NULL);CHKERRQ(ierr);
    if (ncols) {
      lbw[0] = PetscMax(lbw[0], r - cols[0]);
      lbw[1] = PetscMax(lbw[1], cols[ncols-1] - r);
    }
    ierr = MatRestoreRow(A, r, &ncols, &cols, NULL);CHKERRQ(ierr);
  }
  ierr = MPI_Allreduce(lbw, gbw, 2, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject) A));CHKERRQ(ierr);
  *bw = 2*PetscMax(gbw[0], gbw[1]) + 1;
  PetscFunctionReturn(0);
}

PetscErrorCode MatGetDiagonalHermitian_Normal(Mat N,Vec v)
{
  Mat_Normal        *Na = (Mat_Normal*)N->data;
  Mat               A   = Na->A;
  PetscErrorCode    ierr;
  PetscInt          i,j,rstart,rend,nnz;
  const PetscInt    *cols;
  PetscScalar       *diag,*work,*values;
  const PetscScalar *mvalues;

  PetscFunctionBegin;
  ierr = PetscMalloc2(A->cmap->N,&diag,A->cmap->N,&work);CHKERRQ(ierr);
  ierr = PetscMemzero(work,A->cmap->N*sizeof(PetscScalar));CHKERRQ(ierr);
  ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr);
  for (i=rstart; i<rend; i++) {
    ierr = MatGetRow(A,i,&nnz,&cols,&mvalues);CHKERRQ(ierr);
    for (j=0; j<nnz; j++) {
      work[cols[j]] += mvalues[j]*PetscConj(mvalues[j]);
    }
    ierr = MatRestoreRow(A,i,&nnz,&cols,&mvalues);CHKERRQ(ierr);
  }
  ierr   = MPIU_Allreduce(work,diag,A->cmap->N,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)N));CHKERRQ(ierr);
  rstart = N->cmap->rstart;
  rend   = N->cmap->rend;
  ierr   = VecGetArray(v,&values);CHKERRQ(ierr);
  ierr   = PetscMemcpy(values,diag+rstart,(rend-rstart)*sizeof(PetscScalar));CHKERRQ(ierr);
  ierr   = VecRestoreArray(v,&values);CHKERRQ(ierr);
  ierr   = PetscFree2(diag,work);CHKERRQ(ierr);
  ierr   = VecScale(v,Na->scale);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

void
PetscSparseMtrx :: times(const FloatMatrix &B, FloatMatrix &answer) const
{
    if ( this->giveNumberOfColumns() != B.giveNumberOfRows() ) {
        OOFEM_ERROR("Dimension mismatch");
    }

#ifdef __PARALLEL_MODE
    if ( emodel->isParallel() ) {
        OOFEM_ERROR("PetscSparseMtrx :: times - Not implemented");
    }
#endif
    // I'm opting to work with a set of vectors, as i think it might be faster and more robust. / Mikael

    int nr = this->giveNumberOfRows();
    int nc = B.giveNumberOfColumns();
    answer.resize(nr, nc);
    double *aptr = answer.givePointer();

#if 0
     // Approach using several vectors. Not sure if it is optimal, but it includes petsc calls which i suspect are inefficient. / Mikael
     // UNTESTED!
     Vec globX, globY;
     VecCreate(PETSC_COMM_SELF, &globY);
     VecSetType(globY, VECSEQ);
     VecSetSizes(globY, PETSC_DECIDE, nr);
     int nrB = B.giveNumberOfRows();
     for (int k = 0; k < nc; k++) {
         double colVals[nrB];
         for (int i = 0; i < nrB; i++) colVals[i] = B(i,k); // B.copyColumn(Bk,k);
         VecCreateSeqWithArray(PETSC_COMM_SELF, nrB, colVals, &globX);
         MatMult(this->mtrx, globX, globY );
                 double *ptr;
                 VecGetArray(globY, &ptr);
                 for (int i = 0; i < nr; i++) *aptr++ = ptr[i]; // answer.setColumn(Ak,k);
                 VecRestoreArray(globY, &ptr);
                 VecDestroy(globX);
     }
     VecDestroy(globY);
#endif

    Mat globB, globC;
    MatCreateSeqDense(PETSC_COMM_SELF, B.giveNumberOfRows(), B.giveNumberOfColumns(), B.givePointer(), & globB);
    MatMatMult(this->mtrx, globB, MAT_INITIAL_MATRIX, PETSC_DEFAULT, & globC);
    const double *vals;
    for ( int r = 0; r < nr; r++ ) {
        MatGetRow(globC, r, NULL, NULL, & vals);
        for ( int i = 0, i2 = r; i < nc; i++, i2 += nr ) {
            aptr [ i2 ] = vals [ i ];
        }
        MatRestoreRow(globC, r, NULL, NULL, & vals);
    }

    MatDestroy(&globB);
    MatDestroy(&globC);
}

PetscErrorCode MatGetRow_SMF(Mat mat,PetscInt row,PetscInt *ncols,const PetscInt **cols,const PetscScalar **vals)
{
  PetscErrorCode   ierr;
  MatSubMatFreeCtx ctx;

  PetscFunctionBegin;
  ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr);
  ierr = MatGetRow(ctx->A,row,ncols,cols,vals);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PetscErrorCode PCGAMGOptProl_Classical_Jacobi(PC pc,const Mat A,Mat *P)
{

  PetscErrorCode    ierr;
  PetscInt          f,s,n,cf,cs,i,idx;
  PetscInt          *coarserows;
  PetscInt          ncols;
  const PetscInt    *pcols;
  const PetscScalar *pvals;
  Mat               Pnew;
  Vec               diag;
  PC_MG             *mg          = (PC_MG*)pc->data;
  PC_GAMG           *pc_gamg     = (PC_GAMG*)mg->innerctx;
  PC_GAMG_Classical *cls         = (PC_GAMG_Classical*)pc_gamg->subctx;

  PetscFunctionBegin;
  if (cls->nsmooths == 0) {
    ierr = PCGAMGTruncateProlongator_Private(pc,P);CHKERRQ(ierr);
    PetscFunctionReturn(0);
  }
  ierr = MatGetOwnershipRange(*P,&s,&f);CHKERRQ(ierr);
  n = f-s;
  ierr = MatGetOwnershipRangeColumn(*P,&cs,&cf);CHKERRQ(ierr);
  ierr = PetscMalloc(sizeof(PetscInt)*n,&coarserows);CHKERRQ(ierr);
  /* identify the rows corresponding to coarse unknowns */
  idx = 0;
  for (i=s;i<f;i++) {
    ierr = MatGetRow(*P,i,&ncols,&pcols,&pvals);CHKERRQ(ierr);
    /* assume, for now, that it's a coarse unknown if it has a single unit entry */
    if (ncols == 1) {
      if (pvals[0] == 1.) {
        coarserows[idx] = i;
        idx++;
      }
    }
    ierr = MatRestoreRow(*P,i,&ncols,&pcols,&pvals);CHKERRQ(ierr);
  }
  ierr = MatGetVecs(A,&diag,0);CHKERRQ(ierr);
  ierr = MatGetDiagonal(A,diag);CHKERRQ(ierr);
  ierr = VecReciprocal(diag);CHKERRQ(ierr);
  for (i=0;i<cls->nsmooths;i++) {
    ierr = MatMatMult(A,*P,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&Pnew);CHKERRQ(ierr);
    ierr = MatZeroRows(Pnew,idx,coarserows,0.,NULL,NULL);CHKERRQ(ierr);
    ierr = MatDiagonalScale(Pnew,diag,0);CHKERRQ(ierr);
    ierr = MatAYPX(Pnew,-1.0,*P,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
    ierr = MatDestroy(P);CHKERRQ(ierr);
    *P  = Pnew;
    Pnew = NULL;
  }
  ierr = VecDestroy(&diag);CHKERRQ(ierr);
  ierr = PetscFree(coarserows);CHKERRQ(ierr);
  ierr = PCGAMGTruncateProlongator_Private(pc,P);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}