1 #include <../src/mat/impls/aij/seq/aij.h> 2 #include <../src/mat/impls/baij/seq/baij.h> 3 #include <../src/mat/impls/sbaij/seq/sbaij.h> 4 5 PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat) 6 { 7 Mat B; 8 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data; 9 Mat_SeqAIJ *b; 10 PetscInt *ai = a->i, *aj = a->j, m = A->rmap->N, n = A->cmap->n, i, j, k, *bi, *bj, *rowlengths, nz, *rowstart, itmp; 11 PetscInt bs = A->rmap->bs, bs2 = bs * bs, mbs = A->rmap->N / bs, diagcnt = 0; 12 MatScalar *av, *bv; 13 const int aconj = PetscDefined(USE_COMPLEX) && A->hermitian == PETSC_BOOL3_TRUE ? 1 : 0; 14 15 PetscFunctionBegin; 16 /* compute rowlengths of newmat */ 17 PetscCall(PetscMalloc2(m, &rowlengths, m + 1, &rowstart)); 18 19 for (i = 0; i < mbs; i++) rowlengths[i * bs] = 0; 20 k = 0; 21 for (i = 0; i < mbs; i++) { 22 nz = ai[i + 1] - ai[i]; 23 if (nz) { 24 rowlengths[k] += nz; /* no. of upper triangular blocks */ 25 if (*aj == i) { 26 aj++; 27 diagcnt++; 28 nz--; 29 } /* skip diagonal */ 30 for (j = 0; j < nz; j++) { /* no. of lower triangular blocks */ 31 rowlengths[(*aj) * bs]++; 32 aj++; 33 } 34 } 35 rowlengths[k] *= bs; 36 for (j = 1; j < bs; j++) rowlengths[k + j] = rowlengths[k]; 37 k += bs; 38 } 39 40 if (reuse != MAT_REUSE_MATRIX) { 41 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B)); 42 PetscCall(MatSetSizes(B, m, n, m, n)); 43 PetscCall(MatSetType(B, MATSEQAIJ)); 44 PetscCall(MatSeqAIJSetPreallocation(B, 0, rowlengths)); 45 PetscCall(MatSetBlockSize(B, A->rmap->bs)); 46 } else B = *newmat; 47 48 b = (Mat_SeqAIJ *)B->data; 49 bi = b->i; 50 bj = b->j; 51 bv = b->a; 52 53 /* set b->i */ 54 bi[0] = 0; 55 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 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); 64 65 /* set b->j and b->a */ 66 aj = a->j; 67 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++; 82 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++; 105 av += bs2; 106 } 107 } 108 PetscCall(PetscFree2(rowlengths, rowstart)); 109 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 110 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 111 112 if (reuse == MAT_INPLACE_MATRIX) { 113 PetscCall(MatHeaderReplace(A, &B)); 114 } else { 115 *newmat = B; 116 } 117 PetscFunctionReturn(PETSC_SUCCESS); 118 } 119 120 #include <../src/mat/impls/aij/seq/aij.h> 121 122 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqSBAIJ_Preallocate(Mat A, PetscInt **nnz) 123 { 124 Mat_SeqAIJ *Aa = (Mat_SeqAIJ *)A->data; 125 PetscInt m, n, bs = A->rmap->bs; 126 const PetscInt *ai = Aa->i, *aj = Aa->j; 127 128 PetscFunctionBegin; 129 PetscCall(MatGetSize(A, &m, &n)); 130 131 if (bs == 1) { 132 const PetscInt *adiag; 133 134 PetscCall(MatGetDiagonalMarkers_SeqAIJ(A, &adiag, NULL)); 135 PetscCall(PetscMalloc1(m, nnz)); 136 for (PetscInt i = 0; i < m; i++) { 137 if (adiag[i] == ai[i + 1]) { 138 (*nnz)[i] = 0; 139 for (PetscInt j = ai[i]; j < ai[i + 1]; j++) (*nnz)[i] += (aj[j] > i); 140 } else (*nnz)[i] = ai[i + 1] - adiag[i]; 141 } 142 } else { 143 PetscHSetIJ ht; 144 PetscHashIJKey key; 145 PetscBool missing; 146 147 PetscCall(PetscHSetIJCreate(&ht)); 148 PetscCall(PetscCalloc1(m / bs, nnz)); 149 for (PetscInt i = 0; i < m; i++) { 150 key.i = i / bs; 151 for (PetscInt k = ai[i]; k < ai[i + 1]; k++) { 152 key.j = aj[k] / bs; 153 if (key.j < key.i) continue; 154 PetscCall(PetscHSetIJQueryAdd(ht, key, &missing)); 155 if (missing) (*nnz)[key.i]++; 156 } 157 } 158 PetscCall(PetscHSetIJDestroy(&ht)); 159 } 160 PetscFunctionReturn(PETSC_SUCCESS); 161 } 162 163 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqSBAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat) 164 { 165 Mat B; 166 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 167 Mat_SeqSBAIJ *b; 168 PetscInt *ai = a->i, *aj, m = A->rmap->N, n = A->cmap->N, i, j, *bi, *bj, *rowlengths, bs = A->rmap->bs; 169 MatScalar *av, *bv; 170 PetscBool miss = PETSC_FALSE; 171 const PetscInt *adiag; 172 173 PetscFunctionBegin; 174 #if !defined(PETSC_USE_COMPLEX) 175 PetscCheck(A->symmetric == PETSC_BOOL3_TRUE, PetscObjectComm((PetscObject)A), PETSC_ERR_USER, "Matrix must be symmetric. Call MatSetOption(mat,MAT_SYMMETRIC,PETSC_TRUE)"); 176 #else 177 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)"); 178 #endif 179 PetscCheck(n == m, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Matrix must be square"); 180 PetscCall(MatGetDiagonalMarkers_SeqAIJ(A, &adiag, NULL)); 181 if (bs == 1) { 182 PetscCall(PetscMalloc1(m, &rowlengths)); 183 for (i = 0; i < m; i++) { 184 if (adiag[i] == ai[i + 1]) { /* missing diagonal */ 185 rowlengths[i] = (ai[i + 1] - ai[i]) + 1; /* allocate some extra space */ 186 miss = PETSC_TRUE; 187 } else { 188 rowlengths[i] = (ai[i + 1] - adiag[i]); 189 } 190 } 191 } else PetscCall(MatConvert_SeqAIJ_SeqSBAIJ_Preallocate(A, &rowlengths)); 192 if (reuse != MAT_REUSE_MATRIX) { 193 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B)); 194 PetscCall(MatSetSizes(B, m, n, m, n)); 195 PetscCall(MatSetType(B, MATSEQSBAIJ)); 196 PetscCall(MatSeqSBAIJSetPreallocation(B, bs, 0, rowlengths)); 197 } else B = *newmat; 198 199 if (bs == 1 && !miss) { 200 b = (Mat_SeqSBAIJ *)B->data; 201 bi = b->i; 202 bj = b->j; 203 bv = b->a; 204 205 bi[0] = 0; 206 for (i = 0; i < m; i++) { 207 aj = a->j + adiag[i]; 208 av = a->a + adiag[i]; 209 for (j = 0; j < rowlengths[i]; j++) { 210 *bj = *aj; 211 bj++; 212 aj++; 213 *bv = *av; 214 bv++; 215 av++; 216 } 217 bi[i + 1] = bi[i] + rowlengths[i]; 218 b->ilen[i] = rowlengths[i]; 219 } 220 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 221 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 222 } else { 223 PetscCall(MatSetOption(B, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE)); 224 /* reuse may not be equal to MAT_REUSE_MATRIX, but the basic converter will reallocate or replace newmat if this value is not used */ 225 /* if reuse is equal to MAT_INITIAL_MATRIX, it has been appropriately preallocated before */ 226 /* MAT_INPLACE_MATRIX, it will be replaced with MatHeaderReplace below */ 227 PetscCall(MatConvert_Basic(A, newtype, MAT_REUSE_MATRIX, &B)); 228 } 229 PetscCall(PetscFree(rowlengths)); 230 if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &B)); 231 else *newmat = B; 232 PetscFunctionReturn(PETSC_SUCCESS); 233 } 234 235 PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat) 236 { 237 Mat B; 238 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data; 239 Mat_SeqBAIJ *b; 240 PetscInt *ai = a->i, *aj = a->j, m = A->rmap->N, n = A->cmap->n, *bi, *bj, *browlengths, nz, *browstart, itmp; 241 PetscInt bs = A->rmap->bs, bs2 = bs * bs, mbs = m / bs; 242 MatScalar *av, *bv; 243 const int aconj = PetscDefined(USE_COMPLEX) && A->hermitian == PETSC_BOOL3_TRUE ? 1 : 0; 244 245 PetscFunctionBegin; 246 /* compute browlengths of newmat */ 247 PetscCall(PetscMalloc2(mbs, &browlengths, mbs, &browstart)); 248 for (PetscInt i = 0; i < mbs; i++) browlengths[i] = 0; 249 for (PetscInt i = 0; i < mbs; i++) { 250 nz = ai[i + 1] - ai[i]; 251 aj++; /* skip diagonal */ 252 for (PetscInt k = 1; k < nz; k++) { /* no. of lower triangular blocks */ 253 browlengths[*aj]++; 254 aj++; 255 } 256 browlengths[i] += nz; /* no. of upper triangular blocks */ 257 } 258 259 if (reuse != MAT_REUSE_MATRIX) { 260 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B)); 261 PetscCall(MatSetSizes(B, m, n, m, n)); 262 PetscCall(MatSetType(B, MATSEQBAIJ)); 263 PetscCall(MatSeqBAIJSetPreallocation(B, bs, 0, browlengths)); 264 } else B = *newmat; 265 266 b = (Mat_SeqBAIJ *)B->data; 267 bi = b->i; 268 bj = b->j; 269 bv = b->a; 270 271 /* set b->i */ 272 bi[0] = 0; 273 for (PetscInt i = 0; i < mbs; i++) { 274 b->ilen[i] = browlengths[i]; 275 bi[i + 1] = bi[i] + browlengths[i]; 276 browstart[i] = bi[i]; 277 } 278 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); 279 280 /* set b->j and b->a */ 281 aj = a->j; 282 av = a->a; 283 for (PetscInt i = 0; i < mbs; i++) { 284 /* diagonal block */ 285 *(bj + browstart[i]) = *aj; 286 aj++; 287 288 itmp = bs2 * browstart[i]; 289 for (PetscInt k = 0; k < bs2; k++) { 290 *(bv + itmp + k) = *av; 291 av++; 292 } 293 browstart[i]++; 294 295 nz = ai[i + 1] - ai[i] - 1; 296 while (nz--) { 297 /* lower triangular blocks - transpose blocks of A */ 298 *(bj + browstart[*aj]) = i; /* block col index */ 299 300 itmp = bs2 * browstart[*aj]; /* row index */ 301 for (PetscInt col = 0; col < bs; col++) { 302 PetscInt k = col; 303 for (PetscInt row = 0; row < bs; row++) { 304 bv[itmp + col * bs + row] = aconj ? PetscConj(av[k]) : av[k]; 305 k += bs; 306 } 307 } 308 browstart[*aj]++; 309 310 /* upper triangular blocks */ 311 *(bj + browstart[i]) = *aj; 312 aj++; 313 314 itmp = bs2 * browstart[i]; 315 for (PetscInt k = 0; k < bs2; k++) bv[itmp + k] = av[k]; 316 av += bs2; 317 browstart[i]++; 318 } 319 } 320 PetscCall(PetscFree2(browlengths, browstart)); 321 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 322 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 323 324 if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &B)); 325 else *newmat = B; 326 PetscFunctionReturn(PETSC_SUCCESS); 327 } 328 329 PETSC_INTERN PetscErrorCode MatConvert_SeqBAIJ_SeqSBAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat) 330 { 331 Mat B; 332 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data; 333 Mat_SeqSBAIJ *b; 334 PetscInt *ai = a->i, *aj, m = A->rmap->N, n = A->cmap->n, k, *bi, *bj, *browlengths; 335 PetscInt bs = A->rmap->bs, bs2 = bs * bs, mbs = m / bs; 336 MatScalar *av, *bv; 337 const PetscInt *adiag; 338 339 PetscFunctionBegin; 340 PetscCheck(A->symmetric == PETSC_BOOL3_TRUE, PetscObjectComm((PetscObject)A), PETSC_ERR_USER, "Matrix must be symmetric. Call MatSetOption(mat,MAT_SYMMETRIC,PETSC_TRUE)"); 341 PetscCheck(n == m, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Matrix must be square"); 342 PetscCall(MatGetDiagonalMarkers_SeqBAIJ(A, &adiag, NULL)); 343 PetscCall(PetscMalloc1(mbs, &browlengths)); 344 for (PetscInt i = 0; i < mbs; i++) browlengths[i] = ai[i + 1] - adiag[i]; 345 346 if (reuse != MAT_REUSE_MATRIX) { 347 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B)); 348 PetscCall(MatSetSizes(B, m, n, m, n)); 349 PetscCall(MatSetType(B, MATSEQSBAIJ)); 350 PetscCall(MatSeqSBAIJSetPreallocation(B, bs, 0, browlengths)); 351 } else B = *newmat; 352 353 b = (Mat_SeqSBAIJ *)B->data; 354 bi = b->i; 355 bj = b->j; 356 bv = b->a; 357 358 bi[0] = 0; 359 for (PetscInt i = 0; i < mbs; i++) { 360 aj = a->j + adiag[i]; 361 av = a->a + (adiag[i]) * bs2; 362 for (PetscInt j = 0; j < browlengths[i]; j++) { 363 *bj = *aj; 364 bj++; 365 aj++; 366 for (k = 0; k < bs2; k++) { 367 *bv = *av; 368 bv++; 369 av++; 370 } 371 } 372 bi[i + 1] = bi[i] + browlengths[i]; 373 b->ilen[i] = browlengths[i]; 374 } 375 PetscCall(PetscFree(browlengths)); 376 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 377 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 378 379 if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &B)); 380 else *newmat = B; 381 PetscFunctionReturn(PETSC_SUCCESS); 382 } 383