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