1 2 #include <../src/mat/impls/baij/seq/baij.h> 3 4 #undef __FUNCT__ 5 #define __FUNCT__ "MatConvert_SeqBAIJ_SeqAIJ" 6 PETSC_EXTERN PetscErrorCode MatConvert_SeqBAIJ_SeqAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 7 { 8 Mat B; 9 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 10 PetscErrorCode ierr; 11 PetscInt bs = A->rmap->bs,*ai = a->i,*aj = a->j,n = A->rmap->N/bs,i,j,k; 12 PetscInt *rowlengths,*rows,*cols,maxlen = 0,ncols; 13 MatScalar *aa = a->a; 14 15 PetscFunctionBegin; 16 ierr = PetscMalloc1(n*bs,&rowlengths);CHKERRQ(ierr); 17 for (i=0; i<n; i++) { 18 maxlen = PetscMax(maxlen,(ai[i+1] - ai[i])); 19 for (j=0; j<bs; j++) { 20 rowlengths[i*bs+j] = bs*(ai[i+1] - ai[i]); 21 } 22 } 23 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 24 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 25 ierr = MatSetType(B,MATSEQAIJ);CHKERRQ(ierr); 26 ierr = MatSeqAIJSetPreallocation(B,0,rowlengths);CHKERRQ(ierr); 27 ierr = MatSetOption(B,MAT_ROW_ORIENTED,PETSC_FALSE);CHKERRQ(ierr); 28 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 29 30 ierr = PetscMalloc1(bs,&rows);CHKERRQ(ierr); 31 ierr = PetscMalloc1(bs*maxlen,&cols);CHKERRQ(ierr); 32 for (i=0; i<n; i++) { 33 for (j=0; j<bs; j++) { 34 rows[j] = i*bs+j; 35 } 36 ncols = ai[i+1] - ai[i]; 37 for (k=0; k<ncols; k++) { 38 for (j=0; j<bs; j++) { 39 cols[k*bs+j] = bs*(*aj) + j; 40 } 41 aj++; 42 } 43 ierr = MatSetValues(B,bs,rows,bs*ncols,cols,aa,INSERT_VALUES);CHKERRQ(ierr); 44 aa += ncols*bs*bs; 45 } 46 ierr = PetscFree(cols);CHKERRQ(ierr); 47 ierr = PetscFree(rows);CHKERRQ(ierr); 48 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 49 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 50 51 B->rmap->bs = A->rmap->bs; 52 53 if (reuse == MAT_REUSE_MATRIX) { 54 ierr = MatHeaderReplace(A,B);CHKERRQ(ierr); 55 } else { 56 *newmat = B; 57 } 58 PetscFunctionReturn(0); 59 } 60 61 #include <../src/mat/impls/aij/seq/aij.h> 62 63 #undef __FUNCT__ 64 #define __FUNCT__ "MatConvert_SeqAIJ_SeqBAIJ" 65 PETSC_EXTERN PetscErrorCode MatConvert_SeqAIJ_SeqBAIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 66 { 67 Mat B; 68 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 69 Mat_SeqBAIJ *b; 70 PetscErrorCode ierr; 71 PetscInt *ai=a->i,m=A->rmap->N,n=A->cmap->N,i,*rowlengths; 72 73 PetscFunctionBegin; 74 if (n != m) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix must be square"); 75 if (A->rmap->bs > 1) { 76 ierr = MatConvert_Basic(A,newtype,reuse,newmat);CHKERRQ(ierr); 77 PetscFunctionReturn(0); 78 } 79 80 ierr = PetscMalloc1(m,&rowlengths);CHKERRQ(ierr); 81 for (i=0; i<m; i++) { 82 rowlengths[i] = ai[i+1] - ai[i]; 83 } 84 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 85 ierr = MatSetSizes(B,m,n,m,n);CHKERRQ(ierr); 86 ierr = MatSetType(B,MATSEQBAIJ);CHKERRQ(ierr); 87 ierr = MatSeqBAIJSetPreallocation_SeqBAIJ(B,1,0,rowlengths);CHKERRQ(ierr); 88 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 89 90 ierr = MatSetOption(B,MAT_ROW_ORIENTED,PETSC_TRUE);CHKERRQ(ierr); 91 92 b = (Mat_SeqBAIJ*)(B->data); 93 94 ierr = PetscMemcpy(b->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 95 ierr = PetscMemcpy(b->ilen,a->ilen,m*sizeof(PetscInt));CHKERRQ(ierr); 96 ierr = PetscMemcpy(b->j,a->j,a->nz*sizeof(PetscInt));CHKERRQ(ierr); 97 ierr = PetscMemcpy(b->a,a->a,a->nz*sizeof(MatScalar));CHKERRQ(ierr); 98 99 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 100 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 101 102 if (reuse == MAT_REUSE_MATRIX) { 103 ierr = MatHeaderReplace(A,B);CHKERRQ(ierr); 104 } else { 105 *newmat = B; 106 } 107 PetscFunctionReturn(0); 108 } 109