prom_petsc.hh

//  $Id: prom_petsc.hh,v 1.34 2004/05/07 20:55:04 adams Exp $
//  Author:             Mark F. Adams
// Copyright (c) 2000 by Mark F. Adams 
//  Filename:           prom_petsc.hh
//
#ifndef __PROM_PETSC_H__
#define __PROM_PETSC_H__

#if defined(PROM_USE_C_PETSC)
extern "C" {
#endif

#include "src/vec/impls/mpi/pvecimpl.h"
  //#include "src/ksp/kspimpl.h"
#include "src/ksp/pc/impls/mg/mgimpl.h"
#include "src/ksp/ksp/impls/cg/cgctx.h"
#include "petscksp.h"
#include "petscmg.h"
#include "petscpc.h"
#include "petscmat.h"
#include "src/mat/impls/baij/mpi/mpibaij.h"
#include "src/mat/impls/sbaij/mpi/mpisbaij.h"
#include "src/mat/impls/aij/mpi/mpiaij.h"

#if defined(PROM_USE_C_PETSC)
}
#endif

typedef IS PromIS;

#include "prom_base.hh"

#define PromISGetSize ISGetLocalSize
#define PromISGetIndices ISGetIndices
#define PromISRestoreIndices ISRestoreIndices
#define PromISCreateGeneral ISCreateGeneral
#define PromISDestroy ISDestroy

#undef __FUNCT__
#define __FUNCT__ "PromComm" 

class PromComm : public PromComm_base
{
public:
  PromComm( MPI_Comm comm ) : mpi_comm_( comm ) {
    assert(mpi_comm_!=MPI_COMM_NULL); assert(mpi_comm_!=MPI_COMM_WORLD);
  }
  MPI_Comm operator()() const { return mpi_comm_; }
  MPI_Comm MPIComm() const { return mpi_comm_; }
  int mype() const { 
    int tt; assert(mpi_comm_ != MPI_COMM_NULL);
    MPI_Comm_rank( mpi_comm_, &tt); 
    return tt; 
  }
  int np()const{ 
    assert(mpi_comm_!=MPI_COMM_NULL); int tt; 
    MPI_Comm_size(mpi_comm_,&tt); 
    return tt;
  }
  int getNewTag() const { 
    int tag; PetscCommGetNewTag( mpi_comm_, &tag ); return tag;
  }
  int getNewTag(MPI_Comm comm) const { 
    int tag; PetscCommGetNewTag( comm, &tag ); return tag;
  }
  const MPI_Comm mpi_comm_;
};

#undef __FUNCT__
#define __FUNCT__ "PromMap" 

class PromMap : public PromMap_base
{
  friend class PromMatrix;
public:
  PromMap( int nloc, const int ElementSize, const PromComm &comm ) 
    : PromMap_base(nloc,ElementSize*nloc,comm), bsize_(ElementSize), 
      id_leq_(NULL) {}
  PromMap( int nloc, const PromComm &comm, const int nloceq ) // hack for VBR
    : PromMap_base(nloc,nloceq,comm), bsize_(-1), id_leq_(NULL) {}
  PromMap( int nloc, const int *ElementSizeList, const PromComm &comm ) 
    :PromMap_base(nloc, sumeqs(nloc,ElementSizeList),comm), bsize_(-1),
     id_leq_(NULL) {
    SetLidLeq( ElementSizeList );
  }
protected:
  int SetLidLeq( const int *ElementSizeList ) {
    int ii,jj,tt,s2,r2,nloc=nlocal(); 
    // determine if constant ndf
    if( nloc > 0 ) {
      for(ii=1,tt=ElementSizeList[0];ii<nloc;ii++) {
        if(ElementSizeList[ii] != tt) break;
      }
      if( ii != nloc ) s2 = tt = 0; // broke
      else s2 = tt; // high / low
    }
    else{ s2 = 0; tt = 999; }
    MPI_Allreduce(&tt, &jj, 1, MPI_INT, MPI_MIN, comm_() );
    MPI_Allreduce(&s2, &r2, 1, MPI_INT, MPI_MAX, comm_() );
    const bool constndf = (jj==r2 && jj>0); const int ndf = jj;
    if( constndf ) bsize_ = ndf;
    else{
      PetscMalloc(sizeof(int)*(nloc+1),&id_leq_); 
      for(ii=0,id_leq_[0] = 0 ; ii < nloc ; ii++ ){
        id_leq_[ii+1] = id_leq_[ii] + ElementSizeList[ii];
        assert(ElementSizeList[ii]>0);
      }
    }
    return 0;
  }

