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