1 2 /* 3 Defines the basic matrix operations for the SBAIJ (compressed row) 4 matrix storage format. 5 */ 6 #include <../src/mat/impls/baij/seq/baij.h> /*I "petscmat.h" I*/ 7 #include <../src/mat/impls/sbaij/seq/sbaij.h> 8 #include <petscblaslapack.h> 9 10 #include <../src/mat/impls/sbaij/seq/relax.h> 11 #define USESHORT 12 #include <../src/mat/impls/sbaij/seq/relax.h> 13 14 #if defined(PETSC_HAVE_ELEMENTAL) 15 PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_Elemental(Mat, MatType, MatReuse, Mat *); 16 #endif 17 #if defined(PETSC_HAVE_SCALAPACK) 18 PETSC_INTERN PetscErrorCode MatConvert_SBAIJ_ScaLAPACK(Mat, MatType, MatReuse, Mat *); 19 #endif 20 PETSC_INTERN PetscErrorCode MatConvert_MPISBAIJ_Basic(Mat, MatType, MatReuse, Mat *); 21 22 /* 23 Checks for missing diagonals 24 */ 25 PetscErrorCode MatMissingDiagonal_SeqSBAIJ(Mat A, PetscBool *missing, PetscInt *dd) { 26 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data; 27 PetscInt *diag, *ii = a->i, i; 28 29 PetscFunctionBegin; 30 PetscCall(MatMarkDiagonal_SeqSBAIJ(A)); 31 *missing = PETSC_FALSE; 32 if (A->rmap->n > 0 && !ii) { 33 *missing = PETSC_TRUE; 34 if (dd) *dd = 0; 35 PetscCall(PetscInfo(A, "Matrix has no entries therefore is missing diagonal\n")); 36 } else { 37 diag = a->diag; 38 for (i = 0; i < a->mbs; i++) { 39 if (diag[i] >= ii[i + 1]) { 40 *missing = PETSC_TRUE; 41 if (dd) *dd = i; 42 break; 43 } 44 } 45 } 46 PetscFunctionReturn(0); 47 } 48 49 PetscErrorCode MatMarkDiagonal_SeqSBAIJ(Mat A) { 50 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data; 51 PetscInt i, j; 52 53 PetscFunctionBegin; 54 if (!a->diag) { 55 PetscCall(PetscMalloc1(a->mbs, &a->diag)); 56 a->free_diag = PETSC_TRUE; 57 } 58 for (i = 0; i < a->mbs; i++) { 59 a->diag[i] = a->i[i + 1]; 60 for (j = a->i[i]; j < a->i[i + 1]; j++) { 61 if (a->j[j] == i) { 62 a->diag[i] = j; 63 break; 64 } 65 } 66 } 67 PetscFunctionReturn(0); 68 } 69 70 static PetscErrorCode MatGetRowIJ_SeqSBAIJ(Mat A, PetscInt oshift, PetscBool symmetric, PetscBool blockcompressed, PetscInt *nn, const PetscInt *inia[], const PetscInt *inja[], PetscBool *done) { 71 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data; 72 PetscInt i, j, n = a->mbs, nz = a->i[n], *tia, *tja, bs = A->rmap->bs, k, l, cnt; 73 PetscInt **ia = (PetscInt **)inia, **ja = (PetscInt **)inja; 74 75 PetscFunctionBegin; 76 *nn = n; 77 if (!ia) PetscFunctionReturn(0); 78 if (symmetric) { 79 PetscCall(MatToSymmetricIJ_SeqAIJ(n, a->i, a->j, PETSC_FALSE, 0, 0, &tia, &tja)); 80 nz = tia[n]; 81 } else { 82 tia = a->i; 83 tja = a->j; 84 } 85 86 if (!blockcompressed && bs > 1) { 87 (*nn) *= bs; 88 /* malloc & create the natural set of indices */ 89 PetscCall(PetscMalloc1((n + 1) * bs, ia)); 90 if (n) { 91 (*ia)[0] = oshift; 92 for (j = 1; j < bs; j++) (*ia)[j] = (tia[1] - tia[0]) * bs + (*ia)[j - 1]; 93 } 94 95 for (i = 1; i < n; i++) { 96 (*ia)[i * bs] = (tia[i] - tia[i - 1]) * bs + (*ia)[i * bs - 1]; 97 for (j = 1; j < bs; j++) (*ia)[i * bs + j] = (tia[i + 1] - tia[i]) * bs + (*ia)[i * bs + j - 1]; 98 } 99 if (n) (*ia)[n * bs] = (tia[n] - tia[n - 1]) * bs + (*ia)[n * bs - 1]; 100 101 if (inja) { 102 PetscCall(PetscMalloc1(nz * bs * bs, ja)); 103 cnt = 0; 104 for (i = 0; i < n; i++) { 105 for (j = 0; j < bs; j++) { 106 for (k = tia[i]; k < tia[i + 1]; k++) { 107 for (l = 0; l < bs; l++) (*ja)[cnt++] = bs * tja[k] + l; 108 } 109 } 110 } 111 } 112 113 if (symmetric) { /* deallocate memory allocated in MatToSymmetricIJ_SeqAIJ() */ 114 PetscCall(PetscFree(tia)); 115 PetscCall(PetscFree(tja)); 116 } 117 } else if (oshift == 1) { 118 if (symmetric) { 119 nz = tia[A->rmap->n / bs]; 120 /* add 1 to i and j indices */ 121 for (i = 0; i < A->rmap->n / bs + 1; i++) tia[i] = tia[i] + 1; 122 *ia = tia; 123 if (ja) { 124 for (i = 0; i < nz; i++) tja[i] = tja[i] + 1; 125 *ja = tja; 126 } 127 } else { 128 nz = a->i[A->rmap->n / bs]; 129 /* malloc space and add 1 to i and j indices */ 130 PetscCall(PetscMalloc1(A->rmap->n / bs + 1, ia)); 131 for (i = 0; i < A->rmap->n / bs + 1; i++) (*ia)[i] = a->i[i] + 1; 132 if (ja) { 133 PetscCall(PetscMalloc1(nz, ja)); 134 for (i = 0; i < nz; i++) (*ja)[i] = a->j[i] + 1; 135 } 136 } 137 } else { 138 *ia = tia; 139 if (ja) *ja = tja; 140 } 141 PetscFunctionReturn(0); 142 } 143 144 static PetscErrorCode MatRestoreRowIJ_SeqSBAIJ(Mat A, PetscInt oshift, PetscBool symmetric, PetscBool blockcompressed, PetscInt *nn, const PetscInt *ia[], const PetscInt *ja[], PetscBool *done) { 145 PetscFunctionBegin; 146 if (!ia) PetscFunctionReturn(0); 147 if ((!blockcompressed && A->rmap->bs > 1) || (symmetric || oshift == 1)) { 148 PetscCall(PetscFree(*ia)); 149 if (ja) PetscCall(PetscFree(*ja)); 150 } 151 PetscFunctionReturn(0); 152 } 153 154 PetscErrorCode MatDestroy_SeqSBAIJ(Mat A) { 155 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data; 156 157 PetscFunctionBegin; 158 #if defined(PETSC_USE_LOG) 159 PetscLogObjectState((PetscObject)A, "Rows=%" PetscInt_FMT ", NZ=%" PetscInt_FMT, A->rmap->N, a->nz); 160 #endif 161 PetscCall(MatSeqXAIJFreeAIJ(A, &a->a, &a->j, &a->i)); 162 if (a->free_diag) PetscCall(PetscFree(a->diag)); 163 PetscCall(ISDestroy(&a->row)); 164 PetscCall(ISDestroy(&a->col)); 165 PetscCall(ISDestroy(&a->icol)); 166 PetscCall(PetscFree(a->idiag)); 167 PetscCall(PetscFree(a->inode.size)); 168 if (a->free_imax_ilen) PetscCall(PetscFree2(a->imax, a->ilen)); 169 PetscCall(PetscFree(a->solve_work)); 170 PetscCall(PetscFree(a->sor_work)); 171 PetscCall(PetscFree(a->solves_work)); 172 PetscCall(PetscFree(a->mult_work)); 173 PetscCall(PetscFree(a->saved_values)); 174 if (a->free_jshort) PetscCall(PetscFree(a->jshort)); 175 PetscCall(PetscFree(a->inew)); 176 PetscCall(MatDestroy(&a->parent)); 177 PetscCall(PetscFree(A->data)); 178 179 PetscCall(PetscObjectChangeTypeName((PetscObject)A, NULL)); 180 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSeqSBAIJGetArray_C", NULL)); 181 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSeqSBAIJRestoreArray_C", NULL)); 182 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatStoreValues_C", NULL)); 183 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatRetrieveValues_C", NULL)); 184 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSeqSBAIJSetColumnIndices_C", NULL)); 185 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqsbaij_seqaij_C", NULL)); 186 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqsbaij_seqbaij_C", NULL)); 187 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSeqSBAIJSetPreallocation_C", NULL)); 188 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSeqSBAIJSetPreallocationCSR_C", NULL)); 189 #if defined(PETSC_HAVE_ELEMENTAL) 190 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqsbaij_elemental_C", NULL)); 191 #endif 192 #if defined(PETSC_HAVE_SCALAPACK) 193 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqsbaij_scalapack_C", NULL)); 194 #endif 195 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatFactorGetSolverType_C", NULL)); 196 PetscFunctionReturn(0); 197 } 198 199 PetscErrorCode MatSetOption_SeqSBAIJ(Mat A, MatOption op, PetscBool flg) { 200 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data; 201 #if defined(PETSC_USE_COMPLEX) 202 PetscInt bs; 203 #endif 204 205 PetscFunctionBegin; 206 #if defined(PETSC_USE_COMPLEX) 207 PetscCall(MatGetBlockSize(A, &bs)); 208 #endif 209 switch (op) { 210 case MAT_ROW_ORIENTED: a->roworiented = flg; break; 211 case MAT_KEEP_NONZERO_PATTERN: a->keepnonzeropattern = flg; break; 212 case MAT_NEW_NONZERO_LOCATIONS: a->nonew = (flg ? 0 : 1); break; 213 case MAT_NEW_NONZERO_LOCATION_ERR: a->nonew = (flg ? -1 : 0); break; 214 case MAT_NEW_NONZERO_ALLOCATION_ERR: a->nonew = (flg ? -2 : 0); break; 215 case MAT_UNUSED_NONZERO_LOCATION_ERR: a->nounused = (flg ? -1 : 0); break; 216 case MAT_FORCE_DIAGONAL_ENTRIES: 217 case MAT_IGNORE_OFF_PROC_ENTRIES: 218 case MAT_USE_HASH_TABLE: 219 case MAT_SORTED_FULL: PetscCall(PetscInfo(A, "Option %s ignored\n", MatOptions[op])); break; 220 case MAT_HERMITIAN: 221 #if defined(PETSC_USE_COMPLEX) 222 if (flg) { /* disable transpose ops */ 223 PetscCheck(bs <= 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for Hermitian with block size greater than 1"); 224 A->ops->multtranspose = NULL; 225 A->ops->multtransposeadd = NULL; 226 A->symmetric = PETSC_BOOL3_FALSE; 227 } 228 #endif 229 break; 230 case MAT_SYMMETRIC: 231 case MAT_SPD: 232 #if defined(PETSC_USE_COMPLEX) 233 if (flg) { /* An hermitian and symmetric matrix has zero imaginary part (restore back transpose ops) */ 234 A->ops->multtranspose = A->ops->mult; 235 A->ops->multtransposeadd = A->ops->multadd; 236 } 237 #endif 238 break; 239 /* These options are handled directly by MatSetOption() */ 240 case MAT_STRUCTURALLY_SYMMETRIC: 241 case MAT_SYMMETRY_ETERNAL: 242 case MAT_STRUCTURAL_SYMMETRY_ETERNAL: 243 case MAT_STRUCTURE_ONLY: 244 case MAT_SPD_ETERNAL: 245 /* These options are handled directly by MatSetOption() */ 246 break; 247 case MAT_IGNORE_LOWER_TRIANGULAR: a->ignore_ltriangular = flg; break; 248 case MAT_ERROR_LOWER_TRIANGULAR: a->ignore_ltriangular = flg; break; 249 case MAT_GETROW_UPPERTRIANGULAR: a->getrow_utriangular = flg; break; 250 case MAT_SUBMAT_SINGLEIS: break; 251 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "unknown option %d", op); 252 } 253 PetscFunctionReturn(0); 254 } 255 256 PetscErrorCode MatGetRow_SeqSBAIJ(Mat A, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v) { 257 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data; 258 259 PetscFunctionBegin; 260 PetscCheck(!A || a->getrow_utriangular, PETSC_COMM_SELF, PETSC_ERR_SUP, "MatGetRow is not supported for SBAIJ matrix format. Getting the upper triangular part of row, run with -mat_getrow_uppertriangular, call MatSetOption(mat,MAT_GETROW_UPPERTRIANGULAR,PETSC_TRUE) or MatGetRowUpperTriangular()"); 261 262 /* Get the upper triangular part of the row */ 263 PetscCall(MatGetRow_SeqBAIJ_private(A, row, nz, idx, v, a->i, a->j, a->a)); 264 PetscFunctionReturn(0); 265 } 266 267 PetscErrorCode MatRestoreRow_SeqSBAIJ(Mat A, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v) { 268 PetscFunctionBegin; 269 if (nz) *nz = 0; 270 if (idx) PetscCall(PetscFree(*idx)); 271 if (v) PetscCall(PetscFree(*v)); 272 PetscFunctionReturn(0); 273 } 274 275 PetscErrorCode MatGetRowUpperTriangular_SeqSBAIJ(Mat A) { 276 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data; 277 278 PetscFunctionBegin; 279 a->getrow_utriangular = PETSC_TRUE; 280 PetscFunctionReturn(0); 281 } 282 283 PetscErrorCode MatRestoreRowUpperTriangular_SeqSBAIJ(Mat A) { 284 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data; 285 286 PetscFunctionBegin; 287 a->getrow_utriangular = PETSC_FALSE; 288 PetscFunctionReturn(0); 289 } 290 291 PetscErrorCode MatTranspose_SeqSBAIJ(Mat A, MatReuse reuse, Mat *B) { 292 PetscFunctionBegin; 293 if (reuse == MAT_REUSE_MATRIX) PetscCall(MatTransposeCheckNonzeroState_Private(A, *B)); 294 if (reuse == MAT_INITIAL_MATRIX) { 295 PetscCall(MatDuplicate(A, MAT_COPY_VALUES, B)); 296 } else if (reuse == MAT_REUSE_MATRIX) { 297 PetscCall(MatCopy(A, *B, SAME_NONZERO_PATTERN)); 298 } 299 PetscFunctionReturn(0); 300 } 301 302 PetscErrorCode MatView_SeqSBAIJ_ASCII(Mat A, PetscViewer viewer) { 303 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data; 304 PetscInt i, j, bs = A->rmap->bs, k, l, bs2 = a->bs2; 305 PetscViewerFormat format; 306 PetscInt *diag; 307 308 PetscFunctionBegin; 309 PetscCall(PetscViewerGetFormat(viewer, &format)); 310 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 311 PetscCall(PetscViewerASCIIPrintf(viewer, " block size is %" PetscInt_FMT "\n", bs)); 312 } else if (format == PETSC_VIEWER_ASCII_MATLAB) { 313 Mat aij; 314 const char *matname; 315 316 if (A->factortype && bs > 1) { 317 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Warning: matrix is factored with bs>1. MatView() with PETSC_VIEWER_ASCII_MATLAB is not supported and ignored!\n")); 318 PetscFunctionReturn(0); 319 } 320 PetscCall(MatConvert(A, MATSEQAIJ, MAT_INITIAL_MATRIX, &aij)); 321 PetscCall(PetscObjectGetName((PetscObject)A, &matname)); 322 PetscCall(PetscObjectSetName((PetscObject)aij, matname)); 323 PetscCall(MatView(aij, viewer)); 324 PetscCall(MatDestroy(&aij)); 325 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 326 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE)); 327 for (i = 0; i < a->mbs; i++) { 328 for (j = 0; j < bs; j++) { 329 PetscCall(PetscViewerASCIIPrintf(viewer, "row %" PetscInt_FMT ":", i * bs + j)); 330 for (k = a->i[i]; k < a->i[i + 1]; k++) { 331 for (l = 0; l < bs; l++) { 332 #if defined(PETSC_USE_COMPLEX) 333 if (PetscImaginaryPart(a->a[bs2 * k + l * bs + j]) > 0.0 && PetscRealPart(a->a[bs2 * k + l * bs + j]) != 0.0) { 334 PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g + %g i) ", bs * a->j[k] + l, (double)PetscRealPart(a->a[bs2 * k + l * bs + j]), (double)PetscImaginaryPart(a->a[bs2 * k + l * bs + j]))); 335 } else if (PetscImaginaryPart(a->a[bs2 * k + l * bs + j]) < 0.0 && PetscRealPart(a->a[bs2 * k + l * bs + j]) != 0.0) { 336 PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g - %g i) ", bs * a->j[k] + l, (double)PetscRealPart(a->a[bs2 * k + l * bs + j]), -(double)PetscImaginaryPart(a->a[bs2 * k + l * bs + j]))); 337 } else if (PetscRealPart(a->a[bs2 * k + l * bs + j]) != 0.0) { 338 PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", bs * a->j[k] + l, (double)PetscRealPart(a->a[bs2 * k + l * bs + j]))); 339 } 340 #else 341 if (a->a[bs2 * k + l * bs + j] != 0.0) PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", bs * a->j[k] + l, (double)a->a[bs2 * k + l * bs + j])); 342 #endif 343 } 344 } 345 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 346 } 347 } 348 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE)); 349 } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) { 350 PetscFunctionReturn(0); 351 } else { 352 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE)); 353 if (A->factortype) { /* for factored matrix */ 354 PetscCheck(bs <= 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "matrix is factored with bs>1. Not implemented yet"); 355 356 diag = a->diag; 357 for (i = 0; i < a->mbs; i++) { /* for row block i */ 358 PetscCall(PetscViewerASCIIPrintf(viewer, "row %" PetscInt_FMT ":", i)); 359 /* diagonal entry */ 360 #if defined(PETSC_USE_COMPLEX) 361 if (PetscImaginaryPart(a->a[diag[i]]) > 0.0) { 362 PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g + %g i) ", a->j[diag[i]], (double)PetscRealPart(1.0 / a->a[diag[i]]), (double)PetscImaginaryPart(1.0 / a->a[diag[i]]))); 363 } else if (PetscImaginaryPart(a->a[diag[i]]) < 0.0) { 364 PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g - %g i) ", a->j[diag[i]], (double)PetscRealPart(1.0 / a->a[diag[i]]), -(double)PetscImaginaryPart(1.0 / a->a[diag[i]]))); 365 } else { 366 PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", a->j[diag[i]], (double)PetscRealPart(1.0 / a->a[diag[i]]))); 367 } 368 #else 369 PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", a->j[diag[i]], (double)(1.0 / a->a[diag[i]]))); 370 #endif 371 /* off-diagonal entries */ 372 for (k = a->i[i]; k < a->i[i + 1] - 1; k++) { 373 #if defined(PETSC_USE_COMPLEX) 374 if (PetscImaginaryPart(a->a[k]) > 0.0) { 375 PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g + %g i) ", bs * a->j[k], (double)PetscRealPart(a->a[k]), (double)PetscImaginaryPart(a->a[k]))); 376 } else if (PetscImaginaryPart(a->a[k]) < 0.0) { 377 PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g - %g i) ", bs * a->j[k], (double)PetscRealPart(a->a[k]), -(double)PetscImaginaryPart(a->a[k]))); 378 } else { 379 PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", bs * a->j[k], (double)PetscRealPart(a->a[k]))); 380 } 381 #else 382 PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", a->j[k], (double)a->a[k])); 383 #endif 384 } 385 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 386 } 387 388 } else { /* for non-factored matrix */ 389 for (i = 0; i < a->mbs; i++) { /* for row block i */ 390 for (j = 0; j < bs; j++) { /* for row bs*i + j */ 391 PetscCall(PetscViewerASCIIPrintf(viewer, "row %" PetscInt_FMT ":", i * bs + j)); 392 for (k = a->i[i]; k < a->i[i + 1]; k++) { /* for column block */ 393 for (l = 0; l < bs; l++) { /* for column */ 394 #if defined(PETSC_USE_COMPLEX) 395 if (PetscImaginaryPart(a->a[bs2 * k + l * bs + j]) > 0.0) { 396 PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g + %g i) ", bs * a->j[k] + l, (double)PetscRealPart(a->a[bs2 * k + l * bs + j]), (double)PetscImaginaryPart(a->a[bs2 * k + l * bs + j]))); 397 } else if (PetscImaginaryPart(a->a[bs2 * k + l * bs + j]) < 0.0) { 398 PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g - %g i) ", bs * a->j[k] + l, (double)PetscRealPart(a->a[bs2 * k + l * bs + j]), -(double)PetscImaginaryPart(a->a[bs2 * k + l * bs + j]))); 399 } else { 400 PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", bs * a->j[k] + l, (double)PetscRealPart(a->a[bs2 * k + l * bs + j]))); 401 } 402 #else 403 PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", bs * a->j[k] + l, (double)a->a[bs2 * k + l * bs + j])); 404 #endif 405 } 406 } 407 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 408 } 409 } 410 } 411 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE)); 412 } 413 PetscCall(PetscViewerFlush(viewer)); 414 PetscFunctionReturn(0); 415 } 416 417 #include <petscdraw.h> 418 static PetscErrorCode MatView_SeqSBAIJ_Draw_Zoom(PetscDraw draw, void *Aa) { 419 Mat A = (Mat)Aa; 420 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data; 421 PetscInt row, i, j, k, l, mbs = a->mbs, color, bs = A->rmap->bs, bs2 = a->bs2; 422 PetscReal xl, yl, xr, yr, x_l, x_r, y_l, y_r; 423 MatScalar *aa; 424 PetscViewer viewer; 425 426 PetscFunctionBegin; 427 PetscCall(PetscObjectQuery((PetscObject)A, "Zoomviewer", (PetscObject *)&viewer)); 428 PetscCall(PetscDrawGetCoordinates(draw, &xl, &yl, &xr, &yr)); 429 430 /* loop over matrix elements drawing boxes */ 431 432 PetscDrawCollectiveBegin(draw); 433 PetscCall(PetscDrawString(draw, .3 * (xl + xr), .3 * (yl + yr), PETSC_DRAW_BLACK, "symmetric")); 434 /* Blue for negative, Cyan for zero and Red for positive */ 435 color = PETSC_DRAW_BLUE; 436 for (i = 0, row = 0; i < mbs; i++, row += bs) { 437 for (j = a->i[i]; j < a->i[i + 1]; j++) { 438 y_l = A->rmap->N - row - 1.0; 439 y_r = y_l + 1.0; 440 x_l = a->j[j] * bs; 441 x_r = x_l + 1.0; 442 aa = a->a + j * bs2; 443 for (k = 0; k < bs; k++) { 444 for (l = 0; l < bs; l++) { 445 if (PetscRealPart(*aa++) >= 0.) continue; 446 PetscCall(PetscDrawRectangle(draw, x_l + k, y_l - l, x_r + k, y_r - l, color, color, color, color)); 447 } 448 } 449 } 450 } 451 color = PETSC_DRAW_CYAN; 452 for (i = 0, row = 0; i < mbs; i++, row += bs) { 453 for (j = a->i[i]; j < a->i[i + 1]; j++) { 454 y_l = A->rmap->N - row - 1.0; 455 y_r = y_l + 1.0; 456 x_l = a->j[j] * bs; 457 x_r = x_l + 1.0; 458 aa = a->a + j * bs2; 459 for (k = 0; k < bs; k++) { 460 for (l = 0; l < bs; l++) { 461 if (PetscRealPart(*aa++) != 0.) continue; 462 PetscCall(PetscDrawRectangle(draw, x_l + k, y_l - l, x_r + k, y_r - l, color, color, color, color)); 463 } 464 } 465 } 466 } 467 color = PETSC_DRAW_RED; 468 for (i = 0, row = 0; i < mbs; i++, row += bs) { 469 for (j = a->i[i]; j < a->i[i + 1]; j++) { 470 y_l = A->rmap->N - row - 1.0; 471 y_r = y_l + 1.0; 472 x_l = a->j[j] * bs; 473 x_r = x_l + 1.0; 474 aa = a->a + j * bs2; 475 for (k = 0; k < bs; k++) { 476 for (l = 0; l < bs; l++) { 477 if (PetscRealPart(*aa++) <= 0.) continue; 478 PetscCall(PetscDrawRectangle(draw, x_l + k, y_l - l, x_r + k, y_r - l, color, color, color, color)); 479 } 480 } 481 } 482 } 483 PetscDrawCollectiveEnd(draw); 484 PetscFunctionReturn(0); 485 } 486 487 static PetscErrorCode MatView_SeqSBAIJ_Draw(Mat A, PetscViewer viewer) { 488 PetscReal xl, yl, xr, yr, w, h; 489 PetscDraw draw; 490 PetscBool isnull; 491 492 PetscFunctionBegin; 493 PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw)); 494 PetscCall(PetscDrawIsNull(draw, &isnull)); 495 if (isnull) PetscFunctionReturn(0); 496 497 xr = A->rmap->N; 498 yr = A->rmap->N; 499 h = yr / 10.0; 500 w = xr / 10.0; 501 xr += w; 502 yr += h; 503 xl = -w; 504 yl = -h; 505 PetscCall(PetscDrawSetCoordinates(draw, xl, yl, xr, yr)); 506 PetscCall(PetscObjectCompose((PetscObject)A, "Zoomviewer", (PetscObject)viewer)); 507 PetscCall(PetscDrawZoom(draw, MatView_SeqSBAIJ_Draw_Zoom, A)); 508 PetscCall(PetscObjectCompose((PetscObject)A, "Zoomviewer", NULL)); 509 PetscCall(PetscDrawSave(draw)); 510 PetscFunctionReturn(0); 511 } 512 513 /* Used for both MPIBAIJ and MPISBAIJ matrices */ 514 #define MatView_SeqSBAIJ_Binary MatView_SeqBAIJ_Binary 515 516 PetscErrorCode MatView_SeqSBAIJ(Mat A, PetscViewer viewer) { 517 PetscBool iascii, isbinary, isdraw; 518 519 PetscFunctionBegin; 520 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 521 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary)); 522 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw)); 523 if (iascii) { 524 PetscCall(MatView_SeqSBAIJ_ASCII(A, viewer)); 525 } else if (isbinary) { 526 PetscCall(MatView_SeqSBAIJ_Binary(A, viewer)); 527 } else if (isdraw) { 528 PetscCall(MatView_SeqSBAIJ_Draw(A, viewer)); 529 } else { 530 Mat B; 531 const char *matname; 532 PetscCall(MatConvert(A, MATSEQAIJ, MAT_INITIAL_MATRIX, &B)); 533 PetscCall(PetscObjectGetName((PetscObject)A, &matname)); 534 PetscCall(PetscObjectSetName((PetscObject)B, matname)); 535 PetscCall(MatView(B, viewer)); 536 PetscCall(MatDestroy(&B)); 537 } 538 PetscFunctionReturn(0); 539 } 540 541 PetscErrorCode MatGetValues_SeqSBAIJ(Mat A, PetscInt m, const PetscInt im[], PetscInt n, const PetscInt in[], PetscScalar v[]) { 542 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data; 543 PetscInt *rp, k, low, high, t, row, nrow, i, col, l, *aj = a->j; 544 PetscInt *ai = a->i, *ailen = a->ilen; 545 PetscInt brow, bcol, ridx, cidx, bs = A->rmap->bs, bs2 = a->bs2; 546 MatScalar *ap, *aa = a->a; 547 548 PetscFunctionBegin; 549 for (k = 0; k < m; k++) { /* loop over rows */ 550 row = im[k]; 551 brow = row / bs; 552 if (row < 0) { 553 v += n; 554 continue; 555 } /* negative row */ 556 PetscCheck(row < A->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, row, A->rmap->N - 1); 557 rp = aj + ai[brow]; 558 ap = aa + bs2 * ai[brow]; 559 nrow = ailen[brow]; 560 for (l = 0; l < n; l++) { /* loop over columns */ 561 if (in[l] < 0) { 562 v++; 563 continue; 564 } /* negative column */ 565 PetscCheck(in[l] < A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT, in[l], A->cmap->n - 1); 566 col = in[l]; 567 bcol = col / bs; 568 cidx = col % bs; 569 ridx = row % bs; 570 high = nrow; 571 low = 0; /* assume unsorted */ 572 while (high - low > 5) { 573 t = (low + high) / 2; 574 if (rp[t] > bcol) high = t; 575 else low = t; 576 } 577 for (i = low; i < high; i++) { 578 if (rp[i] > bcol) break; 579 if (rp[i] == bcol) { 580 *v++ = ap[bs2 * i + bs * cidx + ridx]; 581 goto finished; 582 } 583 } 584 *v++ = 0.0; 585 finished:; 586 } 587 } 588 PetscFunctionReturn(0); 589 } 590 591 PetscErrorCode MatPermute_SeqSBAIJ(Mat A, IS rowp, IS colp, Mat *B) { 592 Mat C; 593 594 PetscFunctionBegin; 595 PetscCall(MatConvert(A, MATSEQBAIJ, MAT_INITIAL_MATRIX, &C)); 596 PetscCall(MatPermute(C, rowp, colp, B)); 597 PetscCall(MatDestroy(&C)); 598 if (rowp == colp) PetscCall(MatConvert(*B, MATSEQSBAIJ, MAT_INPLACE_MATRIX, B)); 599 PetscFunctionReturn(0); 600 } 601 602 PetscErrorCode MatSetValuesBlocked_SeqSBAIJ(Mat A, PetscInt m, const PetscInt im[], PetscInt n, const PetscInt in[], const PetscScalar v[], InsertMode is) { 603 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data; 604 PetscInt *rp, k, low, high, t, ii, jj, row, nrow, i, col, l, rmax, N, lastcol = -1; 605 PetscInt *imax = a->imax, *ai = a->i, *ailen = a->ilen; 606 PetscInt *aj = a->j, nonew = a->nonew, bs2 = a->bs2, bs = A->rmap->bs, stepval; 607 PetscBool roworiented = a->roworiented; 608 const PetscScalar *value = v; 609 MatScalar *ap, *aa = a->a, *bap; 610 611 PetscFunctionBegin; 612 if (roworiented) stepval = (n - 1) * bs; 613 else stepval = (m - 1) * bs; 614 615 for (k = 0; k < m; k++) { /* loop over added rows */ 616 row = im[k]; 617 if (row < 0) continue; 618 PetscCheck(row < a->mbs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Block index row too large %" PetscInt_FMT " max %" PetscInt_FMT, row, a->mbs - 1); 619 rp = aj + ai[row]; 620 ap = aa + bs2 * ai[row]; 621 rmax = imax[row]; 622 nrow = ailen[row]; 623 low = 0; 624 high = nrow; 625 for (l = 0; l < n; l++) { /* loop over added columns */ 626 if (in[l] < 0) continue; 627 col = in[l]; 628 PetscCheck(col < a->nbs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Block index column too large %" PetscInt_FMT " max %" PetscInt_FMT, col, a->nbs - 1); 629 if (col < row) { 630 if (a->ignore_ltriangular) continue; /* ignore lower triangular block */ 631 else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Lower triangular value cannot be set for sbaij format. Ignoring these values, run with -mat_ignore_lower_triangular or call MatSetOption(mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE)"); 632 } 633 if (roworiented) value = v + k * (stepval + bs) * bs + l * bs; 634 else value = v + l * (stepval + bs) * bs + k * bs; 635 636 if (col <= lastcol) low = 0; 637 else high = nrow; 638 639 lastcol = col; 640 while (high - low > 7) { 641 t = (low + high) / 2; 642 if (rp[t] > col) high = t; 643 else low = t; 644 } 645 for (i = low; i < high; i++) { 646 if (rp[i] > col) break; 647 if (rp[i] == col) { 648 bap = ap + bs2 * i; 649 if (roworiented) { 650 if (is == ADD_VALUES) { 651 for (ii = 0; ii < bs; ii++, value += stepval) { 652 for (jj = ii; jj < bs2; jj += bs) bap[jj] += *value++; 653 } 654 } else { 655 for (ii = 0; ii < bs; ii++, value += stepval) { 656 for (jj = ii; jj < bs2; jj += bs) bap[jj] = *value++; 657 } 658 } 659 } else { 660 if (is == ADD_VALUES) { 661 for (ii = 0; ii < bs; ii++, value += stepval) { 662 for (jj = 0; jj < bs; jj++) *bap++ += *value++; 663 } 664 } else { 665 for (ii = 0; ii < bs; ii++, value += stepval) { 666 for (jj = 0; jj < bs; jj++) *bap++ = *value++; 667 } 668 } 669 } 670 goto noinsert2; 671 } 672 } 673 if (nonew == 1) goto noinsert2; 674 PetscCheck(nonew != -1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Inserting a new block index nonzero block (%" PetscInt_FMT ", %" PetscInt_FMT ") in the matrix", row, col); 675 MatSeqXAIJReallocateAIJ(A, a->mbs, bs2, nrow, row, col, rmax, aa, ai, aj, rp, ap, imax, nonew, MatScalar); 676 N = nrow++ - 1; 677 high++; 678 /* shift up all the later entries in this row */ 679 PetscCall(PetscArraymove(rp + i + 1, rp + i, N - i + 1)); 680 PetscCall(PetscArraymove(ap + bs2 * (i + 1), ap + bs2 * i, bs2 * (N - i + 1))); 681 PetscCall(PetscArrayzero(ap + bs2 * i, bs2)); 682 rp[i] = col; 683 bap = ap + bs2 * i; 684 if (roworiented) { 685 for (ii = 0; ii < bs; ii++, value += stepval) { 686 for (jj = ii; jj < bs2; jj += bs) bap[jj] = *value++; 687 } 688 } else { 689 for (ii = 0; ii < bs; ii++, value += stepval) { 690 for (jj = 0; jj < bs; jj++) *bap++ = *value++; 691 } 692 } 693 noinsert2:; 694 low = i; 695 } 696 ailen[row] = nrow; 697 } 698 PetscFunctionReturn(0); 699 } 700 701 /* 702 This is not yet used 703 */ 704 PetscErrorCode MatAssemblyEnd_SeqSBAIJ_SeqAIJ_Inode(Mat A) { 705 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data; 706 const PetscInt *ai = a->i, *aj = a->j, *cols; 707 PetscInt i = 0, j, blk_size, m = A->rmap->n, node_count = 0, nzx, nzy, *ns, row, nz, cnt, cnt2, *counts; 708 PetscBool flag; 709 710 PetscFunctionBegin; 711 PetscCall(PetscMalloc1(m, &ns)); 712 while (i < m) { 713 nzx = ai[i + 1] - ai[i]; /* Number of nonzeros */ 714 /* Limits the number of elements in a node to 'a->inode.limit' */ 715 for (j = i + 1, blk_size = 1; j < m && blk_size < a->inode.limit; ++j, ++blk_size) { 716 nzy = ai[j + 1] - ai[j]; 717 if (nzy != (nzx - j + i)) break; 718 PetscCall(PetscArraycmp(aj + ai[i] + j - i, aj + ai[j], nzy, &flag)); 719 if (!flag) break; 720 } 721 ns[node_count++] = blk_size; 722 723 i = j; 724 } 725 if (!a->inode.size && m && node_count > .9 * m) { 726 PetscCall(PetscFree(ns)); 727 PetscCall(PetscInfo(A, "Found %" PetscInt_FMT " nodes out of %" PetscInt_FMT " rows. Not using Inode routines\n", node_count, m)); 728 } else { 729 a->inode.node_count = node_count; 730 731 PetscCall(PetscMalloc1(node_count, &a->inode.size)); 732 PetscCall(PetscArraycpy(a->inode.size, ns, node_count)); 733 PetscCall(PetscFree(ns)); 734 PetscCall(PetscInfo(A, "Found %" PetscInt_FMT " nodes of %" PetscInt_FMT ". Limit used: %" PetscInt_FMT ". Using Inode routines\n", node_count, m, a->inode.limit)); 735 736 /* count collections of adjacent columns in each inode */ 737 row = 0; 738 cnt = 0; 739 for (i = 0; i < node_count; i++) { 740 cols = aj + ai[row] + a->inode.size[i]; 741 nz = ai[row + 1] - ai[row] - a->inode.size[i]; 742 for (j = 1; j < nz; j++) { 743 if (cols[j] != cols[j - 1] + 1) cnt++; 744 } 745 cnt++; 746 row += a->inode.size[i]; 747 } 748 PetscCall(PetscMalloc1(2 * cnt, &counts)); 749 cnt = 0; 750 row = 0; 751 for (i = 0; i < node_count; i++) { 752 cols = aj + ai[row] + a->inode.size[i]; 753 counts[2 * cnt] = cols[0]; 754 nz = ai[row + 1] - ai[row] - a->inode.size[i]; 755 cnt2 = 1; 756 for (j = 1; j < nz; j++) { 757 if (cols[j] != cols[j - 1] + 1) { 758 counts[2 * (cnt++) + 1] = cnt2; 759 counts[2 * cnt] = cols[j]; 760 cnt2 = 1; 761 } else cnt2++; 762 } 763 counts[2 * (cnt++) + 1] = cnt2; 764 row += a->inode.size[i]; 765 } 766 PetscCall(PetscIntView(2 * cnt, counts, NULL)); 767 } 768 PetscFunctionReturn(0); 769 } 770 771 PetscErrorCode MatAssemblyEnd_SeqSBAIJ(Mat A, MatAssemblyType mode) { 772 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data; 773 PetscInt fshift = 0, i, *ai = a->i, *aj = a->j, *imax = a->imax; 774 PetscInt m = A->rmap->N, *ip, N, *ailen = a->ilen; 775 PetscInt mbs = a->mbs, bs2 = a->bs2, rmax = 0; 776 MatScalar *aa = a->a, *ap; 777 778 PetscFunctionBegin; 779 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 780 781 if (m) rmax = ailen[0]; 782 for (i = 1; i < mbs; i++) { 783 /* move each row back by the amount of empty slots (fshift) before it*/ 784 fshift += imax[i - 1] - ailen[i - 1]; 785 rmax = PetscMax(rmax, ailen[i]); 786 if (fshift) { 787 ip = aj + ai[i]; 788 ap = aa + bs2 * ai[i]; 789 N = ailen[i]; 790 PetscCall(PetscArraymove(ip - fshift, ip, N)); 791 PetscCall(PetscArraymove(ap - bs2 * fshift, ap, bs2 * N)); 792 } 793 ai[i] = ai[i - 1] + ailen[i - 1]; 794 } 795 if (mbs) { 796 fshift += imax[mbs - 1] - ailen[mbs - 1]; 797 ai[mbs] = ai[mbs - 1] + ailen[mbs - 1]; 798 } 799 /* reset ilen and imax for each row */ 800 for (i = 0; i < mbs; i++) ailen[i] = imax[i] = ai[i + 1] - ai[i]; 801 a->nz = ai[mbs]; 802 803 /* diagonals may have moved, reset it */ 804 if (a->diag) PetscCall(PetscArraycpy(a->diag, ai, mbs)); 805 PetscCheck(!fshift || a->nounused != -1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unused space detected in matrix: %" PetscInt_FMT " X %" PetscInt_FMT " block size %" PetscInt_FMT ", %" PetscInt_FMT " unneeded", m, A->cmap->n, A->rmap->bs, fshift * bs2); 806 807 PetscCall(PetscInfo(A, "Matrix size: %" PetscInt_FMT " X %" PetscInt_FMT ", block size %" PetscInt_FMT "; storage space: %" PetscInt_FMT " unneeded, %" PetscInt_FMT " used\n", m, A->rmap->N, A->rmap->bs, fshift * bs2, a->nz * bs2)); 808 PetscCall(PetscInfo(A, "Number of mallocs during MatSetValues is %" PetscInt_FMT "\n", a->reallocs)); 809 PetscCall(PetscInfo(A, "Most nonzeros blocks in any row is %" PetscInt_FMT "\n", rmax)); 810 811 A->info.mallocs += a->reallocs; 812 a->reallocs = 0; 813 A->info.nz_unneeded = (PetscReal)fshift * bs2; 814 a->idiagvalid = PETSC_FALSE; 815 a->rmax = rmax; 816 817 if (A->cmap->n < 65536 && A->cmap->bs == 1) { 818 if (a->jshort && a->free_jshort) { 819 /* when matrix data structure is changed, previous jshort must be replaced */ 820 PetscCall(PetscFree(a->jshort)); 821 } 822 PetscCall(PetscMalloc1(a->i[A->rmap->n], &a->jshort)); 823 for (i = 0; i < a->i[A->rmap->n]; i++) a->jshort[i] = a->j[i]; 824 A->ops->mult = MatMult_SeqSBAIJ_1_ushort; 825 A->ops->sor = MatSOR_SeqSBAIJ_ushort; 826 a->free_jshort = PETSC_TRUE; 827 } 828 PetscFunctionReturn(0); 829 } 830 831 /* 832 This function returns an array of flags which indicate the locations of contiguous 833 blocks that should be zeroed. for eg: if bs = 3 and is = [0,1,2,3,5,6,7,8,9] 834 then the resulting sizes = [3,1,1,3,1] corresponding to sets [(0,1,2),(3),(5),(6,7,8),(9)] 835 Assume: sizes should be long enough to hold all the values. 836 */ 837 PetscErrorCode MatZeroRows_SeqSBAIJ_Check_Blocks(PetscInt idx[], PetscInt n, PetscInt bs, PetscInt sizes[], PetscInt *bs_max) { 838 PetscInt i, j, k, row; 839 PetscBool flg; 840 841 PetscFunctionBegin; 842 for (i = 0, j = 0; i < n; j++) { 843 row = idx[i]; 844 if (row % bs != 0) { /* Not the beginning of a block */ 845 sizes[j] = 1; 846 i++; 847 } else if (i + bs > n) { /* Beginning of a block, but complete block doesn't exist (at idx end) */ 848 sizes[j] = 1; /* Also makes sure at least 'bs' values exist for next else */ 849 i++; 850 } else { /* Beginning of the block, so check if the complete block exists */ 851 flg = PETSC_TRUE; 852 for (k = 1; k < bs; k++) { 853 if (row + k != idx[i + k]) { /* break in the block */ 854 flg = PETSC_FALSE; 855 break; 856 } 857 } 858 if (flg) { /* No break in the bs */ 859 sizes[j] = bs; 860 i += bs; 861 } else { 862 sizes[j] = 1; 863 i++; 864 } 865 } 866 } 867 *bs_max = j; 868 PetscFunctionReturn(0); 869 } 870 871 /* Only add/insert a(i,j) with i<=j (blocks). 872 Any a(i,j) with i>j input by user is ingored. 873 */ 874 875 PetscErrorCode MatSetValues_SeqSBAIJ(Mat A, PetscInt m, const PetscInt im[], PetscInt n, const PetscInt in[], const PetscScalar v[], InsertMode is) { 876 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data; 877 PetscInt *rp, k, low, high, t, ii, row, nrow, i, col, l, rmax, N, lastcol = -1; 878 PetscInt *imax = a->imax, *ai = a->i, *ailen = a->ilen, roworiented = a->roworiented; 879 PetscInt *aj = a->j, nonew = a->nonew, bs = A->rmap->bs, brow, bcol; 880 PetscInt ridx, cidx, bs2 = a->bs2; 881 MatScalar *ap, value, *aa = a->a, *bap; 882 883 PetscFunctionBegin; 884 for (k = 0; k < m; k++) { /* loop over added rows */ 885 row = im[k]; /* row number */ 886 brow = row / bs; /* block row number */ 887 if (row < 0) continue; 888 PetscCheck(row < A->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, row, A->rmap->N - 1); 889 rp = aj + ai[brow]; /*ptr to beginning of column value of the row block*/ 890 ap = aa + bs2 * ai[brow]; /*ptr to beginning of element value of the row block*/ 891 rmax = imax[brow]; /* maximum space allocated for this row */ 892 nrow = ailen[brow]; /* actual length of this row */ 893 low = 0; 894 high = nrow; 895 for (l = 0; l < n; l++) { /* loop over added columns */ 896 if (in[l] < 0) continue; 897 PetscCheck(in[l] < A->cmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT, in[l], A->cmap->N - 1); 898 col = in[l]; 899 bcol = col / bs; /* block col number */ 900 901 if (brow > bcol) { 902 if (a->ignore_ltriangular) continue; /* ignore lower triangular values */ 903 else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Lower triangular value cannot be set for sbaij format. Ignoring these values, run with -mat_ignore_lower_triangular or call MatSetOption(mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE)"); 904 } 905 906 ridx = row % bs; 907 cidx = col % bs; /*row and col index inside the block */ 908 if ((brow == bcol && ridx <= cidx) || (brow < bcol)) { 909 /* element value a(k,l) */ 910 if (roworiented) value = v[l + k * n]; 911 else value = v[k + l * m]; 912 913 /* move pointer bap to a(k,l) quickly and add/insert value */ 914 if (col <= lastcol) low = 0; 915 else high = nrow; 916 917 lastcol = col; 918 while (high - low > 7) { 919 t = (low + high) / 2; 920 if (rp[t] > bcol) high = t; 921 else low = t; 922 } 923 for (i = low; i < high; i++) { 924 if (rp[i] > bcol) break; 925 if (rp[i] == bcol) { 926 bap = ap + bs2 * i + bs * cidx + ridx; 927 if (is == ADD_VALUES) *bap += value; 928 else *bap = value; 929 /* for diag block, add/insert its symmetric element a(cidx,ridx) */ 930 if (brow == bcol && ridx < cidx) { 931 bap = ap + bs2 * i + bs * ridx + cidx; 932 if (is == ADD_VALUES) *bap += value; 933 else *bap = value; 934 } 935 goto noinsert1; 936 } 937 } 938 939 if (nonew == 1) goto noinsert1; 940 PetscCheck(nonew != -1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Inserting a new nonzero (%" PetscInt_FMT ", %" PetscInt_FMT ") in the matrix", row, col); 941 MatSeqXAIJReallocateAIJ(A, a->mbs, bs2, nrow, brow, bcol, rmax, aa, ai, aj, rp, ap, imax, nonew, MatScalar); 942 943 N = nrow++ - 1; 944 high++; 945 /* shift up all the later entries in this row */ 946 PetscCall(PetscArraymove(rp + i + 1, rp + i, N - i + 1)); 947 PetscCall(PetscArraymove(ap + bs2 * (i + 1), ap + bs2 * i, bs2 * (N - i + 1))); 948 PetscCall(PetscArrayzero(ap + bs2 * i, bs2)); 949 rp[i] = bcol; 950 ap[bs2 * i + bs * cidx + ridx] = value; 951 /* for diag block, add/insert its symmetric element a(cidx,ridx) */ 952 if (brow == bcol && ridx < cidx) ap[bs2 * i + bs * ridx + cidx] = value; 953 A->nonzerostate++; 954 noinsert1:; 955 low = i; 956 } 957 } /* end of loop over added columns */ 958 ailen[brow] = nrow; 959 } /* end of loop over added rows */ 960 PetscFunctionReturn(0); 961 } 962 963 PetscErrorCode MatICCFactor_SeqSBAIJ(Mat inA, IS row, const MatFactorInfo *info) { 964 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)inA->data; 965 Mat outA; 966 PetscBool row_identity; 967 968 PetscFunctionBegin; 969 PetscCheck(info->levels == 0, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only levels=0 is supported for in-place icc"); 970 PetscCall(ISIdentity(row, &row_identity)); 971 PetscCheck(row_identity, PETSC_COMM_SELF, PETSC_ERR_SUP, "Matrix reordering is not supported"); 972 PetscCheck(inA->rmap->bs == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Matrix block size %" PetscInt_FMT " is not supported", inA->rmap->bs); /* Need to replace MatCholeskyFactorSymbolic_SeqSBAIJ_MSR()! */ 973 974 outA = inA; 975 inA->factortype = MAT_FACTOR_ICC; 976 PetscCall(PetscFree(inA->solvertype)); 977 PetscCall(PetscStrallocpy(MATSOLVERPETSC, &inA->solvertype)); 978 979 PetscCall(MatMarkDiagonal_SeqSBAIJ(inA)); 980 PetscCall(MatSeqSBAIJSetNumericFactorization_inplace(inA, row_identity)); 981 982 PetscCall(PetscObjectReference((PetscObject)row)); 983 PetscCall(ISDestroy(&a->row)); 984 a->row = row; 985 PetscCall(PetscObjectReference((PetscObject)row)); 986 PetscCall(ISDestroy(&a->col)); 987 a->col = row; 988 989 /* Create the invert permutation so that it can be used in MatCholeskyFactorNumeric() */ 990 if (a->icol) PetscCall(ISInvertPermutation(row, PETSC_DECIDE, &a->icol)); 991 992 if (!a->solve_work) { PetscCall(PetscMalloc1(inA->rmap->N + inA->rmap->bs, &a->solve_work)); } 993 994 PetscCall(MatCholeskyFactorNumeric(outA, inA, info)); 995 PetscFunctionReturn(0); 996 } 997 998 PetscErrorCode MatSeqSBAIJSetColumnIndices_SeqSBAIJ(Mat mat, PetscInt *indices) { 999 Mat_SeqSBAIJ *baij = (Mat_SeqSBAIJ *)mat->data; 1000 PetscInt i, nz, n; 1001 1002 PetscFunctionBegin; 1003 nz = baij->maxnz; 1004 n = mat->cmap->n; 1005 for (i = 0; i < nz; i++) baij->j[i] = indices[i]; 1006 1007 baij->nz = nz; 1008 for (i = 0; i < n; i++) baij->ilen[i] = baij->imax[i]; 1009 1010 PetscCall(MatSetOption(mat, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE)); 1011 PetscFunctionReturn(0); 1012 } 1013 1014 /*@ 1015 MatSeqSBAIJSetColumnIndices - Set the column indices for all the rows 1016 in a `MATSEQSBAIJ` matrix. 1017 1018 Input Parameters: 1019 + mat - the `MATSEQSBAIJ` matrix 1020 - indices - the column indices 1021 1022 Level: advanced 1023 1024 Notes: 1025 This can be called if you have precomputed the nonzero structure of the 1026 matrix and want to provide it to the matrix object to improve the performance 1027 of the `MatSetValues()` operation. 1028 1029 You MUST have set the correct numbers of nonzeros per row in the call to 1030 `MatCreateSeqSBAIJ()`, and the columns indices MUST be sorted. 1031 1032 MUST be called before any calls to MatSetValues() 1033 1034 .seealso: `MATSEQSBAIJ`, `MatCreateSeqSBAIJ` 1035 @*/ 1036 PetscErrorCode MatSeqSBAIJSetColumnIndices(Mat mat, PetscInt *indices) { 1037 PetscFunctionBegin; 1038 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 1039 PetscValidIntPointer(indices, 2); 1040 PetscUseMethod(mat, "MatSeqSBAIJSetColumnIndices_C", (Mat, PetscInt *), (mat, indices)); 1041 PetscFunctionReturn(0); 1042 } 1043 1044 PetscErrorCode MatCopy_SeqSBAIJ(Mat A, Mat B, MatStructure str) { 1045 PetscBool isbaij; 1046 1047 PetscFunctionBegin; 1048 PetscCall(PetscObjectTypeCompareAny((PetscObject)B, &isbaij, MATSEQSBAIJ, MATMPISBAIJ, "")); 1049 PetscCheck(isbaij, PetscObjectComm((PetscObject)B), PETSC_ERR_SUP, "Not for matrix type %s", ((PetscObject)B)->type_name); 1050 /* If the two matrices have the same copy implementation and nonzero pattern, use fast copy. */ 1051 if (str == SAME_NONZERO_PATTERN && A->ops->copy == B->ops->copy) { 1052 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data; 1053 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ *)B->data; 1054 1055 PetscCheck(a->i[a->mbs] == b->i[b->mbs], PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of nonzeros in two matrices are different"); 1056 PetscCheck(a->mbs == b->mbs, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of rows in two matrices are different"); 1057 PetscCheck(a->bs2 == b->bs2, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Different block size"); 1058 PetscCall(PetscArraycpy(b->a, a->a, a->bs2 * a->i[a->mbs])); 1059 PetscCall(PetscObjectStateIncrease((PetscObject)B)); 1060 } else { 1061 PetscCall(MatGetRowUpperTriangular(A)); 1062 PetscCall(MatCopy_Basic(A, B, str)); 1063 PetscCall(MatRestoreRowUpperTriangular(A)); 1064 } 1065 PetscFunctionReturn(0); 1066 } 1067 1068 PetscErrorCode MatSetUp_SeqSBAIJ(Mat A) { 1069 PetscFunctionBegin; 1070 PetscCall(MatSeqSBAIJSetPreallocation(A, A->rmap->bs, PETSC_DEFAULT, NULL)); 1071 PetscFunctionReturn(0); 1072 } 1073 1074 static PetscErrorCode MatSeqSBAIJGetArray_SeqSBAIJ(Mat A, PetscScalar *array[]) { 1075 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data; 1076 1077 PetscFunctionBegin; 1078 *array = a->a; 1079 PetscFunctionReturn(0); 1080 } 1081 1082 static PetscErrorCode MatSeqSBAIJRestoreArray_SeqSBAIJ(Mat A, PetscScalar *array[]) { 1083 PetscFunctionBegin; 1084 *array = NULL; 1085 PetscFunctionReturn(0); 1086 } 1087 1088 PetscErrorCode MatAXPYGetPreallocation_SeqSBAIJ(Mat Y, Mat X, PetscInt *nnz) { 1089 PetscInt bs = Y->rmap->bs, mbs = Y->rmap->N / bs; 1090 Mat_SeqSBAIJ *x = (Mat_SeqSBAIJ *)X->data; 1091 Mat_SeqSBAIJ *y = (Mat_SeqSBAIJ *)Y->data; 1092 1093 PetscFunctionBegin; 1094 /* Set the number of nonzeros in the new matrix */ 1095 PetscCall(MatAXPYGetPreallocation_SeqX_private(mbs, x->i, x->j, y->i, y->j, nnz)); 1096 PetscFunctionReturn(0); 1097 } 1098 1099 PetscErrorCode MatAXPY_SeqSBAIJ(Mat Y, PetscScalar a, Mat X, MatStructure str) { 1100 Mat_SeqSBAIJ *x = (Mat_SeqSBAIJ *)X->data, *y = (Mat_SeqSBAIJ *)Y->data; 1101 PetscInt bs = Y->rmap->bs, bs2 = bs * bs; 1102 PetscBLASInt one = 1; 1103 1104 PetscFunctionBegin; 1105 if (str == UNKNOWN_NONZERO_PATTERN || (PetscDefined(USE_DEBUG) && str == SAME_NONZERO_PATTERN)) { 1106 PetscBool e = x->nz == y->nz && x->mbs == y->mbs ? PETSC_TRUE : PETSC_FALSE; 1107 if (e) { 1108 PetscCall(PetscArraycmp(x->i, y->i, x->mbs + 1, &e)); 1109 if (e) { 1110 PetscCall(PetscArraycmp(x->j, y->j, x->i[x->mbs], &e)); 1111 if (e) str = SAME_NONZERO_PATTERN; 1112 } 1113 } 1114 if (!e) PetscCheck(str != SAME_NONZERO_PATTERN, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "MatStructure is not SAME_NONZERO_PATTERN"); 1115 } 1116 if (str == SAME_NONZERO_PATTERN) { 1117 PetscScalar alpha = a; 1118 PetscBLASInt bnz; 1119 PetscCall(PetscBLASIntCast(x->nz * bs2, &bnz)); 1120 PetscCallBLAS("BLASaxpy", BLASaxpy_(&bnz, &alpha, x->a, &one, y->a, &one)); 1121 PetscCall(PetscObjectStateIncrease((PetscObject)Y)); 1122 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 1123 PetscCall(MatSetOption(X, MAT_GETROW_UPPERTRIANGULAR, PETSC_TRUE)); 1124 PetscCall(MatAXPY_Basic(Y, a, X, str)); 1125 PetscCall(MatSetOption(X, MAT_GETROW_UPPERTRIANGULAR, PETSC_FALSE)); 1126 } else { 1127 Mat B; 1128 PetscInt *nnz; 1129 PetscCheck(bs == X->rmap->bs, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Matrices must have same block size"); 1130 PetscCall(MatGetRowUpperTriangular(X)); 1131 PetscCall(MatGetRowUpperTriangular(Y)); 1132 PetscCall(PetscMalloc1(Y->rmap->N, &nnz)); 1133 PetscCall(MatCreate(PetscObjectComm((PetscObject)Y), &B)); 1134 PetscCall(PetscObjectSetName((PetscObject)B, ((PetscObject)Y)->name)); 1135 PetscCall(MatSetSizes(B, Y->rmap->n, Y->cmap->n, Y->rmap->N, Y->cmap->N)); 1136 PetscCall(MatSetBlockSizesFromMats(B, Y, Y)); 1137 PetscCall(MatSetType(B, ((PetscObject)Y)->type_name)); 1138 PetscCall(MatAXPYGetPreallocation_SeqSBAIJ(Y, X, nnz)); 1139 PetscCall(MatSeqSBAIJSetPreallocation(B, bs, 0, nnz)); 1140 1141 PetscCall(MatAXPY_BasicWithPreallocation(B, Y, a, X, str)); 1142 1143 PetscCall(MatHeaderMerge(Y, &B)); 1144 PetscCall(PetscFree(nnz)); 1145 PetscCall(MatRestoreRowUpperTriangular(X)); 1146 PetscCall(MatRestoreRowUpperTriangular(Y)); 1147 } 1148 PetscFunctionReturn(0); 1149 } 1150 1151 PetscErrorCode MatIsSymmetric_SeqSBAIJ(Mat A, PetscReal tol, PetscBool *flg) { 1152 PetscFunctionBegin; 1153 *flg = PETSC_TRUE; 1154 PetscFunctionReturn(0); 1155 } 1156 1157 PetscErrorCode MatIsStructurallySymmetric_SeqSBAIJ(Mat A, PetscBool *flg) { 1158 PetscFunctionBegin; 1159 *flg = PETSC_TRUE; 1160 PetscFunctionReturn(0); 1161 } 1162 1163 PetscErrorCode MatIsHermitian_SeqSBAIJ(Mat A, PetscReal tol, PetscBool *flg) { 1164 PetscFunctionBegin; 1165 *flg = PETSC_FALSE; 1166 PetscFunctionReturn(0); 1167 } 1168 1169 PetscErrorCode MatConjugate_SeqSBAIJ(Mat A) { 1170 #if defined(PETSC_USE_COMPLEX) 1171 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data; 1172 PetscInt i, nz = a->bs2 * a->i[a->mbs]; 1173 MatScalar *aa = a->a; 1174 1175 PetscFunctionBegin; 1176 for (i = 0; i < nz; i++) aa[i] = PetscConj(aa[i]); 1177 #else 1178 PetscFunctionBegin; 1179 #endif 1180 PetscFunctionReturn(0); 1181 } 1182 1183 PetscErrorCode MatRealPart_SeqSBAIJ(Mat A) { 1184 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data; 1185 PetscInt i, nz = a->bs2 * a->i[a->mbs]; 1186 MatScalar *aa = a->a; 1187 1188 PetscFunctionBegin; 1189 for (i = 0; i < nz; i++) aa[i] = PetscRealPart(aa[i]); 1190 PetscFunctionReturn(0); 1191 } 1192 1193 PetscErrorCode MatImaginaryPart_SeqSBAIJ(Mat A) { 1194 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data; 1195 PetscInt i, nz = a->bs2 * a->i[a->mbs]; 1196 MatScalar *aa = a->a; 1197 1198 PetscFunctionBegin; 1199 for (i = 0; i < nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1200 PetscFunctionReturn(0); 1201 } 1202 1203 PetscErrorCode MatZeroRowsColumns_SeqSBAIJ(Mat A, PetscInt is_n, const PetscInt is_idx[], PetscScalar diag, Vec x, Vec b) { 1204 Mat_SeqSBAIJ *baij = (Mat_SeqSBAIJ *)A->data; 1205 PetscInt i, j, k, count; 1206 PetscInt bs = A->rmap->bs, bs2 = baij->bs2, row, col; 1207 PetscScalar zero = 0.0; 1208 MatScalar *aa; 1209 const PetscScalar *xx; 1210 PetscScalar *bb; 1211 PetscBool *zeroed, vecs = PETSC_FALSE; 1212 1213 PetscFunctionBegin; 1214 /* fix right hand side if needed */ 1215 if (x && b) { 1216 PetscCall(VecGetArrayRead(x, &xx)); 1217 PetscCall(VecGetArray(b, &bb)); 1218 vecs = PETSC_TRUE; 1219 } 1220 1221 /* zero the columns */ 1222 PetscCall(PetscCalloc1(A->rmap->n, &zeroed)); 1223 for (i = 0; i < is_n; i++) { 1224 PetscCheck(is_idx[i] >= 0 && is_idx[i] < A->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "row %" PetscInt_FMT " out of range", is_idx[i]); 1225 zeroed[is_idx[i]] = PETSC_TRUE; 1226 } 1227 if (vecs) { 1228 for (i = 0; i < A->rmap->N; i++) { 1229 row = i / bs; 1230 for (j = baij->i[row]; j < baij->i[row + 1]; j++) { 1231 for (k = 0; k < bs; k++) { 1232 col = bs * baij->j[j] + k; 1233 if (col <= i) continue; 1234 aa = ((MatScalar *)(baij->a)) + j * bs2 + (i % bs) + bs * k; 1235 if (!zeroed[i] && zeroed[col]) bb[i] -= aa[0] * xx[col]; 1236 if (zeroed[i] && !zeroed[col]) bb[col] -= aa[0] * xx[i]; 1237 } 1238 } 1239 } 1240 for (i = 0; i < is_n; i++) bb[is_idx[i]] = diag * xx[is_idx[i]]; 1241 } 1242 1243 for (i = 0; i < A->rmap->N; i++) { 1244 if (!zeroed[i]) { 1245 row = i / bs; 1246 for (j = baij->i[row]; j < baij->i[row + 1]; j++) { 1247 for (k = 0; k < bs; k++) { 1248 col = bs * baij->j[j] + k; 1249 if (zeroed[col]) { 1250 aa = ((MatScalar *)(baij->a)) + j * bs2 + (i % bs) + bs * k; 1251 aa[0] = 0.0; 1252 } 1253 } 1254 } 1255 } 1256 } 1257 PetscCall(PetscFree(zeroed)); 1258 if (vecs) { 1259 PetscCall(VecRestoreArrayRead(x, &xx)); 1260 PetscCall(VecRestoreArray(b, &bb)); 1261 } 1262 1263 /* zero the rows */ 1264 for (i = 0; i < is_n; i++) { 1265 row = is_idx[i]; 1266 count = (baij->i[row / bs + 1] - baij->i[row / bs]) * bs; 1267 aa = ((MatScalar *)(baij->a)) + baij->i[row / bs] * bs2 + (row % bs); 1268 for (k = 0; k < count; k++) { 1269 aa[0] = zero; 1270 aa += bs; 1271 } 1272 if (diag != 0.0) PetscUseTypeMethod(A, setvalues, 1, &row, 1, &row, &diag, INSERT_VALUES); 1273 } 1274 PetscCall(MatAssemblyEnd_SeqSBAIJ(A, MAT_FINAL_ASSEMBLY)); 1275 PetscFunctionReturn(0); 1276 } 1277 1278 PetscErrorCode MatShift_SeqSBAIJ(Mat Y, PetscScalar a) { 1279 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)Y->data; 1280 1281 PetscFunctionBegin; 1282 if (!Y->preallocated || !aij->nz) PetscCall(MatSeqSBAIJSetPreallocation(Y, Y->rmap->bs, 1, NULL)); 1283 PetscCall(MatShift_Basic(Y, a)); 1284 PetscFunctionReturn(0); 1285 } 1286 1287 /* -------------------------------------------------------------------*/ 1288 static struct _MatOps MatOps_Values = {MatSetValues_SeqSBAIJ, 1289 MatGetRow_SeqSBAIJ, 1290 MatRestoreRow_SeqSBAIJ, 1291 MatMult_SeqSBAIJ_N, 1292 /* 4*/ MatMultAdd_SeqSBAIJ_N, 1293 MatMult_SeqSBAIJ_N, /* transpose versions are same as non-transpose versions */ 1294 MatMultAdd_SeqSBAIJ_N, 1295 NULL, 1296 NULL, 1297 NULL, 1298 /* 10*/ NULL, 1299 NULL, 1300 MatCholeskyFactor_SeqSBAIJ, 1301 MatSOR_SeqSBAIJ, 1302 MatTranspose_SeqSBAIJ, 1303 /* 15*/ MatGetInfo_SeqSBAIJ, 1304 MatEqual_SeqSBAIJ, 1305 MatGetDiagonal_SeqSBAIJ, 1306 MatDiagonalScale_SeqSBAIJ, 1307 MatNorm_SeqSBAIJ, 1308 /* 20*/ NULL, 1309 MatAssemblyEnd_SeqSBAIJ, 1310 MatSetOption_SeqSBAIJ, 1311 MatZeroEntries_SeqSBAIJ, 1312 /* 24*/ NULL, 1313 NULL, 1314 NULL, 1315 NULL, 1316 NULL, 1317 /* 29*/ MatSetUp_SeqSBAIJ, 1318 NULL, 1319 NULL, 1320 NULL, 1321 NULL, 1322 /* 34*/ MatDuplicate_SeqSBAIJ, 1323 NULL, 1324 NULL, 1325 NULL, 1326 MatICCFactor_SeqSBAIJ, 1327 /* 39*/ MatAXPY_SeqSBAIJ, 1328 MatCreateSubMatrices_SeqSBAIJ, 1329 MatIncreaseOverlap_SeqSBAIJ, 1330 MatGetValues_SeqSBAIJ, 1331 MatCopy_SeqSBAIJ, 1332 /* 44*/ NULL, 1333 MatScale_SeqSBAIJ, 1334 MatShift_SeqSBAIJ, 1335 NULL, 1336 MatZeroRowsColumns_SeqSBAIJ, 1337 /* 49*/ NULL, 1338 MatGetRowIJ_SeqSBAIJ, 1339 MatRestoreRowIJ_SeqSBAIJ, 1340 NULL, 1341 NULL, 1342 /* 54*/ NULL, 1343 NULL, 1344 NULL, 1345 MatPermute_SeqSBAIJ, 1346 MatSetValuesBlocked_SeqSBAIJ, 1347 /* 59*/ MatCreateSubMatrix_SeqSBAIJ, 1348 NULL, 1349 NULL, 1350 NULL, 1351 NULL, 1352 /* 64*/ NULL, 1353 NULL, 1354 NULL, 1355 NULL, 1356 NULL, 1357 /* 69*/ MatGetRowMaxAbs_SeqSBAIJ, 1358 NULL, 1359 MatConvert_MPISBAIJ_Basic, 1360 NULL, 1361 NULL, 1362 /* 74*/ NULL, 1363 NULL, 1364 NULL, 1365 NULL, 1366 NULL, 1367 /* 79*/ NULL, 1368 NULL, 1369 NULL, 1370 MatGetInertia_SeqSBAIJ, 1371 MatLoad_SeqSBAIJ, 1372 /* 84*/ MatIsSymmetric_SeqSBAIJ, 1373 MatIsHermitian_SeqSBAIJ, 1374 MatIsStructurallySymmetric_SeqSBAIJ, 1375 NULL, 1376 NULL, 1377 /* 89*/ NULL, 1378 NULL, 1379 NULL, 1380 NULL, 1381 NULL, 1382 /* 94*/ NULL, 1383 NULL, 1384 NULL, 1385 NULL, 1386 NULL, 1387 /* 99*/ NULL, 1388 NULL, 1389 NULL, 1390 MatConjugate_SeqSBAIJ, 1391 NULL, 1392 /*104*/ NULL, 1393 MatRealPart_SeqSBAIJ, 1394 MatImaginaryPart_SeqSBAIJ, 1395 MatGetRowUpperTriangular_SeqSBAIJ, 1396 MatRestoreRowUpperTriangular_SeqSBAIJ, 1397 /*109*/ NULL, 1398 NULL, 1399 NULL, 1400 NULL, 1401 MatMissingDiagonal_SeqSBAIJ, 1402 /*114*/ NULL, 1403 NULL, 1404 NULL, 1405 NULL, 1406 NULL, 1407 /*119*/ NULL, 1408 NULL, 1409 NULL, 1410 NULL, 1411 NULL, 1412 /*124*/ NULL, 1413 NULL, 1414 NULL, 1415 NULL, 1416 NULL, 1417 /*129*/ NULL, 1418 NULL, 1419 NULL, 1420 NULL, 1421 NULL, 1422 /*134*/ NULL, 1423 NULL, 1424 NULL, 1425 NULL, 1426 NULL, 1427 /*139*/ MatSetBlockSizes_Default, 1428 NULL, 1429 NULL, 1430 NULL, 1431 NULL, 1432 /*144*/ MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ, 1433 NULL, 1434 NULL, 1435 NULL, 1436 NULL, 1437 NULL, 1438 /*150*/ NULL}; 1439 1440 PetscErrorCode MatStoreValues_SeqSBAIJ(Mat mat) { 1441 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data; 1442 PetscInt nz = aij->i[mat->rmap->N] * mat->rmap->bs * aij->bs2; 1443 1444 PetscFunctionBegin; 1445 PetscCheck(aij->nonew == 1, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1446 1447 /* allocate space for values if not already there */ 1448 if (!aij->saved_values) PetscCall(PetscMalloc1(nz + 1, &aij->saved_values)); 1449 1450 /* copy values over */ 1451 PetscCall(PetscArraycpy(aij->saved_values, aij->a, nz)); 1452 PetscFunctionReturn(0); 1453 } 1454 1455 PetscErrorCode MatRetrieveValues_SeqSBAIJ(Mat mat) { 1456 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data; 1457 PetscInt nz = aij->i[mat->rmap->N] * mat->rmap->bs * aij->bs2; 1458 1459 PetscFunctionBegin; 1460 PetscCheck(aij->nonew == 1, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1461 PetscCheck(aij->saved_values, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call MatStoreValues(A);first"); 1462 1463 /* copy values over */ 1464 PetscCall(PetscArraycpy(aij->a, aij->saved_values, nz)); 1465 PetscFunctionReturn(0); 1466 } 1467 1468 static PetscErrorCode MatSeqSBAIJSetPreallocation_SeqSBAIJ(Mat B, PetscInt bs, PetscInt nz, PetscInt *nnz) { 1469 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ *)B->data; 1470 PetscInt i, mbs, nbs, bs2; 1471 PetscBool skipallocation = PETSC_FALSE, flg = PETSC_FALSE, realalloc = PETSC_FALSE; 1472 1473 PetscFunctionBegin; 1474 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 1475 1476 PetscCall(MatSetBlockSize(B, PetscAbs(bs))); 1477 PetscCall(PetscLayoutSetUp(B->rmap)); 1478 PetscCall(PetscLayoutSetUp(B->cmap)); 1479 PetscCheck(B->rmap->N <= B->cmap->N, PETSC_COMM_SELF, PETSC_ERR_SUP, "SEQSBAIJ matrix cannot have more rows %" PetscInt_FMT " than columns %" PetscInt_FMT, B->rmap->N, B->cmap->N); 1480 PetscCall(PetscLayoutGetBlockSize(B->rmap, &bs)); 1481 1482 B->preallocated = PETSC_TRUE; 1483 1484 mbs = B->rmap->N / bs; 1485 nbs = B->cmap->n / bs; 1486 bs2 = bs * bs; 1487 1488 PetscCheck(mbs * bs == B->rmap->N && nbs * bs == B->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number rows, cols must be divisible by blocksize"); 1489 1490 if (nz == MAT_SKIP_ALLOCATION) { 1491 skipallocation = PETSC_TRUE; 1492 nz = 0; 1493 } 1494 1495 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 3; 1496 PetscCheck(nz >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "nz cannot be less than 0: value %" PetscInt_FMT, nz); 1497 if (nnz) { 1498 for (i = 0; i < mbs; i++) { 1499 PetscCheck(nnz[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "nnz cannot be less than 0: local row %" PetscInt_FMT " value %" PetscInt_FMT, i, nnz[i]); 1500 PetscCheck(nnz[i] <= nbs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "nnz cannot be greater than block row length: local row %" PetscInt_FMT " value %" PetscInt_FMT " block rowlength %" PetscInt_FMT, i, nnz[i], nbs); 1501 } 1502 } 1503 1504 B->ops->mult = MatMult_SeqSBAIJ_N; 1505 B->ops->multadd = MatMultAdd_SeqSBAIJ_N; 1506 B->ops->multtranspose = MatMult_SeqSBAIJ_N; 1507 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_N; 1508 1509 PetscCall(PetscOptionsGetBool(((PetscObject)B)->options, ((PetscObject)B)->prefix, "-mat_no_unroll", &flg, NULL)); 1510 if (!flg) { 1511 switch (bs) { 1512 case 1: 1513 B->ops->mult = MatMult_SeqSBAIJ_1; 1514 B->ops->multadd = MatMultAdd_SeqSBAIJ_1; 1515 B->ops->multtranspose = MatMult_SeqSBAIJ_1; 1516 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_1; 1517 break; 1518 case 2: 1519 B->ops->mult = MatMult_SeqSBAIJ_2; 1520 B->ops->multadd = MatMultAdd_SeqSBAIJ_2; 1521 B->ops->multtranspose = MatMult_SeqSBAIJ_2; 1522 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_2; 1523 break; 1524 case 3: 1525 B->ops->mult = MatMult_SeqSBAIJ_3; 1526 B->ops->multadd = MatMultAdd_SeqSBAIJ_3; 1527 B->ops->multtranspose = MatMult_SeqSBAIJ_3; 1528 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_3; 1529 break; 1530 case 4: 1531 B->ops->mult = MatMult_SeqSBAIJ_4; 1532 B->ops->multadd = MatMultAdd_SeqSBAIJ_4; 1533 B->ops->multtranspose = MatMult_SeqSBAIJ_4; 1534 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_4; 1535 break; 1536 case 5: 1537 B->ops->mult = MatMult_SeqSBAIJ_5; 1538 B->ops->multadd = MatMultAdd_SeqSBAIJ_5; 1539 B->ops->multtranspose = MatMult_SeqSBAIJ_5; 1540 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_5; 1541 break; 1542 case 6: 1543 B->ops->mult = MatMult_SeqSBAIJ_6; 1544 B->ops->multadd = MatMultAdd_SeqSBAIJ_6; 1545 B->ops->multtranspose = MatMult_SeqSBAIJ_6; 1546 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_6; 1547 break; 1548 case 7: 1549 B->ops->mult = MatMult_SeqSBAIJ_7; 1550 B->ops->multadd = MatMultAdd_SeqSBAIJ_7; 1551 B->ops->multtranspose = MatMult_SeqSBAIJ_7; 1552 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_7; 1553 break; 1554 } 1555 } 1556 1557 b->mbs = mbs; 1558 b->nbs = nbs; 1559 if (!skipallocation) { 1560 if (!b->imax) { 1561 PetscCall(PetscMalloc2(mbs, &b->imax, mbs, &b->ilen)); 1562 1563 b->free_imax_ilen = PETSC_TRUE; 1564 } 1565 if (!nnz) { 1566 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 1567 else if (nz <= 0) nz = 1; 1568 nz = PetscMin(nbs, nz); 1569 for (i = 0; i < mbs; i++) b->imax[i] = nz; 1570 PetscCall(PetscIntMultError(nz, mbs, &nz)); 1571 } else { 1572 PetscInt64 nz64 = 0; 1573 for (i = 0; i < mbs; i++) { 1574 b->imax[i] = nnz[i]; 1575 nz64 += nnz[i]; 1576 } 1577 PetscCall(PetscIntCast(nz64, &nz)); 1578 } 1579 /* b->ilen will count nonzeros in each block row so far. */ 1580 for (i = 0; i < mbs; i++) b->ilen[i] = 0; 1581 /* nz=(nz+mbs)/2; */ /* total diagonal and superdiagonal nonzero blocks */ 1582 1583 /* allocate the matrix space */ 1584 PetscCall(MatSeqXAIJFreeAIJ(B, &b->a, &b->j, &b->i)); 1585 PetscCall(PetscMalloc3(bs2 * nz, &b->a, nz, &b->j, B->rmap->N + 1, &b->i)); 1586 PetscCall(PetscArrayzero(b->a, nz * bs2)); 1587 PetscCall(PetscArrayzero(b->j, nz)); 1588 1589 b->singlemalloc = PETSC_TRUE; 1590 1591 /* pointer to beginning of each row */ 1592 b->i[0] = 0; 1593 for (i = 1; i < mbs + 1; i++) b->i[i] = b->i[i - 1] + b->imax[i - 1]; 1594 1595 b->free_a = PETSC_TRUE; 1596 b->free_ij = PETSC_TRUE; 1597 } else { 1598 b->free_a = PETSC_FALSE; 1599 b->free_ij = PETSC_FALSE; 1600 } 1601 1602 b->bs2 = bs2; 1603 b->nz = 0; 1604 b->maxnz = nz; 1605 b->inew = NULL; 1606 b->jnew = NULL; 1607 b->anew = NULL; 1608 b->a2anew = NULL; 1609 b->permute = PETSC_FALSE; 1610 1611 B->was_assembled = PETSC_FALSE; 1612 B->assembled = PETSC_FALSE; 1613 if (realalloc) PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE)); 1614 PetscFunctionReturn(0); 1615 } 1616 1617 PetscErrorCode MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ(Mat B, PetscInt bs, const PetscInt ii[], const PetscInt jj[], const PetscScalar V[]) { 1618 PetscInt i, j, m, nz, anz, nz_max = 0, *nnz; 1619 PetscScalar *values = NULL; 1620 PetscBool roworiented = ((Mat_SeqSBAIJ *)B->data)->roworiented; 1621 1622 PetscFunctionBegin; 1623 PetscCheck(bs >= 1, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "Invalid block size specified, must be positive but it is %" PetscInt_FMT, bs); 1624 PetscCall(PetscLayoutSetBlockSize(B->rmap, bs)); 1625 PetscCall(PetscLayoutSetBlockSize(B->cmap, bs)); 1626 PetscCall(PetscLayoutSetUp(B->rmap)); 1627 PetscCall(PetscLayoutSetUp(B->cmap)); 1628 PetscCall(PetscLayoutGetBlockSize(B->rmap, &bs)); 1629 m = B->rmap->n / bs; 1630 1631 PetscCheck(!ii[0], PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "ii[0] must be 0 but it is %" PetscInt_FMT, ii[0]); 1632 PetscCall(PetscMalloc1(m + 1, &nnz)); 1633 for (i = 0; i < m; i++) { 1634 nz = ii[i + 1] - ii[i]; 1635 PetscCheck(nz >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row %" PetscInt_FMT " has a negative number of columns %" PetscInt_FMT, i, nz); 1636 anz = 0; 1637 for (j = 0; j < nz; j++) { 1638 /* count only values on the diagonal or above */ 1639 if (jj[ii[i] + j] >= i) { 1640 anz = nz - j; 1641 break; 1642 } 1643 } 1644 nz_max = PetscMax(nz_max, anz); 1645 nnz[i] = anz; 1646 } 1647 PetscCall(MatSeqSBAIJSetPreallocation(B, bs, 0, nnz)); 1648 PetscCall(PetscFree(nnz)); 1649 1650 values = (PetscScalar *)V; 1651 if (!values) PetscCall(PetscCalloc1(bs * bs * nz_max, &values)); 1652 for (i = 0; i < m; i++) { 1653 PetscInt ncols = ii[i + 1] - ii[i]; 1654 const PetscInt *icols = jj + ii[i]; 1655 if (!roworiented || bs == 1) { 1656 const PetscScalar *svals = values + (V ? (bs * bs * ii[i]) : 0); 1657 PetscCall(MatSetValuesBlocked_SeqSBAIJ(B, 1, &i, ncols, icols, svals, INSERT_VALUES)); 1658 } else { 1659 for (j = 0; j < ncols; j++) { 1660 const PetscScalar *svals = values + (V ? (bs * bs * (ii[i] + j)) : 0); 1661 PetscCall(MatSetValuesBlocked_SeqSBAIJ(B, 1, &i, 1, &icols[j], svals, INSERT_VALUES)); 1662 } 1663 } 1664 } 1665 if (!V) PetscCall(PetscFree(values)); 1666 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 1667 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 1668 PetscCall(MatSetOption(B, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE)); 1669 PetscFunctionReturn(0); 1670 } 1671 1672 /* 1673 This is used to set the numeric factorization for both Cholesky and ICC symbolic factorization 1674 */ 1675 PetscErrorCode MatSeqSBAIJSetNumericFactorization_inplace(Mat B, PetscBool natural) { 1676 PetscBool flg = PETSC_FALSE; 1677 PetscInt bs = B->rmap->bs; 1678 1679 PetscFunctionBegin; 1680 PetscCall(PetscOptionsGetBool(((PetscObject)B)->options, ((PetscObject)B)->prefix, "-mat_no_unroll", &flg, NULL)); 1681 if (flg) bs = 8; 1682 1683 if (!natural) { 1684 switch (bs) { 1685 case 1: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_inplace; break; 1686 case 2: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2; break; 1687 case 3: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3; break; 1688 case 4: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4; break; 1689 case 5: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5; break; 1690 case 6: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6; break; 1691 case 7: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7; break; 1692 default: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N; break; 1693 } 1694 } else { 1695 switch (bs) { 1696 case 1: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering_inplace; break; 1697 case 2: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering; break; 1698 case 3: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3_NaturalOrdering; break; 1699 case 4: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4_NaturalOrdering; break; 1700 case 5: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5_NaturalOrdering; break; 1701 case 6: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6_NaturalOrdering; break; 1702 case 7: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7_NaturalOrdering; break; 1703 default: B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N_NaturalOrdering; break; 1704 } 1705 } 1706 PetscFunctionReturn(0); 1707 } 1708 1709 PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqAIJ(Mat, MatType, MatReuse, Mat *); 1710 PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat, MatType, MatReuse, Mat *); 1711 static PetscErrorCode MatFactorGetSolverType_petsc(Mat A, MatSolverType *type) { 1712 PetscFunctionBegin; 1713 *type = MATSOLVERPETSC; 1714 PetscFunctionReturn(0); 1715 } 1716 1717 PETSC_INTERN PetscErrorCode MatGetFactor_seqsbaij_petsc(Mat A, MatFactorType ftype, Mat *B) { 1718 PetscInt n = A->rmap->n; 1719 1720 PetscFunctionBegin; 1721 #if defined(PETSC_USE_COMPLEX) 1722 PetscCheck(A->hermitian != PETSC_BOOL3_TRUE || A->symmetric == PETSC_BOOL3_TRUE || (ftype != MAT_FACTOR_CHOLESKY && ftype != MAT_FACTOR_ICC), PETSC_COMM_SELF, PETSC_ERR_SUP, "Hermitian CHOLESKY or ICC Factor is not supported"); 1723 #endif 1724 1725 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), B)); 1726 PetscCall(MatSetSizes(*B, n, n, n, n)); 1727 if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) { 1728 PetscCall(MatSetType(*B, MATSEQSBAIJ)); 1729 PetscCall(MatSeqSBAIJSetPreallocation(*B, A->rmap->bs, MAT_SKIP_ALLOCATION, NULL)); 1730 1731 (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJ; 1732 (*B)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqSBAIJ; 1733 PetscCall(PetscStrallocpy(MATORDERINGNATURAL, (char **)&(*B)->preferredordering[MAT_FACTOR_CHOLESKY])); 1734 PetscCall(PetscStrallocpy(MATORDERINGNATURAL, (char **)&(*B)->preferredordering[MAT_FACTOR_ICC])); 1735 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Factor type not supported"); 1736 1737 (*B)->factortype = ftype; 1738 (*B)->canuseordering = PETSC_TRUE; 1739 PetscCall(PetscFree((*B)->solvertype)); 1740 PetscCall(PetscStrallocpy(MATSOLVERPETSC, &(*B)->solvertype)); 1741 PetscCall(PetscObjectComposeFunction((PetscObject)*B, "MatFactorGetSolverType_C", MatFactorGetSolverType_petsc)); 1742 PetscFunctionReturn(0); 1743 } 1744 1745 /*@C 1746 MatSeqSBAIJGetArray - gives access to the array where the data for a `MATSEQSBAIJ` matrix is stored 1747 1748 Not Collective 1749 1750 Input Parameter: 1751 . mat - a `MATSEQSBAIJ` matrix 1752 1753 Output Parameter: 1754 . array - pointer to the data 1755 1756 Level: intermediate 1757 1758 .seealso: `MATSEQSBAIJ`, `MatSeqSBAIJRestoreArray()`, `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArray()` 1759 @*/ 1760 PetscErrorCode MatSeqSBAIJGetArray(Mat A, PetscScalar **array) { 1761 PetscFunctionBegin; 1762 PetscUseMethod(A, "MatSeqSBAIJGetArray_C", (Mat, PetscScalar **), (A, array)); 1763 PetscFunctionReturn(0); 1764 } 1765 1766 /*@C 1767 MatSeqSBAIJRestoreArray - returns access to the array where the data for a `MATSEQSBAIJ` matrix is stored obtained by `MatSeqSBAIJGetArray()` 1768 1769 Not Collective 1770 1771 Input Parameters: 1772 + mat - a MATSEQSBAIJ matrix 1773 - array - pointer to the data 1774 1775 Level: intermediate 1776 1777 .seealso: `MATSEQSBAIJ`, `MatSeqSBAIJGetArray()`, `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArray()` 1778 @*/ 1779 PetscErrorCode MatSeqSBAIJRestoreArray(Mat A, PetscScalar **array) { 1780 PetscFunctionBegin; 1781 PetscUseMethod(A, "MatSeqSBAIJRestoreArray_C", (Mat, PetscScalar **), (A, array)); 1782 PetscFunctionReturn(0); 1783 } 1784 1785 /*MC 1786 MATSEQSBAIJ - MATSEQSBAIJ = "seqsbaij" - A matrix type to be used for sequential symmetric block sparse matrices, 1787 based on block compressed sparse row format. Only the upper triangular portion of the matrix is stored. 1788 1789 For complex numbers by default this matrix is symmetric, NOT Hermitian symmetric. To make it Hermitian symmetric you 1790 can call `MatSetOption`(`Mat`, `MAT_HERMITIAN`). 1791 1792 Options Database Keys: 1793 . -mat_type seqsbaij - sets the matrix type to "seqsbaij" during a call to `MatSetFromOptions()` 1794 1795 Notes: 1796 By default if you insert values into the lower triangular part of the matrix they are simply ignored (since they are not 1797 stored and it is assumed they symmetric to the upper triangular). If you call `MatSetOption`(`Mat`,`MAT_IGNORE_LOWER_TRIANGULAR`,`PETSC_FALSE`) or use 1798 the options database -mat_ignore_lower_triangular false it will generate an error if you try to set a value in the lower triangular portion. 1799 1800 The number of rows in the matrix must be less than or equal to the number of columns 1801 1802 Level: beginner 1803 1804 .seealso: `MATSEQSBAIJ`, `MatCreateSeqSBAIJ()`, `MatType`, `MATMPISBAIJ` 1805 M*/ 1806 PETSC_EXTERN PetscErrorCode MatCreate_SeqSBAIJ(Mat B) { 1807 Mat_SeqSBAIJ *b; 1808 PetscMPIInt size; 1809 PetscBool no_unroll = PETSC_FALSE, no_inode = PETSC_FALSE; 1810 1811 PetscFunctionBegin; 1812 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B), &size)); 1813 PetscCheck(size <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Comm must be of size 1"); 1814 1815 PetscCall(PetscNew(&b)); 1816 B->data = (void *)b; 1817 PetscCall(PetscMemcpy(B->ops, &MatOps_Values, sizeof(struct _MatOps))); 1818 1819 B->ops->destroy = MatDestroy_SeqSBAIJ; 1820 B->ops->view = MatView_SeqSBAIJ; 1821 b->row = NULL; 1822 b->icol = NULL; 1823 b->reallocs = 0; 1824 b->saved_values = NULL; 1825 b->inode.limit = 5; 1826 b->inode.max_limit = 5; 1827 1828 b->roworiented = PETSC_TRUE; 1829 b->nonew = 0; 1830 b->diag = NULL; 1831 b->solve_work = NULL; 1832 b->mult_work = NULL; 1833 B->spptr = NULL; 1834 B->info.nz_unneeded = (PetscReal)b->maxnz * b->bs2; 1835 b->keepnonzeropattern = PETSC_FALSE; 1836 1837 b->inew = NULL; 1838 b->jnew = NULL; 1839 b->anew = NULL; 1840 b->a2anew = NULL; 1841 b->permute = PETSC_FALSE; 1842 1843 b->ignore_ltriangular = PETSC_TRUE; 1844 1845 PetscCall(PetscOptionsGetBool(((PetscObject)B)->options, ((PetscObject)B)->prefix, "-mat_ignore_lower_triangular", &b->ignore_ltriangular, NULL)); 1846 1847 b->getrow_utriangular = PETSC_FALSE; 1848 1849 PetscCall(PetscOptionsGetBool(((PetscObject)B)->options, ((PetscObject)B)->prefix, "-mat_getrow_uppertriangular", &b->getrow_utriangular, NULL)); 1850 1851 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqSBAIJGetArray_C", MatSeqSBAIJGetArray_SeqSBAIJ)); 1852 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqSBAIJRestoreArray_C", MatSeqSBAIJRestoreArray_SeqSBAIJ)); 1853 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatStoreValues_C", MatStoreValues_SeqSBAIJ)); 1854 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatRetrieveValues_C", MatRetrieveValues_SeqSBAIJ)); 1855 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqSBAIJSetColumnIndices_C", MatSeqSBAIJSetColumnIndices_SeqSBAIJ)); 1856 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqsbaij_seqaij_C", MatConvert_SeqSBAIJ_SeqAIJ)); 1857 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqsbaij_seqbaij_C", MatConvert_SeqSBAIJ_SeqBAIJ)); 1858 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqSBAIJSetPreallocation_C", MatSeqSBAIJSetPreallocation_SeqSBAIJ)); 1859 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqSBAIJSetPreallocationCSR_C", MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ)); 1860 #if defined(PETSC_HAVE_ELEMENTAL) 1861 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqsbaij_elemental_C", MatConvert_SeqSBAIJ_Elemental)); 1862 #endif 1863 #if defined(PETSC_HAVE_SCALAPACK) 1864 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqsbaij_scalapack_C", MatConvert_SBAIJ_ScaLAPACK)); 1865 #endif 1866 1867 B->symmetry_eternal = PETSC_TRUE; 1868 B->structural_symmetry_eternal = PETSC_TRUE; 1869 B->symmetric = PETSC_BOOL3_TRUE; 1870 B->structurally_symmetric = PETSC_BOOL3_TRUE; 1871 #if defined(PETSC_USE_COMPLEX) 1872 B->hermitian = PETSC_BOOL3_FALSE; 1873 #else 1874 B->hermitian = PETSC_BOOL3_TRUE; 1875 #endif 1876 1877 PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATSEQSBAIJ)); 1878 1879 PetscOptionsBegin(PetscObjectComm((PetscObject)B), ((PetscObject)B)->prefix, "Options for SEQSBAIJ matrix", "Mat"); 1880 PetscCall(PetscOptionsBool("-mat_no_unroll", "Do not optimize for inodes (slower)", NULL, no_unroll, &no_unroll, NULL)); 1881 if (no_unroll) PetscCall(PetscInfo(B, "Not using Inode routines due to -mat_no_unroll\n")); 1882 PetscCall(PetscOptionsBool("-mat_no_inode", "Do not optimize for inodes (slower)", NULL, no_inode, &no_inode, NULL)); 1883 if (no_inode) PetscCall(PetscInfo(B, "Not using Inode routines due to -mat_no_inode\n")); 1884 PetscCall(PetscOptionsInt("-mat_inode_limit", "Do not use inodes larger then this value", NULL, b->inode.limit, &b->inode.limit, NULL)); 1885 PetscOptionsEnd(); 1886 b->inode.use = (PetscBool)(!(no_unroll || no_inode)); 1887 if (b->inode.limit > b->inode.max_limit) b->inode.limit = b->inode.max_limit; 1888 PetscFunctionReturn(0); 1889 } 1890 1891 /*@C 1892 MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block 1893 compressed row) `MATSEQSBAIJ` format. For good matrix assembly performance the 1894 user should preallocate the matrix storage by setting the parameter nz 1895 (or the array nnz). By setting these parameters accurately, performance 1896 during matrix assembly can be increased by more than a factor of 50. 1897 1898 Collective on Mat 1899 1900 Input Parameters: 1901 + B - the symmetric matrix 1902 . bs - size of block, the blocks are ALWAYS square. One can use `MatSetBlockSizes()` to set a different row and column blocksize but the row 1903 blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with `MatCreateVecs()` 1904 . nz - number of block nonzeros per block row (same for all rows) 1905 - nnz - array containing the number of block nonzeros in the upper triangular plus 1906 diagonal portion of each block (possibly different for each block row) or NULL 1907 1908 Options Database Keys: 1909 + -mat_no_unroll - uses code that does not unroll the loops in the 1910 block calculations (much slower) 1911 - -mat_block_size - size of the blocks to use (only works if a negative bs is passed in 1912 1913 Level: intermediate 1914 1915 Notes: 1916 Specify the preallocated storage with either nz or nnz (not both). 1917 Set nz = `PETSC_DEFAULT` and nnz = NULL for PETSc to control dynamic memory 1918 allocation. See [Sparse Matrices](sec_matsparse) for details. 1919 1920 You can call `MatGetInfo()` to get information on how effective the preallocation was; 1921 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 1922 You can also run with the option -info and look for messages with the string 1923 malloc in them to see if additional memory allocation was needed. 1924 1925 If the nnz parameter is given then the nz parameter is ignored 1926 1927 .seealso: [Sparse Matrices](sec_matsparse), `MATSEQSBAIJ`, `MatCreate()`, `MatCreateSeqAIJ()`, `MatSetValues()`, `MatCreateSBAIJ()` 1928 @*/ 1929 PetscErrorCode MatSeqSBAIJSetPreallocation(Mat B, PetscInt bs, PetscInt nz, const PetscInt nnz[]) { 1930 PetscFunctionBegin; 1931 PetscValidHeaderSpecific(B, MAT_CLASSID, 1); 1932 PetscValidType(B, 1); 1933 PetscValidLogicalCollectiveInt(B, bs, 2); 1934 PetscTryMethod(B, "MatSeqSBAIJSetPreallocation_C", (Mat, PetscInt, PetscInt, const PetscInt[]), (B, bs, nz, nnz)); 1935 PetscFunctionReturn(0); 1936 } 1937 1938 /*@C 1939 MatSeqSBAIJSetPreallocationCSR - Creates a sparse parallel matrix in `MATSEQSBAIJ` format using the given nonzero structure and (optional) numerical values 1940 1941 Input Parameters: 1942 + B - the matrix 1943 . bs - size of block, the blocks are ALWAYS square. 1944 . i - the indices into j for the start of each local row (starts with zero) 1945 . j - the column indices for each local row (starts with zero) these must be sorted for each row 1946 - v - optional values in the matrix 1947 1948 Level: advanced 1949 1950 Notes: 1951 The order of the entries in values is specified by the `MatOption` `MAT_ROW_ORIENTED`. For example, C programs 1952 may want to use the default `MAT_ROW_ORIENTED` = `PETSC_TRUE` and use an array v[nnz][bs][bs] where the second index is 1953 over rows within a block and the last index is over columns within a block row. Fortran programs will likely set 1954 `MAT_ROW_ORIENTED` = `PETSC_FALSE` and use a Fortran array v(bs,bs,nnz) in which the first index is over rows within a 1955 block column and the second index is over columns within a block. 1956 1957 Any entries below the diagonal are ignored 1958 1959 Though this routine has Preallocation() in the name it also sets the exact nonzero locations of the matrix entries 1960 and usually the numerical values as well 1961 1962 .seealso: `MATSEQSBAIJ`, `MatCreate()`, `MatCreateSeqSBAIJ()`, `MatSetValuesBlocked()`, `MatSeqSBAIJSetPreallocation()`, `MATSEQSBAIJ` 1963 @*/ 1964 PetscErrorCode MatSeqSBAIJSetPreallocationCSR(Mat B, PetscInt bs, const PetscInt i[], const PetscInt j[], const PetscScalar v[]) { 1965 PetscFunctionBegin; 1966 PetscValidHeaderSpecific(B, MAT_CLASSID, 1); 1967 PetscValidType(B, 1); 1968 PetscValidLogicalCollectiveInt(B, bs, 2); 1969 PetscTryMethod(B, "MatSeqSBAIJSetPreallocationCSR_C", (Mat, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[]), (B, bs, i, j, v)); 1970 PetscFunctionReturn(0); 1971 } 1972 1973 /*@C 1974 MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block 1975 compressed row) `MATSEQSBAIJ` format. For good matrix assembly performance the 1976 user should preallocate the matrix storage by setting the parameter nz 1977 (or the array nnz). By setting these parameters accurately, performance 1978 during matrix assembly can be increased by more than a factor of 50. 1979 1980 Collective 1981 1982 Input Parameters: 1983 + comm - MPI communicator, set to `PETSC_COMM_SELF` 1984 . bs - size of block, the blocks are ALWAYS square. One can use `MatSetBlockSizes()` to set a different row and column blocksize but the row 1985 blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs() 1986 . m - number of rows, or number of columns 1987 . nz - number of block nonzeros per block row (same for all rows) 1988 - nnz - array containing the number of block nonzeros in the upper triangular plus 1989 diagonal portion of each block (possibly different for each block row) or NULL 1990 1991 Output Parameter: 1992 . A - the symmetric matrix 1993 1994 Options Database Keys: 1995 + -mat_no_unroll - uses code that does not unroll the loops in the 1996 block calculations (much slower) 1997 - -mat_block_size - size of the blocks to use 1998 1999 Level: intermediate 2000 2001 It is recommended that one use the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`, 2002 MatXXXXSetPreallocation() paradigm instead of this routine directly. 2003 [MatXXXXSetPreallocation() is, for example, `MatSeqAIJSetPreallocation()`] 2004 2005 Notes: 2006 The number of rows and columns must be divisible by blocksize. 2007 This matrix type does not support complex Hermitian operation. 2008 2009 Specify the preallocated storage with either nz or nnz (not both). 2010 Set nz = `PETSC_DEFAULT` and nnz = NULL for PETSc to control dynamic memory 2011 allocation. See [Sparse Matrices](sec_matsparse) for details. 2012 2013 If the nnz parameter is given then the nz parameter is ignored 2014 2015 .seealso: [Sparse Matrices](sec_matsparse), `MATSEQSBAIJ`, `MatCreate()`, `MatCreateSeqAIJ()`, `MatSetValues()`, `MatCreateSBAIJ()` 2016 @*/ 2017 PetscErrorCode MatCreateSeqSBAIJ(MPI_Comm comm, PetscInt bs, PetscInt m, PetscInt n, PetscInt nz, const PetscInt nnz[], Mat *A) { 2018 PetscFunctionBegin; 2019 PetscCall(MatCreate(comm, A)); 2020 PetscCall(MatSetSizes(*A, m, n, m, n)); 2021 PetscCall(MatSetType(*A, MATSEQSBAIJ)); 2022 PetscCall(MatSeqSBAIJSetPreallocation(*A, bs, nz, (PetscInt *)nnz)); 2023 PetscFunctionReturn(0); 2024 } 2025 2026 PetscErrorCode MatDuplicate_SeqSBAIJ(Mat A, MatDuplicateOption cpvalues, Mat *B) { 2027 Mat C; 2028 Mat_SeqSBAIJ *c, *a = (Mat_SeqSBAIJ *)A->data; 2029 PetscInt i, mbs = a->mbs, nz = a->nz, bs2 = a->bs2; 2030 2031 PetscFunctionBegin; 2032 PetscCheck(a->i[mbs] == nz, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Corrupt matrix"); 2033 2034 *B = NULL; 2035 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &C)); 2036 PetscCall(MatSetSizes(C, A->rmap->N, A->cmap->n, A->rmap->N, A->cmap->n)); 2037 PetscCall(MatSetBlockSizesFromMats(C, A, A)); 2038 PetscCall(MatSetType(C, MATSEQSBAIJ)); 2039 c = (Mat_SeqSBAIJ *)C->data; 2040 2041 C->preallocated = PETSC_TRUE; 2042 C->factortype = A->factortype; 2043 c->row = NULL; 2044 c->icol = NULL; 2045 c->saved_values = NULL; 2046 c->keepnonzeropattern = a->keepnonzeropattern; 2047 C->assembled = PETSC_TRUE; 2048 2049 PetscCall(PetscLayoutReference(A->rmap, &C->rmap)); 2050 PetscCall(PetscLayoutReference(A->cmap, &C->cmap)); 2051 c->bs2 = a->bs2; 2052 c->mbs = a->mbs; 2053 c->nbs = a->nbs; 2054 2055 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2056 c->imax = a->imax; 2057 c->ilen = a->ilen; 2058 c->free_imax_ilen = PETSC_FALSE; 2059 } else { 2060 PetscCall(PetscMalloc2((mbs + 1), &c->imax, (mbs + 1), &c->ilen)); 2061 for (i = 0; i < mbs; i++) { 2062 c->imax[i] = a->imax[i]; 2063 c->ilen[i] = a->ilen[i]; 2064 } 2065 c->free_imax_ilen = PETSC_TRUE; 2066 } 2067 2068 /* allocate the matrix space */ 2069 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2070 PetscCall(PetscMalloc1(bs2 * nz, &c->a)); 2071 c->i = a->i; 2072 c->j = a->j; 2073 c->singlemalloc = PETSC_FALSE; 2074 c->free_a = PETSC_TRUE; 2075 c->free_ij = PETSC_FALSE; 2076 c->parent = A; 2077 PetscCall(PetscObjectReference((PetscObject)A)); 2078 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE)); 2079 PetscCall(MatSetOption(C, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE)); 2080 } else { 2081 PetscCall(PetscMalloc3(bs2 * nz, &c->a, nz, &c->j, mbs + 1, &c->i)); 2082 PetscCall(PetscArraycpy(c->i, a->i, mbs + 1)); 2083 c->singlemalloc = PETSC_TRUE; 2084 c->free_a = PETSC_TRUE; 2085 c->free_ij = PETSC_TRUE; 2086 } 2087 if (mbs > 0) { 2088 if (cpvalues != MAT_SHARE_NONZERO_PATTERN) PetscCall(PetscArraycpy(c->j, a->j, nz)); 2089 if (cpvalues == MAT_COPY_VALUES) { 2090 PetscCall(PetscArraycpy(c->a, a->a, bs2 * nz)); 2091 } else { 2092 PetscCall(PetscArrayzero(c->a, bs2 * nz)); 2093 } 2094 if (a->jshort) { 2095 /* cannot share jshort, it is reallocated in MatAssemblyEnd_SeqSBAIJ() */ 2096 /* if the parent matrix is reassembled, this child matrix will never notice */ 2097 PetscCall(PetscMalloc1(nz, &c->jshort)); 2098 PetscCall(PetscArraycpy(c->jshort, a->jshort, nz)); 2099 2100 c->free_jshort = PETSC_TRUE; 2101 } 2102 } 2103 2104 c->roworiented = a->roworiented; 2105 c->nonew = a->nonew; 2106 2107 if (a->diag) { 2108 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2109 c->diag = a->diag; 2110 c->free_diag = PETSC_FALSE; 2111 } else { 2112 PetscCall(PetscMalloc1(mbs, &c->diag)); 2113 for (i = 0; i < mbs; i++) c->diag[i] = a->diag[i]; 2114 c->free_diag = PETSC_TRUE; 2115 } 2116 } 2117 c->nz = a->nz; 2118 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 2119 c->solve_work = NULL; 2120 c->mult_work = NULL; 2121 2122 *B = C; 2123 PetscCall(PetscFunctionListDuplicate(((PetscObject)A)->qlist, &((PetscObject)C)->qlist)); 2124 PetscFunctionReturn(0); 2125 } 2126 2127 /* Used for both SeqBAIJ and SeqSBAIJ matrices */ 2128 #define MatLoad_SeqSBAIJ_Binary MatLoad_SeqBAIJ_Binary 2129 2130 PetscErrorCode MatLoad_SeqSBAIJ(Mat mat, PetscViewer viewer) { 2131 PetscBool isbinary; 2132 2133 PetscFunctionBegin; 2134 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary)); 2135 PetscCheck(isbinary, PetscObjectComm((PetscObject)viewer), PETSC_ERR_SUP, "Viewer type %s not yet supported for reading %s matrices", ((PetscObject)viewer)->type_name, ((PetscObject)mat)->type_name); 2136 PetscCall(MatLoad_SeqSBAIJ_Binary(mat, viewer)); 2137 PetscFunctionReturn(0); 2138 } 2139 2140 /*@ 2141 MatCreateSeqSBAIJWithArrays - Creates an sequential `MATSEQSBAIJ` matrix using matrix elements 2142 (upper triangular entries in CSR format) provided by the user. 2143 2144 Collective 2145 2146 Input Parameters: 2147 + comm - must be an MPI communicator of size 1 2148 . bs - size of block 2149 . m - number of rows 2150 . n - number of columns 2151 . i - row indices; that is i[0] = 0, i[row] = i[row-1] + number of block elements in that row block row of the matrix 2152 . j - column indices 2153 - a - matrix values 2154 2155 Output Parameter: 2156 . mat - the matrix 2157 2158 Level: advanced 2159 2160 Notes: 2161 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2162 once the matrix is destroyed 2163 2164 You cannot set new nonzero locations into this matrix, that will generate an error. 2165 2166 The i and j indices are 0 based 2167 2168 When block size is greater than 1 the matrix values must be stored using the SBAIJ storage format (see the SBAIJ source code to determine this). For block size of 1 2169 it is the regular CSR format excluding the lower triangular elements. 2170 2171 .seealso: `MATSEQSBAIJ`, `MatCreate()`, `MatCreateSBAIJ()`, `MatCreateSeqSBAIJ()` 2172 @*/ 2173 PetscErrorCode MatCreateSeqSBAIJWithArrays(MPI_Comm comm, PetscInt bs, PetscInt m, PetscInt n, PetscInt i[], PetscInt j[], PetscScalar a[], Mat *mat) { 2174 PetscInt ii; 2175 Mat_SeqSBAIJ *sbaij; 2176 2177 PetscFunctionBegin; 2178 PetscCheck(bs == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "block size %" PetscInt_FMT " > 1 is not supported yet", bs); 2179 PetscCheck(m == 0 || i[0] == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "i (row indices) must start with 0"); 2180 2181 PetscCall(MatCreate(comm, mat)); 2182 PetscCall(MatSetSizes(*mat, m, n, m, n)); 2183 PetscCall(MatSetType(*mat, MATSEQSBAIJ)); 2184 PetscCall(MatSeqSBAIJSetPreallocation(*mat, bs, MAT_SKIP_ALLOCATION, NULL)); 2185 sbaij = (Mat_SeqSBAIJ *)(*mat)->data; 2186 PetscCall(PetscMalloc2(m, &sbaij->imax, m, &sbaij->ilen)); 2187 2188 sbaij->i = i; 2189 sbaij->j = j; 2190 sbaij->a = a; 2191 2192 sbaij->singlemalloc = PETSC_FALSE; 2193 sbaij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2194 sbaij->free_a = PETSC_FALSE; 2195 sbaij->free_ij = PETSC_FALSE; 2196 sbaij->free_imax_ilen = PETSC_TRUE; 2197 2198 for (ii = 0; ii < m; ii++) { 2199 sbaij->ilen[ii] = sbaij->imax[ii] = i[ii + 1] - i[ii]; 2200 PetscCheck(i[ii + 1] >= i[ii], PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative row length in i (row indices) row = %" PetscInt_FMT " length = %" PetscInt_FMT, ii, i[ii + 1] - i[ii]); 2201 } 2202 if (PetscDefined(USE_DEBUG)) { 2203 for (ii = 0; ii < sbaij->i[m]; ii++) { 2204 PetscCheck(j[ii] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative column index at location = %" PetscInt_FMT " index = %" PetscInt_FMT, ii, j[ii]); 2205 PetscCheck(j[ii] < n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column index too large at location = %" PetscInt_FMT " index = %" PetscInt_FMT, ii, j[ii]); 2206 } 2207 } 2208 2209 PetscCall(MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY)); 2210 PetscCall(MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY)); 2211 PetscFunctionReturn(0); 2212 } 2213 2214 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ(MPI_Comm comm, Mat inmat, PetscInt n, MatReuse scall, Mat *outmat) { 2215 PetscFunctionBegin; 2216 PetscCall(MatCreateMPIMatConcatenateSeqMat_MPISBAIJ(comm, inmat, n, scall, outmat)); 2217 PetscFunctionReturn(0); 2218 } 2219