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