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 #if defined(PETSC_USE_COMPLEX) 1031 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data; 1032 PetscInt i, nz = a->bs2 * a->i[a->mbs]; 1033 MatScalar *aa = a->a; 1034 1035 PetscFunctionBegin; 1036 for (i = 0; i < nz; i++) aa[i] = PetscConj(aa[i]); 1037 #else 1038 PetscFunctionBegin; 1039 #endif 1040 PetscFunctionReturn(PETSC_SUCCESS); 1041 } 1042 1043 static PetscErrorCode MatRealPart_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] = PetscRealPart(aa[i]); 1051 PetscFunctionReturn(PETSC_SUCCESS); 1052 } 1053 1054 static PetscErrorCode MatImaginaryPart_SeqSBAIJ(Mat A) 1055 { 1056 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data; 1057 PetscInt i, nz = a->bs2 * a->i[a->mbs]; 1058 MatScalar *aa = a->a; 1059 1060 PetscFunctionBegin; 1061 for (i = 0; i < nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1062 PetscFunctionReturn(PETSC_SUCCESS); 1063 } 1064 1065 static PetscErrorCode MatZeroRowsColumns_SeqSBAIJ(Mat A, PetscInt is_n, const PetscInt is_idx[], PetscScalar diag, Vec x, Vec b) 1066 { 1067 Mat_SeqSBAIJ *baij = (Mat_SeqSBAIJ *)A->data; 1068 PetscInt i, j, k, count; 1069 PetscInt bs = A->rmap->bs, bs2 = baij->bs2, row, col; 1070 PetscScalar zero = 0.0; 1071 MatScalar *aa; 1072 const PetscScalar *xx; 1073 PetscScalar *bb; 1074 PetscBool *zeroed, vecs = PETSC_FALSE; 1075 1076 PetscFunctionBegin; 1077 /* fix right-hand side if needed */ 1078 if (x && b) { 1079 PetscCall(VecGetArrayRead(x, &xx)); 1080 PetscCall(VecGetArray(b, &bb)); 1081 vecs = PETSC_TRUE; 1082 } 1083 1084 /* zero the columns */ 1085 PetscCall(PetscCalloc1(A->rmap->n, &zeroed)); 1086 for (i = 0; i < is_n; i++) { 1087 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]); 1088 zeroed[is_idx[i]] = PETSC_TRUE; 1089 } 1090 if (vecs) { 1091 for (i = 0; i < A->rmap->N; i++) { 1092 row = i / bs; 1093 for (j = baij->i[row]; j < baij->i[row + 1]; j++) { 1094 for (k = 0; k < bs; k++) { 1095 col = bs * baij->j[j] + k; 1096 if (col <= i) continue; 1097 aa = baij->a + j * bs2 + (i % bs) + bs * k; 1098 if (!zeroed[i] && zeroed[col]) bb[i] -= aa[0] * xx[col]; 1099 if (zeroed[i] && !zeroed[col]) bb[col] -= aa[0] * xx[i]; 1100 } 1101 } 1102 } 1103 for (i = 0; i < is_n; i++) bb[is_idx[i]] = diag * xx[is_idx[i]]; 1104 } 1105 1106 for (i = 0; i < A->rmap->N; i++) { 1107 if (!zeroed[i]) { 1108 row = i / bs; 1109 for (j = baij->i[row]; j < baij->i[row + 1]; j++) { 1110 for (k = 0; k < bs; k++) { 1111 col = bs * baij->j[j] + k; 1112 if (zeroed[col]) { 1113 aa = baij->a + j * bs2 + (i % bs) + bs * k; 1114 aa[0] = 0.0; 1115 } 1116 } 1117 } 1118 } 1119 } 1120 PetscCall(PetscFree(zeroed)); 1121 if (vecs) { 1122 PetscCall(VecRestoreArrayRead(x, &xx)); 1123 PetscCall(VecRestoreArray(b, &bb)); 1124 } 1125 1126 /* zero the rows */ 1127 for (i = 0; i < is_n; i++) { 1128 row = is_idx[i]; 1129 count = (baij->i[row / bs + 1] - baij->i[row / bs]) * bs; 1130 aa = baij->a + baij->i[row / bs] * bs2 + (row % bs); 1131 for (k = 0; k < count; k++) { 1132 aa[0] = zero; 1133 aa += bs; 1134 } 1135 if (diag != 0.0) PetscUseTypeMethod(A, setvalues, 1, &row, 1, &row, &diag, INSERT_VALUES); 1136 } 1137 PetscCall(MatAssemblyEnd_SeqSBAIJ(A, MAT_FINAL_ASSEMBLY)); 1138 PetscFunctionReturn(PETSC_SUCCESS); 1139 } 1140 1141 static PetscErrorCode MatShift_SeqSBAIJ(Mat Y, PetscScalar a) 1142 { 1143 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)Y->data; 1144 1145 PetscFunctionBegin; 1146 if (!Y->preallocated || !aij->nz) PetscCall(MatSeqSBAIJSetPreallocation(Y, Y->rmap->bs, 1, NULL)); 1147 PetscCall(MatShift_Basic(Y, a)); 1148 PetscFunctionReturn(PETSC_SUCCESS); 1149 } 1150 1151 PetscErrorCode MatEliminateZeros_SeqSBAIJ(Mat A, PetscBool keep) 1152 { 1153 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data; 1154 PetscInt fshift = 0, fshift_prev = 0, i, *ai = a->i, *aj = a->j, *imax = a->imax, j, k; 1155 PetscInt m = A->rmap->N, *ailen = a->ilen; 1156 PetscInt mbs = a->mbs, bs2 = a->bs2, rmax = 0; 1157 MatScalar *aa = a->a, *ap; 1158 PetscBool zero; 1159 1160 PetscFunctionBegin; 1161 PetscCheck(A->assembled, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cannot eliminate zeros for unassembled matrix"); 1162 if (m) rmax = ailen[0]; 1163 for (i = 1; i <= mbs; i++) { 1164 for (k = ai[i - 1]; k < ai[i]; k++) { 1165 zero = PETSC_TRUE; 1166 ap = aa + bs2 * k; 1167 for (j = 0; j < bs2 && zero; j++) { 1168 if (ap[j] != 0.0) zero = PETSC_FALSE; 1169 } 1170 if (zero && (aj[k] != i - 1 || !keep)) fshift++; 1171 else { 1172 if (zero && aj[k] == i - 1) PetscCall(PetscInfo(A, "Keep the diagonal block at row %" PetscInt_FMT "\n", i - 1)); 1173 aj[k - fshift] = aj[k]; 1174 PetscCall(PetscArraymove(ap - bs2 * fshift, ap, bs2)); 1175 } 1176 } 1177 ai[i - 1] -= fshift_prev; 1178 fshift_prev = fshift; 1179 ailen[i - 1] = imax[i - 1] = ai[i] - fshift - ai[i - 1]; 1180 a->nonzerorowcnt += ((ai[i] - fshift - ai[i - 1]) > 0); 1181 rmax = PetscMax(rmax, ailen[i - 1]); 1182 } 1183 if (fshift) { 1184 if (mbs) { 1185 ai[mbs] -= fshift; 1186 a->nz = ai[mbs]; 1187 } 1188 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)); 1189 A->nonzerostate++; 1190 A->info.nz_unneeded += (PetscReal)fshift; 1191 a->rmax = rmax; 1192 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 1193 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 1194 } 1195 PetscFunctionReturn(PETSC_SUCCESS); 1196 } 1197 1198 static struct _MatOps MatOps_Values = {MatSetValues_SeqSBAIJ, 1199 MatGetRow_SeqSBAIJ, 1200 MatRestoreRow_SeqSBAIJ, 1201 MatMult_SeqSBAIJ_N, 1202 /* 4*/ MatMultAdd_SeqSBAIJ_N, 1203 MatMult_SeqSBAIJ_N, /* transpose versions are same as non-transpose versions */ 1204 MatMultAdd_SeqSBAIJ_N, 1205 NULL, 1206 NULL, 1207 NULL, 1208 /* 10*/ NULL, 1209 NULL, 1210 MatCholeskyFactor_SeqSBAIJ, 1211 MatSOR_SeqSBAIJ, 1212 MatTranspose_SeqSBAIJ, 1213 /* 15*/ MatGetInfo_SeqSBAIJ, 1214 MatEqual_SeqSBAIJ, 1215 MatGetDiagonal_SeqSBAIJ, 1216 MatDiagonalScale_SeqSBAIJ, 1217 MatNorm_SeqSBAIJ, 1218 /* 20*/ NULL, 1219 MatAssemblyEnd_SeqSBAIJ, 1220 MatSetOption_SeqSBAIJ, 1221 MatZeroEntries_SeqSBAIJ, 1222 /* 24*/ NULL, 1223 NULL, 1224 NULL, 1225 NULL, 1226 NULL, 1227 /* 29*/ MatSetUp_Seq_Hash, 1228 NULL, 1229 NULL, 1230 NULL, 1231 NULL, 1232 /* 34*/ MatDuplicate_SeqSBAIJ, 1233 NULL, 1234 NULL, 1235 NULL, 1236 MatICCFactor_SeqSBAIJ, 1237 /* 39*/ MatAXPY_SeqSBAIJ, 1238 MatCreateSubMatrices_SeqSBAIJ, 1239 MatIncreaseOverlap_SeqSBAIJ, 1240 MatGetValues_SeqSBAIJ, 1241 MatCopy_SeqSBAIJ, 1242 /* 44*/ NULL, 1243 MatScale_SeqSBAIJ, 1244 MatShift_SeqSBAIJ, 1245 NULL, 1246 MatZeroRowsColumns_SeqSBAIJ, 1247 /* 49*/ NULL, 1248 MatGetRowIJ_SeqSBAIJ, 1249 MatRestoreRowIJ_SeqSBAIJ, 1250 NULL, 1251 NULL, 1252 /* 54*/ NULL, 1253 NULL, 1254 NULL, 1255 MatPermute_SeqSBAIJ, 1256 MatSetValuesBlocked_SeqSBAIJ, 1257 /* 59*/ MatCreateSubMatrix_SeqSBAIJ, 1258 NULL, 1259 NULL, 1260 NULL, 1261 NULL, 1262 /* 64*/ NULL, 1263 NULL, 1264 NULL, 1265 NULL, 1266 MatGetRowMaxAbs_SeqSBAIJ, 1267 /* 69*/ NULL, 1268 MatConvert_MPISBAIJ_Basic, 1269 NULL, 1270 NULL, 1271 NULL, 1272 /* 74*/ NULL, 1273 NULL, 1274 NULL, 1275 MatGetInertia_SeqSBAIJ, 1276 MatLoad_SeqSBAIJ, 1277 /* 79*/ NULL, 1278 NULL, 1279 MatIsStructurallySymmetric_SeqSBAIJ, 1280 NULL, 1281 NULL, 1282 /* 84*/ NULL, 1283 NULL, 1284 NULL, 1285 NULL, 1286 NULL, 1287 /* 89*/ NULL, 1288 NULL, 1289 NULL, 1290 NULL, 1291 MatConjugate_SeqSBAIJ, 1292 /* 94*/ NULL, 1293 NULL, 1294 MatRealPart_SeqSBAIJ, 1295 MatImaginaryPart_SeqSBAIJ, 1296 MatGetRowUpperTriangular_SeqSBAIJ, 1297 /* 99*/ MatRestoreRowUpperTriangular_SeqSBAIJ, 1298 NULL, 1299 NULL, 1300 NULL, 1301 NULL, 1302 /*104*/ NULL, 1303 NULL, 1304 NULL, 1305 NULL, 1306 NULL, 1307 /*109*/ NULL, 1308 NULL, 1309 NULL, 1310 NULL, 1311 NULL, 1312 /*114*/ NULL, 1313 NULL, 1314 NULL, 1315 NULL, 1316 NULL, 1317 /*119*/ NULL, 1318 NULL, 1319 NULL, 1320 NULL, 1321 NULL, 1322 /*124*/ NULL, 1323 NULL, 1324 MatSetBlockSizes_Default, 1325 NULL, 1326 NULL, 1327 /*129*/ NULL, 1328 MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ, 1329 NULL, 1330 NULL, 1331 NULL, 1332 /*134*/ NULL, 1333 NULL, 1334 MatEliminateZeros_SeqSBAIJ, 1335 NULL, 1336 NULL, 1337 /*139*/ NULL, 1338 NULL, 1339 MatCopyHashToXAIJ_Seq_Hash, 1340 NULL, 1341 NULL}; 1342 1343 static PetscErrorCode MatStoreValues_SeqSBAIJ(Mat mat) 1344 { 1345 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data; 1346 PetscInt nz = aij->i[mat->rmap->N] * mat->rmap->bs * aij->bs2; 1347 1348 PetscFunctionBegin; 1349 PetscCheck(aij->nonew == 1, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1350 1351 /* allocate space for values if not already there */ 1352 if (!aij->saved_values) PetscCall(PetscMalloc1(nz + 1, &aij->saved_values)); 1353 1354 /* copy values over */ 1355 PetscCall(PetscArraycpy(aij->saved_values, aij->a, nz)); 1356 PetscFunctionReturn(PETSC_SUCCESS); 1357 } 1358 1359 static PetscErrorCode MatRetrieveValues_SeqSBAIJ(Mat mat) 1360 { 1361 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data; 1362 PetscInt nz = aij->i[mat->rmap->N] * mat->rmap->bs * aij->bs2; 1363 1364 PetscFunctionBegin; 1365 PetscCheck(aij->nonew == 1, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1366 PetscCheck(aij->saved_values, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call MatStoreValues(A);first"); 1367 1368 /* copy values over */ 1369 PetscCall(PetscArraycpy(aij->a, aij->saved_values, nz)); 1370 PetscFunctionReturn(PETSC_SUCCESS); 1371 } 1372 1373 static PetscErrorCode MatSeqSBAIJSetPreallocation_SeqSBAIJ(Mat B, PetscInt bs, PetscInt nz, const PetscInt nnz[]) 1374 { 1375 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ *)B->data; 1376 PetscInt i, mbs, nbs, bs2; 1377 PetscBool skipallocation = PETSC_FALSE, flg = PETSC_FALSE, realalloc = PETSC_FALSE; 1378 1379 PetscFunctionBegin; 1380 if (B->hash_active) { 1381 PetscInt bs; 1382 B->ops[0] = b->cops; 1383 PetscCall(PetscHMapIJVDestroy(&b->ht)); 1384 PetscCall(MatGetBlockSize(B, &bs)); 1385 if (bs > 1) PetscCall(PetscHSetIJDestroy(&b->bht)); 1386 PetscCall(PetscFree(b->dnz)); 1387 PetscCall(PetscFree(b->bdnz)); 1388 B->hash_active = PETSC_FALSE; 1389 } 1390 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 1391 1392 PetscCall(MatSetBlockSize(B, bs)); 1393 PetscCall(PetscLayoutSetUp(B->rmap)); 1394 PetscCall(PetscLayoutSetUp(B->cmap)); 1395 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); 1396 PetscCall(PetscLayoutGetBlockSize(B->rmap, &bs)); 1397 1398 B->preallocated = PETSC_TRUE; 1399 1400 mbs = B->rmap->N / bs; 1401 nbs = B->cmap->n / bs; 1402 bs2 = bs * bs; 1403 1404 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"); 1405 1406 if (nz == MAT_SKIP_ALLOCATION) { 1407 skipallocation = PETSC_TRUE; 1408 nz = 0; 1409 } 1410 1411 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 3; 1412 PetscCheck(nz >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "nz cannot be less than 0: value %" PetscInt_FMT, nz); 1413 if (nnz) { 1414 for (i = 0; i < mbs; i++) { 1415 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]); 1416 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); 1417 } 1418 } 1419 1420 B->ops->mult = MatMult_SeqSBAIJ_N; 1421 B->ops->multadd = MatMultAdd_SeqSBAIJ_N; 1422 B->ops->multtranspose = MatMult_SeqSBAIJ_N; 1423 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_N; 1424 1425 PetscCall(PetscOptionsGetBool(((PetscObject)B)->options, ((PetscObject)B)->prefix, "-mat_no_unroll", &flg, NULL)); 1426 if (!flg) { 1427 switch (bs) { 1428 case 1: 1429 B->ops->mult = MatMult_SeqSBAIJ_1; 1430 B->ops->multadd = MatMultAdd_SeqSBAIJ_1; 1431 B->ops->multtranspose = MatMult_SeqSBAIJ_1; 1432 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_1; 1433 break; 1434 case 2: 1435 B->ops->mult = MatMult_SeqSBAIJ_2; 1436 B->ops->multadd = MatMultAdd_SeqSBAIJ_2; 1437 B->ops->multtranspose = MatMult_SeqSBAIJ_2; 1438 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_2; 1439 break; 1440 case 3: 1441 B->ops->mult = MatMult_SeqSBAIJ_3; 1442 B->ops->multadd = MatMultAdd_SeqSBAIJ_3; 1443 B->ops->multtranspose = MatMult_SeqSBAIJ_3; 1444 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_3; 1445 break; 1446 case 4: 1447 B->ops->mult = MatMult_SeqSBAIJ_4; 1448 B->ops->multadd = MatMultAdd_SeqSBAIJ_4; 1449 B->ops->multtranspose = MatMult_SeqSBAIJ_4; 1450 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_4; 1451 break; 1452 case 5: 1453 B->ops->mult = MatMult_SeqSBAIJ_5; 1454 B->ops->multadd = MatMultAdd_SeqSBAIJ_5; 1455 B->ops->multtranspose = MatMult_SeqSBAIJ_5; 1456 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_5; 1457 break; 1458 case 6: 1459 B->ops->mult = MatMult_SeqSBAIJ_6; 1460 B->ops->multadd = MatMultAdd_SeqSBAIJ_6; 1461 B->ops->multtranspose = MatMult_SeqSBAIJ_6; 1462 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_6; 1463 break; 1464 case 7: 1465 B->ops->mult = MatMult_SeqSBAIJ_7; 1466 B->ops->multadd = MatMultAdd_SeqSBAIJ_7; 1467 B->ops->multtranspose = MatMult_SeqSBAIJ_7; 1468 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_7; 1469 break; 1470 } 1471 } 1472 1473 b->mbs = mbs; 1474 b->nbs = nbs; 1475 if (!skipallocation) { 1476 if (!b->imax) { 1477 PetscCall(PetscMalloc2(mbs, &b->imax, mbs, &b->ilen)); 1478 1479 b->free_imax_ilen = PETSC_TRUE; 1480 } 1481 if (!nnz) { 1482 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 1483 else if (nz <= 0) nz = 1; 1484 nz = PetscMin(nbs, nz); 1485 for (i = 0; i < mbs; i++) b->imax[i] = nz; 1486 PetscCall(PetscIntMultError(nz, mbs, &nz)); 1487 } else { 1488 PetscInt64 nz64 = 0; 1489 for (i = 0; i < mbs; i++) { 1490 b->imax[i] = nnz[i]; 1491 nz64 += nnz[i]; 1492 } 1493 PetscCall(PetscIntCast(nz64, &nz)); 1494 } 1495 /* b->ilen will count nonzeros in each block row so far. */ 1496 for (i = 0; i < mbs; i++) b->ilen[i] = 0; 1497 /* nz=(nz+mbs)/2; */ /* total diagonal and superdiagonal nonzero blocks */ 1498 1499 /* allocate the matrix space */ 1500 PetscCall(MatSeqXAIJFreeAIJ(B, &b->a, &b->j, &b->i)); 1501 PetscCall(PetscShmgetAllocateArray(bs2 * nz, sizeof(PetscScalar), (void **)&b->a)); 1502 PetscCall(PetscShmgetAllocateArray(nz, sizeof(PetscInt), (void **)&b->j)); 1503 PetscCall(PetscShmgetAllocateArray(B->rmap->n + 1, sizeof(PetscInt), (void **)&b->i)); 1504 b->free_a = PETSC_TRUE; 1505 b->free_ij = PETSC_TRUE; 1506 PetscCall(PetscArrayzero(b->a, nz * bs2)); 1507 PetscCall(PetscArrayzero(b->j, nz)); 1508 b->free_a = PETSC_TRUE; 1509 b->free_ij = PETSC_TRUE; 1510 1511 /* pointer to beginning of each row */ 1512 b->i[0] = 0; 1513 for (i = 1; i < mbs + 1; i++) b->i[i] = b->i[i - 1] + b->imax[i - 1]; 1514 1515 } else { 1516 b->free_a = PETSC_FALSE; 1517 b->free_ij = PETSC_FALSE; 1518 } 1519 1520 b->bs2 = bs2; 1521 b->nz = 0; 1522 b->maxnz = nz; 1523 b->inew = NULL; 1524 b->jnew = NULL; 1525 b->anew = NULL; 1526 b->a2anew = NULL; 1527 b->permute = PETSC_FALSE; 1528 1529 B->was_assembled = PETSC_FALSE; 1530 B->assembled = PETSC_FALSE; 1531 if (realalloc) PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE)); 1532 PetscFunctionReturn(PETSC_SUCCESS); 1533 } 1534 1535 static PetscErrorCode MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ(Mat B, PetscInt bs, const PetscInt ii[], const PetscInt jj[], const PetscScalar V[]) 1536 { 1537 PetscInt i, j, m, nz, anz, nz_max = 0, *nnz; 1538 PetscScalar *values = NULL; 1539 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ *)B->data; 1540 PetscBool roworiented = b->roworiented; 1541 PetscBool ilw = b->ignore_ltriangular; 1542 1543 PetscFunctionBegin; 1544 PetscCheck(bs >= 1, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "Invalid block size specified, must be positive but it is %" PetscInt_FMT, bs); 1545 PetscCall(PetscLayoutSetBlockSize(B->rmap, bs)); 1546 PetscCall(PetscLayoutSetBlockSize(B->cmap, bs)); 1547 PetscCall(PetscLayoutSetUp(B->rmap)); 1548 PetscCall(PetscLayoutSetUp(B->cmap)); 1549 PetscCall(PetscLayoutGetBlockSize(B->rmap, &bs)); 1550 m = B->rmap->n / bs; 1551 1552 PetscCheck(!ii[0], PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "ii[0] must be 0 but it is %" PetscInt_FMT, ii[0]); 1553 PetscCall(PetscMalloc1(m + 1, &nnz)); 1554 for (i = 0; i < m; i++) { 1555 nz = ii[i + 1] - ii[i]; 1556 PetscCheck(nz >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row %" PetscInt_FMT " has a negative number of columns %" PetscInt_FMT, i, nz); 1557 PetscCheckSorted(nz, jj + ii[i]); 1558 anz = 0; 1559 for (j = 0; j < nz; j++) { 1560 /* count only values on the diagonal or above */ 1561 if (jj[ii[i] + j] >= i) { 1562 anz = nz - j; 1563 break; 1564 } 1565 } 1566 nz_max = PetscMax(nz_max, nz); 1567 nnz[i] = anz; 1568 } 1569 PetscCall(MatSeqSBAIJSetPreallocation(B, bs, 0, nnz)); 1570 PetscCall(PetscFree(nnz)); 1571 1572 values = (PetscScalar *)V; 1573 if (!values) PetscCall(PetscCalloc1(bs * bs * nz_max, &values)); 1574 b->ignore_ltriangular = PETSC_TRUE; 1575 for (i = 0; i < m; i++) { 1576 PetscInt ncols = ii[i + 1] - ii[i]; 1577 const PetscInt *icols = jj + ii[i]; 1578 1579 if (!roworiented || bs == 1) { 1580 const PetscScalar *svals = values + (V ? (bs * bs * ii[i]) : 0); 1581 PetscCall(MatSetValuesBlocked_SeqSBAIJ(B, 1, &i, ncols, icols, svals, INSERT_VALUES)); 1582 } else { 1583 for (j = 0; j < ncols; j++) { 1584 const PetscScalar *svals = values + (V ? (bs * bs * (ii[i] + j)) : 0); 1585 PetscCall(MatSetValuesBlocked_SeqSBAIJ(B, 1, &i, 1, &icols[j], svals, INSERT_VALUES)); 1586 } 1587 } 1588 } 1589 if (!V) PetscCall(PetscFree(values)); 1590 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 1591 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 1592 PetscCall(MatSetOption(B, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE)); 1593 b->ignore_ltriangular = ilw; 1594 PetscFunctionReturn(PETSC_SUCCESS); 1595 } 1596 1597 /* 1598 This is used to set the numeric factorization for both Cholesky and ICC symbolic factorization 1599 */ 1600 PetscErrorCode MatSeqSBAIJSetNumericFactorization_inplace(Mat B, PetscBool natural) 1601 { 1602 PetscBool flg = PETSC_FALSE; 1603 PetscInt bs = B->rmap->bs; 1604 1605 PetscFunctionBegin; 1606 PetscCall(PetscOptionsGetBool(((PetscObject)B)->options, ((PetscObject)B)->prefix, "-mat_no_unroll", &flg, NULL)); 1607 if (flg) bs = 8; 1608 1609 if (!natural) { 1610 switch (bs) { 1611 case 1: 1612 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_inplace; 1613 break; 1614 case 2: 1615 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2; 1616 break; 1617 case 3: 1618 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3; 1619 break; 1620 case 4: 1621 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4; 1622 break; 1623 case 5: 1624 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5; 1625 break; 1626 case 6: 1627 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6; 1628 break; 1629 case 7: 1630 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7; 1631 break; 1632 default: 1633 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N; 1634 break; 1635 } 1636 } else { 1637 switch (bs) { 1638 case 1: 1639 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering_inplace; 1640 break; 1641 case 2: 1642 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering; 1643 break; 1644 case 3: 1645 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3_NaturalOrdering; 1646 break; 1647 case 4: 1648 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4_NaturalOrdering; 1649 break; 1650 case 5: 1651 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5_NaturalOrdering; 1652 break; 1653 case 6: 1654 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6_NaturalOrdering; 1655 break; 1656 case 7: 1657 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7_NaturalOrdering; 1658 break; 1659 default: 1660 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N_NaturalOrdering; 1661 break; 1662 } 1663 } 1664 PetscFunctionReturn(PETSC_SUCCESS); 1665 } 1666 1667 PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqAIJ(Mat, MatType, MatReuse, Mat *); 1668 PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat, MatType, MatReuse, Mat *); 1669 static PetscErrorCode MatFactorGetSolverType_petsc(Mat A, MatSolverType *type) 1670 { 1671 PetscFunctionBegin; 1672 *type = MATSOLVERPETSC; 1673 PetscFunctionReturn(PETSC_SUCCESS); 1674 } 1675 1676 PETSC_INTERN PetscErrorCode MatGetFactor_seqsbaij_petsc(Mat A, MatFactorType ftype, Mat *B) 1677 { 1678 PetscInt n = A->rmap->n; 1679 1680 PetscFunctionBegin; 1681 #if defined(PETSC_USE_COMPLEX) 1682 if ((ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) && A->hermitian == PETSC_BOOL3_TRUE && A->symmetric != PETSC_BOOL3_TRUE) { 1683 PetscCall(PetscInfo(A, "Hermitian MAT_FACTOR_CHOLESKY or MAT_FACTOR_ICC are not supported. Use MAT_FACTOR_LU instead.\n")); 1684 *B = NULL; 1685 PetscFunctionReturn(PETSC_SUCCESS); 1686 } 1687 #endif 1688 1689 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), B)); 1690 PetscCall(MatSetSizes(*B, n, n, n, n)); 1691 PetscCheck(ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC, PETSC_COMM_SELF, PETSC_ERR_SUP, "Factor type not supported"); 1692 PetscCall(MatSetType(*B, MATSEQSBAIJ)); 1693 PetscCall(MatSeqSBAIJSetPreallocation(*B, A->rmap->bs, MAT_SKIP_ALLOCATION, NULL)); 1694 1695 (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJ; 1696 (*B)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqSBAIJ; 1697 PetscCall(PetscStrallocpy(MATORDERINGNATURAL, (char **)&(*B)->preferredordering[MAT_FACTOR_CHOLESKY])); 1698 PetscCall(PetscStrallocpy(MATORDERINGNATURAL, (char **)&(*B)->preferredordering[MAT_FACTOR_ICC])); 1699 1700 (*B)->factortype = ftype; 1701 (*B)->canuseordering = PETSC_TRUE; 1702 PetscCall(PetscFree((*B)->solvertype)); 1703 PetscCall(PetscStrallocpy(MATSOLVERPETSC, &(*B)->solvertype)); 1704 PetscCall(PetscObjectComposeFunction((PetscObject)*B, "MatFactorGetSolverType_C", MatFactorGetSolverType_petsc)); 1705 PetscFunctionReturn(PETSC_SUCCESS); 1706 } 1707 1708 /*@C 1709 MatSeqSBAIJGetArray - gives access to the array where the numerical data for a `MATSEQSBAIJ` matrix is stored 1710 1711 Not Collective 1712 1713 Input Parameter: 1714 . A - a `MATSEQSBAIJ` matrix 1715 1716 Output Parameter: 1717 . array - pointer to the data 1718 1719 Level: intermediate 1720 1721 .seealso: [](ch_matrices), `Mat`, `MATSEQSBAIJ`, `MatSeqSBAIJRestoreArray()`, `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArray()` 1722 @*/ 1723 PetscErrorCode MatSeqSBAIJGetArray(Mat A, PetscScalar *array[]) 1724 { 1725 PetscFunctionBegin; 1726 PetscUseMethod(A, "MatSeqSBAIJGetArray_C", (Mat, PetscScalar **), (A, array)); 1727 PetscFunctionReturn(PETSC_SUCCESS); 1728 } 1729 1730 /*@C 1731 MatSeqSBAIJRestoreArray - returns access to the array where the numerical data for a `MATSEQSBAIJ` matrix is stored obtained by `MatSeqSBAIJGetArray()` 1732 1733 Not Collective 1734 1735 Input Parameters: 1736 + A - a `MATSEQSBAIJ` matrix 1737 - array - pointer to the data 1738 1739 Level: intermediate 1740 1741 .seealso: [](ch_matrices), `Mat`, `MATSEQSBAIJ`, `MatSeqSBAIJGetArray()`, `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArray()` 1742 @*/ 1743 PetscErrorCode MatSeqSBAIJRestoreArray(Mat A, PetscScalar *array[]) 1744 { 1745 PetscFunctionBegin; 1746 PetscUseMethod(A, "MatSeqSBAIJRestoreArray_C", (Mat, PetscScalar **), (A, array)); 1747 PetscFunctionReturn(PETSC_SUCCESS); 1748 } 1749 1750 /*MC 1751 MATSEQSBAIJ - MATSEQSBAIJ = "seqsbaij" - A matrix type to be used for sequential symmetric block sparse matrices, 1752 based on block compressed sparse row format. Only the upper triangular portion of the matrix is stored. 1753 1754 For complex numbers by default this matrix is symmetric, NOT Hermitian symmetric. To make it Hermitian symmetric you 1755 can call `MatSetOption`(`Mat`, `MAT_HERMITIAN`). 1756 1757 Options Database Key: 1758 . -mat_type seqsbaij - sets the matrix type to "seqsbaij" during a call to `MatSetFromOptions()` 1759 1760 Level: beginner 1761 1762 Notes: 1763 By default if you insert values into the lower triangular part of the matrix they are simply ignored (since they are not 1764 stored and it is assumed they symmetric to the upper triangular). If you call `MatSetOption`(`Mat`,`MAT_IGNORE_LOWER_TRIANGULAR`,`PETSC_FALSE`) or use 1765 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. 1766 1767 The number of rows in the matrix must be less than or equal to the number of columns 1768 1769 .seealso: [](ch_matrices), `Mat`, `MATSEQSBAIJ`, `MatCreateSeqSBAIJ()`, `MatType`, `MATMPISBAIJ` 1770 M*/ 1771 PETSC_EXTERN PetscErrorCode MatCreate_SeqSBAIJ(Mat B) 1772 { 1773 Mat_SeqSBAIJ *b; 1774 PetscMPIInt size; 1775 PetscBool no_unroll = PETSC_FALSE, no_inode = PETSC_FALSE; 1776 1777 PetscFunctionBegin; 1778 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B), &size)); 1779 PetscCheck(size <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Comm must be of size 1"); 1780 1781 PetscCall(PetscNew(&b)); 1782 B->data = (void *)b; 1783 B->ops[0] = MatOps_Values; 1784 1785 B->ops->destroy = MatDestroy_SeqSBAIJ; 1786 B->ops->view = MatView_SeqSBAIJ; 1787 b->row = NULL; 1788 b->icol = NULL; 1789 b->reallocs = 0; 1790 b->saved_values = NULL; 1791 b->inode.limit = 5; 1792 b->inode.max_limit = 5; 1793 1794 b->roworiented = PETSC_TRUE; 1795 b->nonew = 0; 1796 b->diag = NULL; 1797 b->solve_work = NULL; 1798 b->mult_work = NULL; 1799 B->spptr = NULL; 1800 B->info.nz_unneeded = (PetscReal)b->maxnz * b->bs2; 1801 b->keepnonzeropattern = PETSC_FALSE; 1802 1803 b->inew = NULL; 1804 b->jnew = NULL; 1805 b->anew = NULL; 1806 b->a2anew = NULL; 1807 b->permute = PETSC_FALSE; 1808 1809 b->ignore_ltriangular = PETSC_TRUE; 1810 1811 PetscCall(PetscOptionsGetBool(((PetscObject)B)->options, ((PetscObject)B)->prefix, "-mat_ignore_lower_triangular", &b->ignore_ltriangular, NULL)); 1812 1813 b->getrow_utriangular = PETSC_FALSE; 1814 1815 PetscCall(PetscOptionsGetBool(((PetscObject)B)->options, ((PetscObject)B)->prefix, "-mat_getrow_uppertriangular", &b->getrow_utriangular, NULL)); 1816 1817 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqSBAIJGetArray_C", MatSeqSBAIJGetArray_SeqSBAIJ)); 1818 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqSBAIJRestoreArray_C", MatSeqSBAIJRestoreArray_SeqSBAIJ)); 1819 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatStoreValues_C", MatStoreValues_SeqSBAIJ)); 1820 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatRetrieveValues_C", MatRetrieveValues_SeqSBAIJ)); 1821 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqSBAIJSetColumnIndices_C", MatSeqSBAIJSetColumnIndices_SeqSBAIJ)); 1822 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqsbaij_seqaij_C", MatConvert_SeqSBAIJ_SeqAIJ)); 1823 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqsbaij_seqbaij_C", MatConvert_SeqSBAIJ_SeqBAIJ)); 1824 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqSBAIJSetPreallocation_C", MatSeqSBAIJSetPreallocation_SeqSBAIJ)); 1825 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqSBAIJSetPreallocationCSR_C", MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ)); 1826 #if defined(PETSC_HAVE_ELEMENTAL) 1827 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqsbaij_elemental_C", MatConvert_SeqSBAIJ_Elemental)); 1828 #endif 1829 #if defined(PETSC_HAVE_SCALAPACK) && (defined(PETSC_USE_REAL_SINGLE) || defined(PETSC_USE_REAL_DOUBLE)) 1830 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqsbaij_scalapack_C", MatConvert_SBAIJ_ScaLAPACK)); 1831 #endif 1832 1833 B->symmetry_eternal = PETSC_TRUE; 1834 B->structural_symmetry_eternal = PETSC_TRUE; 1835 B->symmetric = PETSC_BOOL3_TRUE; 1836 B->structurally_symmetric = PETSC_BOOL3_TRUE; 1837 #if !defined(PETSC_USE_COMPLEX) 1838 B->hermitian = PETSC_BOOL3_TRUE; 1839 #endif 1840 1841 PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATSEQSBAIJ)); 1842 1843 PetscOptionsBegin(PetscObjectComm((PetscObject)B), ((PetscObject)B)->prefix, "Options for SEQSBAIJ matrix", "Mat"); 1844 PetscCall(PetscOptionsBool("-mat_no_unroll", "Do not optimize for inodes (slower)", NULL, no_unroll, &no_unroll, NULL)); 1845 if (no_unroll) PetscCall(PetscInfo(B, "Not using Inode routines due to -mat_no_unroll\n")); 1846 PetscCall(PetscOptionsBool("-mat_no_inode", "Do not optimize for inodes (slower)", NULL, no_inode, &no_inode, NULL)); 1847 if (no_inode) PetscCall(PetscInfo(B, "Not using Inode routines due to -mat_no_inode\n")); 1848 PetscCall(PetscOptionsInt("-mat_inode_limit", "Do not use inodes larger than this value", NULL, b->inode.limit, &b->inode.limit, NULL)); 1849 PetscOptionsEnd(); 1850 b->inode.use = (PetscBool)(!(no_unroll || no_inode)); 1851 if (b->inode.limit > b->inode.max_limit) b->inode.limit = b->inode.max_limit; 1852 PetscFunctionReturn(PETSC_SUCCESS); 1853 } 1854 1855 /*@ 1856 MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block 1857 compressed row) `MATSEQSBAIJ` format. For good matrix assembly performance the 1858 user should preallocate the matrix storage by setting the parameter `nz` 1859 (or the array `nnz`). 1860 1861 Collective 1862 1863 Input Parameters: 1864 + B - the symmetric matrix 1865 . bs - size of block, the blocks are ALWAYS square. One can use `MatSetBlockSizes()` to set a different row and column blocksize but the row 1866 blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with `MatCreateVecs()` 1867 . nz - number of block nonzeros per block row (same for all rows) 1868 - nnz - array containing the number of block nonzeros in the upper triangular plus 1869 diagonal portion of each block (possibly different for each block row) or `NULL` 1870 1871 Options Database Keys: 1872 + -mat_no_unroll - uses code that does not unroll the loops in the block calculations (much slower) 1873 - -mat_block_size - size of the blocks to use (only works if a negative bs is passed in 1874 1875 Level: intermediate 1876 1877 Notes: 1878 Specify the preallocated storage with either `nz` or `nnz` (not both). 1879 Set `nz` = `PETSC_DEFAULT` and `nnz` = `NULL` for PETSc to control dynamic memory 1880 allocation. See [Sparse Matrices](sec_matsparse) for details. 1881 1882 You can call `MatGetInfo()` to get information on how effective the preallocation was; 1883 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 1884 You can also run with the option `-info` and look for messages with the string 1885 malloc in them to see if additional memory allocation was needed. 1886 1887 If the `nnz` parameter is given then the `nz` parameter is ignored 1888 1889 .seealso: [](ch_matrices), `Mat`, [Sparse Matrices](sec_matsparse), `MATSEQSBAIJ`, `MatCreate()`, `MatCreateSeqAIJ()`, `MatSetValues()`, `MatCreateSBAIJ()` 1890 @*/ 1891 PetscErrorCode MatSeqSBAIJSetPreallocation(Mat B, PetscInt bs, PetscInt nz, const PetscInt nnz[]) 1892 { 1893 PetscFunctionBegin; 1894 PetscValidHeaderSpecific(B, MAT_CLASSID, 1); 1895 PetscValidType(B, 1); 1896 PetscValidLogicalCollectiveInt(B, bs, 2); 1897 PetscTryMethod(B, "MatSeqSBAIJSetPreallocation_C", (Mat, PetscInt, PetscInt, const PetscInt[]), (B, bs, nz, nnz)); 1898 PetscFunctionReturn(PETSC_SUCCESS); 1899 } 1900 1901 /*@C 1902 MatSeqSBAIJSetPreallocationCSR - Creates a sparse parallel matrix in `MATSEQSBAIJ` format using the given nonzero structure and (optional) numerical values 1903 1904 Input Parameters: 1905 + B - the matrix 1906 . bs - size of block, the blocks are ALWAYS square. 1907 . i - the indices into `j` for the start of each local row (indices start with zero) 1908 . j - the column indices for each local row (indices start with zero) these must be sorted for each row 1909 - v - optional values in the matrix, use `NULL` if not provided 1910 1911 Level: advanced 1912 1913 Notes: 1914 The `i`,`j`,`v` values are COPIED with this routine; to avoid the copy use `MatCreateSeqSBAIJWithArrays()` 1915 1916 The order of the entries in values is specified by the `MatOption` `MAT_ROW_ORIENTED`. For example, C programs 1917 may want to use the default `MAT_ROW_ORIENTED` = `PETSC_TRUE` and use an array v[nnz][bs][bs] where the second index is 1918 over rows within a block and the last index is over columns within a block row. Fortran programs will likely set 1919 `MAT_ROW_ORIENTED` = `PETSC_FALSE` and use a Fortran array v(bs,bs,nnz) in which the first index is over rows within a 1920 block column and the second index is over columns within a block. 1921 1922 Any entries provided that lie below the diagonal are ignored 1923 1924 Though this routine has Preallocation() in the name it also sets the exact nonzero locations of the matrix entries 1925 and usually the numerical values as well 1926 1927 .seealso: [](ch_matrices), `Mat`, `MATSEQSBAIJ`, `MatCreate()`, `MatCreateSeqSBAIJ()`, `MatSetValuesBlocked()`, `MatSeqSBAIJSetPreallocation()` 1928 @*/ 1929 PetscErrorCode MatSeqSBAIJSetPreallocationCSR(Mat B, PetscInt bs, const PetscInt i[], const PetscInt j[], const PetscScalar v[]) 1930 { 1931 PetscFunctionBegin; 1932 PetscValidHeaderSpecific(B, MAT_CLASSID, 1); 1933 PetscValidType(B, 1); 1934 PetscValidLogicalCollectiveInt(B, bs, 2); 1935 PetscTryMethod(B, "MatSeqSBAIJSetPreallocationCSR_C", (Mat, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[]), (B, bs, i, j, v)); 1936 PetscFunctionReturn(PETSC_SUCCESS); 1937 } 1938 1939 /*@ 1940 MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in (block 1941 compressed row) `MATSEQSBAIJ` format. For good matrix assembly performance the 1942 user should preallocate the matrix storage by setting the parameter `nz` 1943 (or the array `nnz`). 1944 1945 Collective 1946 1947 Input Parameters: 1948 + comm - MPI communicator, set to `PETSC_COMM_SELF` 1949 . bs - size of block, the blocks are ALWAYS square. One can use `MatSetBlockSizes()` to set a different row and column blocksize but the row 1950 blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs() 1951 . m - number of rows 1952 . n - number of columns 1953 . nz - number of block nonzeros per block row (same for all rows) 1954 - nnz - array containing the number of block nonzeros in the upper triangular plus 1955 diagonal portion of each block (possibly different for each block row) or `NULL` 1956 1957 Output Parameter: 1958 . A - the symmetric matrix 1959 1960 Options Database Keys: 1961 + -mat_no_unroll - uses code that does not unroll the loops in the block calculations (much slower) 1962 - -mat_block_size - size of the blocks to use 1963 1964 Level: intermediate 1965 1966 Notes: 1967 It is recommended that one use `MatCreateFromOptions()` or the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`, 1968 MatXXXXSetPreallocation() paradigm instead of this routine directly. 1969 [MatXXXXSetPreallocation() is, for example, `MatSeqAIJSetPreallocation()`] 1970 1971 The number of rows and columns must be divisible by blocksize. 1972 This matrix type does not support complex Hermitian operation. 1973 1974 Specify the preallocated storage with either `nz` or `nnz` (not both). 1975 Set `nz` = `PETSC_DEFAULT` and `nnz` = `NULL` for PETSc to control dynamic memory 1976 allocation. See [Sparse Matrices](sec_matsparse) for details. 1977 1978 If the `nnz` parameter is given then the `nz` parameter is ignored 1979 1980 .seealso: [](ch_matrices), `Mat`, [Sparse Matrices](sec_matsparse), `MATSEQSBAIJ`, `MatCreate()`, `MatCreateSeqAIJ()`, `MatSetValues()`, `MatCreateSBAIJ()` 1981 @*/ 1982 PetscErrorCode MatCreateSeqSBAIJ(MPI_Comm comm, PetscInt bs, PetscInt m, PetscInt n, PetscInt nz, const PetscInt nnz[], Mat *A) 1983 { 1984 PetscFunctionBegin; 1985 PetscCall(MatCreate(comm, A)); 1986 PetscCall(MatSetSizes(*A, m, n, m, n)); 1987 PetscCall(MatSetType(*A, MATSEQSBAIJ)); 1988 PetscCall(MatSeqSBAIJSetPreallocation(*A, bs, nz, (PetscInt *)nnz)); 1989 PetscFunctionReturn(PETSC_SUCCESS); 1990 } 1991 1992 PetscErrorCode MatDuplicate_SeqSBAIJ(Mat A, MatDuplicateOption cpvalues, Mat *B) 1993 { 1994 Mat C; 1995 Mat_SeqSBAIJ *c, *a = (Mat_SeqSBAIJ *)A->data; 1996 PetscInt i, mbs = a->mbs, nz = a->nz, bs2 = a->bs2; 1997 1998 PetscFunctionBegin; 1999 PetscCheck(A->assembled, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Cannot duplicate unassembled matrix"); 2000 PetscCheck(a->i[mbs] == nz, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Corrupt matrix"); 2001 2002 *B = NULL; 2003 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &C)); 2004 PetscCall(MatSetSizes(C, A->rmap->N, A->cmap->n, A->rmap->N, A->cmap->n)); 2005 PetscCall(MatSetBlockSizesFromMats(C, A, A)); 2006 PetscCall(MatSetType(C, MATSEQSBAIJ)); 2007 c = (Mat_SeqSBAIJ *)C->data; 2008 2009 C->preallocated = PETSC_TRUE; 2010 C->factortype = A->factortype; 2011 c->row = NULL; 2012 c->icol = NULL; 2013 c->saved_values = NULL; 2014 c->keepnonzeropattern = a->keepnonzeropattern; 2015 C->assembled = PETSC_TRUE; 2016 2017 PetscCall(PetscLayoutReference(A->rmap, &C->rmap)); 2018 PetscCall(PetscLayoutReference(A->cmap, &C->cmap)); 2019 c->bs2 = a->bs2; 2020 c->mbs = a->mbs; 2021 c->nbs = a->nbs; 2022 2023 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2024 c->imax = a->imax; 2025 c->ilen = a->ilen; 2026 c->free_imax_ilen = PETSC_FALSE; 2027 } else { 2028 PetscCall(PetscMalloc2(mbs + 1, &c->imax, mbs + 1, &c->ilen)); 2029 for (i = 0; i < mbs; i++) { 2030 c->imax[i] = a->imax[i]; 2031 c->ilen[i] = a->ilen[i]; 2032 } 2033 c->free_imax_ilen = PETSC_TRUE; 2034 } 2035 2036 /* allocate the matrix space */ 2037 PetscCall(PetscShmgetAllocateArray(bs2 * nz, sizeof(PetscScalar), (void **)&c->a)); 2038 c->free_a = PETSC_TRUE; 2039 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2040 PetscCall(PetscArrayzero(c->a, bs2 * nz)); 2041 c->i = a->i; 2042 c->j = a->j; 2043 c->free_ij = PETSC_FALSE; 2044 c->parent = A; 2045 PetscCall(PetscObjectReference((PetscObject)A)); 2046 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE)); 2047 PetscCall(MatSetOption(C, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE)); 2048 } else { 2049 PetscCall(PetscShmgetAllocateArray(nz, sizeof(PetscInt), (void **)&c->j)); 2050 PetscCall(PetscShmgetAllocateArray(mbs + 1, sizeof(PetscInt), (void **)&c->i)); 2051 PetscCall(PetscArraycpy(c->i, a->i, mbs + 1)); 2052 c->free_ij = PETSC_TRUE; 2053 } 2054 if (mbs > 0) { 2055 if (cpvalues != MAT_SHARE_NONZERO_PATTERN) PetscCall(PetscArraycpy(c->j, a->j, nz)); 2056 if (cpvalues == MAT_COPY_VALUES) { 2057 PetscCall(PetscArraycpy(c->a, a->a, bs2 * nz)); 2058 } else { 2059 PetscCall(PetscArrayzero(c->a, bs2 * nz)); 2060 } 2061 if (a->jshort) { 2062 /* cannot share jshort, it is reallocated in MatAssemblyEnd_SeqSBAIJ() */ 2063 /* if the parent matrix is reassembled, this child matrix will never notice */ 2064 PetscCall(PetscMalloc1(nz, &c->jshort)); 2065 PetscCall(PetscArraycpy(c->jshort, a->jshort, nz)); 2066 2067 c->free_jshort = PETSC_TRUE; 2068 } 2069 } 2070 2071 c->roworiented = a->roworiented; 2072 c->nonew = a->nonew; 2073 c->nz = a->nz; 2074 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 2075 c->solve_work = NULL; 2076 c->mult_work = NULL; 2077 2078 *B = C; 2079 PetscCall(PetscFunctionListDuplicate(((PetscObject)A)->qlist, &((PetscObject)C)->qlist)); 2080 PetscFunctionReturn(PETSC_SUCCESS); 2081 } 2082 2083 /* Used for both SeqBAIJ and SeqSBAIJ matrices */ 2084 #define MatLoad_SeqSBAIJ_Binary MatLoad_SeqBAIJ_Binary 2085 2086 PetscErrorCode MatLoad_SeqSBAIJ(Mat mat, PetscViewer viewer) 2087 { 2088 PetscBool isbinary; 2089 2090 PetscFunctionBegin; 2091 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary)); 2092 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); 2093 PetscCall(MatLoad_SeqSBAIJ_Binary(mat, viewer)); 2094 PetscFunctionReturn(PETSC_SUCCESS); 2095 } 2096 2097 /*@ 2098 MatCreateSeqSBAIJWithArrays - Creates an sequential `MATSEQSBAIJ` matrix using matrix elements 2099 (upper triangular entries in CSR format) provided by the user. 2100 2101 Collective 2102 2103 Input Parameters: 2104 + comm - must be an MPI communicator of size 1 2105 . bs - size of block 2106 . m - number of rows 2107 . n - number of columns 2108 . 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 2109 . j - column indices 2110 - a - matrix values 2111 2112 Output Parameter: 2113 . mat - the matrix 2114 2115 Level: advanced 2116 2117 Notes: 2118 The `i`, `j`, and `a` arrays are not copied by this routine, the user must free these arrays 2119 once the matrix is destroyed 2120 2121 You cannot set new nonzero locations into this matrix, that will generate an error. 2122 2123 The `i` and `j` indices are 0 based 2124 2125 When block size is greater than 1 the matrix values must be stored using the `MATSBAIJ` storage format. For block size of 1 2126 it is the regular CSR format excluding the lower triangular elements. 2127 2128 .seealso: [](ch_matrices), `Mat`, `MATSEQSBAIJ`, `MatCreate()`, `MatCreateSBAIJ()`, `MatCreateSeqSBAIJ()` 2129 @*/ 2130 PetscErrorCode MatCreateSeqSBAIJWithArrays(MPI_Comm comm, PetscInt bs, PetscInt m, PetscInt n, PetscInt i[], PetscInt j[], PetscScalar a[], Mat *mat) 2131 { 2132 PetscInt ii; 2133 Mat_SeqSBAIJ *sbaij; 2134 2135 PetscFunctionBegin; 2136 PetscCheck(bs == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "block size %" PetscInt_FMT " > 1 is not supported yet", bs); 2137 PetscCheck(m == 0 || i[0] == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "i (row indices) must start with 0"); 2138 2139 PetscCall(MatCreate(comm, mat)); 2140 PetscCall(MatSetSizes(*mat, m, n, m, n)); 2141 PetscCall(MatSetType(*mat, MATSEQSBAIJ)); 2142 PetscCall(MatSeqSBAIJSetPreallocation(*mat, bs, MAT_SKIP_ALLOCATION, NULL)); 2143 sbaij = (Mat_SeqSBAIJ *)(*mat)->data; 2144 PetscCall(PetscMalloc2(m, &sbaij->imax, m, &sbaij->ilen)); 2145 2146 sbaij->i = i; 2147 sbaij->j = j; 2148 sbaij->a = a; 2149 2150 sbaij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2151 sbaij->free_a = PETSC_FALSE; 2152 sbaij->free_ij = PETSC_FALSE; 2153 sbaij->free_imax_ilen = PETSC_TRUE; 2154 2155 for (ii = 0; ii < m; ii++) { 2156 sbaij->ilen[ii] = sbaij->imax[ii] = i[ii + 1] - i[ii]; 2157 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]); 2158 } 2159 if (PetscDefined(USE_DEBUG)) { 2160 for (ii = 0; ii < sbaij->i[m]; ii++) { 2161 PetscCheck(j[ii] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative column index at location = %" PetscInt_FMT " index = %" PetscInt_FMT, ii, j[ii]); 2162 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]); 2163 } 2164 } 2165 2166 PetscCall(MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY)); 2167 PetscCall(MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY)); 2168 PetscFunctionReturn(PETSC_SUCCESS); 2169 } 2170 2171 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ(MPI_Comm comm, Mat inmat, PetscInt n, MatReuse scall, Mat *outmat) 2172 { 2173 PetscFunctionBegin; 2174 PetscCall(MatCreateMPIMatConcatenateSeqMat_MPISBAIJ(comm, inmat, n, scall, outmat)); 2175 PetscFunctionReturn(PETSC_SUCCESS); 2176 } 2177