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