  static int sumeqs(const int nloc, const int *ElementSizeList ){
    int ii,cc;
    for(ii = cc = 0 ; ii < nloc ; ii++ )cc += ElementSizeList[ii];
    return cc;
  }
public:
  virtual ~PromMap() {
    if( id_leq_ != NULL ) PetscFree( id_leq_ ); 
  }
  //
  int nlocal() const {
    assert(this!=NULL);
    int mype = comm_.mype();
    return (proc_gnode_[mype+1] - proc_gnode_[mype]);
  }
  int nglobal() const {
    assert(this!=NULL); 
    int np = comm_.np();
    return (proc_gnode_[np]); 
  } 
  int nlocalEq() const {
    assert(this!=NULL); 
    int mype = comm_.mype();
    return (proc_geq_[mype+1] - proc_geq_[mype]);
  }
  int nglobalEq() const {
    assert(this!=NULL); 
    int np = comm_.np();
    return proc_geq_[np]; 
  }
  bool isConstNDF() const { assert(this!=NULL && bsize_ != 0); return (bsize_>0); }
  int constNDFSize() const { assert(this!=NULL && bsize_>0); return bsize_; }
  int locID_ndf(const int i)const{ 
    assert(this); // this will fail while VBR hack is not fixed
    return (bsize_>0) ? bsize_ : id_leq_[i+1] - id_leq_[i];
  }
private:
  int *id_leq_;
  short int bsize_;
};

#undef __FUNCT__
#define __FUNCT__ "PromVector" 


