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