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