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