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