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: `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 /* -------------------------------------------------------------------*/ 1341 static struct _MatOps MatOps_Values = {MatSetValues_SeqSBAIJ, 1342 MatGetRow_SeqSBAIJ, 1343 MatRestoreRow_SeqSBAIJ, 1344 MatMult_SeqSBAIJ_N, 1345 /* 4*/ MatMultAdd_SeqSBAIJ_N, 1346 MatMult_SeqSBAIJ_N, /* transpose versions are same as non-transpose versions */ 1347 MatMultAdd_SeqSBAIJ_N, 1348 NULL, 1349 NULL, 1350 NULL, 1351 /* 10*/ NULL, 1352 NULL, 1353 MatCholeskyFactor_SeqSBAIJ, 1354 MatSOR_SeqSBAIJ, 1355 MatTranspose_SeqSBAIJ, 1356 /* 15*/ MatGetInfo_SeqSBAIJ, 1357 MatEqual_SeqSBAIJ, 1358 MatGetDiagonal_SeqSBAIJ, 1359 MatDiagonalScale_SeqSBAIJ, 1360 MatNorm_SeqSBAIJ, 1361 /* 20*/ NULL, 1362 MatAssemblyEnd_SeqSBAIJ, 1363 MatSetOption_SeqSBAIJ, 1364 MatZeroEntries_SeqSBAIJ, 1365 /* 24*/ NULL, 1366 NULL, 1367 NULL, 1368 NULL, 1369 NULL, 1370 /* 29*/ MatSetUp_Seq_Hash, 1371 NULL, 1372 NULL, 1373 NULL, 1374 NULL, 1375 /* 34*/ MatDuplicate_SeqSBAIJ, 1376 NULL, 1377 NULL, 1378 NULL, 1379 MatICCFactor_SeqSBAIJ, 1380 /* 39*/ MatAXPY_SeqSBAIJ, 1381 MatCreateSubMatrices_SeqSBAIJ, 1382 MatIncreaseOverlap_SeqSBAIJ, 1383 MatGetValues_SeqSBAIJ, 1384 MatCopy_SeqSBAIJ, 1385 /* 44*/ NULL, 1386 MatScale_SeqSBAIJ, 1387 MatShift_SeqSBAIJ, 1388 NULL, 1389 MatZeroRowsColumns_SeqSBAIJ, 1390 /* 49*/ NULL, 1391 MatGetRowIJ_SeqSBAIJ, 1392 MatRestoreRowIJ_SeqSBAIJ, 1393 NULL, 1394 NULL, 1395 /* 54*/ NULL, 1396 NULL, 1397 NULL, 1398 MatPermute_SeqSBAIJ, 1399 MatSetValuesBlocked_SeqSBAIJ, 1400 /* 59*/ MatCreateSubMatrix_SeqSBAIJ, 1401 NULL, 1402 NULL, 1403 NULL, 1404 NULL, 1405 /* 64*/ NULL, 1406 NULL, 1407 NULL, 1408 NULL, 1409 NULL, 1410 /* 69*/ MatGetRowMaxAbs_SeqSBAIJ, 1411 NULL, 1412 MatConvert_MPISBAIJ_Basic, 1413 NULL, 1414 NULL, 1415 /* 74*/ NULL, 1416 NULL, 1417 NULL, 1418 NULL, 1419 NULL, 1420 /* 79*/ NULL, 1421 NULL, 1422 NULL, 1423 MatGetInertia_SeqSBAIJ, 1424 MatLoad_SeqSBAIJ, 1425 /* 84*/ MatIsSymmetric_SeqSBAIJ, 1426 MatIsHermitian_SeqSBAIJ, 1427 MatIsStructurallySymmetric_SeqSBAIJ, 1428 NULL, 1429 NULL, 1430 /* 89*/ NULL, 1431 NULL, 1432 NULL, 1433 NULL, 1434 NULL, 1435 /* 94*/ NULL, 1436 NULL, 1437 NULL, 1438 NULL, 1439 NULL, 1440 /* 99*/ NULL, 1441 NULL, 1442 NULL, 1443 MatConjugate_SeqSBAIJ, 1444 NULL, 1445 /*104*/ NULL, 1446 MatRealPart_SeqSBAIJ, 1447 MatImaginaryPart_SeqSBAIJ, 1448 MatGetRowUpperTriangular_SeqSBAIJ, 1449 MatRestoreRowUpperTriangular_SeqSBAIJ, 1450 /*109*/ NULL, 1451 NULL, 1452 NULL, 1453 NULL, 1454 MatMissingDiagonal_SeqSBAIJ, 1455 /*114*/ NULL, 1456 NULL, 1457 NULL, 1458 NULL, 1459 NULL, 1460 /*119*/ NULL, 1461 NULL, 1462 NULL, 1463 NULL, 1464 NULL, 1465 /*124*/ NULL, 1466 NULL, 1467 NULL, 1468 NULL, 1469 NULL, 1470 /*129*/ NULL, 1471 NULL, 1472 NULL, 1473 NULL, 1474 NULL, 1475 /*134*/ NULL, 1476 NULL, 1477 NULL, 1478 NULL, 1479 NULL, 1480 /*139*/ MatSetBlockSizes_Default, 1481 NULL, 1482 NULL, 1483 NULL, 1484 NULL, 1485 /*144*/ MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ, 1486 NULL, 1487 NULL, 1488 NULL, 1489 NULL, 1490 NULL, 1491 /*150*/ NULL, 1492 NULL}; 1493 1494 PetscErrorCode MatStoreValues_SeqSBAIJ(Mat mat) 1495 { 1496 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data; 1497 PetscInt nz = aij->i[mat->rmap->N] * mat->rmap->bs * aij->bs2; 1498 1499 PetscFunctionBegin; 1500 PetscCheck(aij->nonew == 1, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1501 1502 /* allocate space for values if not already there */ 1503 if (!aij->saved_values) PetscCall(PetscMalloc1(nz + 1, &aij->saved_values)); 1504 1505 /* copy values over */ 1506 PetscCall(PetscArraycpy(aij->saved_values, aij->a, nz)); 1507 PetscFunctionReturn(PETSC_SUCCESS); 1508 } 1509 1510 PetscErrorCode MatRetrieveValues_SeqSBAIJ(Mat mat) 1511 { 1512 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data; 1513 PetscInt nz = aij->i[mat->rmap->N] * mat->rmap->bs * aij->bs2; 1514 1515 PetscFunctionBegin; 1516 PetscCheck(aij->nonew == 1, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1517 PetscCheck(aij->saved_values, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call MatStoreValues(A);first"); 1518 1519 /* copy values over */ 1520 PetscCall(PetscArraycpy(aij->a, aij->saved_values, nz)); 1521 PetscFunctionReturn(PETSC_SUCCESS); 1522 } 1523 1524 static PetscErrorCode MatSeqSBAIJSetPreallocation_SeqSBAIJ(Mat B, PetscInt bs, PetscInt nz, PetscInt *nnz) 1525 { 1526 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ *)B->data; 1527 PetscInt i, mbs, nbs, bs2; 1528 PetscBool skipallocation = PETSC_FALSE, flg = PETSC_FALSE, realalloc = PETSC_FALSE; 1529 1530 PetscFunctionBegin; 1531 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 1532 1533 PetscCall(MatSetBlockSize(B, PetscAbs(bs))); 1534 PetscCall(PetscLayoutSetUp(B->rmap)); 1535 PetscCall(PetscLayoutSetUp(B->cmap)); 1536 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); 1537 PetscCall(PetscLayoutGetBlockSize(B->rmap, &bs)); 1538 1539 B->preallocated = PETSC_TRUE; 1540 1541 mbs = B->rmap->N / bs; 1542 nbs = B->cmap->n / bs; 1543 bs2 = bs * bs; 1544 1545 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"); 1546 1547 if (nz == MAT_SKIP_ALLOCATION) { 1548 skipallocation = PETSC_TRUE; 1549 nz = 0; 1550 } 1551 1552 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 3; 1553 PetscCheck(nz >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "nz cannot be less than 0: value %" PetscInt_FMT, nz); 1554 if (nnz) { 1555 for (i = 0; i < mbs; i++) { 1556 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]); 1557 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); 1558 } 1559 } 1560 1561 B->ops->mult = MatMult_SeqSBAIJ_N; 1562 B->ops->multadd = MatMultAdd_SeqSBAIJ_N; 1563 B->ops->multtranspose = MatMult_SeqSBAIJ_N; 1564 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_N; 1565 1566 PetscCall(PetscOptionsGetBool(((PetscObject)B)->options, ((PetscObject)B)->prefix, "-mat_no_unroll", &flg, NULL)); 1567 if (!flg) { 1568 switch (bs) { 1569 case 1: 1570 B->ops->mult = MatMult_SeqSBAIJ_1; 1571 B->ops->multadd = MatMultAdd_SeqSBAIJ_1; 1572 B->ops->multtranspose = MatMult_SeqSBAIJ_1; 1573 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_1; 1574 break; 1575 case 2: 1576 B->ops->mult = MatMult_SeqSBAIJ_2; 1577 B->ops->multadd = MatMultAdd_SeqSBAIJ_2; 1578 B->ops->multtranspose = MatMult_SeqSBAIJ_2; 1579 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_2; 1580 break; 1581 case 3: 1582 B->ops->mult = MatMult_SeqSBAIJ_3; 1583 B->ops->multadd = MatMultAdd_SeqSBAIJ_3; 1584 B->ops->multtranspose = MatMult_SeqSBAIJ_3; 1585 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_3; 1586 break; 1587 case 4: 1588 B->ops->mult = MatMult_SeqSBAIJ_4; 1589 B->ops->multadd = MatMultAdd_SeqSBAIJ_4; 1590 B->ops->multtranspose = MatMult_SeqSBAIJ_4; 1591 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_4; 1592 break; 1593 case 5: 1594 B->ops->mult = MatMult_SeqSBAIJ_5; 1595 B->ops->multadd = MatMultAdd_SeqSBAIJ_5; 1596 B->ops->multtranspose = MatMult_SeqSBAIJ_5; 1597 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_5; 1598 break; 1599 case 6: 1600 B->ops->mult = MatMult_SeqSBAIJ_6; 1601 B->ops->multadd = MatMultAdd_SeqSBAIJ_6; 1602 B->ops->multtranspose = MatMult_SeqSBAIJ_6; 1603 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_6; 1604 break; 1605 case 7: 1606 B->ops->mult = MatMult_SeqSBAIJ_7; 1607 B->ops->multadd = MatMultAdd_SeqSBAIJ_7; 1608 B->ops->multtranspose = MatMult_SeqSBAIJ_7; 1609 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_7; 1610 break; 1611 } 1612 } 1613 1614 b->mbs = mbs; 1615 b->nbs = nbs; 1616 if (!skipallocation) { 1617 if (!b->imax) { 1618 PetscCall(PetscMalloc2(mbs, &b->imax, mbs, &b->ilen)); 1619 1620 b->free_imax_ilen = PETSC_TRUE; 1621 } 1622 if (!nnz) { 1623 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 1624 else if (nz <= 0) nz = 1; 1625 nz = PetscMin(nbs, nz); 1626 for (i = 0; i < mbs; i++) b->imax[i] = nz; 1627 PetscCall(PetscIntMultError(nz, mbs, &nz)); 1628 } else { 1629 PetscInt64 nz64 = 0; 1630 for (i = 0; i < mbs; i++) { 1631 b->imax[i] = nnz[i]; 1632 nz64 += nnz[i]; 1633 } 1634 PetscCall(PetscIntCast(nz64, &nz)); 1635 } 1636 /* b->ilen will count nonzeros in each block row so far. */ 1637 for (i = 0; i < mbs; i++) b->ilen[i] = 0; 1638 /* nz=(nz+mbs)/2; */ /* total diagonal and superdiagonal nonzero blocks */ 1639 1640 /* allocate the matrix space */ 1641 PetscCall(MatSeqXAIJFreeAIJ(B, &b->a, &b->j, &b->i)); 1642 PetscCall(PetscMalloc3(bs2 * nz, &b->a, nz, &b->j, B->rmap->N + 1, &b->i)); 1643 PetscCall(PetscArrayzero(b->a, nz * bs2)); 1644 PetscCall(PetscArrayzero(b->j, nz)); 1645 1646 b->singlemalloc = PETSC_TRUE; 1647 1648 /* pointer to beginning of each row */ 1649 b->i[0] = 0; 1650 for (i = 1; i < mbs + 1; i++) b->i[i] = b->i[i - 1] + b->imax[i - 1]; 1651 1652 b->free_a = PETSC_TRUE; 1653 b->free_ij = PETSC_TRUE; 1654 } else { 1655 b->free_a = PETSC_FALSE; 1656 b->free_ij = PETSC_FALSE; 1657 } 1658 1659 b->bs2 = bs2; 1660 b->nz = 0; 1661 b->maxnz = nz; 1662 b->inew = NULL; 1663 b->jnew = NULL; 1664 b->anew = NULL; 1665 b->a2anew = NULL; 1666 b->permute = PETSC_FALSE; 1667 1668 B->was_assembled = PETSC_FALSE; 1669 B->assembled = PETSC_FALSE; 1670 if (realalloc) PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE)); 1671 PetscFunctionReturn(PETSC_SUCCESS); 1672 } 1673 1674 PetscErrorCode MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ(Mat B, PetscInt bs, const PetscInt ii[], const PetscInt jj[], const PetscScalar V[]) 1675 { 1676 PetscInt i, j, m, nz, anz, nz_max = 0, *nnz; 1677 PetscScalar *values = NULL; 1678 PetscBool roworiented = ((Mat_SeqSBAIJ *)B->data)->roworiented; 1679 1680 PetscFunctionBegin; 1681 PetscCheck(bs >= 1, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "Invalid block size specified, must be positive but it is %" PetscInt_FMT, bs); 1682 PetscCall(PetscLayoutSetBlockSize(B->rmap, bs)); 1683 PetscCall(PetscLayoutSetBlockSize(B->cmap, bs)); 1684 PetscCall(PetscLayoutSetUp(B->rmap)); 1685 PetscCall(PetscLayoutSetUp(B->cmap)); 1686 PetscCall(PetscLayoutGetBlockSize(B->rmap, &bs)); 1687 m = B->rmap->n / bs; 1688 1689 PetscCheck(!ii[0], PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "ii[0] must be 0 but it is %" PetscInt_FMT, ii[0]); 1690 PetscCall(PetscMalloc1(m + 1, &nnz)); 1691 for (i = 0; i < m; i++) { 1692 nz = ii[i + 1] - ii[i]; 1693 PetscCheck(nz >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row %" PetscInt_FMT " has a negative number of columns %" PetscInt_FMT, i, nz); 1694 anz = 0; 1695 for (j = 0; j < nz; j++) { 1696 /* count only values on the diagonal or above */ 1697 if (jj[ii[i] + j] >= i) { 1698 anz = nz - j; 1699 break; 1700 } 1701 } 1702 nz_max = PetscMax(nz_max, anz); 1703 nnz[i] = anz; 1704 } 1705 PetscCall(MatSeqSBAIJSetPreallocation(B, bs, 0, nnz)); 1706 PetscCall(PetscFree(nnz)); 1707 1708 values = (PetscScalar *)V; 1709 if (!values) PetscCall(PetscCalloc1(bs * bs * nz_max, &values)); 1710 for (i = 0; i < m; i++) { 1711 PetscInt ncols = ii[i + 1] - ii[i]; 1712 const PetscInt *icols = jj + ii[i]; 1713 if (!roworiented || bs == 1) { 1714 const PetscScalar *svals = values + (V ? (bs * bs * ii[i]) : 0); 1715 PetscCall(MatSetValuesBlocked_SeqSBAIJ(B, 1, &i, ncols, icols, svals, INSERT_VALUES)); 1716 } else { 1717 for (j = 0; j < ncols; j++) { 1718 const PetscScalar *svals = values + (V ? (bs * bs * (ii[i] + j)) : 0); 1719 PetscCall(MatSetValuesBlocked_SeqSBAIJ(B, 1, &i, 1, &icols[j], svals, INSERT_VALUES)); 1720 } 1721 } 1722 } 1723 if (!V) PetscCall(PetscFree(values)); 1724 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 1725 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 1726 PetscCall(MatSetOption(B, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE)); 1727 PetscFunctionReturn(PETSC_SUCCESS); 1728 } 1729 1730 /* 1731 This is used to set the numeric factorization for both Cholesky and ICC symbolic factorization 1732 */ 1733 PetscErrorCode MatSeqSBAIJSetNumericFactorization_inplace(Mat B, PetscBool natural) 1734 { 1735 PetscBool flg = PETSC_FALSE; 1736 PetscInt bs = B->rmap->bs; 1737 1738 PetscFunctionBegin; 1739 PetscCall(PetscOptionsGetBool(((PetscObject)B)->options, ((PetscObject)B)->prefix, "-mat_no_unroll", &flg, NULL)); 1740 if (flg) bs = 8; 1741 1742 if (!natural) { 1743 switch (bs) { 1744 case 1: 1745 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_inplace; 1746 break; 1747 case 2: 1748 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2; 1749 break; 1750 case 3: 1751 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3; 1752 break; 1753 case 4: 1754 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4; 1755 break; 1756 case 5: 1757 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5; 1758 break; 1759 case 6: 1760 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6; 1761 break; 1762 case 7: 1763 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7; 1764 break; 1765 default: 1766 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N; 1767 break; 1768 } 1769 } else { 1770 switch (bs) { 1771 case 1: 1772 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering_inplace; 1773 break; 1774 case 2: 1775 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering; 1776 break; 1777 case 3: 1778 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3_NaturalOrdering; 1779 break; 1780 case 4: 1781 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4_NaturalOrdering; 1782 break; 1783 case 5: 1784 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5_NaturalOrdering; 1785 break; 1786 case 6: 1787 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6_NaturalOrdering; 1788 break; 1789 case 7: 1790 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7_NaturalOrdering; 1791 break; 1792 default: 1793 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N_NaturalOrdering; 1794 break; 1795 } 1796 } 1797 PetscFunctionReturn(PETSC_SUCCESS); 1798 } 1799 1800 PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqAIJ(Mat, MatType, MatReuse, Mat *); 1801 PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat, MatType, MatReuse, Mat *); 1802 static PetscErrorCode MatFactorGetSolverType_petsc(Mat A, MatSolverType *type) 1803 { 1804 PetscFunctionBegin; 1805 *type = MATSOLVERPETSC; 1806 PetscFunctionReturn(PETSC_SUCCESS); 1807 } 1808 1809 PETSC_INTERN PetscErrorCode MatGetFactor_seqsbaij_petsc(Mat A, MatFactorType ftype, Mat *B) 1810 { 1811 PetscInt n = A->rmap->n; 1812 1813 PetscFunctionBegin; 1814 #if defined(PETSC_USE_COMPLEX) 1815 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"); 1816 #endif 1817 1818 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), B)); 1819 PetscCall(MatSetSizes(*B, n, n, n, n)); 1820 if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) { 1821 PetscCall(MatSetType(*B, MATSEQSBAIJ)); 1822 PetscCall(MatSeqSBAIJSetPreallocation(*B, A->rmap->bs, MAT_SKIP_ALLOCATION, NULL)); 1823 1824 (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJ; 1825 (*B)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqSBAIJ; 1826 PetscCall(PetscStrallocpy(MATORDERINGNATURAL, (char **)&(*B)->preferredordering[MAT_FACTOR_CHOLESKY])); 1827 PetscCall(PetscStrallocpy(MATORDERINGNATURAL, (char **)&(*B)->preferredordering[MAT_FACTOR_ICC])); 1828 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Factor type not supported"); 1829 1830 (*B)->factortype = ftype; 1831 (*B)->canuseordering = PETSC_TRUE; 1832 PetscCall(PetscFree((*B)->solvertype)); 1833 PetscCall(PetscStrallocpy(MATSOLVERPETSC, &(*B)->solvertype)); 1834 PetscCall(PetscObjectComposeFunction((PetscObject)*B, "MatFactorGetSolverType_C", MatFactorGetSolverType_petsc)); 1835 PetscFunctionReturn(PETSC_SUCCESS); 1836 } 1837 1838 /*@C 1839 MatSeqSBAIJGetArray - gives access to the array where the data for a `MATSEQSBAIJ` matrix is stored 1840 1841 Not Collective 1842 1843 Input Parameter: 1844 . mat - a `MATSEQSBAIJ` matrix 1845 1846 Output Parameter: 1847 . array - pointer to the data 1848 1849 Level: intermediate 1850 1851 .seealso: `MATSEQSBAIJ`, `MatSeqSBAIJRestoreArray()`, `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArray()` 1852 @*/ 1853 PetscErrorCode MatSeqSBAIJGetArray(Mat A, PetscScalar **array) 1854 { 1855 PetscFunctionBegin; 1856 PetscUseMethod(A, "MatSeqSBAIJGetArray_C", (Mat, PetscScalar **), (A, array)); 1857 PetscFunctionReturn(PETSC_SUCCESS); 1858 } 1859 1860 /*@C 1861 MatSeqSBAIJRestoreArray - returns access to the array where the data for a `MATSEQSBAIJ` matrix is stored obtained by `MatSeqSBAIJGetArray()` 1862 1863 Not Collective 1864 1865 Input Parameters: 1866 + mat - a MATSEQSBAIJ matrix 1867 - array - pointer to the data 1868 1869 Level: intermediate 1870 1871 .seealso: `MATSEQSBAIJ`, `MatSeqSBAIJGetArray()`, `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArray()` 1872 @*/ 1873 PetscErrorCode MatSeqSBAIJRestoreArray(Mat A, PetscScalar **array) 1874 { 1875 PetscFunctionBegin; 1876 PetscUseMethod(A, "MatSeqSBAIJRestoreArray_C", (Mat, PetscScalar **), (A, array)); 1877 PetscFunctionReturn(PETSC_SUCCESS); 1878 } 1879 1880 /*MC 1881 MATSEQSBAIJ - MATSEQSBAIJ = "seqsbaij" - A matrix type to be used for sequential symmetric block sparse matrices, 1882 based on block compressed sparse row format. Only the upper triangular portion of the matrix is stored. 1883 1884 For complex numbers by default this matrix is symmetric, NOT Hermitian symmetric. To make it Hermitian symmetric you 1885 can call `MatSetOption`(`Mat`, `MAT_HERMITIAN`). 1886 1887 Options Database Keys: 1888 . -mat_type seqsbaij - sets the matrix type to "seqsbaij" during a call to `MatSetFromOptions()` 1889 1890 Notes: 1891 By default if you insert values into the lower triangular part of the matrix they are simply ignored (since they are not 1892 stored and it is assumed they symmetric to the upper triangular). If you call `MatSetOption`(`Mat`,`MAT_IGNORE_LOWER_TRIANGULAR`,`PETSC_FALSE`) or use 1893 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. 1894 1895 The number of rows in the matrix must be less than or equal to the number of columns 1896 1897 Level: beginner 1898 1899 .seealso: `MATSEQSBAIJ`, `MatCreateSeqSBAIJ()`, `MatType`, `MATMPISBAIJ` 1900 M*/ 1901 PETSC_EXTERN PetscErrorCode MatCreate_SeqSBAIJ(Mat B) 1902 { 1903 Mat_SeqSBAIJ *b; 1904 PetscMPIInt size; 1905 PetscBool no_unroll = PETSC_FALSE, no_inode = PETSC_FALSE; 1906 1907 PetscFunctionBegin; 1908 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B), &size)); 1909 PetscCheck(size <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Comm must be of size 1"); 1910 1911 PetscCall(PetscNew(&b)); 1912 B->data = (void *)b; 1913 PetscCall(PetscMemcpy(B->ops, &MatOps_Values, sizeof(struct _MatOps))); 1914 1915 B->ops->destroy = MatDestroy_SeqSBAIJ; 1916 B->ops->view = MatView_SeqSBAIJ; 1917 b->row = NULL; 1918 b->icol = NULL; 1919 b->reallocs = 0; 1920 b->saved_values = NULL; 1921 b->inode.limit = 5; 1922 b->inode.max_limit = 5; 1923 1924 b->roworiented = PETSC_TRUE; 1925 b->nonew = 0; 1926 b->diag = NULL; 1927 b->solve_work = NULL; 1928 b->mult_work = NULL; 1929 B->spptr = NULL; 1930 B->info.nz_unneeded = (PetscReal)b->maxnz * b->bs2; 1931 b->keepnonzeropattern = PETSC_FALSE; 1932 1933 b->inew = NULL; 1934 b->jnew = NULL; 1935 b->anew = NULL; 1936 b->a2anew = NULL; 1937 b->permute = PETSC_FALSE; 1938 1939 b->ignore_ltriangular = PETSC_TRUE; 1940 1941 PetscCall(PetscOptionsGetBool(((PetscObject)B)->options, ((PetscObject)B)->prefix, "-mat_ignore_lower_triangular", &b->ignore_ltriangular, NULL)); 1942 1943 b->getrow_utriangular = PETSC_FALSE; 1944 1945 PetscCall(PetscOptionsGetBool(((PetscObject)B)->options, ((PetscObject)B)->prefix, "-mat_getrow_uppertriangular", &b->getrow_utriangular, NULL)); 1946 1947 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqSBAIJGetArray_C", MatSeqSBAIJGetArray_SeqSBAIJ)); 1948 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqSBAIJRestoreArray_C", MatSeqSBAIJRestoreArray_SeqSBAIJ)); 1949 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatStoreValues_C", MatStoreValues_SeqSBAIJ)); 1950 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatRetrieveValues_C", MatRetrieveValues_SeqSBAIJ)); 1951 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqSBAIJSetColumnIndices_C", MatSeqSBAIJSetColumnIndices_SeqSBAIJ)); 1952 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqsbaij_seqaij_C", MatConvert_SeqSBAIJ_SeqAIJ)); 1953 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqsbaij_seqbaij_C", MatConvert_SeqSBAIJ_SeqBAIJ)); 1954 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqSBAIJSetPreallocation_C", MatSeqSBAIJSetPreallocation_SeqSBAIJ)); 1955 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqSBAIJSetPreallocationCSR_C", MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ)); 1956 #if defined(PETSC_HAVE_ELEMENTAL) 1957 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqsbaij_elemental_C", MatConvert_SeqSBAIJ_Elemental)); 1958 #endif 1959 #if defined(PETSC_HAVE_SCALAPACK) 1960 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqsbaij_scalapack_C", MatConvert_SBAIJ_ScaLAPACK)); 1961 #endif 1962 1963 B->symmetry_eternal = PETSC_TRUE; 1964 B->structural_symmetry_eternal = PETSC_TRUE; 1965 B->symmetric = PETSC_BOOL3_TRUE; 1966 B->structurally_symmetric = PETSC_BOOL3_TRUE; 1967 #if defined(PETSC_USE_COMPLEX) 1968 B->hermitian = PETSC_BOOL3_FALSE; 1969 #else 1970 B->hermitian = PETSC_BOOL3_TRUE; 1971 #endif 1972 1973 PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATSEQSBAIJ)); 1974 1975 PetscOptionsBegin(PetscObjectComm((PetscObject)B), ((PetscObject)B)->prefix, "Options for SEQSBAIJ matrix", "Mat"); 1976 PetscCall(PetscOptionsBool("-mat_no_unroll", "Do not optimize for inodes (slower)", NULL, no_unroll, &no_unroll, NULL)); 1977 if (no_unroll) PetscCall(PetscInfo(B, "Not using Inode routines due to -mat_no_unroll\n")); 1978 PetscCall(PetscOptionsBool("-mat_no_inode", "Do not optimize for inodes (slower)", NULL, no_inode, &no_inode, NULL)); 1979 if (no_inode) PetscCall(PetscInfo(B, "Not using Inode routines due to -mat_no_inode\n")); 1980 PetscCall(PetscOptionsInt("-mat_inode_limit", "Do not use inodes larger then this value", NULL, b->inode.limit, &b->inode.limit, NULL)); 1981 PetscOptionsEnd(); 1982 b->inode.use = (PetscBool)(!(no_unroll || no_inode)); 1983 if (b->inode.limit > b->inode.max_limit) b->inode.limit = b->inode.max_limit; 1984 PetscFunctionReturn(PETSC_SUCCESS); 1985 } 1986 1987 /*@C 1988 MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block 1989 compressed row) `MATSEQSBAIJ` format. For good matrix assembly performance the 1990 user should preallocate the matrix storage by setting the parameter nz 1991 (or the array nnz). By setting these parameters accurately, performance 1992 during matrix assembly can be increased by more than a factor of 50. 1993 1994 Collective 1995 1996 Input Parameters: 1997 + B - the symmetric matrix 1998 . bs - size of block, the blocks are ALWAYS square. One can use `MatSetBlockSizes()` to set a different row and column blocksize but the row 1999 blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with `MatCreateVecs()` 2000 . nz - number of block nonzeros per block row (same for all rows) 2001 - nnz - array containing the number of block nonzeros in the upper triangular plus 2002 diagonal portion of each block (possibly different for each block row) or NULL 2003 2004 Options Database Keys: 2005 + -mat_no_unroll - uses code that does not unroll the loops in the 2006 block calculations (much slower) 2007 - -mat_block_size - size of the blocks to use (only works if a negative bs is passed in 2008 2009 Level: intermediate 2010 2011 Notes: 2012 Specify the preallocated storage with either nz or nnz (not both). 2013 Set nz = `PETSC_DEFAULT` and nnz = NULL for PETSc to control dynamic memory 2014 allocation. See [Sparse Matrices](sec_matsparse) for details. 2015 2016 You can call `MatGetInfo()` to get information on how effective the preallocation was; 2017 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 2018 You can also run with the option -info and look for messages with the string 2019 malloc in them to see if additional memory allocation was needed. 2020 2021 If the nnz parameter is given then the nz parameter is ignored 2022 2023 .seealso: [Sparse Matrices](sec_matsparse), `MATSEQSBAIJ`, `MatCreate()`, `MatCreateSeqAIJ()`, `MatSetValues()`, `MatCreateSBAIJ()` 2024 @*/ 2025 PetscErrorCode MatSeqSBAIJSetPreallocation(Mat B, PetscInt bs, PetscInt nz, const PetscInt nnz[]) 2026 { 2027 PetscFunctionBegin; 2028 PetscValidHeaderSpecific(B, MAT_CLASSID, 1); 2029 PetscValidType(B, 1); 2030 PetscValidLogicalCollectiveInt(B, bs, 2); 2031 PetscTryMethod(B, "MatSeqSBAIJSetPreallocation_C", (Mat, PetscInt, PetscInt, const PetscInt[]), (B, bs, nz, nnz)); 2032 PetscFunctionReturn(PETSC_SUCCESS); 2033 } 2034 2035 /*@C 2036 MatSeqSBAIJSetPreallocationCSR - Creates a sparse parallel matrix in `MATSEQSBAIJ` format using the given nonzero structure and (optional) numerical values 2037 2038 Input Parameters: 2039 + B - the matrix 2040 . bs - size of block, the blocks are ALWAYS square. 2041 . i - the indices into j for the start of each local row (starts with zero) 2042 . j - the column indices for each local row (starts with zero) these must be sorted for each row 2043 - v - optional values in the matrix 2044 2045 Level: advanced 2046 2047 Notes: 2048 The order of the entries in values is specified by the `MatOption` `MAT_ROW_ORIENTED`. For example, C programs 2049 may want to use the default `MAT_ROW_ORIENTED` = `PETSC_TRUE` and use an array v[nnz][bs][bs] where the second index is 2050 over rows within a block and the last index is over columns within a block row. Fortran programs will likely set 2051 `MAT_ROW_ORIENTED` = `PETSC_FALSE` and use a Fortran array v(bs,bs,nnz) in which the first index is over rows within a 2052 block column and the second index is over columns within a block. 2053 2054 Any entries below the diagonal are ignored 2055 2056 Though this routine has Preallocation() in the name it also sets the exact nonzero locations of the matrix entries 2057 and usually the numerical values as well 2058 2059 .seealso: `MATSEQSBAIJ`, `MatCreate()`, `MatCreateSeqSBAIJ()`, `MatSetValuesBlocked()`, `MatSeqSBAIJSetPreallocation()`, `MATSEQSBAIJ` 2060 @*/ 2061 PetscErrorCode MatSeqSBAIJSetPreallocationCSR(Mat B, PetscInt bs, const PetscInt i[], const PetscInt j[], const PetscScalar v[]) 2062 { 2063 PetscFunctionBegin; 2064 PetscValidHeaderSpecific(B, MAT_CLASSID, 1); 2065 PetscValidType(B, 1); 2066 PetscValidLogicalCollectiveInt(B, bs, 2); 2067 PetscTryMethod(B, "MatSeqSBAIJSetPreallocationCSR_C", (Mat, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[]), (B, bs, i, j, v)); 2068 PetscFunctionReturn(PETSC_SUCCESS); 2069 } 2070 2071 /*@C 2072 MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block 2073 compressed row) `MATSEQSBAIJ` format. For good matrix assembly performance the 2074 user should preallocate the matrix storage by setting the parameter nz 2075 (or the array nnz). By setting these parameters accurately, performance 2076 during matrix assembly can be increased by more than a factor of 50. 2077 2078 Collective 2079 2080 Input Parameters: 2081 + comm - MPI communicator, set to `PETSC_COMM_SELF` 2082 . bs - size of block, the blocks are ALWAYS square. One can use `MatSetBlockSizes()` to set a different row and column blocksize but the row 2083 blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs() 2084 . m - number of rows, or number of columns 2085 . nz - number of block nonzeros per block row (same for all rows) 2086 - nnz - array containing the number of block nonzeros in the upper triangular plus 2087 diagonal portion of each block (possibly different for each block row) or NULL 2088 2089 Output Parameter: 2090 . A - the symmetric matrix 2091 2092 Options Database Keys: 2093 + -mat_no_unroll - uses code that does not unroll the loops in the 2094 block calculations (much slower) 2095 - -mat_block_size - size of the blocks to use 2096 2097 Level: intermediate 2098 2099 It is recommended that one use the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`, 2100 MatXXXXSetPreallocation() paradigm instead of this routine directly. 2101 [MatXXXXSetPreallocation() is, for example, `MatSeqAIJSetPreallocation()`] 2102 2103 Notes: 2104 The number of rows and columns must be divisible by blocksize. 2105 This matrix type does not support complex Hermitian operation. 2106 2107 Specify the preallocated storage with either nz or nnz (not both). 2108 Set nz = `PETSC_DEFAULT` and nnz = NULL for PETSc to control dynamic memory 2109 allocation. See [Sparse Matrices](sec_matsparse) for details. 2110 2111 If the nnz parameter is given then the nz parameter is ignored 2112 2113 .seealso: [Sparse Matrices](sec_matsparse), `MATSEQSBAIJ`, `MatCreate()`, `MatCreateSeqAIJ()`, `MatSetValues()`, `MatCreateSBAIJ()` 2114 @*/ 2115 PetscErrorCode MatCreateSeqSBAIJ(MPI_Comm comm, PetscInt bs, PetscInt m, PetscInt n, PetscInt nz, const PetscInt nnz[], Mat *A) 2116 { 2117 PetscFunctionBegin; 2118 PetscCall(MatCreate(comm, A)); 2119 PetscCall(MatSetSizes(*A, m, n, m, n)); 2120 PetscCall(MatSetType(*A, MATSEQSBAIJ)); 2121 PetscCall(MatSeqSBAIJSetPreallocation(*A, bs, nz, (PetscInt *)nnz)); 2122 PetscFunctionReturn(PETSC_SUCCESS); 2123 } 2124 2125 PetscErrorCode MatDuplicate_SeqSBAIJ(Mat A, MatDuplicateOption cpvalues, Mat *B) 2126 { 2127 Mat C; 2128 Mat_SeqSBAIJ *c, *a = (Mat_SeqSBAIJ *)A->data; 2129 PetscInt i, mbs = a->mbs, nz = a->nz, bs2 = a->bs2; 2130 2131 PetscFunctionBegin; 2132 PetscCheck(a->i[mbs] == nz, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Corrupt matrix"); 2133 2134 *B = NULL; 2135 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &C)); 2136 PetscCall(MatSetSizes(C, A->rmap->N, A->cmap->n, A->rmap->N, A->cmap->n)); 2137 PetscCall(MatSetBlockSizesFromMats(C, A, A)); 2138 PetscCall(MatSetType(C, MATSEQSBAIJ)); 2139 c = (Mat_SeqSBAIJ *)C->data; 2140 2141 C->preallocated = PETSC_TRUE; 2142 C->factortype = A->factortype; 2143 c->row = NULL; 2144 c->icol = NULL; 2145 c->saved_values = NULL; 2146 c->keepnonzeropattern = a->keepnonzeropattern; 2147 C->assembled = PETSC_TRUE; 2148 2149 PetscCall(PetscLayoutReference(A->rmap, &C->rmap)); 2150 PetscCall(PetscLayoutReference(A->cmap, &C->cmap)); 2151 c->bs2 = a->bs2; 2152 c->mbs = a->mbs; 2153 c->nbs = a->nbs; 2154 2155 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2156 c->imax = a->imax; 2157 c->ilen = a->ilen; 2158 c->free_imax_ilen = PETSC_FALSE; 2159 } else { 2160 PetscCall(PetscMalloc2((mbs + 1), &c->imax, (mbs + 1), &c->ilen)); 2161 for (i = 0; i < mbs; i++) { 2162 c->imax[i] = a->imax[i]; 2163 c->ilen[i] = a->ilen[i]; 2164 } 2165 c->free_imax_ilen = PETSC_TRUE; 2166 } 2167 2168 /* allocate the matrix space */ 2169 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2170 PetscCall(PetscMalloc1(bs2 * nz, &c->a)); 2171 c->i = a->i; 2172 c->j = a->j; 2173 c->singlemalloc = PETSC_FALSE; 2174 c->free_a = PETSC_TRUE; 2175 c->free_ij = PETSC_FALSE; 2176 c->parent = A; 2177 PetscCall(PetscObjectReference((PetscObject)A)); 2178 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE)); 2179 PetscCall(MatSetOption(C, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE)); 2180 } else { 2181 PetscCall(PetscMalloc3(bs2 * nz, &c->a, nz, &c->j, mbs + 1, &c->i)); 2182 PetscCall(PetscArraycpy(c->i, a->i, mbs + 1)); 2183 c->singlemalloc = PETSC_TRUE; 2184 c->free_a = PETSC_TRUE; 2185 c->free_ij = PETSC_TRUE; 2186 } 2187 if (mbs > 0) { 2188 if (cpvalues != MAT_SHARE_NONZERO_PATTERN) PetscCall(PetscArraycpy(c->j, a->j, nz)); 2189 if (cpvalues == MAT_COPY_VALUES) { 2190 PetscCall(PetscArraycpy(c->a, a->a, bs2 * nz)); 2191 } else { 2192 PetscCall(PetscArrayzero(c->a, bs2 * nz)); 2193 } 2194 if (a->jshort) { 2195 /* cannot share jshort, it is reallocated in MatAssemblyEnd_SeqSBAIJ() */ 2196 /* if the parent matrix is reassembled, this child matrix will never notice */ 2197 PetscCall(PetscMalloc1(nz, &c->jshort)); 2198 PetscCall(PetscArraycpy(c->jshort, a->jshort, nz)); 2199 2200 c->free_jshort = PETSC_TRUE; 2201 } 2202 } 2203 2204 c->roworiented = a->roworiented; 2205 c->nonew = a->nonew; 2206 2207 if (a->diag) { 2208 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2209 c->diag = a->diag; 2210 c->free_diag = PETSC_FALSE; 2211 } else { 2212 PetscCall(PetscMalloc1(mbs, &c->diag)); 2213 for (i = 0; i < mbs; i++) c->diag[i] = a->diag[i]; 2214 c->free_diag = PETSC_TRUE; 2215 } 2216 } 2217 c->nz = a->nz; 2218 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 2219 c->solve_work = NULL; 2220 c->mult_work = NULL; 2221 2222 *B = C; 2223 PetscCall(PetscFunctionListDuplicate(((PetscObject)A)->qlist, &((PetscObject)C)->qlist)); 2224 PetscFunctionReturn(PETSC_SUCCESS); 2225 } 2226 2227 /* Used for both SeqBAIJ and SeqSBAIJ matrices */ 2228 #define MatLoad_SeqSBAIJ_Binary MatLoad_SeqBAIJ_Binary 2229 2230 PetscErrorCode MatLoad_SeqSBAIJ(Mat mat, PetscViewer viewer) 2231 { 2232 PetscBool isbinary; 2233 2234 PetscFunctionBegin; 2235 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary)); 2236 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); 2237 PetscCall(MatLoad_SeqSBAIJ_Binary(mat, viewer)); 2238 PetscFunctionReturn(PETSC_SUCCESS); 2239 } 2240 2241 /*@ 2242 MatCreateSeqSBAIJWithArrays - Creates an sequential `MATSEQSBAIJ` matrix using matrix elements 2243 (upper triangular entries in CSR format) provided by the user. 2244 2245 Collective 2246 2247 Input Parameters: 2248 + comm - must be an MPI communicator of size 1 2249 . bs - size of block 2250 . m - number of rows 2251 . n - number of columns 2252 . 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 2253 . j - column indices 2254 - a - matrix values 2255 2256 Output Parameter: 2257 . mat - the matrix 2258 2259 Level: advanced 2260 2261 Notes: 2262 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2263 once the matrix is destroyed 2264 2265 You cannot set new nonzero locations into this matrix, that will generate an error. 2266 2267 The i and j indices are 0 based 2268 2269 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 2270 it is the regular CSR format excluding the lower triangular elements. 2271 2272 .seealso: `MATSEQSBAIJ`, `MatCreate()`, `MatCreateSBAIJ()`, `MatCreateSeqSBAIJ()` 2273 @*/ 2274 PetscErrorCode MatCreateSeqSBAIJWithArrays(MPI_Comm comm, PetscInt bs, PetscInt m, PetscInt n, PetscInt i[], PetscInt j[], PetscScalar a[], Mat *mat) 2275 { 2276 PetscInt ii; 2277 Mat_SeqSBAIJ *sbaij; 2278 2279 PetscFunctionBegin; 2280 PetscCheck(bs == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "block size %" PetscInt_FMT " > 1 is not supported yet", bs); 2281 PetscCheck(m == 0 || i[0] == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "i (row indices) must start with 0"); 2282 2283 PetscCall(MatCreate(comm, mat)); 2284 PetscCall(MatSetSizes(*mat, m, n, m, n)); 2285 PetscCall(MatSetType(*mat, MATSEQSBAIJ)); 2286 PetscCall(MatSeqSBAIJSetPreallocation(*mat, bs, MAT_SKIP_ALLOCATION, NULL)); 2287 sbaij = (Mat_SeqSBAIJ *)(*mat)->data; 2288 PetscCall(PetscMalloc2(m, &sbaij->imax, m, &sbaij->ilen)); 2289 2290 sbaij->i = i; 2291 sbaij->j = j; 2292 sbaij->a = a; 2293 2294 sbaij->singlemalloc = PETSC_FALSE; 2295 sbaij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2296 sbaij->free_a = PETSC_FALSE; 2297 sbaij->free_ij = PETSC_FALSE; 2298 sbaij->free_imax_ilen = PETSC_TRUE; 2299 2300 for (ii = 0; ii < m; ii++) { 2301 sbaij->ilen[ii] = sbaij->imax[ii] = i[ii + 1] - i[ii]; 2302 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]); 2303 } 2304 if (PetscDefined(USE_DEBUG)) { 2305 for (ii = 0; ii < sbaij->i[m]; ii++) { 2306 PetscCheck(j[ii] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative column index at location = %" PetscInt_FMT " index = %" PetscInt_FMT, ii, j[ii]); 2307 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]); 2308 } 2309 } 2310 2311 PetscCall(MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY)); 2312 PetscCall(MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY)); 2313 PetscFunctionReturn(PETSC_SUCCESS); 2314 } 2315 2316 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ(MPI_Comm comm, Mat inmat, PetscInt n, MatReuse scall, Mat *outmat) 2317 { 2318 PetscFunctionBegin; 2319 PetscCall(MatCreateMPIMatConcatenateSeqMat_MPISBAIJ(comm, inmat, n, scall, outmat)); 2320 PetscFunctionReturn(PETSC_SUCCESS); 2321 } 2322