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