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