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