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,*ai=a->i,*aj=a->j,m=A->m,n=A->n,i,j,k,*bi,*bj, 15 *rowlengths,nz,*rowstart,itmp; 16 int bs=a->bs,bs2=bs*bs,mbs=A->m/bs; 17 PetscScalar *av,*bv; 18 19 PetscFunctionBegin; 20 21 /* compute rowlengths of newmat */ 22 ierr = PetscMalloc((2*m+1)*sizeof(int),&rowlengths);CHKERRQ(ierr); 23 rowstart = rowlengths + m; 24 25 for (i=0; i<mbs; i++) rowlengths[i*bs] = 0; 26 aj = a->j; 27 k = 0; 28 for (i=0; i<mbs; i++) { 29 nz = ai[i+1] - ai[i]; 30 aj++; /* skip diagonal */ 31 for (j=1; j<nz; j++) { /* no. of lower triangular blocks */ 32 rowlengths[(*aj)*bs]++; aj++; 33 } 34 rowlengths[k] += nz; /* no. of upper triangular blocks */ 35 rowlengths[k] *= bs; 36 for (j=1; j<bs; j++) { 37 rowlengths[k+j] = rowlengths[k]; 38 } 39 k += bs; 40 /* printf(" rowlengths[%d]: %d\n",i, rowlengths[i]); */ 41 } 42 43 ierr = MatCreateSeqAIJ(PETSC_COMM_SELF,m,n,0,rowlengths,&B);CHKERRQ(ierr); 44 ierr = MatSetOption(B,MAT_COLUMN_ORIENTED);CHKERRQ(ierr); 45 ierr = MatSetOption(B,MAT_ROWS_SORTED);CHKERRQ(ierr); 46 ierr = MatSetOption(B,MAT_COLUMNS_SORTED);CHKERRQ(ierr); 47 48 b = (Mat_SeqAIJ*)(B->data); 49 bi = b->i; 50 bj = b->j; 51 bv = b->a; 52 53 /* set b->i */ 54 bi[0] = 0; rowstart[0] = 0; 55 for (i=0; i<mbs; i++){ 56 for (j=0; j<bs; j++){ 57 b->ilen[i*bs+j] = rowlengths[i*bs]; 58 rowstart[i*bs+j+1] = rowstart[i*bs+j] + rowlengths[i*bs]; 59 } 60 bi[i+1] = bi[i] + rowlengths[i*bs]/bs; 61 } 62 if (bi[mbs] != 2*a->nz - mbs) SETERRQ2(1,"bi[mbs]: %d != 2*a->nz-mbs: %d\n",bi[mbs],2*a->nz - mbs); 63 64 /* set b->j and b->a */ 65 aj = a->j; av = a->a; 66 for (i=0; i<mbs; i++) { 67 /* diagonal block */ 68 for (j=0; j<bs; j++){ /* row i*bs+j */ 69 itmp = i*bs+j; 70 for (k=0; k<bs; k++){ /* col i*bs+k */ 71 *(bj + rowstart[itmp]) = (*aj)*bs+k; 72 *(bv + rowstart[itmp]) = *(av+k*bs+j); 73 rowstart[itmp]++; 74 } 75 } 76 aj++; av += bs2; 77 78 nz = ai[i+1] - ai[i] -1; 79 while (nz--){ 80 /* lower triangular blocks */ 81 for (j=0; j<bs; j++){ /* row (*aj)*bs+j */ 82 itmp = (*aj)*bs+j; 83 for (k=0; k<bs; k++){ /* col i*bs+k */ 84 *(bj + rowstart[itmp]) = i*bs+k; 85 *(bv + rowstart[itmp]) = *(av+k*bs+j); 86 rowstart[itmp]++; 87 } 88 } 89 /* upper triangular blocks */ 90 for (j=0; j<bs; j++){ /* row i*bs+j */ 91 itmp = i*bs+j; 92 for (k=0; k<bs; k++){ /* col (*aj)*bs+k */ 93 *(bj + rowstart[itmp]) = (*aj)*bs+k; 94 *(bv + rowstart[itmp]) = *(av+k*bs+j); 95 rowstart[itmp]++; 96 } 97 } 98 aj++; av += bs2; 99 } 100 } 101 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 102 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 103 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 104 105 /* Fake support for "inplace" convert. */ 106 if (*newmat == A) { 107 ierr = MatDestroy(A);CHKERRQ(ierr); 108 } 109 *newmat = B; 110 PetscFunctionReturn(0); 111 } 112 #undef __FUNCT__ 113 #define __FUNCT__ "MatConvert_SeqAIJ_SeqSBAIJ" 114 int MatConvert_SeqAIJ_SeqSBAIJ(Mat A,const MatType newtype,Mat *newmat) { 115 Mat B; 116 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 117 Mat_SeqSBAIJ *b; 118 int ierr,*ai=a->i,*aj,m=A->M,n=A->N,i,j, 119 *bi,*bj,*rowlengths; 120 PetscScalar *av,*bv; 121 122 PetscFunctionBegin; 123 if (n != m) SETERRQ(PETSC_ERR_ARG_WRONG,"Matrix must be square"); 124 if (!a->diag){ 125 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 126 } 127 128 ierr = PetscMalloc(m*sizeof(int),&rowlengths);CHKERRQ(ierr); 129 for (i=0; i<m; i++) { 130 rowlengths[i] = ai[i+1] - a->diag[i]; 131 } 132 ierr = MatCreateSeqSBAIJ(PETSC_COMM_SELF,1,m,n,0,rowlengths,&B);CHKERRQ(ierr); 133 134 ierr = MatSetOption(B,MAT_ROW_ORIENTED);CHKERRQ(ierr); 135 ierr = MatSetOption(B,MAT_ROWS_SORTED);CHKERRQ(ierr); 136 ierr = MatSetOption(B,MAT_COLUMNS_SORTED);CHKERRQ(ierr); 137 138 b = (Mat_SeqSBAIJ*)(B->data); 139 bi = b->i; 140 bj = b->j; 141 bv = b->a; 142 143 bi[0] = 0; 144 for (i=0; i<m; i++) { 145 aj = a->j + a->diag[i]; 146 av = a->a + a->diag[i]; 147 for (j=0; j<rowlengths[i]; j++){ 148 *bj = *aj; bj++; aj++; 149 *bv = *av; bv++; av++; 150 } 151 bi[i+1] = bi[i] + rowlengths[i]; 152 b->ilen[i] = rowlengths[i]; 153 } 154 155 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 156 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 157 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 158 159 /* Fake support for "inplace" convert. */ 160 if (*newmat == A) { 161 ierr = MatDestroy(A);CHKERRQ(ierr); 162 } 163 *newmat = B; 164 165 PetscFunctionReturn(0); 166 } 167 EXTERN_C_END 168 169