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