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