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