class PromVector : public PromVector_base
{
  friend class PromMatrix;
  friend class PromPC;
  friend class PromZRP;
  friend class PromSmoother;
  friend class Prometheus;
  friend class pFEAP;
  friend class PromPCGaussSeidel;
  friend class PromPCKKT_ASM;
public:
  PromVector(const PromMap &map, const int fact = 1 ) : 
    vec_(NULL), PromVector_base(map)
  {
    PetscFunctionBegin;
    VecCreateMPI(map.comm_(),map.nlocalEq()*fact,map.nglobalEq()*fact,&vec_);
    assert(vec_ != NULL);
#if defined(PETSC_USE_BOPT_g)
    double b=1.e206;
    VecSet(&b,vec_);
#endif
    PetscStackPop; /* since function returns void cannot use PetscFunctionReturn(); */
  }
  PromVector( const PromVector_base * const vec_b ) : 
    vec_(NULL), PromVector_base( vec_b->map_ ) { 
    const PromVector *vec = dynamic_cast<const PromVector*>(vec_b);
    assert(vec!=NULL);
    PetscFunctionBegin;
    int ierr = VecDuplicate( vec->vec_, &vec_ );       assert(ierr==0);
    PetscStackPop; /* since function returns void cannot use PetscFunctionReturn(); */
  }
  virtual ~PromVector(){ 
    //if( malled_ ) 
    VecDestroy(vec_); vec_ = NULL;
  }
  virtual int GetArray( double **arr ) const {
    assert(this!=NULL);
    return VecGetArray( vec_, arr );
  }
  virtual int RestoreArray( double **arr ) const {
    assert(this!=NULL);
    return VecRestoreArray( vec_, arr );
  }
  virtual double operator[](int i) const { 
    double *arr, tmp; assert(this!=NULL); 
    int ierr = VecGetArray( vec_, &arr ); CHKERRQ(ierr);
    tmp = arr[i];
    ierr = VecRestoreArray( vec_, &arr ); CHKERRQ(ierr);
    return tmp;
  }
  virtual int Set( double val ){
    assert(this!=NULL);
    return VecSet( &val, vec_ );
  }
  virtual int Scale( double val ) {
    assert(this!=NULL);
    return VecScale( &val, vec_ );
  }
  virtual int Sqrt(){
    assert(this!=NULL); // todo!
    int ierr = VecSqrt( vec_ ); CHKERRQ(ierr);
    return 0;
  }
  virtual int Dot(const PromVector_base * const b, double * const dot)const{
    assert(this!=NULL);
    const PromVector *vec = dynamic_cast<const PromVector*>(b);
    assert(vec!=NULL);
    return VecDot( vec_, vec->vec_, dot );
  }
  virtual int SetValues( const int nrow, int *rowp, const double *vals,
                         const int add_type ){
    assert(this!=NULL);
    return VecSetValues( vec_, nrow, rowp, vals, (InsertMode)add_type );
  }
  virtual int Norm2( double * const norm ) const {
    assert(this!=NULL); PetscScalar nn;
    int ierr = VecNorm( vec_, NORM_2, &nn ); *norm = (double)nn;
    return ierr;
  } 
  virtual int NormMax( double * const norm ) const {
    assert(this!=NULL); PetscScalar nn;
    int ierr = VecNorm( vec_, NORM_MAX, &nn ); *norm = (double)nn;
    return ierr;
  } 
  virtual int Copy( const PromVector_base * const pvec ) {
    const PromVector *vec = dynamic_cast<const PromVector*>(pvec);
    assert(vec!=NULL);
    int ierr = VecCopy( vec->vec_, vec_ ); CHKERRQ(ierr); 
    return 0;
  }
  virtual int AXPBY( double a, double b, const PromVector_base *const X ) {
    assert(this!=NULL);
    const PromVector *vec = dynamic_cast<const PromVector*>(X);
    assert(vec!=NULL);
    return VecAXPBY( &a, &b, vec->vec_, vec_ ); 
  }
  virtual int AYPX( double a, const PromVector_base *const X ) {
    const PromVector *vec = dynamic_cast<const PromVector*>(X);
    assert(vec!=NULL);
    int ierr =  VecAYPX( &a, vec->vec_, vec_ ); CHKERRQ(ierr); 
    return 0;
  }
  virtual int AXPY( double a, const PromVector_base *const X ) {
    const PromVector *vec = dynamic_cast<const PromVector*>(X);
    assert(vec!=NULL);
    int ierr = VecAXPY( &a, vec->vec_, vec_ ); CHKERRQ(ierr); 
    return 0;
  }
  virtual int WAXPY( double alpha, const PromVector_base * const X, 
                     const PromVector_base * const Y ) {
    assert(this!=NULL);
    const PromVector *x = dynamic_cast<const PromVector*>(X);assert(x!=NULL);
    const PromVector *y = dynamic_cast<const PromVector*>(Y);assert(y!=NULL);
    return VecWAXPY( &alpha, x->vec_, y->vec_, vec_ ); 
  }
  virtual int MAXPY(const int nv, double alpha[], const PromVector_base *const*Xarr){
    assert(this!=NULL); assert(nv<1000); assert(Xarr!=NULL);
    Vec varr[1000]; 
    for(int i=0;i<nv;i++){
      const PromVector *vv = dynamic_cast<const PromVector*>(Xarr[i]);
      assert(vv!=NULL); varr[i] = vv->vec_;
    }
    return VecMAXPY( nv, alpha, vec_, varr ); 
  }
  virtual int MDot(const int nv, const PromVector_base *const*Xarr, double arr[])const{
    assert(this!=NULL); assert(nv<1000); assert(Xarr!=NULL);
    Vec varr[1000]; 
    for(int i=0;i<nv;i++){
      const PromVector *vv = dynamic_cast<const PromVector*>(Xarr[i]);
      assert(vv!=NULL); varr[i] = vv->vec_;
    }
    return VecMDot( nv, vec_, varr, arr );
  }
  virtual int Assembly(){
    assert(this!=NULL);
    int ierr = VecAssemblyBegin( vec_ ); CHKERRQ(ierr);
    return VecAssemblyEnd( vec_ );
  }
  virtual int PointwiseMult( const PromVector_base * const X,
                             const PromVector_base * const Y){
    assert(this!=NULL);  
    const PromVector *x = dynamic_cast<const PromVector*>(X);assert(x!=NULL);
    const PromVector *y = dynamic_cast<const PromVector*>(Y);assert(y!=NULL);
    return VecPointwiseMult( x->vec_, y->vec_, vec_ );
  } 
  virtual int PointwiseDiv( const PromVector_base *const X, 
                    const PromVector_base *const Y )  {
    assert(this!=NULL);  
    const PromVector *x = dynamic_cast<const PromVector*>(X);assert(x!=NULL);
    const PromVector *y = dynamic_cast<const PromVector*>(Y);assert(y!=NULL);
    return VecPointwiseDivide( x->vec_, y->vec_, vec_ );
  } 
  virtual int Reciprocal()  {
    assert(this!=NULL);  
    return VecReciprocal( vec_ );
  } 
  virtual int getN() const { return vec_->N; }
  virtual int getn() const { return vec_->n; }
  virtual int my0() const { return map_.my0Eq(); }
protected:
  Vec vec_;
};

#undef __FUNCT__
#define __FUNCT__ "PromMatrix" 


