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