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