1 /*$Id: aijsbaij.c,v 1.9 2001/08/07 03:02:55 balay Exp $*/ 2 3 #include "src/mat/impls/aij/seq/aij.h" 4 #include "src/mat/impls/sbaij/seq/sbaij.h" 5 6 EXTERN_C_BEGIN 7 #undef __FUNCT__ 8 #define __FUNCT__ "MatConvert_SeqSBAI_SeqAIJ" 9 int MatConvert_SeqSBAIJ_SeqAIJ(Mat A,const MatType newtype,Mat *newmat) 10 { 11 Mat B; 12 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 13 Mat_SeqAIJ *b; 14 int ierr,bs = a->bs,*ai=a->i,*aj=a->j,m=A->M/bs,*bi,*bj, 15 i,*rowlengths,nz,*rowstart; 16 PetscScalar *av,*bv; 17 18 PetscFunctionBegin; 19 /* compute rowlengths of newmat */ 20 ierr = PetscMalloc(m*bs*sizeof(int),&rowlengths);CHKERRQ(ierr); 21 for (i=0; i<m; i++) rowlengths[i] = 0; 22 aj = a->j; 23 for (i=0; i<m; i++) { 24 nz = ai[i+1] - ai[i]; 25 rowlengths[i] += nz; /* upper triangular part */ 26 aj++; nz--; /* skip diagonal */ 27 while (nz--) { rowlengths[*aj++]++; } /* lower triangular part */ 28 } 29 30 ierr = MatCreateSeqAIJ(PETSC_COMM_SELF,m,m,0,rowlengths,&B);CHKERRQ(ierr); 31 ierr = MatSetOption(B,MAT_ROW_ORIENTED);CHKERRQ(ierr); 32 ierr = MatSetOption(B,MAT_ROWS_SORTED);CHKERRQ(ierr); 33 ierr = MatSetOption(B,MAT_COLUMNS_SORTED);CHKERRQ(ierr); 34 35 b = (Mat_SeqAIJ*)(B->data); 36 bi = b->i; 37 bj = b->j; 38 bv = b->a; 39 40 /* set b->i */ 41 rowstart = rowlengths; /* rowstart renames rowlengths for code understanding */ 42 bi[0] = 0; 43 for (i=0; i<m; i++){ 44 b->ilen[i] = rowlengths[i]; 45 bi[i+1] = bi[i] + rowlengths[i]; 46 rowstart[i] = bi[i]; 47 } 48 if (bi[m] != 2*a->nz - m) SETERRQ2(1,"bi[m]: %d != 2*a->nz-m: %d\n",bi[m],2*a->nz - m); 49 50 /* set b->j and b->a */ 51 aj = a->j; av = a->a; 52 for (i=0; i<m; i++) { 53 /* diagonal */ 54 *(bj + rowstart[i]) = *aj; aj++; 55 *(bv + rowstart[i]) = *av; av++; 56 rowstart[i]++; 57 nz = ai[i+1] - ai[i] - 1; 58 while (nz--){ 59 /* lower triangular part */ 60 *(bj + rowstart[*aj]) = i; 61 *(bv + rowstart[*aj]) = *av; 62 rowstart[*aj]++; 63 /* upper triangular part */ 64 *(bj + rowstart[i]) = *aj; aj++; 65 *(bv + rowstart[i]) = *av; av++; 66 rowstart[i]++; 67 } 68 } 69 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 70 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 71 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 72 73 /* Fake support for "inplace" convert. */ 74 if (*newmat == A) { 75 ierr = MatDestroy(A);CHKERRQ(ierr); 76 } 77 *newmat = B; 78 79 PetscFunctionReturn(0); 80 } 81 #undef __FUNCT__ 82 #define __FUNCT__ "MatConvert_SeqAIJ_SeqSBAIJ" 83 int MatConvert_SeqAIJ_SeqSBAIJ(Mat A,const MatType newtype,Mat *newmat) { 84 Mat B; 85 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 86 Mat_SeqSBAIJ *b; 87 int ierr,*ai=a->i,*aj,m=A->M,n=A->N,i,j, 88 *bi,*bj,*rowlengths; 89 PetscScalar *av,*bv; 90 91 PetscFunctionBegin; 92 if (n != m) SETERRQ(PETSC_ERR_ARG_WRONG,"Matrix must be square"); 93 if (!a->diag){ 94 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 95 } 96 97 ierr = PetscMalloc(m*sizeof(int),&rowlengths);CHKERRQ(ierr); 98 for (i=0; i<m; i++) { 99 rowlengths[i] = ai[i+1] - a->diag[i]; 100 } 101 ierr = MatCreateSeqSBAIJ(PETSC_COMM_SELF,1,m,n,0,rowlengths,&B);CHKERRQ(ierr); 102 103 ierr = MatSetOption(B,MAT_ROW_ORIENTED);CHKERRQ(ierr); 104 ierr = MatSetOption(B,MAT_ROWS_SORTED);CHKERRQ(ierr); 105 ierr = MatSetOption(B,MAT_COLUMNS_SORTED);CHKERRQ(ierr); 106 107 b = (Mat_SeqSBAIJ*)(B->data); 108 bi = b->i; 109 bj = b->j; 110 bv = b->a; 111 112 bi[0] = 0; 113 for (i=0; i<m; i++) { 114 aj = a->j + a->diag[i]; 115 av = a->a + a->diag[i]; 116 for (j=0; j<rowlengths[i]; j++){ 117 *bj = *aj; bj++; aj++; 118 *bv = *av; bv++; av++; 119 } 120 bi[i+1] = bi[i] + rowlengths[i]; 121 b->ilen[i] = rowlengths[i]; 122 } 123 124 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 125 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 126 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 127 128 /* Fake support for "inplace" convert. */ 129 if (*newmat == A) { 130 ierr = MatDestroy(A);CHKERRQ(ierr); 131 } 132 *newmat = B; 133 134 PetscFunctionReturn(0); 135 } 136 EXTERN_C_END 137 138