class PromMatrix : public PromMatrix_base
{
  friend class PromPC;
  friend class PromMG;
  friend class Prometheus;
  friend class PromContext;
  friend class PromPCGaussSeidel;
  friend class PromPCKKT_ASM;
  friend class PromCRMatrix;
public:
  PromMatrix( const PromGrid *grid, int nA, int *pnA,int nB,int *pnB, int fact=1); 
  PromMatrix( const PromMap &map, int nA, int *pnA, int nB, int *pnB, 
              int lcols = -1, int gcols = -1, int fact=1);  
  virtual int SetUp( int nA,int *pnA,int nB,int *pnB, int lcols, int gcols, int fact=1); 
  PromMatrix( const PromMatrix_base* B ) : 
    mat_(NULL), full_mat_(NULL), varray_(NULL), varr_sz_(0), 
    array_(NULL), ghostGID_BID_(NULL), ghostGID_ndf_(NULL), 
    ghostBID_GID_(NULL), PromMatrix_base( B->map_ ) {
    const PromMatrix *A = dynamic_cast<const PromMatrix*>(B); assert(A!=NULL);
    MatDuplicate( A->mat_, MAT_COPY_VALUES, &mat_ );
  }
  PromMatrix( const PromMap &cmap, Mat a ) : 
    mat_(a), varray_(NULL), full_mat_(NULL), //malled_(TRUE),
    varr_sz_(0), array_(NULL), ghostGID_BID_(NULL), ghostGID_ndf_(NULL),
    PromMatrix_base( cmap ), ghostBID_GID_(NULL) {}
  virtual int Transpose( const PromMap &cmap, PromMatrix_base **out ) const {
    assert(this!=NULL); assert(0); // not used
    int ierr; Mat aa;
    ierr = MatTranspose( mat_, &aa ); CHKERRQ(ierr); 
    PromMatrix *trans = new PromMatrix( cmap, aa );
    *out = trans;
    return 0;
  }
  virtual ~PromMatrix() { 
    //if(malled_) {
    MatDestroy(mat_); mat_ = NULL;
    if( ghostGID_BID_ != NULL ) { 
      delete ghostGID_BID_; delete ghostGID_ndf_; delete ghostBID_GID_;
    }
    if( full_mat_ != NULL )  MatDestroy( full_mat_ );
    if( varray_ != NULL )  PetscFree( varray_ );
    //}
  }
  virtual int ConvertToSymm();
  virtual int ConvertFromSymm();
  virtual int ConvertToSuperLU() {
    if( isSuperLU() ) return 0;
#if defined(PROM_USE_SUPERLU)
    Mat M;
    int ierr = MatConvert( mat_, MATSUPERLU_DIST, &M); CHKERRQ(ierr); 
    MatDestroy(mat_); mat_ = M;
#endif
    return 0;
  }
  virtual bool isSuperLU()const{ 
    return FALSE; 
  }
  int GetLocalNodeRow( const int brow, const int len, int &ncols, int adjacs[],
                       double **Aij = NULL ) const ;
private:
  int GetNodeScalarRow_PETSc( const int row, const int len, int &ncols, 
                              int Brp[], double **Bap ) const;
  int GetLocal_constBlk_NodeRow( const int brow, const int len, int &ncols, 
                                 int adjacs[], double **Aij = NULL ) const ; 
  int GetSymmetricLocal_constBlk_NodeRow( const int brow, const int len, int &ncols,
                                          int adjacs[], double **Aij = NULL ) const ; 
  int GetLocalVBRNodeRow( const int brow, const int len, int &ncols, 
                          int adjacs[], double **Aij = NULL ) const ;
public:
  int GetDiagonal( PromVector_base *work ) const {
    assert(this!=NULL && work != NULL);
    const PromVector *vec = dynamic_cast<const PromVector*>(work);
    assert(vec!=NULL);
    return MatGetDiagonal( mat_, vec->vec_ ); 
  } 
  virtual int MultTranspose( const PromVector_base * const X,
                             PromVector_base * const Y ) const {
    assert(this!=NULL);
    const PromVector *x = dynamic_cast<const PromVector*>(X);assert(x!=NULL);
    const PromVector *y = dynamic_cast<const PromVector*>(Y);assert(y!=NULL);
    return MatMultTranspose( mat_, x->vec_, y->vec_ );
  }
  virtual int MultTransposeAdd( const PromVector_base * const W,
                                const PromVector_base * const X, 
                                PromVector_base * const Y ) const {
    const PromVector *x = dynamic_cast<const PromVector*>(X);assert(x!=NULL);
    const PromVector *y = dynamic_cast<const PromVector*>(Y);assert(y!=NULL);
    const PromVector *w = dynamic_cast<const PromVector*>(W);assert(w!=NULL);
    return MatMultTransposeAdd( mat_, w->vec_, x->vec_, y->vec_ );
  }
  virtual int Mult( const PromVector_base * const X, 
                    PromVector_base * const Y ) const {
    assert(this!=NULL);
    const PromVector *x = dynamic_cast<const PromVector*>(X);assert(x!=NULL);
    const PromVector *y = dynamic_cast<const PromVector*>(Y);assert(y!=NULL);
    return MatMult( mat_, x->vec_, y->vec_ );
  }
  virtual int KSPMult( const PromVector_base * const x, 
                       PromVector_base * const b ) const {
    return Mult(x,b);
  }
  virtual int MultAdd( const PromVector_base * const X, 
                       const PromVector_base * const Y, 
                       PromVector_base * const W )const{
    assert(this!=NULL);
    const PromVector *x = dynamic_cast<const PromVector*>(X);assert(x!=NULL);
    const PromVector *y = dynamic_cast<const PromVector*>(Y);assert(y!=NULL);
    const PromVector *w = dynamic_cast<const PromVector*>(W);assert(w!=NULL);
    return MatMultAdd( mat_, x->vec_, y->vec_, w->vec_ );
  }
  virtual int GetOwnershipRange( int *first, int *end ) const {
    assert(this!=NULL);
    int ierr = MatGetOwnershipRange( mat_, first, end );
    return ierr;
  }
  int GetRow(int row, int *ncols, int **colpp=NULL,double **arrdat=NULL)const{
    assert(this!=NULL);
    if ( colpp == NULL && arrdat == NULL ) {
      return MatGetRow( mat_, row, ncols, PETSC_NULL, PETSC_NULL );
    }
    else if( arrdat == NULL && colpp != NULL ) {
      return MatGetRow( mat_, row, ncols, colpp, PETSC_NULL );
    }
    else if( colpp == NULL && arrdat != NULL ) {
      return MatGetRow( mat_, row, ncols, PETSC_NULL, arrdat );
    }
    else {
      return MatGetRow( mat_, row, ncols, colpp, arrdat );
    }
  }
  int RestoreRow(int row,int*ncols,int**colpp=NULL,double**arrdat=NULL)const{
    assert(this!=NULL);
    return MatRestoreRow( mat_, row, ncols, colpp, arrdat );
  }
  int SetValues( const int nrow, int *pgeqi, const int ncol, int *pcol,
                 double *vals, const int add_type );
  int SetValuesBlocked( const int nrow, int *pgidi, const int ncol,
                        int *pcol, double *vals, const int add_type);
  int ZeroRows( PromIS is, double diag = 0.0 ){
    assert(this!=NULL);
    dirty_data_ = TRUE;
    return MatZeroRows( mat_, is, diag == 0.0 ? PETSC_NULL : &diag );
  }
  virtual int ZeroEntries(){
    assert(this!=NULL); int ierr;
    if( isSymmetricMatrix() ) {
      ierr = ConvertFromSymm(); CHKERRQ(ierr);
    }
    else {
      ierr = MatZeroEntries( mat_ ); CHKERRQ(ierr);
    }
    dirty_data_ = TRUE; // indefinite_ = FALSE; DRhoA_ = -2.0; -- do for coarse grids w/ helpmholtz!
    return 0;
  }
  int Assembly(){
    assert(this!=NULL);
    int ierr = MatAssemblyBegin( mat_, MAT_FINAL_ASSEMBLY ); CHKERRQ(ierr);
    ierr = MatAssemblyEnd( mat_, MAT_FINAL_ASSEMBLY ); CHKERRQ(ierr);
    return 0;
  } 
  int FlushAssembly(){
    assert(this!=NULL);
    int ierr = MatAssemblyBegin( mat_, MAT_FLUSH_ASSEMBLY ); CHKERRQ(ierr);
    return MatAssemblyEnd( mat_, MAT_FLUSH_ASSEMBLY );
  }
  int Shift( double val ){
    assert(this!=NULL);
    return MatShift( &val, mat_ );
  }
  int Scale(const PromVector_base *const L, 
            const PromVector_base *const R ){
    assert(this!=NULL);
    const PromVector *l = dynamic_cast<const PromVector*>(L);
    const PromVector *r = dynamic_cast<const PromVector*>(R);
    return MatDiagonalScale( mat_,  (l != NULL) ? l->vec_ : NULL,
                             (r != NULL) ? r->vec_ : NULL );
  }
  int Scale( double val ) {
    assert(this!=NULL);
    if( val == 0.0 ) return ZeroEntries();
    else return MatScale( &val, mat_ );
  }
  int NormFrob( double *const norm ) const {
    assert(this!=NULL); PetscScalar nn;
    int ierr = MatNorm( mat_, NORM_FROBENIUS, &nn ); CHKERRQ(ierr);
    *norm = (double)nn;
    return 0;
  }  
  int NormInf( double *const norm ) const {
    assert(this!=NULL); PetscScalar nn;
    if( mat_->stash.size == 1 ) {
      int ierr = MatNorm( mat_, NORM_INFINITY, &nn ); CHKERRQ(ierr);
    }
    else {
      int ierr = MatNorm( mat_, NORM_FROBENIUS, &nn ); CHKERRQ(ierr);
    }
    *norm = (double)nn/sqrt((double)mat_->N);
    return 0;
  }
  int AXPY( double val, const PromMatrix_base * const A ) {
    assert(this!=NULL);
    const PromMatrix *AA = dynamic_cast<const PromMatrix*>(A);assert(AA!=NULL);
    if( AA->dirty_data_ == FALSE ) {
      AA->dirty_data_ = TRUE;
      AA->SetDRho( 0.0 ); // this is collective so set here
    }
    return MatAXPY( &val, AA->mat_, mat_, SAME_NONZERO_PATTERN );
  }
  int SetOption( int opt ) {
    assert(this!=NULL);
    return MatSetOption( mat_, (MatOption)opt );
  }
  virtual int getN() const { return mat_->N; }
  virtual int getn() const { return mat_->n; }
  virtual int getM() const { return mat_->M; }
  virtual int getm() const { return mat_->m; }
  int getMNodes()const{ return map_.nglobal(); }
  int getLocalNNZ() const {
    MatInfo info; 
    MatGetInfo( mat_, MAT_LOCAL, &info );
    return (int)(info.nz_used+.05);
  }
  int getGlobalNNZ() const {
    MatInfo info; 
    MatGetInfo( mat_, MAT_GLOBAL_SUM, &info );
    return (int)(info.nz_used+.05);
  }
  MPI_Comm MPIComm() const { return map_.comm_(); }
  virtual bool isSymmetricValued()const{
    assert(this != NULL && mat_ != NULL);
    if( isSymmetricMatrix() ) return TRUE;
    //else return FALSE; // diabled
    PetscTruth flg; 
    int ierr = MatIsSymmetric(mat_, &flg); // not implemented for BAIJ!!!
    return (flg==PETSC_TRUE);
  }
  virtual bool isSymmetricMatrix()const{
    assert(this != NULL && mat_ != NULL);
    int ierr, size; MatType type;
    MPI_Comm_size( mat_->comm, &size );
    ierr = MatGetType( mat_, &type); CHKERRQ(ierr);
    if(size==1) return (strcmp(type, "seqsbaij")==0);
    else return (strcmp(type, "mpisbaij")==0);
  }
  int GetGArray( int **garr, int *nghosts ) const {
    int ierr, bs, size, *out = NULL;
    MPI_Comm_size( mat_->comm, &size );
    if( size == 1 ) *nghosts = 0;
    else {
      if( isSymmetricMatrix() ) {
        SETERRQ(-1,"GetGArray for symmetric matrix not implemented ");
        //Mat_MPISBAIJ *baij = (Mat_MPISBAIJ*)mat_->data; // upper left matix
        //Mat_SeqSBAIJ *B = (Mat_SeqSBAIJ*)baij->B->data;
        //out = baij->garray;     
        //*nghosts = (out==NULL) ? 0 : B->nbs;
      }
      else { 
        ierr = MatGetBlockSize( mat_, &bs); CHKERRQ(ierr);
        if( bs == 1 ) { 
          Mat_MPIAIJ *aij = (Mat_MPIAIJ*)mat_->data; 
          out = aij->garray;
          *nghosts = (out==NULL) ? 0 : aij->B->n;
        }
        else { // blocked matrix
          Mat_MPIBAIJ *baij = (Mat_MPIBAIJ*)mat_->data; // upper left matix
          Mat_SeqBAIJ *B = (Mat_SeqBAIJ*)baij->B->data;
          out = baij->garray;     
          *nghosts = (out==NULL) ? 0 : B->nbs;
        }
      }
    }
    if(garr) *garr = out;
    return 0;
  }
  int RestoreGArray( int **arr ) const { return 0; }
  bool isOK() const { assert(this!=NULL); return (mat_ != NULL); }
  virtual bool isAssembled()const{
    assert(this!=NULL); return (mat_ && mat_->assembled);
  }
  int Print( char *str=NULL, FILE *ifile = stderr ) const{
    int rank,idi,idj1,kk,first_row,last_row,ncols,*colpp;
    MatGetOwnershipRange(mat_, &first_row, &last_row);
    double *arrdat; 
    FILE *file; 
    if(str == NULL ) file = ifile;
    else {
      char path[128]; MPI_Comm_rank(PETSC_COMM_WORLD, &rank); 
      sprintf(path,"%s_%d", str, rank );
      file = fopen(path,"w"); assert(file!=NULL);
    }
    for( idi = first_row ; idi < last_row ; idi++ ){ 
      MatGetRow( mat_, idi, &ncols, &colpp, &arrdat );
      for ( kk = 0 ; kk < ncols ; kk++ ){
        if( arrdat[kk] == 0.0 ) continue; 
        idj1 = colpp[kk] + 1;
        fprintf(file,"%d %d %.14e\n",idi + 1, idj1, arrdat[kk]);
      }
      MatRestoreRow(mat_,idi,&ncols, &colpp,&arrdat);
    }
    //fprintf( file,"%d %d %.14e\n",last_row, last_row, 0.0 );
    fclose(file);
    return 0;
  }

protected:
  Mat mat_;
  Mat full_mat_;
  mutable int *array_;
public:
  // vbr
  static int getLidWithLBid(const int blid, const int bs, 
                            const int nloc, const int lid_leq[]);
  mutable double *varray_;
  PromTable<int> *ghostGID_BID_;
  PromTable<int> *ghostGID_ndf_; 
  PromTable<int> *ghostBID_GID_;
  mutable int varr_sz_;
  //const bool malled_;
};

