Commits

Debojyoti Ghosh  committed 1d7ab32

MATTAIJ: Extended the TAIJ type to handle operations of the form [S \otimes I + T \otimes A] where S and T can be rectangular matrices; If NULL passed as S or T to MatCreateTAIJ, then set them to zero.

  • Participants
  • Parent commits 39c5290
  • Branches debo/mat-taij, jed/mat-taij

Comments (0)

Files changed (3)

File include/petscmat.h

 PETSC_EXTERN PetscErrorCode MatMAIJRedimension(Mat,PetscInt,Mat*);
 PETSC_EXTERN PetscErrorCode MatMAIJGetAIJ(Mat,Mat*);
 
-PETSC_EXTERN PetscErrorCode MatCreateTAIJ(Mat,PetscInt,const PetscScalar[],Mat*);
+PETSC_EXTERN PetscErrorCode MatCreateTAIJ(Mat,PetscInt,PetscInt,const PetscScalar[],const PetscScalar[],Mat*);
 PETSC_EXTERN PetscErrorCode MatTAIJGetAIJ(Mat,Mat*);
 PETSC_EXTERN PetscErrorCode MatTAIJGetS(Mat,const PetscScalar**);
+PETSC_EXTERN PetscErrorCode MatTAIJGetT(Mat,const PetscScalar**);
 
 PETSC_EXTERN PetscErrorCode MatComputeExplicitOperator(Mat,Mat*);
 

File src/mat/impls/taij/taij.c

   Defines the basic matrix operations for the TAIJ  matrix storage format.
   This format is used to evaluate matrices of the form:
 
-    [S \otimes I + I \otimes A]
+    [S \otimes I + T \otimes A]
 
   where
-    S is a dense (m \times m) matrix
+    S is a dense (p \times q) matrix
+    T is a dense (p \times q) matrix
     A is an AIJ  (n \times n) matrix
     I is the identity matrix
 
-  The resulting matrix is (mn \times mn)
+  The resulting matrix is (np \times nq)
 
   We provide:
      MatMult()
      MatMultAdd()
      MatInvertBlockDiagonal()
   and
-     MatCreateTAIJ(Mat,dof,Mat*)
+     MatCreateTAIJ(Mat,p,q,,Mat*)
 
   This single directory handles both the sequential and parallel codes
 */
 #undef __FUNCT__
 #define __FUNCT__ "MatTAIJGetS"
 /*@C
-   MatTAIJGetSIJ - Get the S matrix describing the shift action of the TAIJ matrix
+   MatTAIJGetS - Get the S matrix describing the shift action of the TAIJ matrix
 
    Not Collective, the entire S is stored and returned independently on all processes
 
 
    Level: advanced
 
-   Notes: The reference count on the SIJ matrix is not increased so you should not destroy it.
+   Notes: The reference count on the S matrix is not increased so you should not destroy it.
 
 .seealso: MatCreateTAIJ()
 @*/
 }
 
 #undef __FUNCT__
