1 2 #include "src/mat/impls/aij/seq/aij.h" 3 #include "src/mat/impls/baij/seq/baij.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 PetscErrorCode 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 PetscErrorCode ierr; 15 int *ai=a->i,*aj=a->j,m=A->m,n=A->n,i,j,k,*bi,*bj, 16 *rowlengths,nz,*rowstart,itmp; 17 int bs=A->bs,bs2=bs*bs,mbs=A->m/bs; 18 PetscScalar *av,*bv; 19 20 PetscFunctionBegin; 21 22 /* compute rowlengths of newmat */ 23 ierr = PetscMalloc((2*m+1)*sizeof(int),&rowlengths);CHKERRQ(ierr); 24 rowstart = rowlengths + m; 25 26 for (i=0; i<mbs; i++) rowlengths[i*bs] = 0; 27 aj = a->j; 28 k = 0; 29 for (i=0; i<mbs; i++) { 30 nz = ai[i+1] - ai[i]; 31 aj++; /* skip diagonal */ 32 for (j=1; j<nz; j++) { /* no. of lower triangular blocks */ 33 rowlengths[(*aj)*bs]++; aj++; 34 } 35 rowlengths[k] += nz; /* no. of upper triangular blocks */ 36 rowlengths[k] *= bs; 37 for (j=1; j<bs; j++) { 38 rowlengths[k+j] = rowlengths[k]; 39 } 40 k += bs; 41 /* printf(" rowlengths[%d]: %d\n",i, rowlengths[i]); */ 42 } 43 44 ierr = MatCreate(A->comm,m,n,m,n,&B);CHKERRQ(ierr); 45 ierr = MatSetType(B,newtype);CHKERRQ(ierr); 46 ierr = MatSeqAIJSetPreallocation(B,0,rowlengths);CHKERRQ(ierr); 47 ierr = MatSetOption(B,MAT_COLUMN_ORIENTED);CHKERRQ(ierr); 48 ierr = MatSetOption(B,MAT_ROWS_SORTED);CHKERRQ(ierr); 49 ierr = MatSetOption(B,MAT_COLUMNS_SORTED);CHKERRQ(ierr); 50 B->bs = A->bs; 51 52 b = (Mat_SeqAIJ*)(B->data); 53 bi = b->i; 54 bj = b->j; 55 bv = b->a; 56 57 /* set b->i */ 58 bi[0] = 0; rowstart[0] = 0; 59 for (i=0; i<mbs; i++){ 60 for (j=0; j<bs; j++){ 61 b->ilen[i*bs+j] = rowlengths[i*bs]; 62 rowstart[i*bs+j+1] = rowstart[i*bs+j] + rowlengths[i*bs]; 63 } 64 bi[i+1] = bi[i] + rowlengths[i*bs]/bs; 65 } 66 if (bi[mbs] != 2*a->nz - mbs) SETERRQ2(PETSC_ERR_PLIB,"bi[mbs]: %D != 2*a->nz-mbs: %D\n",bi[mbs],2*a->nz - mbs); 67 68 /* set b->j and b->a */ 69 aj = a->j; av = a->a; 70 for (i=0; i<mbs; i++) { 71 /* diagonal block */ 72 for (j=0; j<bs; j++){ /* row i*bs+j */ 73 itmp = i*bs+j; 74 for (k=0; k<bs; k++){ /* col i*bs+k */ 75 *(bj + rowstart[itmp]) = (*aj)*bs+k; 76 *(bv + rowstart[itmp]) = *(av+k*bs+j); 77 rowstart[itmp]++; 78 } 79 } 80 aj++; av += bs2; 81 82 nz = ai[i+1] - ai[i] -1; 83 while (nz--){ 84 /* lower triangular blocks */ 85 for (j=0; j<bs; j++){ /* row (*aj)*bs+j */ 86 itmp = (*aj)*bs+j; 87 for (k=0; k<bs; k++){ /* col i*bs+k */ 88 *(bj + rowstart[itmp]) = i*bs+k; 89 *(bv + rowstart[itmp]) = *(av+k*bs+j); 90 rowstart[itmp]++; 91 } 92 } 93 /* upper triangular blocks */ 94 for (j=0; j<bs; j++){ /* row i*bs+j */ 95 itmp = i*bs+j; 96 for (k=0; k<bs; k++){ /* col (*aj)*bs+k */ 97 *(bj + rowstart[itmp]) = (*aj)*bs+k; 98 *(bv + rowstart[itmp]) = *(av+k*bs+j); 99 rowstart[itmp]++; 100 } 101 } 102 aj++; av += bs2; 103 } 104 } 105 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 106 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 107 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 108 109 /* Fake support for "inplace" convert. */ 110 if (*newmat == A) { 111 ierr = MatDestroy(A);CHKERRQ(ierr); 112 } 113 *newmat = B; 114 PetscFunctionReturn(0); 115 } 116 #undef __FUNCT__ 117 #define __FUNCT__ "MatConvert_SeqAIJ_SeqSBAIJ" 118 PetscErrorCode MatConvert_SeqAIJ_SeqSBAIJ(Mat A,const MatType newtype,Mat *newmat) { 119 Mat B; 120 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 121 Mat_SeqSBAIJ *b; 122 PetscErrorCode ierr; 123 int *ai=a->i,*aj,m=A->M,n=A->N,i,j, 124 *bi,*bj,*rowlengths; 125 PetscScalar *av,*bv; 126 127 PetscFunctionBegin; 128 if (n != m) SETERRQ(PETSC_ERR_ARG_WRONG,"Matrix must be square"); 129 ierr = MatMissingDiagonal_SeqAIJ(A);CHKERRQ(ierr); /* check for missing diagonals, then mark diag */ 130 131 ierr = PetscMalloc(m*sizeof(int),&rowlengths);CHKERRQ(ierr); 132 for (i=0; i<m; i++) { 133 rowlengths[i] = ai[i+1] - a->diag[i]; 134 } 135 ierr = MatCreate(A->comm,m,n,m,n,&B);CHKERRQ(ierr); 136 ierr = MatSetType(B,newtype);CHKERRQ(ierr); 137 ierr = MatSeqSBAIJSetPreallocation(B,1,0,rowlengths);CHKERRQ(ierr); 138 139 ierr = MatSetOption(B,MAT_ROW_ORIENTED);CHKERRQ(ierr); 140 ierr = MatSetOption(B,MAT_ROWS_SORTED);CHKERRQ(ierr); 141 ierr = MatSetOption(B,MAT_COLUMNS_SORTED);CHKERRQ(ierr); 142 143 b = (Mat_SeqSBAIJ*)(B->data); 144 bi = b->i; 145 bj = b->j; 146 bv = b->a; 147 148 bi[0] = 0; 149 for (i=0; i<m; i++) { 150 aj = a->j + a->diag[i]; 151 av = a->a + a->diag[i]; 152 for (j=0; j<rowlengths[i]; j++){ 153 *bj = *aj; bj++; aj++; 154 *bv = *av; bv++; av++; 155 } 156 bi[i+1] = bi[i] + rowlengths[i]; 157 b->ilen[i] = rowlengths[i]; 158 } 159 160 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 161 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 162 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 163 164 /* Fake support for "inplace" convert. */ 165 if (*newmat == A) { 166 ierr = MatDestroy(A);CHKERRQ(ierr); 167 } 168 *newmat = B; 169 170 PetscFunctionReturn(0); 171 } 172 EXTERN_C_END 173 174 EXTERN_C_BEGIN 175 #undef __FUNCT__ 176 #define __FUNCT__ "MatConvert_SeqSBAI_SeqBAIJ" 177 PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat A,const MatType newtype,Mat *newmat) 178 { 179 Mat B; 180 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 181 Mat_SeqBAIJ *b; 182 PetscErrorCode ierr; 183 int *ai=a->i,*aj=a->j,m=A->m,n=A->n,i,k,*bi,*bj, 184 *browlengths,nz,*browstart,itmp; 185 int bs=A->bs,bs2=bs*bs,mbs=m/bs; 186 PetscScalar *av,*bv; 187 188 PetscFunctionBegin; 189 /* compute browlengths of newmat */ 190 ierr = PetscMalloc(2*mbs*sizeof(int),&browlengths);CHKERRQ(ierr); 191 browstart = browlengths + mbs; 192 for (i=0; i<mbs; i++) browlengths[i] = 0; 193 aj = a->j; 194 for (i=0; i<mbs; i++) { 195 nz = ai[i+1] - ai[i]; 196 aj++; /* skip diagonal */ 197 for (k=1; k<nz; k++) { /* no. of lower triangular blocks */ 198 browlengths[*aj]++; aj++; 199 } 200 browlengths[i] += nz; /* no. of upper triangular blocks */ 201 } 202 203 ierr = MatCreate(A->comm,m,n,m,n,&B);CHKERRQ(ierr); 204 ierr = MatSetType(B,newtype);CHKERRQ(ierr); 205 ierr = MatSeqBAIJSetPreallocation(B,bs,0,browlengths);CHKERRQ(ierr); 206 ierr = MatSetOption(B,MAT_ROW_ORIENTED);CHKERRQ(ierr); 207 ierr = MatSetOption(B,MAT_ROWS_SORTED);CHKERRQ(ierr); 208 ierr = MatSetOption(B,MAT_COLUMNS_SORTED);CHKERRQ(ierr); 209 210 b = (Mat_SeqBAIJ*)(B->data); 211 bi = b->i; 212 bj = b->j; 213 bv = b->a; 214 215 /* set b->i */ 216 bi[0] = 0; 217 for (i=0; i<mbs; i++){ 218 b->ilen[i] = browlengths[i]; 219 bi[i+1] = bi[i] + browlengths[i]; 220 browstart[i] = bi[i]; 221 } 222 if (bi[mbs] != 2*a->nz - mbs) SETERRQ2(PETSC_ERR_PLIB,"bi[mbs]: %D != 2*a->nz - mbs: %D\n",bi[mbs],2*a->nz - mbs); 223 224 /* set b->j and b->a */ 225 aj = a->j; av = a->a; 226 for (i=0; i<mbs; i++) { 227 /* diagonal block */ 228 *(bj + browstart[i]) = *aj; aj++; 229 itmp = bs2*browstart[i]; 230 for (k=0; k<bs2; k++){ 231 *(bv + itmp + k) = *av; av++; 232 } 233 browstart[i]++; 234 235 nz = ai[i+1] - ai[i] -1; 236 while (nz--){ 237 /* lower triangular blocks */ 238 *(bj + browstart[*aj]) = i; 239 itmp = bs2*browstart[*aj]; 240 for (k=0; k<bs2; k++){ 241 *(bv + itmp + k) = *(av + k); 242 } 243 browstart[*aj]++; 244 245 /* upper triangular blocks */ 246 *(bj + browstart[i]) = *aj; aj++; 247 itmp = bs2*browstart[i]; 248 for (k=0; k<bs2; k++){ 249 *(bv + itmp + k) = *av; av++; 250 } 251 browstart[i]++; 252 } 253 } 254 ierr = PetscFree(browlengths);CHKERRQ(ierr); 255 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 256 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 257 258 /* Fake support for "inplace" convert. */ 259 if (*newmat == A) { 260 ierr = MatDestroy(A);CHKERRQ(ierr); 261 } 262 *newmat = B; 263 PetscFunctionReturn(0); 264 } 265 #undef __FUNCT__ 266 #define __FUNCT__ "MatConvert_SeqBAIJ_SeqSBAIJ" 267 PetscErrorCode MatConvert_SeqBAIJ_SeqSBAIJ(Mat A,const MatType newtype,Mat *newmat) 268 { 269 Mat B; 270 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 271 Mat_SeqSBAIJ *b; 272 PetscErrorCode ierr; 273 int *ai=a->i,*aj,m=A->m,n=A->n,i,j,k, 274 *bi,*bj,*browlengths; 275 int bs=A->bs,bs2=bs*bs,mbs=m/bs; 276 PetscScalar *av,*bv; 277 278 PetscFunctionBegin; 279 if (n != m) SETERRQ(PETSC_ERR_ARG_WRONG,"Matrix must be square"); 280 ierr = MatMissingDiagonal_SeqBAIJ(A);CHKERRQ(ierr); /* check for missing diagonals, then mark diag */ 281 282 ierr = PetscMalloc(mbs*sizeof(int),&browlengths);CHKERRQ(ierr); 283 for (i=0; i<mbs; i++) { 284 browlengths[i] = ai[i+1] - a->diag[i]; 285 } 286 287 ierr = MatCreate(A->comm,m,n,m,n,&B);CHKERRQ(ierr); 288 ierr = MatSetType(B,newtype);CHKERRQ(ierr); 289 ierr = MatSeqSBAIJSetPreallocation(B,bs,0,browlengths);CHKERRQ(ierr); 290 ierr = MatSetOption(B,MAT_ROW_ORIENTED);CHKERRQ(ierr); 291 ierr = MatSetOption(B,MAT_ROWS_SORTED);CHKERRQ(ierr); 292 ierr = MatSetOption(B,MAT_COLUMNS_SORTED);CHKERRQ(ierr); 293 294 b = (Mat_SeqSBAIJ*)(B->data); 295 bi = b->i; 296 bj = b->j; 297 bv = b->a; 298 299 bi[0] = 0; 300 for (i=0; i<mbs; i++) { 301 aj = a->j + a->diag[i]; 302 av = a->a + (a->diag[i])*bs2; 303 for (j=0; j<browlengths[i]; j++){ 304 *bj = *aj; bj++; aj++; 305 for (k=0; k<bs2; k++){ 306 *bv = *av; bv++; av++; 307 } 308 } 309 bi[i+1] = bi[i] + browlengths[i]; 310 b->ilen[i] = browlengths[i]; 311 } 312 ierr = PetscFree(browlengths);CHKERRQ(ierr); 313 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 314 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 315 316 /* Fake support for "inplace" convert. */ 317 if (*newmat == A) { 318 ierr = MatDestroy(A);CHKERRQ(ierr); 319 } 320 *newmat = B; 321 322 PetscFunctionReturn(0); 323 } 324 EXTERN_C_END 325