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