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