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