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