1 2 #include "src/mat/impls/baij/seq/baij.h" 3 4 EXTERN_C_BEGIN 5 #undef __FUNCT__ 6 #define __FUNCT__ "MatConvert_SeqBAI_SeqAIJ" 7 PetscErrorCode MatConvert_SeqBAIJ_SeqAIJ(Mat A,const MatType newtype,Mat *newmat) 8 { 9 Mat B; 10 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 11 PetscErrorCode ierr; 12 PetscInt bs = A->bs,*ai = a->i,*aj = a->j,n = A->M/bs,i,j,k; 13 PetscInt *rowlengths,*rows,*cols,maxlen = 0,ncols; 14 PetscScalar *aa = a->a; 15 16 PetscFunctionBegin; 17 ierr = PetscMalloc(n*bs*sizeof(int),&rowlengths);CHKERRQ(ierr); 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 ierr = MatCreate(A->comm,A->m,A->n,A->m,A->n,newmat);CHKERRQ(ierr); 25 B = *newmat; 26 ierr = MatSetType(B,newtype);CHKERRQ(ierr); 27 ierr = MatSeqAIJSetPreallocation(B,0,rowlengths);CHKERRQ(ierr); 28 ierr = MatSetOption(B,MAT_COLUMN_ORIENTED);CHKERRQ(ierr); 29 ierr = MatSetOption(B,MAT_ROWS_SORTED);CHKERRQ(ierr); 30 ierr = MatSetOption(B,MAT_COLUMNS_SORTED);CHKERRQ(ierr); 31 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 32 33 ierr = PetscMalloc(bs*sizeof(int),&rows);CHKERRQ(ierr); 34 ierr = PetscMalloc(bs*maxlen*sizeof(int),&cols);CHKERRQ(ierr); 35 for (i=0; i<n; i++) { 36 for (j=0; j<bs; j++) { 37 rows[j] = i*bs+j; 38 } 39 ncols = ai[i+1] - ai[i]; 40 for (k=0; k<ncols; k++) { 41 for (j=0; j<bs; j++) { 42 cols[k*bs+j] = bs*(*aj) + j; 43 } 44 aj++; 45 } 46 ierr = MatSetValues(B,bs,rows,bs*ncols,cols,aa,INSERT_VALUES);CHKERRQ(ierr); 47 aa += ncols*bs*bs; 48 } 49 ierr = PetscFree(cols);CHKERRQ(ierr); 50 ierr = PetscFree(rows);CHKERRQ(ierr); 51 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 52 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 53 54 B->bs = A->bs; 55 56 /* Fake support for "inplace" convert. */ 57 if (B == A) { 58 ierr = MatDestroy(A);CHKERRQ(ierr); 59 } 60 PetscFunctionReturn(0); 61 } 62 EXTERN_C_END 63 64 #include "src/mat/impls/aij/seq/aij.h" 65 66 EXTERN_C_BEGIN 67 #undef __FUNCT__ 68 #define __FUNCT__ "MatConvert_SeqAIJ_SeqBAIJ" 69 PetscErrorCode MatConvert_SeqAIJ_SeqBAIJ(Mat A,const MatType newtype,Mat *newmat) 70 { 71 Mat B; 72 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 73 Mat_SeqBAIJ *b; 74 PetscErrorCode ierr; 75 PetscInt *ai=a->i,m=A->M,n=A->N,i,*rowlengths; 76 77 PetscFunctionBegin; 78 if (n != m) SETERRQ(PETSC_ERR_ARG_WRONG,"Matrix must be square"); 79 80 ierr = PetscMalloc(m*sizeof(int),&rowlengths);CHKERRQ(ierr); 81 for (i=0; i<m; i++) { 82 rowlengths[i] = ai[i+1] - ai[i]; 83 } 84 ierr = MatCreate(A->comm,m,n,m,n,newmat);CHKERRQ(ierr); 85 B = *newmat; 86 ierr = MatSetType(B,newtype);CHKERRQ(ierr); 87 ierr = MatSeqBAIJSetPreallocation(B,1,0,rowlengths);CHKERRQ(ierr); 88 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 89 90 ierr = MatSetOption(B,MAT_ROW_ORIENTED);CHKERRQ(ierr); 91 ierr = MatSetOption(B,MAT_ROWS_SORTED);CHKERRQ(ierr); 92 ierr = MatSetOption(B,MAT_COLUMNS_SORTED);CHKERRQ(ierr); 93 94 b = (Mat_SeqBAIJ*)(B->data); 95 96 ierr = PetscMemcpy(b->i,a->i,(m+1)*sizeof(int));CHKERRQ(ierr); 97 ierr = PetscMemcpy(b->ilen,a->ilen,m*sizeof(int));CHKERRQ(ierr); 98 ierr = PetscMemcpy(b->j,a->j,a->nz*sizeof(int));CHKERRQ(ierr); 99 ierr = PetscMemcpy(b->a,a->a,a->nz*sizeof(MatScalar));CHKERRQ(ierr); 100 101 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 102 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 103 104 /* Fake support for "inplace" convert. */ 105 if (B == A) { 106 ierr = MatDestroy(A);CHKERRQ(ierr); 107 } 108 PetscFunctionReturn(0); 109 } 110 EXTERN_C_END 111