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