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