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_SeqAIJ *b; 8 PetscBool roworiented; 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 if (reuse == MAT_REUSE_MATRIX) { 17 B = *newmat; 18 for (i=0; i<n; i++) { 19 maxlen = PetscMax(maxlen,(ai[i+1] - ai[i])); 20 } 21 } else { 22 ierr = PetscMalloc1(n*bs,&rowlengths);CHKERRQ(ierr); 23 for (i=0; i<n; i++) { 24 maxlen = PetscMax(maxlen,(ai[i+1] - ai[i])); 25 for (j=0; j<bs; j++) { 26 rowlengths[i*bs+j] = bs*(ai[i+1] - ai[i]); 27 } 28 } 29 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 30 ierr = MatSetType(B,MATSEQAIJ);CHKERRQ(ierr); 31 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 32 ierr = MatSetBlockSizes(B,A->rmap->bs,A->cmap->bs);CHKERRQ(ierr); 33 ierr = MatSeqAIJSetPreallocation(B,0,rowlengths);CHKERRQ(ierr); 34 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 35 } 36 b = (Mat_SeqAIJ*)B->data; 37 roworiented = b->roworiented; 38 39 ierr = MatSetOption(B,MAT_ROW_ORIENTED,PETSC_FALSE);CHKERRQ(ierr); 40 ierr = PetscMalloc1(bs,&rows);CHKERRQ(ierr); 41 ierr = PetscMalloc1(bs*maxlen,&cols);CHKERRQ(ierr); 42 for (i=0; i<n; i++) { 43 for (j=0; j<bs; j++) { 44 rows[j] = i*bs+j; 45 } 46 ncols = ai[i+1] - ai[i]; 47 for (k=0; k<ncols; k++) { 48 for (j=0; j<bs; j++) { 49 cols[k*bs+j] = bs*(*aj) + j; 50 } 51 aj++; 52 } 53 ierr = MatSetValues(B,bs,rows,bs*ncols,cols,aa,INSERT_VALUES);CHKERRQ(ierr); 54 aa += ncols*bs*bs; 55 } 56 ierr = PetscFree(cols);CHKERRQ(ierr); 57 ierr = PetscFree(rows);CHKERRQ(ierr); 58 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 59 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 60 ierr = MatSetOption(B,MAT_ROW_ORIENTED,roworiented);CHKERRQ(ierr); 61 62 if (reuse == MAT_INPLACE_MATRIX) { 63 ierr = MatHeaderReplace(A,&B);CHKERRQ(ierr); 64 } else *newmat = B; 65 PetscFunctionReturn(0); 66 } 67 68 #include <../src/mat/impls/aij/seq/aij.h> 69 70 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqBAIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 71 { 72 Mat B; 73 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 74 Mat_SeqBAIJ *b; 75 PetscErrorCode ierr; 76 PetscInt *ai=a->i,m=A->rmap->N,n=A->cmap->N,i,*rowlengths,bs=PetscAbs(A->rmap->bs); 77 78 PetscFunctionBegin; 79 if (reuse != MAT_REUSE_MATRIX) { 80 ierr = PetscMalloc1(m/bs,&rowlengths);CHKERRQ(ierr); 81 for (i=0; i<m/bs; i++) { 82 rowlengths[i] = (ai[i*bs+1] - ai[i*bs])/bs; 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(B,bs,0,rowlengths);CHKERRQ(ierr); 88 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 89 } else B = *newmat; 90 91 if (bs == 1) { 92 b = (Mat_SeqBAIJ*)(B->data); 93 94 ierr = PetscArraycpy(b->i,a->i,m+1);CHKERRQ(ierr); 95 ierr = PetscArraycpy(b->ilen,a->ilen,m);CHKERRQ(ierr); 96 ierr = PetscArraycpy(b->j,a->j,a->nz);CHKERRQ(ierr); 97 ierr = PetscArraycpy(b->a,a->a,a->nz);CHKERRQ(ierr); 98 99 ierr = MatSetOption(B,MAT_ROW_ORIENTED,PETSC_TRUE);CHKERRQ(ierr); 100 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 101 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 102 } else { 103 /* reuse may not be equal to MAT_REUSE_MATRIX, but the basic converter will reallocate or replace newmat if this value is not used */ 104 /* if reuse is equal to MAT_INITIAL_MATRIX, it has been appropriately preallocated before */ 105 /* MAT_INPLACE_MATRIX, it will be replaced with MatHeaderReplace below */ 106 ierr = MatConvert_Basic(A,newtype,MAT_REUSE_MATRIX,&B);CHKERRQ(ierr); 107 } 108 109 if (reuse == MAT_INPLACE_MATRIX) { 110 ierr = MatHeaderReplace(A,&B);CHKERRQ(ierr); 111 } else *newmat = B; 112 PetscFunctionReturn(0); 113 } 114