#undef __FUNCT__
#define __FUNCT__ "PromPC" 


class PromPC : public PromPC_base
{
  friend class PromPCGaussSeidel;
  friend class PromMG;
  friend class PromPCKKT_ASM;
  friend class PromPCKKT_Seg;
  friend class PromPCKKT_Shell;
public :
  PromPC(const PromOptions &opt, PromPerfMonitor &perf,
         const MPI_Comm comm) : 
    pc_(NULL), rwork_(NULL), work_(NULL),
    PromPC_base(opt,perf) {
    PetscFunctionBegin;
    PCCreate( comm, &pc_ );
    PetscStackPop;
  }
  virtual ~PromPC(){
    if( pc_ != NULL ) PCDestroy(pc_); 
    if( work_ != NULL ){ VecDestroy(work_); VecDestroy(rwork_); }
  }
  virtual int SetType( PromPCType type, const PromVector_base * const vB );
  virtual PromPCType getType() const {
    assert( this != NULL );   PetscTruth flag;
    // need to translate
    PetscTypeCompare((PetscObject)pc_, PCILU, &flag);
    if(flag) return PROMPCILU;
    
    PetscTypeCompare((PetscObject)pc_, PCLU, &flag);
    if(flag) return PROMPCLU;
    
    PetscTypeCompare((PetscObject)pc_, PCJACOBI, &flag);
    if(flag) return PROMPCJACOBI;
    
    PetscTypeCompare((PetscObject)pc_, PCASM, &flag);
    if(flag) return PROMPCASM;

    PetscTypeCompare((PetscObject)pc_, PCNONE, &flag);
    if(flag) return PROMPCNONE;
    
    PetscPrintf(PETSC_COMM_WORLD,"[?]%s line %d ERROR: bad pc type???\n",
                __FILE__,__LINE__);
    exit(1);
    return PROMPCILU;
  }
  virtual int SetAbstractOperator_private(const PromMatrix_base * const Ab){
    const PromMatrix *AA=dynamic_cast<const PromMatrix*>(Ab);assert(AA!=NULL);
    assert(pc_->setupcalled != 2); // should be called before setup!
    pc_->pmat = AA->mat_; // hard set of matrix
    pc_->mat = AA->mat_; assert(AA->mat_ != NULL);
    return 0;
  }
protected:
  int SetUp_private(){
#if defined(PROM_TRACE_PC_SETUP)
    PetscPrintf(pc_->comm,
                "\t\t[0]PromPC::SetUp_private this=%p setupcalled=%d\n",
                this,pc_->setupcalled);
#endif
    assert(this!=NULL);      int ierr; // new data, refactor
    ierr = PCSetUp( pc_ );                      CHKERRQ(ierr);
    ierr = PCSetUpOnBlocks( pc_ );              CHKERRQ(ierr);
    assert(pc_->setupcalled == 2);
    return 0;
  }
  virtual int SetOperator_private( const PromMatrix_base * const A ){
    assert(this!=NULL); 
    if( pc_ == NULL ) SETERRQ(-1,"pc_ == NULL");
    const PromMatrix *AA=dynamic_cast<const PromMatrix*>(A);assert(AA!=NULL);
    // ierr = PCSetOperators( pc_, AA->mat_, AA->mat_, SAME_NONZERO_PATTERN );
    pc_->mat  = AA->mat_; // petsc does a crazy reset in PCSetOperators
    pc_->pmat = AA->mat_;
    pc_->flag = SAME_NONZERO_PATTERN;
    if( pc_->setupcalled == 2 ) pc_->setupcalled = 1; // ** soft reset

    return 0;
  }
  virtual int SetNewOperator_private( const PromMatrix_base * const A ){
    assert(this!=NULL); 
    if( pc_ == NULL ) SETERRQ(-1,"pc_ == NULL");
    const PromMatrix *AA=dynamic_cast<const PromMatrix*>(A);assert(AA!=NULL);
    // ierr = PCSetOperators( pc_, AA->mat_, AA->mat_, SAME_NONZERO_PATTERN );
    pc_->mat  = AA->mat_; // petsc does a crazy reset in PCSetOperators
    pc_->pmat = AA->mat_;
    pc_->flag = DIFFERENT_NONZERO_PATTERN;
    if( pc_->setupcalled == 2 ) pc_->setupcalled = 1; // ** soft reset

    return 0;
  }
public:
  int CreateASM( const int nblocks, PromIS *blockISs );

