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_SeqSBAIJ_SeqAIJ" 9 PetscErrorCode MatConvert_SeqSBAIJ_SeqAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 10 { 11 Mat B; 12 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 13 Mat_SeqAIJ *b; 14 PetscErrorCode ierr; 15 PetscInt *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->n,i,j,k,*bi,*bj,*rowlengths,nz,*rowstart,itmp; 16 PetscInt bs =A->rmap->bs,bs2=bs*bs,mbs=A->rmap->N/bs,diagcnt=0; 17 MatScalar *av,*bv; 18 19 PetscFunctionBegin; 20 /* compute rowlengths of newmat */ 21 ierr = PetscMalloc2(m,PetscInt,&rowlengths,m+1,PetscInt,&rowstart);CHKERRQ(ierr); 22 23 for (i=0; i<mbs; i++) rowlengths[i*bs] = 0; 24 aj = a->j; 25 k = 0; 26 for (i=0; i<mbs; i++) { 27 nz = ai[i+1] - ai[i]; 28 if (nz) { 29 rowlengths[k] += nz; /* no. of upper triangular blocks */ 30 if (*aj == i) {aj++;diagcnt++;nz--;} /* skip diagonal */ 31 for (j=0; j<nz; j++) { /* no. of lower triangular blocks */ 32 rowlengths[(*aj)*bs]++; aj++; 33 } 34 } 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 = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 44 ierr = MatSetSizes(B,m,n,m,n);CHKERRQ(ierr); 45 ierr = MatSetType(B,MATSEQAIJ);CHKERRQ(ierr); 46 ierr = MatSeqAIJSetPreallocation(B,0,rowlengths);CHKERRQ(ierr); 47 ierr = MatSetOption(B,MAT_ROW_ORIENTED,PETSC_FALSE);CHKERRQ(ierr); 48 49 B->rmap->bs = A->rmap->bs; 50 51 b = (Mat_SeqAIJ*)(B->data); 52 bi = b->i; 53 bj = b->j; 54 bv = b->a; 55 56 /* set b->i */ 57 bi[0] = 0; rowstart[0] = 0; 58 for (i=0; i<mbs; i++) { 59 for (j=0; j<bs; j++) { 60 b->ilen[i*bs+j] = rowlengths[i*bs]; 61 rowstart[i*bs+j+1] = rowstart[i*bs+j] + rowlengths[i*bs]; 62 } 63 bi[i+1] = bi[i] + rowlengths[i*bs]/bs; 64 } 65 if (bi[mbs] != 2*a->nz - diagcnt) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"bi[mbs]: %D != 2*a->nz-diagcnt: %D\n",bi[mbs],2*a->nz - diagcnt); 66 67 /* set b->j and b->a */ 68 aj = a->j; av = a->a; 69 for (i=0; i<mbs; i++) { 70 nz = ai[i+1] - ai[i]; 71 /* diagonal block */ 72 if (nz && *aj == i) { 73 nz--; 74 for (j=0; j<bs; j++) { /* row i*bs+j */ 75 itmp = i*bs+j; 76 for (k=0; k<bs; k++) { /* col i*bs+k */ 77 *(bj + rowstart[itmp]) = (*aj)*bs+k; 78 *(bv + rowstart[itmp]) = *(av+k*bs+j); 79 rowstart[itmp]++; 80 } 81 } 82 aj++; av += bs2; 83 } 84 85 while (nz--) { 86 /* lower triangular blocks */ 87 for (j=0; j<bs; j++) { /* row (*aj)*bs+j */ 88 itmp = (*aj)*bs+j; 89 for (k=0; k<bs; k++) { /* col i*bs+k */ 90 *(bj + rowstart[itmp]) = i*bs+k; 91 *(bv + rowstart[itmp]) = *(av+j*bs+k); 92 rowstart[itmp]++; 93 } 94 } 95 /* upper triangular blocks */ 96 for (j=0; j<bs; j++) { /* row i*bs+j */ 97 itmp = i*bs+j; 98 for (k=0; k<bs; k++) { /* col (*aj)*bs+k */ 99 *(bj + rowstart[itmp]) = (*aj)*bs+k; 100 *(bv + rowstart[itmp]) = *(av+k*bs+j); 101 rowstart[itmp]++; 102 } 103 } 104 aj++; av += bs2; 105 } 106 } 107 ierr = PetscFree2(rowlengths,rowstart);CHKERRQ(ierr); 108 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 109 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 110 111 if (reuse == MAT_REUSE_MATRIX) { 112 ierr = MatHeaderReplace(A,B);CHKERRQ(ierr); 113 } else { 114 *newmat = B; 115 } 116 PetscFunctionReturn(0); 117 } 118 EXTERN_C_END 119 120 EXTERN_C_BEGIN 121 #undef __FUNCT__ 122 #define __FUNCT__ "MatConvert_SeqAIJ_SeqSBAIJ" 123 PetscErrorCode MatConvert_SeqAIJ_SeqSBAIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 124 { 125 Mat B; 126 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 127 Mat_SeqSBAIJ *b; 128 PetscErrorCode ierr; 129 PetscInt *ai=a->i,*aj,m=A->rmap->N,n=A->cmap->N,i,j,*bi,*bj,*rowlengths; 130 MatScalar *av,*bv; 131 132 PetscFunctionBegin; 133 if (!A->symmetric) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_USER,"Matrix must be symmetric. Call MatSetOption(mat,MAT_SYMMETRIC,PETSC_TRUE)"); 134 if (n != m) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix must be square"); 135 136 ierr = PetscMalloc(m*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 137 for (i=0; i<m; i++) { 138 rowlengths[i] = ai[i+1] - a->diag[i]; 139 } 140 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 141 ierr = MatSetSizes(B,m,n,m,n);CHKERRQ(ierr); 142 ierr = MatSetType(B,MATSEQSBAIJ);CHKERRQ(ierr); 143 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(B,1,0,rowlengths);CHKERRQ(ierr); 144 145 ierr = MatSetOption(B,MAT_ROW_ORIENTED,PETSC_TRUE);CHKERRQ(ierr); 146 147 b = (Mat_SeqSBAIJ*)(B->data); 148 bi = b->i; 149 bj = b->j; 150 bv = b->a; 151 152 bi[0] = 0; 153 for (i=0; i<m; i++) { 154 aj = a->j + a->diag[i]; 155 av = a->a + a->diag[i]; 156 for (j=0; j<rowlengths[i]; j++) { 157 *bj = *aj; bj++; aj++; 158 *bv = *av; bv++; av++; 159 } 160 bi[i+1] = bi[i] + rowlengths[i]; 161 b->ilen[i] = rowlengths[i]; 162 } 163 164 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 165 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 166 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 167 168 if (reuse == MAT_REUSE_MATRIX) { 169 ierr = MatHeaderReplace(A,B);CHKERRQ(ierr); 170 } else { 171 *newmat = B; 172 } 173 PetscFunctionReturn(0); 174 } 175 EXTERN_C_END 176 177 EXTERN_C_BEGIN 178 #undef __FUNCT__ 179 #define __FUNCT__ "MatConvert_SeqSBAIJ_SeqBAIJ" 180 PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 181 { 182 Mat B; 183 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 184 Mat_SeqBAIJ *b; 185 PetscErrorCode ierr; 186 PetscInt *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->n,i,k,*bi,*bj,*browlengths,nz,*browstart,itmp; 187 PetscInt bs =A->rmap->bs,bs2=bs*bs,mbs=m/bs,col,row; 188 MatScalar *av,*bv; 189 190 PetscFunctionBegin; 191 /* compute browlengths of newmat */ 192 ierr = PetscMalloc2(mbs,PetscInt,&browlengths,mbs,PetscInt,&browstart);CHKERRQ(ierr); 193 for (i=0; i<mbs; i++) browlengths[i] = 0; 194 aj = a->j; 195 for (i=0; i<mbs; i++) { 196 nz = ai[i+1] - ai[i]; 197 aj++; /* skip diagonal */ 198 for (k=1; k<nz; k++) { /* no. of lower triangular blocks */ 199 browlengths[*aj]++; aj++; 200 } 201 browlengths[i] += nz; /* no. of upper triangular blocks */ 202 } 203 204 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 205 ierr = MatSetSizes(B,m,n,m,n);CHKERRQ(ierr); 206 ierr = MatSetType(B,MATSEQBAIJ);CHKERRQ(ierr); 207 ierr = MatSeqBAIJSetPreallocation(B,bs,0,browlengths);CHKERRQ(ierr); 208 ierr = MatSetOption(B,MAT_ROW_ORIENTED,PETSC_TRUE);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_COMM_SELF,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 230 itmp = bs2*browstart[i]; 231 for (k=0; k<bs2; k++) { 232 *(bv + itmp + k) = *av; av++; 233 } 234 browstart[i]++; 235 236 nz = ai[i+1] - ai[i] -1; 237 while (nz--) { 238 /* lower triangular blocks - transpose blocks of A */ 239 *(bj + browstart[*aj]) = i; /* block col index */ 240 241 itmp = bs2*browstart[*aj]; /* row index */ 242 for (col=0; col<bs; col++) { 243 k = col; 244 for (row=0; row<bs; row++) { 245 bv[itmp + col*bs+row] = av[k]; k+=bs; 246 } 247 } 248 browstart[*aj]++; 249 250 /* upper triangular blocks */ 251 *(bj + browstart[i]) = *aj; aj++; 252 253 itmp = bs2*browstart[i]; 254 for (k=0; k<bs2; k++) { 255 bv[itmp + k] = av[k]; 256 } 257 av += bs2; 258 browstart[i]++; 259 } 260 } 261 ierr = PetscFree2(browlengths,browstart);CHKERRQ(ierr); 262 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 263 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 264 265 if (reuse == MAT_REUSE_MATRIX) { 266 ierr = MatHeaderReplace(A,B);CHKERRQ(ierr); 267 } else { 268 *newmat = B; 269 } 270 PetscFunctionReturn(0); 271 } 272 EXTERN_C_END 273 274 EXTERN_C_BEGIN 275 #undef __FUNCT__ 276 #define __FUNCT__ "MatConvert_SeqBAIJ_SeqSBAIJ" 277 PetscErrorCode MatConvert_SeqBAIJ_SeqSBAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 278 { 279 Mat B; 280 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 281 Mat_SeqSBAIJ *b; 282 PetscErrorCode ierr; 283 PetscInt *ai=a->i,*aj,m=A->rmap->N,n=A->cmap->n,i,j,k,*bi,*bj,*browlengths; 284 PetscInt bs =A->rmap->bs,bs2=bs*bs,mbs=m/bs,dd; 285 MatScalar *av,*bv; 286 PetscBool flg; 287 288 PetscFunctionBegin; 289 if (!A->symmetric) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_USER,"Matrix must be symmetric. Call MatSetOption(mat,MAT_SYMMETRIC,PETSC_TRUE)"); 290 if (n != m) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix must be square"); 291 ierr = MatMissingDiagonal_SeqBAIJ(A,&flg,&dd);CHKERRQ(ierr); /* check for missing diagonals, then mark diag */ 292 if (flg) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal %D",dd); 293 294 ierr = PetscMalloc(mbs*sizeof(PetscInt),&browlengths);CHKERRQ(ierr); 295 for (i=0; i<mbs; i++) { 296 browlengths[i] = ai[i+1] - a->diag[i]; 297 } 298 299 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 300 ierr = MatSetSizes(B,m,n,m,n);CHKERRQ(ierr); 301 ierr = MatSetType(B,MATSEQSBAIJ);CHKERRQ(ierr); 302 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(B,bs,0,browlengths);CHKERRQ(ierr); 303 ierr = MatSetOption(B,MAT_ROW_ORIENTED,PETSC_TRUE);CHKERRQ(ierr); 304 305 b = (Mat_SeqSBAIJ*)(B->data); 306 bi = b->i; 307 bj = b->j; 308 bv = b->a; 309 310 bi[0] = 0; 311 for (i=0; i<mbs; i++) { 312 aj = a->j + a->diag[i]; 313 av = a->a + (a->diag[i])*bs2; 314 for (j=0; j<browlengths[i]; j++) { 315 *bj = *aj; bj++; aj++; 316 for (k=0; k<bs2; k++) { 317 *bv = *av; bv++; av++; 318 } 319 } 320 bi[i+1] = bi[i] + browlengths[i]; 321 b->ilen[i] = browlengths[i]; 322 } 323 ierr = PetscFree(browlengths);CHKERRQ(ierr); 324 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 325 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 326 327 if (reuse == MAT_REUSE_MATRIX) { 328 ierr = MatHeaderReplace(A,B);CHKERRQ(ierr); 329 } else { 330 *newmat = B; 331 } 332 PetscFunctionReturn(0); 333 } 334 EXTERN_C_END 335