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