  virtual const char * getPCString() const {
    assert(this!=NULL); assert(pc_!=NULL);
    PetscFunctionBegin;
    PCType pctype;
    PCGetType( pc_, &pctype );
    if( kkt_ == NULL ){ 
      PetscTruth isasm;   
      PetscTypeCompare((PetscObject)pc_,PCASM,&isasm); 
      if( isasm ){
        int nlocal, first; KSP *psles; PC subpc;
        int ierr = PCASMGetSubKSP( pc_, &nlocal, &first, &psles );
        if(ierr) return  "ERROR-1";
        KSPGetPC( psles[0], &subpc ); PCType subpctype;
        PCGetType( subpc, &subpctype );
        PetscTypeCompare((PetscObject)pc_,PCILU,&isasm); 
        if(isasm) pctype = subpctype;
      }
    }
    PetscFunctionReturn(pctype);
  }
  virtual const char * getSubPCString() const {
    assert(this!=NULL); assert(pc_!=NULL);
    if( kkt_ != NULL ) return kkt_->getPCString(); 
    PetscTruth isasm;   PCType subpctype;
    PetscTypeCompare((PetscObject)pc_,PCASM,&isasm); 
    if( isasm ){
      int nlocal, first; KSP *psles; PC subpc;
      int ierr = PCASMGetSubKSP( pc_, &nlocal, &first, &psles );
      if(ierr) return "ERROR-2";
      KSPGetPC( psles[0], &subpc );
      PCGetType( subpc, &subpctype );
      return subpctype;
    }
    else return PromPCStrings[PROMPCNONE]; 
  }
  virtual int getNumBlocks() const {
    assert(this!=NULL); 
    PetscFunctionBegin;
    if( pc_ == NULL ) SETERRQ(-1,"pc_ == NULL");
    PetscTruth isasm; int ierr;
    ierr = PetscTypeCompare( (PetscObject)pc_, PCASM, &isasm ); CHKERRQ(ierr);
    if( isasm ) {
      int nlocal, first; KSP *psles; 
      ierr = PCASMGetSubKSP( pc_, &nlocal, &first, &psles );  CHKERRQ(ierr);
      PetscFunctionReturn(nlocal);     
    }
    else if( kkt_ != NULL ) PetscFunctionReturn(kkt_->getNumBlocks()); 
    PetscFunctionReturn(0);
  }
  virtual int Apply_private(const PromVector_base * const bbB, 
                            PromVector_base * const xxB,const bool zerox=TRUE);
protected:
  PC  pc_;
  Vec work_;
  Vec rwork_;
}; // PromPC

