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