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