// want access to PETSc's block data -- hack of all hacks!!!
typedef struct { 
  int        n,n_local,n_local_true;
  int        is_flg;            /* flg set to 1 if the IS created in pcsetup */
  int        overlap;           /* overlap requested by user */
  KSP       *ksp;             /* linear solvers for each block */
  VecScatter *scat;             /* mapping to subregion */
  Vec        *x,*y;
  IS         *is;               /* index set that defines each subdomain */
  Mat        *mat,*pmat;        /* mat is not currently used */
  PCASMType  type;        /* use reduced interpolation, restriction or both */
  int        same_local_solves; /* flag indicating whether all local solvers are same */
  int        inplace;    /* indicates that the sub-matrices are deleted after 
                                PCSetUpOnBlocks() is done. Similar to inplace 
                                1  factor_ization in the case of LU and ILU */
}PETSc_ASM;

class PromGSNode;
int PromDoSchwarzBlocks(Mat AA, const int nb,const int blocks[],PETSc_ASM *osm,
                        const PromVector *const bb, PromVector *const xx,
                        const int ndf, const int my0,
                        const int my0Eq, const bool havgh, const bool pre_s, 
                        const PromTable<PromGSNode*> *gid_ghost, 
                        const double ghostx[], int &numops );

#endif // __PROM_PETSC_H__

---