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