2: #include <../src/mat/impls/baij/seq/baij.h>

  4: EXTERN_C_BEGIN
  7: PetscErrorCode  MatConvert_SeqBAIJ_SeqAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
  8: {
  9:   Mat            B;
 10:   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ*)A->data;
 12:   PetscInt       bs = A->rmap->bs,*ai = a->i,*aj = a->j,n = A->rmap->N/bs,i,j,k;
 13:   PetscInt       *rowlengths,*rows,*cols,maxlen = 0,ncols;
 14:   MatScalar      *aa = a->a;

 17:   PetscMalloc(n*bs*sizeof(PetscInt),&rowlengths);
 18:   for (i=0; i<n; i++) {
 19:     maxlen = PetscMax(maxlen,(ai[i+1] - ai[i]));
 20:     for (j=0; j<bs; j++) {
 21:       rowlengths[i*bs+j] = bs*(ai[i+1] - ai[i]);
 22:     }
 23:   }
 24:   MatCreate(((PetscObject)A)->comm,&B);
 25:   MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
 26:   MatSetType(B,MATSEQAIJ);
 27:   MatSeqAIJSetPreallocation(B,0,rowlengths);
 28:   MatSetOption(B,MAT_ROW_ORIENTED,PETSC_FALSE);
 29:   PetscFree(rowlengths);

 31:   PetscMalloc(bs*sizeof(PetscInt),&rows);
 32:   PetscMalloc(bs*maxlen*sizeof(PetscInt),&cols);
 33:   for (i=0; i<n; i++) {
 34:     for (j=0; j<bs; j++) {
 35:       rows[j] = i*bs+j;
 36:     }
 37:     ncols = ai[i+1] - ai[i];
 38:     for (k=0; k<ncols; k++) {
 39:       for (j=0; j<bs; j++) {
 40:         cols[k*bs+j] = bs*(*aj) + j;
 41:       }
 42:       aj++;
 43:     }
 44:     MatSetValues(B,bs,rows,bs*ncols,cols,aa,INSERT_VALUES);
 45:     aa   += ncols*bs*bs;
 46:   }
 47:   PetscFree(cols);
 48:   PetscFree(rows);
 49:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 50:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
 51: 
 52:   B->rmap->bs = A->rmap->bs;

 54:   if (reuse == MAT_REUSE_MATRIX) {
 55:     MatHeaderReplace(A,B);
 56:   } else {
 57:     *newmat = B;
 58:   }
 59:   return(0);
 60: }
 61: EXTERN_C_END

 63: #include <../src/mat/impls/aij/seq/aij.h>

 65: EXTERN_C_BEGIN
 68: PetscErrorCode  MatConvert_SeqAIJ_SeqBAIJ(Mat A,const MatType newtype,MatReuse reuse,Mat *newmat)
 69: {
 70:   Mat            B;
 71:   Mat_SeqAIJ     *a = (Mat_SeqAIJ*)A->data;
 72:   Mat_SeqBAIJ    *b;
 74:   PetscInt       *ai=a->i,m=A->rmap->N,n=A->cmap->N,i,*rowlengths;

 77:   if (n != m) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix must be square");

 79:   PetscMalloc(m*sizeof(PetscInt),&rowlengths);
 80:   for (i=0; i<m; i++) {
 81:     rowlengths[i] = ai[i+1] - ai[i];
 82:   }
 83:   MatCreate(((PetscObject)A)->comm,&B);
 84:   MatSetSizes(B,m,n,m,n);
 85:   MatSetType(B,MATSEQBAIJ);
 86:   MatSeqBAIJSetPreallocation_SeqBAIJ(B,1,0,rowlengths);
 87:   PetscFree(rowlengths);

 89:   MatSetOption(B,MAT_ROW_ORIENTED,PETSC_TRUE);
 90: 
 91:   b  = (Mat_SeqBAIJ*)(B->data);

 93:   PetscMemcpy(b->i,a->i,(m+1)*sizeof(PetscInt));
 94:   PetscMemcpy(b->ilen,a->ilen,m*sizeof(PetscInt));
 95:   PetscMemcpy(b->j,a->j,a->nz*sizeof(PetscInt));
 96:   PetscMemcpy(b->a,a->a,a->nz*sizeof(MatScalar));
 97: 
 98:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 99:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

101:   if (reuse == MAT_REUSE_MATRIX) {
102:     MatHeaderReplace(A,B);
103:   } else {
104:     *newmat = B;
105:   }
106:   return(0);
107: }
108: EXTERN_C_END