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