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 #if defined(PETSC_USE_COMPLEX) 16 const int aconj = A->hermitian ? 1 : 0; 17 #else 18 const int aconj = 0; 19 #endif 20 21 PetscFunctionBegin; 22 /* compute rowlengths of newmat */ 23 ierr = PetscMalloc2(m,&rowlengths,m+1,&rowstart);CHKERRQ(ierr); 24 25 for (i=0; i<mbs; i++) rowlengths[i*bs] = 0; 26 k = 0; 27 for (i=0; i<mbs; i++) { 28 nz = ai[i+1] - ai[i]; 29 if (nz) { 30 rowlengths[k] += nz; /* no. of upper triangular blocks */ 31 if (*aj == i) {aj++;diagcnt++;nz--;} /* skip diagonal */ 32 for (j=0; j<nz; j++) { /* no. of lower triangular blocks */ 33 rowlengths[(*aj)*bs]++; aj++; 34 } 35 } 36 rowlengths[k] *= bs; 37 for (j=1; j<bs; j++) { 38 rowlengths[k+j] = rowlengths[k]; 39 } 40 k += bs; 41 } 42 43 if (reuse != MAT_REUSE_MATRIX) { 44 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 45 ierr = MatSetSizes(B,m,n,m,n);CHKERRQ(ierr); 46 ierr = MatSetType(B,MATSEQAIJ);CHKERRQ(ierr); 47 ierr = MatSeqAIJSetPreallocation(B,0,rowlengths);CHKERRQ(ierr); 48 ierr = MatSetBlockSize(B,A->rmap->bs);CHKERRQ(ierr); 49 } else B = *newmat; 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]) = aconj ? PetscConj(*(av+j*bs+k)) : *(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_INPLACE_MATRIX) { 112 ierr = MatHeaderReplace(A,&B);CHKERRQ(ierr); 113 } else { 114 *newmat = B; 115 } 116 PetscFunctionReturn(0); 117 } 118 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,bs=PetscAbs(A->rmap->bs); 126 MatScalar *av,*bv; 127 PetscBool miss = PETSC_FALSE; 128 129 PetscFunctionBegin; 130 #if !defined(PETSC_USE_COMPLEX) 131 if (!A->symmetric) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Matrix must be symmetric. Call MatSetOption(mat,MAT_SYMMETRIC,PETSC_TRUE)"); 132 #else 133 if (!A->symmetric && !A->hermitian) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Matrix must be either symmetric or hermitian. Call MatSetOption(mat,MAT_SYMMETRIC,PETSC_TRUE) and/or MatSetOption(mat,MAT_HERMITIAN,PETSC_TRUE)"); 134 #endif 135 if (n != m) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix must be square"); 136 137 ierr = PetscMalloc1(m/bs,&rowlengths);CHKERRQ(ierr); 138 for (i=0; i<m/bs; i++) { 139 if (a->diag[i*bs] == ai[i*bs+1]) { /* missing diagonal */ 140 rowlengths[i] = (ai[i*bs+1] - ai[i*bs])/bs; /* allocate some extra space */ 141 miss = PETSC_TRUE; 142 } else { 143 rowlengths[i] = (ai[i*bs+1] - a->diag[i*bs])/bs; 144 } 145 } 146 if (reuse != MAT_REUSE_MATRIX) { 147 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 148 ierr = MatSetSizes(B,m,n,m,n);CHKERRQ(ierr); 149 ierr = MatSetType(B,MATSEQSBAIJ);CHKERRQ(ierr); 150 ierr = MatSeqSBAIJSetPreallocation(B,bs,0,rowlengths);CHKERRQ(ierr); 151 } else B = *newmat; 152 153 if (bs == 1 && !miss) { 154 b = (Mat_SeqSBAIJ*)(B->data); 155 bi = b->i; 156 bj = b->j; 157 bv = b->a; 158 159 bi[0] = 0; 160 for (i=0; i<m; i++) { 161 aj = a->j + a->diag[i]; 162 av = a->a + a->diag[i]; 163 for (j=0; j<rowlengths[i]; j++) { 164 *bj = *aj; bj++; aj++; 165 *bv = *av; bv++; av++; 166 } 167 bi[i+1] = bi[i] + rowlengths[i]; 168 b->ilen[i] = rowlengths[i]; 169 } 170 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 171 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 172 } else { 173 ierr = MatSetOption(B,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE);CHKERRQ(ierr); 174 /* reuse may not be equal to MAT_REUSE_MATRIX, but the basic converter will reallocate or replace newmat if this value is not used */ 175 /* if reuse is equal to MAT_INITIAL_MATRIX, it has been appropriately preallocated before */ 176 /* MAT_INPLACE_MATRIX, it will be replaced with MatHeaderReplace below */ 177 ierr = MatConvert_Basic(A,newtype,MAT_REUSE_MATRIX,&B);CHKERRQ(ierr); 178 } 179 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 180 if (reuse == MAT_INPLACE_MATRIX) { 181 ierr = MatHeaderReplace(A,&B);CHKERRQ(ierr); 182 } else *newmat = B; 183 PetscFunctionReturn(0); 184 } 185 186 PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 187 { 188 Mat B; 189 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 190 Mat_SeqBAIJ *b; 191 PetscErrorCode ierr; 192 PetscInt *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->n,i,k,*bi,*bj,*browlengths,nz,*browstart,itmp; 193 PetscInt bs =A->rmap->bs,bs2=bs*bs,mbs=m/bs,col,row; 194 MatScalar *av,*bv; 195 #if defined(PETSC_USE_COMPLEX) 196 const int aconj = A->hermitian ? 1 : 0; 197 #else 198 const int aconj = 0; 199 #endif 200 201 PetscFunctionBegin; 202 /* compute browlengths of newmat */ 203 ierr = PetscMalloc2(mbs,&browlengths,mbs,&browstart);CHKERRQ(ierr); 204 for (i=0; i<mbs; i++) browlengths[i] = 0; 205 for (i=0; i<mbs; i++) { 206 nz = ai[i+1] - ai[i]; 207 aj++; /* skip diagonal */ 208 for (k=1; k<nz; k++) { /* no. of lower triangular blocks */ 209 browlengths[*aj]++; aj++; 210 } 211 browlengths[i] += nz; /* no. of upper triangular blocks */ 212 } 213 214 if (reuse != MAT_REUSE_MATRIX) { 215 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 216 ierr = MatSetSizes(B,m,n,m,n);CHKERRQ(ierr); 217 ierr = MatSetType(B,MATSEQBAIJ);CHKERRQ(ierr); 218 ierr = MatSeqBAIJSetPreallocation(B,bs,0,browlengths);CHKERRQ(ierr); 219 } else B = *newmat; 220 221 b = (Mat_SeqBAIJ*)(B->data); 222 bi = b->i; 223 bj = b->j; 224 bv = b->a; 225 226 /* set b->i */ 227 bi[0] = 0; 228 for (i=0; i<mbs; i++) { 229 b->ilen[i] = browlengths[i]; 230 bi[i+1] = bi[i] + browlengths[i]; 231 browstart[i] = bi[i]; 232 } 233 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); 234 235 /* set b->j and b->a */ 236 aj = a->j; av = a->a; 237 for (i=0; i<mbs; i++) { 238 /* diagonal block */ 239 *(bj + browstart[i]) = *aj; aj++; 240 241 itmp = bs2*browstart[i]; 242 for (k=0; k<bs2; k++) { 243 *(bv + itmp + k) = *av; av++; 244 } 245 browstart[i]++; 246 247 nz = ai[i+1] - ai[i] -1; 248 while (nz--) { 249 /* lower triangular blocks - transpose blocks of A */ 250 *(bj + browstart[*aj]) = i; /* block col index */ 251 252 itmp = bs2*browstart[*aj]; /* row index */ 253 for (col=0; col<bs; col++) { 254 k = col; 255 for (row=0; row<bs; row++) { 256 bv[itmp + col*bs+row] = aconj ? PetscConj(av[k]) : av[k]; 257 k+=bs; 258 } 259 } 260 browstart[*aj]++; 261 262 /* upper triangular blocks */ 263 *(bj + browstart[i]) = *aj; aj++; 264 265 itmp = bs2*browstart[i]; 266 for (k=0; k<bs2; k++) { 267 bv[itmp + k] = av[k]; 268 } 269 av += bs2; 270 browstart[i]++; 271 } 272 } 273 ierr = PetscFree2(browlengths,browstart);CHKERRQ(ierr); 274 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 275 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 276 277 if (reuse == MAT_INPLACE_MATRIX) { 278 ierr = MatHeaderReplace(A,&B);CHKERRQ(ierr); 279 } else *newmat = B; 280 PetscFunctionReturn(0); 281 } 282 283 PETSC_INTERN PetscErrorCode MatConvert_SeqBAIJ_SeqSBAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 284 { 285 Mat B; 286 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 287 Mat_SeqSBAIJ *b; 288 PetscErrorCode ierr; 289 PetscInt *ai=a->i,*aj,m=A->rmap->N,n=A->cmap->n,i,j,k,*bi,*bj,*browlengths; 290 PetscInt bs =A->rmap->bs,bs2=bs*bs,mbs=m/bs,dd; 291 MatScalar *av,*bv; 292 PetscBool flg; 293 294 PetscFunctionBegin; 295 if (!A->symmetric) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Matrix must be symmetric. Call MatSetOption(mat,MAT_SYMMETRIC,PETSC_TRUE)"); 296 if (n != m) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix must be square"); 297 ierr = MatMissingDiagonal_SeqBAIJ(A,&flg,&dd);CHKERRQ(ierr); /* check for missing diagonals, then mark diag */ 298 if (flg) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal %D",dd); 299 300 ierr = PetscMalloc1(mbs,&browlengths);CHKERRQ(ierr); 301 for (i=0; i<mbs; i++) { 302 browlengths[i] = ai[i+1] - a->diag[i]; 303 } 304 305 if (reuse != MAT_REUSE_MATRIX) { 306 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 307 ierr = MatSetSizes(B,m,n,m,n);CHKERRQ(ierr); 308 ierr = MatSetType(B,MATSEQSBAIJ);CHKERRQ(ierr); 309 ierr = MatSeqSBAIJSetPreallocation(B,bs,0,browlengths);CHKERRQ(ierr); 310 } else B = *newmat; 311 312 b = (Mat_SeqSBAIJ*)(B->data); 313 bi = b->i; 314 bj = b->j; 315 bv = b->a; 316 317 bi[0] = 0; 318 for (i=0; i<mbs; i++) { 319 aj = a->j + a->diag[i]; 320 av = a->a + (a->diag[i])*bs2; 321 for (j=0; j<browlengths[i]; j++) { 322 *bj = *aj; bj++; aj++; 323 for (k=0; k<bs2; k++) { 324 *bv = *av; bv++; av++; 325 } 326 } 327 bi[i+1] = bi[i] + browlengths[i]; 328 b->ilen[i] = browlengths[i]; 329 } 330 ierr = PetscFree(browlengths);CHKERRQ(ierr); 331 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 332 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 333 334 if (reuse == MAT_INPLACE_MATRIX) { 335 ierr = MatHeaderReplace(A,&B);CHKERRQ(ierr); 336 } else *newmat = B; 337 PetscFunctionReturn(0); 338 } 339