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