+#define __FUNCT__ "MatTAIJGetT"
+/*@C
+   MatTAIJGetT - Get the T matrix describing the shift action of the TAIJ matrix
+
+   Not Collective, the entire T is stored and returned independently on all processes
+
+   Input Parameter:
+.  A - the TAIJ matrix
+
+   Output Parameter:
+.  T - the T matrix, in form of a scalar array in column-major format
+
+   Level: advanced
+
+   Notes: The reference count on the T matrix is not increased so you should not destroy it.
+
+.seealso: MatCreateTAIJ()
+@*/
+PetscErrorCode  MatTAIJGetT(Mat A,const PetscScalar **T)
+{
+  Mat_SeqTAIJ *b = (Mat_SeqTAIJ*)A->data;
+  PetscFunctionBegin;
+  *T = b->T;
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
 #define __FUNCT__ "MatDestroy_SeqTAIJ"
 PetscErrorCode MatDestroy_SeqTAIJ(Mat A)
 {
 
   PetscFunctionBegin;
   ierr = MatDestroy(&b->AIJ);CHKERRQ(ierr);
-  ierr = PetscFree(b->S);CHKERRQ(ierr);
-  ierr = PetscFree(b->ibdiag);CHKERRQ(ierr);
+  ierr = PetscFree3(b->S,b->T,b->ibdiag);CHKERRQ(ierr);
   ierr = PetscFree(A->data);CHKERRQ(ierr);
 
   PetscFunctionReturn(0);
   ierr = MatDestroy(&b->A);CHKERRQ(ierr);
   ierr = VecScatterDestroy(&b->ctx);CHKERRQ(ierr);
   ierr = VecDestroy(&b->w);CHKERRQ(ierr);
+  ierr = PetscFree3(b->S,b->T,b->ibdiag);CHKERRQ(ierr);
   ierr = PetscFree(A->data);CHKERRQ(ierr);
 
   PetscFunctionReturn(0);
 /*MC
   MATTAIJ - MATTAIJ = "taij" - A matrix type to be used to evaluate matrices of the following form:
 
-    [S \otimes I + I \otimes A]
+    [S \otimes I + T \otimes A]
 
   where
-    S is a dense (m \times m) matrix
+    S is a dense (p \times q) matrix
+    T is a dense (p \times q) matrix
     A is an AIJ  (n \times n) matrix
     I is the identity matrix
-  The resulting matrix is (mn \times mn)
+  The resulting matrix is (np \times nq)
   
   The matrix type is based on MATSEQAIJ for a sequential matrix A, and MATMPIAIJ for a distributed matrix A. 
-  S is always stored independently on all processes as a PetscScalar array in column-major format.
+  S and T are always stored independently on all processes as a PetscScalar array in column-major format.
 
   Operations provided:
 . MatMult
 
   Level: advanced
 
-.seealso: MatTAIJGetAIJ(), MatCreateTAIJ()
+.seealso: MatTAIJGetAIJ(), MatTAIJGetS(), MatTAIJGetT(), MatCreateTAIJ()
 M*/
 
 #undef __FUNCT__
 {
   Mat_SeqTAIJ       *b = (Mat_SeqTAIJ*)A->data;
   Mat_SeqAIJ        *a = (Mat_SeqAIJ*)b->AIJ->data;
-  const PetscScalar *s = b->S;
+  const PetscScalar *s = b->S, *t = b->T;
   const PetscScalar *x,*v,*bx;
   PetscScalar       *y,*sums;
   PetscErrorCode    ierr;
   const PetscInt    m = b->AIJ->rmap->n,*idx,*ii;
-  PetscInt          n,i,jrow,j,dof = b->dof,k;
+  PetscInt          n,i,jrow,j,l,p=b->p,q=b->q,k;
 
   PetscFunctionBegin;
 
   for (i=0; i<m; i++) {
     jrow = ii[i];
     n    = ii[i+1] - jrow;
-    sums = y + dof*i;
-    bx   = x + dof*i;
-    for (j=0; j<dof; j++) {
-      for (k=0; k<dof; k++) {
-        sums[k] += s[k+j*dof]*bx[j];
+    sums = y + p*i;
+    bx   = x + q*i;
+    for (j=0; j<q; j++) {
+      for (k=0; k<p; k++) {
+        sums[k] += s[k+j*p]*bx[j];
       }
     }
     for (j=0; j<n; j++) {
-      for (k=0; k<dof; k++) {
-        sums[k] += v[jrow+j]*x[dof*idx[jrow+j]+k];
+      for (k=0; k<p; k++) {
+        for (l=0; l<q; l++) {
+          sums[k] += v[jrow+j]*t[k+l*p]*x[q*idx[jrow+j]+l];
+        }
       }
     }
   }
 
-  ierr = PetscLogFlops((2.0*dof*dof-dof)*m+2*dof*a->nz);CHKERRQ(ierr);
+  ierr = PetscLogFlops((2.0*p*q-p)*m+2*p*a->nz);CHKERRQ(ierr);
   ierr = VecRestoreArrayRead(*xx,&x);CHKERRQ(ierr);
   ierr = VecRestoreArray(*zz,&y);CHKERRQ(ierr);
   PetscFunctionReturn(0);
   Mat_SeqAIJ        *a  = (Mat_SeqAIJ*)b->AIJ->data;
   const PetscScalar *s  = b->S;
   const PetscScalar *v  = a->a;
-  const PetscInt    dof = b->dof,dof2=dof*dof,m = b->AIJ->rmap->n,*idx = a->j,*ii = a->i;
+  const PetscInt     p  = b->p, q = b->q, m = b->AIJ->rmap->n, *idx = a->j, *ii = a->i;
   PetscErrorCode    ierr;
-  PetscInt          i,j,*v_pivots;
+  PetscInt          i,j,*v_pivots,dof,dof2;
   PetscScalar       *diag,aval,*v_work;
 
   PetscFunctionBegin;
+  if (p != q) {
+    SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Block size must be square to calculate inverse.");
+  } else {
+    dof  = p;
+    dof2 = dof*dof;
+  }
   if (b->ibdiagvalid) {
     if (values) *values = b->ibdiag;
     PetscFunctionReturn(0);
   Mat            a    = b->AIJ,B;
   Mat_SeqAIJ     *aij = (Mat_SeqAIJ*)a->data;
   PetscErrorCode ierr;
-  PetscInt       m,n,i,ncols,*ilen,nmax = 0,*icols,j,k,ii,dof = b->dof;
-  PetscInt       *cols,*row;
-  PetscScalar    *vals,*s;
+  PetscInt       m,n,i,ncols,*ilen,nmax = 0,*icols,j,k,l,p=b->p,q=b->q;
+  PetscInt       *cols,*srow,*scol;
+  PetscScalar    *vals,*s,*t;
 
   PetscFunctionBegin;
   ierr = MatGetSize(a,&m,&n);CHKERRQ(ierr);
-  ierr = PetscMalloc(dof*m*sizeof(PetscInt),&ilen);CHKERRQ(ierr);
+  ierr = PetscMalloc(p*m*sizeof(PetscInt),&ilen);CHKERRQ(ierr);
   for (i=0; i<m; i++) {
     nmax = PetscMax(nmax,aij->ilen[i]);
-    for (j=0; j<dof; j++) ilen[dof*i+j] = aij->ilen[i] + dof;
+    for (j=0; j<p; j++) ilen[p*i+j] = aij->ilen[i]*q;
   }
-  ierr = MatCreateSeqAIJ(PETSC_COMM_SELF,dof*m,dof*n,0,ilen,&B);CHKERRQ(ierr);
+  ierr = MatCreateSeqAIJ(PETSC_COMM_SELF,m*p,n*q,0,ilen,&B);CHKERRQ(ierr);
   ierr = PetscFree(ilen);CHKERRQ(ierr);
-  ierr = PetscMalloc3(nmax,PetscInt,&icols,dof*dof,PetscScalar,&s,
-                      dof,PetscInt,&row);CHKERRQ(ierr);
-  ii   = 0;
+  ierr = PetscMalloc5(q,PetscInt,&icols,p*q,PetscScalar,&s,
+                      p,PetscInt,&srow,q,PetscInt,&scol,
+                      p*q,PetscScalar,&t);CHKERRQ(ierr);
+
   for (i=0; i<m; i++) {
-    for (j=0; j<dof; j++) {
-      row[j] = i*dof+j;
-      for (k=0; k<dof; k++) {
-        s[j*dof+k] = b->S[j+k*dof];
+    for (j=0; j<p; j++) srow[j] = i*p+j;
+    for (k=0; k<q; k++) scol[k] = i*q+k;
+    for (j=0; j<p; j++) {
+      for (k=0; k<q; k++) {
+        s[j*q+k] = b->S[j+k*p];
       }
     }
-    ierr = MatSetValues_SeqAIJ(B,dof,row,dof,row,s,ADD_VALUES);CHKERRQ(ierr);
+    ierr = MatSetValues_SeqAIJ(B,p,srow,q,scol,s,ADD_VALUES);CHKERRQ(ierr);
     ierr = MatGetRow_SeqAIJ(a,i,&ncols,&cols,&vals);CHKERRQ(ierr);
-    for (j=0; j<dof; j++) {
-      for (k=0; k<ncols; k++) icols[k] = dof*cols[k]+j;
-      ierr = MatSetValues_SeqAIJ(B,1,&ii,ncols,icols,vals,ADD_VALUES);CHKERRQ(ierr);
-      ii++;
+    for (k=0; k<ncols; k++) {
+      for (j=0; j<q; j++) icols[j] = q*cols[k]+j;
+      for (j=0; j<p; j++) {
+        for (l=0; l<q; l++) {
+          t[j*q+l] = b->T[j+l*p]*vals[k];
+        }
+      }
+      ierr = MatSetValues_SeqAIJ(B,p,srow,q,icols,t,ADD_VALUES);CHKERRQ(ierr);
     }
     ierr = MatRestoreRow_SeqAIJ(a,i,&ncols,&cols,&vals);CHKERRQ(ierr);
   }
-  ierr = PetscFree3(icols,s,row);CHKERRQ(ierr);
   ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
   ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+  ierr = PetscFree5(icols,s,srow,scol,t);CHKERRQ(ierr);
 
   if (reuse == MAT_REUSE_MATRIX) {
     ierr = MatHeaderReplace(A,B);CHKERRQ(ierr);
 /*@C
   MatCreateTAIJ - Creates a matrix type to be used for matrices of the following form:
 
-    [S \otimes I + I \otimes A]
+    [S \otimes I + T \otimes A]
 
   where
-    S is a dense (m \times m) matrix
+    S is a dense (p \times q) matrix
+    T is a dense (p \times q) matrix
     A is an AIJ  (n \times n) matrix
     I is the identity matrix
-  The resulting matrix is (mn \times mn)
+  The resulting matrix is (np \times nq)
   
   The matrix type is based on MATSEQAIJ for a sequential matrix A, and MATMPIAIJ for a distributed matrix A. 
   S is always stored independently on all processes as a PetscScalar array in column-major format.
   Input Parameters:
 + A - the AIJ matrix
 + S - the S matrix, stored as a PetscScalar array (column-major)
-- dof - the block size (number of components per node) (dof = m)
++ T - the T matrix, stored as a PetscScalar array (column-major)
+- p - number of rows in S and T
+- q - number of columns in S and T
 
   Output Parameter:
 . taij - the new TAIJ matrix
 
 .seealso: MatTAIJGetAIJ(), MATTAIJ
 @*/
-PetscErrorCode  MatCreateTAIJ(Mat A,PetscInt dof,const PetscScalar S[],Mat *taij)
+PetscErrorCode  MatCreateTAIJ(Mat A,PetscInt p,PetscInt q,const PetscScalar S[],const PetscScalar T[],Mat *taij)
 {
   PetscErrorCode ierr;
   PetscMPIInt    size;
   PetscFunctionBegin;
   ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr);
 
-  if (dof == 1) *taij = A;
-  else {
-    ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr);
-    ierr = MatSetSizes(B,dof*A->rmap->n,dof*A->cmap->n,dof*A->rmap->N,dof*A->cmap->N);CHKERRQ(ierr);
-    ierr = PetscLayoutSetBlockSize(B->rmap,dof);CHKERRQ(ierr);
-    ierr = PetscLayoutSetBlockSize(B->cmap,dof);CHKERRQ(ierr);
-    ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr);
-    ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr);
-
-    B->assembled = PETSC_TRUE;
-
-    ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);CHKERRQ(ierr);
-    if (size == 1) {
-      Mat_SeqTAIJ *b;
-
-      ierr = MatSetType(B,MATSEQTAIJ);CHKERRQ(ierr);
-
-      B->ops->setup   = NULL;
-      B->ops->destroy = MatDestroy_SeqTAIJ;
-      B->ops->view    = MatView_SeqTAIJ;
-      b               = (Mat_SeqTAIJ*)B->data;
-      b->dof          = dof;
-      b->AIJ          = A;
-      PetscMalloc(dof*dof*sizeof(PetscScalar),&b->S);CHKERRQ(ierr);
-      PetscMemcpy(b->S,S,dof*dof*sizeof(PetscScalar));
+  ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr);
+  ierr = MatSetSizes(B,p*A->rmap->n,q*A->cmap->n,p*A->rmap->N,q*A->cmap->N);CHKERRQ(ierr);
+  ierr = PetscLayoutSetBlockSize(B->rmap,p);CHKERRQ(ierr);
+  ierr = PetscLayoutSetBlockSize(B->cmap,q);CHKERRQ(ierr);
+  ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr);
+  ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr);
 
-      B->ops->mult                = MatMult_SeqTAIJ_N;
-      B->ops->multadd             = MatMultAdd_SeqTAIJ_N;
-      B->ops->invertblockdiagonal = MatInvertBlockDiagonal_SeqTAIJ_N;
+  B->assembled = PETSC_TRUE;
+  ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);CHKERRQ(ierr);
 
-      ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqtaij_seqaij_C",MatConvert_SeqTAIJ_SeqAIJ);CHKERRQ(ierr);
-
-    } else {
-      Mat_MPIAIJ  *mpiaij = (Mat_MPIAIJ*)A->data;
-      Mat_MPITAIJ *b;
-      IS          from,to;
-      Vec         gvec;
-
-      ierr = MatSetType(B,MATMPITAIJ);CHKERRQ(ierr);
-
-      B->ops->setup   = NULL;
-      B->ops->destroy = MatDestroy_MPITAIJ;
-      B->ops->view    = MatView_MPITAIJ;
-
-      b      = (Mat_MPITAIJ*)B->data;
-      b->dof = dof;
-      b->A   = A;
-      PetscMalloc(dof*dof*sizeof(PetscScalar),&b->S);CHKERRQ(ierr);
-      PetscMemcpy(b->S,S,dof*dof*sizeof(PetscScalar));
-
-      ierr = MatCreateTAIJ(mpiaij->A,dof,b->S,&b->AIJ);CHKERRQ(ierr); 
-      ierr = MatCreateTAIJ(mpiaij->B,dof,b->S,&b->OAIJ);CHKERRQ(ierr);
-
-      ierr = VecGetSize(mpiaij->lvec,&n);CHKERRQ(ierr);
-      ierr = VecCreate(PETSC_COMM_SELF,&b->w);CHKERRQ(ierr);
-      ierr = VecSetSizes(b->w,n*dof,n*dof);CHKERRQ(ierr);
-      ierr = VecSetBlockSize(b->w,dof);CHKERRQ(ierr);
-      ierr = VecSetType(b->w,VECSEQ);CHKERRQ(ierr);
-
-      /* create two temporary Index sets for build scatter gather */
-      ierr = ISCreateBlock(PetscObjectComm((PetscObject)A),dof,n,mpiaij->garray,PETSC_COPY_VALUES,&from);CHKERRQ(ierr);
-      ierr = ISCreateStride(PETSC_COMM_SELF,n*dof,0,1,&to);CHKERRQ(ierr);
+  if (size == 1) {
+    Mat_SeqTAIJ *b;
 
-      /* create temporary global vector to generate scatter context */
-      ierr = VecCreateMPIWithArray(PetscObjectComm((PetscObject)A),dof,dof*A->cmap->n,dof*A->cmap->N,NULL,&gvec);CHKERRQ(ierr);
+    ierr = MatSetType(B,MATSEQTAIJ);CHKERRQ(ierr);
 
-      /* generate the scatter context */
-      ierr = VecScatterCreate(gvec,from,b->w,to,&b->ctx);CHKERRQ(ierr);
+    B->ops->setup   = NULL;
+    B->ops->destroy = MatDestroy_SeqTAIJ;
+    B->ops->view    = MatView_SeqTAIJ;
+    b               = (Mat_SeqTAIJ*)B->data;
+    b->p            = p;
+    b->q            = q;
+    b->AIJ          = A;
+    PetscMalloc2(p*q,PetscScalar,&b->S,p*q,PetscScalar,&b->T);CHKERRQ(ierr);
+    if (S)  PetscMemcpy (b->S,S,p*q*sizeof(PetscScalar));
+    else    PetscMemzero(b->S,p*q*sizeof(PetscScalar));
+    if (T)  PetscMemcpy(b->T,T,p*q*sizeof(PetscScalar));
+    else    PetscMemzero(b->T,p*q*sizeof(PetscScalar));
 
-      ierr = ISDestroy(&from);CHKERRQ(ierr);
-      ierr = ISDestroy(&to);CHKERRQ(ierr);
-      ierr = VecDestroy(&gvec);CHKERRQ(ierr);
+    B->ops->mult                = MatMult_SeqTAIJ_N;
+    B->ops->multadd             = MatMultAdd_SeqTAIJ_N;
+    B->ops->invertblockdiagonal = MatInvertBlockDiagonal_SeqTAIJ_N;
 
-      B->ops->mult                = MatMult_MPITAIJ_dof;
-      B->ops->multadd             = MatMultAdd_MPITAIJ_dof;
-      B->ops->invertblockdiagonal = MatInvertBlockDiagonal_MPITAIJ_dof;
+    ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqtaij_seqaij_C",MatConvert_SeqTAIJ_SeqAIJ);CHKERRQ(ierr);
 
-      //ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_mpitaij_mpiaij_C",MatConvert_MPITAIJ_MPIAIJ);CHKERRQ(ierr);
-    }
-    B->ops->getsubmatrix = MatGetSubMatrix_TAIJ;
-    ierr  = MatSetUp(B);CHKERRQ(ierr);
-    *taij = B;
-    ierr  = MatViewFromOptions(B,NULL,"-mat_view");CHKERRQ(ierr);
+  } else {
+    Mat_MPIAIJ  *mpiaij = (Mat_MPIAIJ*)A->data;
+    Mat_MPITAIJ *b;
+    IS          from,to;
+    Vec         gvec;
+
+    ierr = MatSetType(B,MATMPITAIJ);CHKERRQ(ierr);
+
+    B->ops->setup   = NULL;
+    B->ops->destroy = MatDestroy_MPITAIJ;
+    B->ops->view    = MatView_MPITAIJ;
+
+    b      = (Mat_MPITAIJ*)B->data;
+    b->p   = p;
+    b->q   = q;
+    b->A   = A;
+    PetscMalloc2(p*q,PetscScalar,&b->S,p*q,PetscScalar,&b->T);CHKERRQ(ierr);
+    if (S)  PetscMemcpy (b->S,S,p*q*sizeof(PetscScalar));
+    else    PetscMemzero(b->S,p*q*sizeof(PetscScalar));
+    if (T)  PetscMemcpy(b->T,T,p*q*sizeof(PetscScalar));
+    else    PetscMemzero(b->T,p*q*sizeof(PetscScalar));
+
+    ierr = MatCreateTAIJ(mpiaij->A,p,q,b->S,b->T,&b->AIJ);CHKERRQ(ierr); 
+    ierr = MatCreateTAIJ(mpiaij->B,p,q,b->S,b->T,&b->OAIJ);CHKERRQ(ierr);
+
+    ierr = VecGetSize(mpiaij->lvec,&n);CHKERRQ(ierr);
+    ierr = VecCreate(PETSC_COMM_SELF,&b->w);CHKERRQ(ierr);
+    ierr = VecSetSizes(b->w,n*q,n*q);CHKERRQ(ierr);
+    ierr = VecSetBlockSize(b->w,q);CHKERRQ(ierr);
+    ierr = VecSetType(b->w,VECSEQ);CHKERRQ(ierr);
+
+    /* create two temporary Index sets for build scatter gather */
+    ierr = ISCreateBlock(PetscObjectComm((PetscObject)A),q,n,mpiaij->garray,PETSC_COPY_VALUES,&from);CHKERRQ(ierr);
+    ierr = ISCreateStride(PETSC_COMM_SELF,n*q,0,1,&to);CHKERRQ(ierr);
+
+    /* create temporary global vector to generate scatter context */
+    ierr = VecCreateMPIWithArray(PetscObjectComm((PetscObject)A),q,q*A->cmap->n,q*A->cmap->N,NULL,&gvec);CHKERRQ(ierr);
+
+    /* generate the scatter context */
+    ierr = VecScatterCreate(gvec,from,b->w,to,&b->ctx);CHKERRQ(ierr);
+
+    ierr = ISDestroy(&from);CHKERRQ(ierr);
+    ierr = ISDestroy(&to);CHKERRQ(ierr);
+    ierr = VecDestroy(&gvec);CHKERRQ(ierr);
+
+    B->ops->mult                = MatMult_MPITAIJ_dof;
+    B->ops->multadd             = MatMultAdd_MPITAIJ_dof;
+    B->ops->invertblockdiagonal = MatInvertBlockDiagonal_MPITAIJ_dof;
+
+    //ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_mpitaij_mpiaij_C",MatConvert_MPITAIJ_MPIAIJ);CHKERRQ(ierr);
   }
+  B->ops->getsubmatrix = MatGetSubMatrix_TAIJ;
+  ierr  = MatSetUp(B);CHKERRQ(ierr);
+  *taij = B;
+  ierr  = MatViewFromOptions(B,NULL,"-mat_view");CHKERRQ(ierr);
+
   PetscFunctionReturn(0);
 }

File src/mat/impls/taij/taij.h

 #include <../src/mat/impls/aij/mpi/mpiaij.h>
 
 #define TAIJHEADER          \
-  PetscInt    dof;          \
+  PetscInt    p,q;          \
   Mat         AIJ;          \
   PetscScalar *S;           \
+  PetscScalar *T;           \
   PetscScalar *ibdiag;      \
   PetscBool   ibdiagvalid;  \
 
 
 typedef struct {
   TAIJHEADER;
-  Mat        OAIJ;    /* representation of interpolation for one component */
+  Mat        OAIJ;    
   Mat        A;
   VecScatter ctx;     /* update ghost points for parallel case */
   Vec        w;       /* work space for ghost values for parallel case */