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