1 /* 2 Defines the basic matrix operations for the AIJ (compressed row) 3 matrix storage format. 4 */ 5 6 #include <../src/mat/impls/aij/seq/aij.h> /*I "petscmat.h" I*/ 7 #include <petscblaslapack.h> 8 #include <petscbt.h> 9 #include <petsc/private/kernels/blocktranspose.h> 10 11 PetscErrorCode MatSeqAIJSetTypeFromOptions(Mat A) 12 { 13 PetscBool flg; 14 char type[256]; 15 16 PetscFunctionBegin; 17 PetscObjectOptionsBegin((PetscObject)A); 18 PetscCall(PetscOptionsFList("-mat_seqaij_type", "Matrix SeqAIJ type", "MatSeqAIJSetType", MatSeqAIJList, "seqaij", type, 256, &flg)); 19 if (flg) PetscCall(MatSeqAIJSetType(A, type)); 20 PetscOptionsEnd(); 21 PetscFunctionReturn(0); 22 } 23 24 PetscErrorCode MatGetColumnReductions_SeqAIJ(Mat A, PetscInt type, PetscReal *reductions) 25 { 26 PetscInt i, m, n; 27 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)A->data; 28 29 PetscFunctionBegin; 30 PetscCall(MatGetSize(A, &m, &n)); 31 PetscCall(PetscArrayzero(reductions, n)); 32 if (type == NORM_2) { 33 for (i = 0; i < aij->i[m]; i++) reductions[aij->j[i]] += PetscAbsScalar(aij->a[i] * aij->a[i]); 34 } else if (type == NORM_1) { 35 for (i = 0; i < aij->i[m]; i++) reductions[aij->j[i]] += PetscAbsScalar(aij->a[i]); 36 } else if (type == NORM_INFINITY) { 37 for (i = 0; i < aij->i[m]; i++) reductions[aij->j[i]] = PetscMax(PetscAbsScalar(aij->a[i]), reductions[aij->j[i]]); 38 } else if (type == REDUCTION_SUM_REALPART || type == REDUCTION_MEAN_REALPART) { 39 for (i = 0; i < aij->i[m]; i++) reductions[aij->j[i]] += PetscRealPart(aij->a[i]); 40 } else if (type == REDUCTION_SUM_IMAGINARYPART || type == REDUCTION_MEAN_IMAGINARYPART) { 41 for (i = 0; i < aij->i[m]; i++) reductions[aij->j[i]] += PetscImaginaryPart(aij->a[i]); 42 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown reduction type"); 43 44 if (type == NORM_2) { 45 for (i = 0; i < n; i++) reductions[i] = PetscSqrtReal(reductions[i]); 46 } else if (type == REDUCTION_MEAN_REALPART || type == REDUCTION_MEAN_IMAGINARYPART) { 47 for (i = 0; i < n; i++) reductions[i] /= m; 48 } 49 PetscFunctionReturn(0); 50 } 51 52 PetscErrorCode MatFindOffBlockDiagonalEntries_SeqAIJ(Mat A, IS *is) 53 { 54 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 55 PetscInt i, m = A->rmap->n, cnt = 0, bs = A->rmap->bs; 56 const PetscInt *jj = a->j, *ii = a->i; 57 PetscInt *rows; 58 59 PetscFunctionBegin; 60 for (i = 0; i < m; i++) { 61 if ((ii[i] != ii[i + 1]) && ((jj[ii[i]] < bs * (i / bs)) || (jj[ii[i + 1] - 1] > bs * ((i + bs) / bs) - 1))) cnt++; 62 } 63 PetscCall(PetscMalloc1(cnt, &rows)); 64 cnt = 0; 65 for (i = 0; i < m; i++) { 66 if ((ii[i] != ii[i + 1]) && ((jj[ii[i]] < bs * (i / bs)) || (jj[ii[i + 1] - 1] > bs * ((i + bs) / bs) - 1))) { 67 rows[cnt] = i; 68 cnt++; 69 } 70 } 71 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, cnt, rows, PETSC_OWN_POINTER, is)); 72 PetscFunctionReturn(0); 73 } 74 75 PetscErrorCode MatFindZeroDiagonals_SeqAIJ_Private(Mat A, PetscInt *nrows, PetscInt **zrows) 76 { 77 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 78 const MatScalar *aa; 79 PetscInt i, m = A->rmap->n, cnt = 0; 80 const PetscInt *ii = a->i, *jj = a->j, *diag; 81 PetscInt *rows; 82 83 PetscFunctionBegin; 84 PetscCall(MatSeqAIJGetArrayRead(A, &aa)); 85 PetscCall(MatMarkDiagonal_SeqAIJ(A)); 86 diag = a->diag; 87 for (i = 0; i < m; i++) { 88 if ((diag[i] >= ii[i + 1]) || (jj[diag[i]] != i) || (aa[diag[i]] == 0.0)) cnt++; 89 } 90 PetscCall(PetscMalloc1(cnt, &rows)); 91 cnt = 0; 92 for (i = 0; i < m; i++) { 93 if ((diag[i] >= ii[i + 1]) || (jj[diag[i]] != i) || (aa[diag[i]] == 0.0)) rows[cnt++] = i; 94 } 95 *nrows = cnt; 96 *zrows = rows; 97 PetscCall(MatSeqAIJRestoreArrayRead(A, &aa)); 98 PetscFunctionReturn(0); 99 } 100 101 PetscErrorCode MatFindZeroDiagonals_SeqAIJ(Mat A, IS *zrows) 102 { 103 PetscInt nrows, *rows; 104 105 PetscFunctionBegin; 106 *zrows = NULL; 107 PetscCall(MatFindZeroDiagonals_SeqAIJ_Private(A, &nrows, &rows)); 108 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)A), nrows, rows, PETSC_OWN_POINTER, zrows)); 109 PetscFunctionReturn(0); 110 } 111 112 PetscErrorCode MatFindNonzeroRows_SeqAIJ(Mat A, IS *keptrows) 113 { 114 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 115 const MatScalar *aa; 116 PetscInt m = A->rmap->n, cnt = 0; 117 const PetscInt *ii; 118 PetscInt n, i, j, *rows; 119 120 PetscFunctionBegin; 121 PetscCall(MatSeqAIJGetArrayRead(A, &aa)); 122 *keptrows = NULL; 123 ii = a->i; 124 for (i = 0; i < m; i++) { 125 n = ii[i + 1] - ii[i]; 126 if (!n) { 127 cnt++; 128 goto ok1; 129 } 130 for (j = ii[i]; j < ii[i + 1]; j++) { 131 if (aa[j] != 0.0) goto ok1; 132 } 133 cnt++; 134 ok1:; 135 } 136 if (!cnt) { 137 PetscCall(MatSeqAIJRestoreArrayRead(A, &aa)); 138 PetscFunctionReturn(0); 139 } 140 PetscCall(PetscMalloc1(A->rmap->n - cnt, &rows)); 141 cnt = 0; 142 for (i = 0; i < m; i++) { 143 n = ii[i + 1] - ii[i]; 144 if (!n) continue; 145 for (j = ii[i]; j < ii[i + 1]; j++) { 146 if (aa[j] != 0.0) { 147 rows[cnt++] = i; 148 break; 149 } 150 } 151 } 152 PetscCall(MatSeqAIJRestoreArrayRead(A, &aa)); 153 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, cnt, rows, PETSC_OWN_POINTER, keptrows)); 154 PetscFunctionReturn(0); 155 } 156 157 PetscErrorCode MatDiagonalSet_SeqAIJ(Mat Y, Vec D, InsertMode is) 158 { 159 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)Y->data; 160 PetscInt i, m = Y->rmap->n; 161 const PetscInt *diag; 162 MatScalar *aa; 163 const PetscScalar *v; 164 PetscBool missing; 165 166 PetscFunctionBegin; 167 if (Y->assembled) { 168 PetscCall(MatMissingDiagonal_SeqAIJ(Y, &missing, NULL)); 169 if (!missing) { 170 diag = aij->diag; 171 PetscCall(VecGetArrayRead(D, &v)); 172 PetscCall(MatSeqAIJGetArray(Y, &aa)); 173 if (is == INSERT_VALUES) { 174 for (i = 0; i < m; i++) aa[diag[i]] = v[i]; 175 } else { 176 for (i = 0; i < m; i++) aa[diag[i]] += v[i]; 177 } 178 PetscCall(MatSeqAIJRestoreArray(Y, &aa)); 179 PetscCall(VecRestoreArrayRead(D, &v)); 180 PetscFunctionReturn(0); 181 } 182 PetscCall(MatSeqAIJInvalidateDiagonal(Y)); 183 } 184 PetscCall(MatDiagonalSet_Default(Y, D, is)); 185 PetscFunctionReturn(0); 186 } 187 188 PetscErrorCode MatGetRowIJ_SeqAIJ(Mat A, PetscInt oshift, PetscBool symmetric, PetscBool inodecompressed, PetscInt *m, const PetscInt *ia[], const PetscInt *ja[], PetscBool *done) 189 { 190 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 191 PetscInt i, ishift; 192 193 PetscFunctionBegin; 194 if (m) *m = A->rmap->n; 195 if (!ia) PetscFunctionReturn(0); 196 ishift = 0; 197 if (symmetric && A->structurally_symmetric != PETSC_BOOL3_TRUE) { 198 PetscCall(MatToSymmetricIJ_SeqAIJ(A->rmap->n, a->i, a->j, PETSC_TRUE, ishift, oshift, (PetscInt **)ia, (PetscInt **)ja)); 199 } else if (oshift == 1) { 200 PetscInt *tia; 201 PetscInt nz = a->i[A->rmap->n]; 202 /* malloc space and add 1 to i and j indices */ 203 PetscCall(PetscMalloc1(A->rmap->n + 1, &tia)); 204 for (i = 0; i < A->rmap->n + 1; i++) tia[i] = a->i[i] + 1; 205 *ia = tia; 206 if (ja) { 207 PetscInt *tja; 208 PetscCall(PetscMalloc1(nz + 1, &tja)); 209 for (i = 0; i < nz; i++) tja[i] = a->j[i] + 1; 210 *ja = tja; 211 } 212 } else { 213 *ia = a->i; 214 if (ja) *ja = a->j; 215 } 216 PetscFunctionReturn(0); 217 } 218 219 PetscErrorCode MatRestoreRowIJ_SeqAIJ(Mat A, PetscInt oshift, PetscBool symmetric, PetscBool inodecompressed, PetscInt *n, const PetscInt *ia[], const PetscInt *ja[], PetscBool *done) 220 { 221 PetscFunctionBegin; 222 if (!ia) PetscFunctionReturn(0); 223 if ((symmetric && A->structurally_symmetric != PETSC_BOOL3_TRUE) || oshift == 1) { 224 PetscCall(PetscFree(*ia)); 225 if (ja) PetscCall(PetscFree(*ja)); 226 } 227 PetscFunctionReturn(0); 228 } 229 230 PetscErrorCode MatGetColumnIJ_SeqAIJ(Mat A, PetscInt oshift, PetscBool symmetric, PetscBool inodecompressed, PetscInt *nn, const PetscInt *ia[], const PetscInt *ja[], PetscBool *done) 231 { 232 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 233 PetscInt i, *collengths, *cia, *cja, n = A->cmap->n, m = A->rmap->n; 234 PetscInt nz = a->i[m], row, *jj, mr, col; 235 236 PetscFunctionBegin; 237 *nn = n; 238 if (!ia) PetscFunctionReturn(0); 239 if (symmetric) { 240 PetscCall(MatToSymmetricIJ_SeqAIJ(A->rmap->n, a->i, a->j, PETSC_TRUE, 0, oshift, (PetscInt **)ia, (PetscInt **)ja)); 241 } else { 242 PetscCall(PetscCalloc1(n, &collengths)); 243 PetscCall(PetscMalloc1(n + 1, &cia)); 244 PetscCall(PetscMalloc1(nz, &cja)); 245 jj = a->j; 246 for (i = 0; i < nz; i++) collengths[jj[i]]++; 247 cia[0] = oshift; 248 for (i = 0; i < n; i++) cia[i + 1] = cia[i] + collengths[i]; 249 PetscCall(PetscArrayzero(collengths, n)); 250 jj = a->j; 251 for (row = 0; row < m; row++) { 252 mr = a->i[row + 1] - a->i[row]; 253 for (i = 0; i < mr; i++) { 254 col = *jj++; 255 256 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 257 } 258 } 259 PetscCall(PetscFree(collengths)); 260 *ia = cia; 261 *ja = cja; 262 } 263 PetscFunctionReturn(0); 264 } 265 266 PetscErrorCode MatRestoreColumnIJ_SeqAIJ(Mat A, PetscInt oshift, PetscBool symmetric, PetscBool inodecompressed, PetscInt *n, const PetscInt *ia[], const PetscInt *ja[], PetscBool *done) 267 { 268 PetscFunctionBegin; 269 if (!ia) PetscFunctionReturn(0); 270 271 PetscCall(PetscFree(*ia)); 272 PetscCall(PetscFree(*ja)); 273 PetscFunctionReturn(0); 274 } 275 276 /* 277 MatGetColumnIJ_SeqAIJ_Color() and MatRestoreColumnIJ_SeqAIJ_Color() are customized from 278 MatGetColumnIJ_SeqAIJ() and MatRestoreColumnIJ_SeqAIJ() by adding an output 279 spidx[], index of a->a, to be used in MatTransposeColoringCreate_SeqAIJ() and MatFDColoringCreate_SeqXAIJ() 280 */ 281 PetscErrorCode MatGetColumnIJ_SeqAIJ_Color(Mat A, PetscInt oshift, PetscBool symmetric, PetscBool inodecompressed, PetscInt *nn, const PetscInt *ia[], const PetscInt *ja[], PetscInt *spidx[], PetscBool *done) 282 { 283 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 284 PetscInt i, *collengths, *cia, *cja, n = A->cmap->n, m = A->rmap->n; 285 PetscInt nz = a->i[m], row, mr, col, tmp; 286 PetscInt *cspidx; 287 const PetscInt *jj; 288 289 PetscFunctionBegin; 290 *nn = n; 291 if (!ia) PetscFunctionReturn(0); 292 293 PetscCall(PetscCalloc1(n, &collengths)); 294 PetscCall(PetscMalloc1(n + 1, &cia)); 295 PetscCall(PetscMalloc1(nz, &cja)); 296 PetscCall(PetscMalloc1(nz, &cspidx)); 297 jj = a->j; 298 for (i = 0; i < nz; i++) collengths[jj[i]]++; 299 cia[0] = oshift; 300 for (i = 0; i < n; i++) cia[i + 1] = cia[i] + collengths[i]; 301 PetscCall(PetscArrayzero(collengths, n)); 302 jj = a->j; 303 for (row = 0; row < m; row++) { 304 mr = a->i[row + 1] - a->i[row]; 305 for (i = 0; i < mr; i++) { 306 col = *jj++; 307 tmp = cia[col] + collengths[col]++ - oshift; 308 cspidx[tmp] = a->i[row] + i; /* index of a->j */ 309 cja[tmp] = row + oshift; 310 } 311 } 312 PetscCall(PetscFree(collengths)); 313 *ia = cia; 314 *ja = cja; 315 *spidx = cspidx; 316 PetscFunctionReturn(0); 317 } 318 319 PetscErrorCode MatRestoreColumnIJ_SeqAIJ_Color(Mat A, PetscInt oshift, PetscBool symmetric, PetscBool inodecompressed, PetscInt *n, const PetscInt *ia[], const PetscInt *ja[], PetscInt *spidx[], PetscBool *done) 320 { 321 PetscFunctionBegin; 322 PetscCall(MatRestoreColumnIJ_SeqAIJ(A, oshift, symmetric, inodecompressed, n, ia, ja, done)); 323 PetscCall(PetscFree(*spidx)); 324 PetscFunctionReturn(0); 325 } 326 327 PetscErrorCode MatSetValuesRow_SeqAIJ(Mat A, PetscInt row, const PetscScalar v[]) 328 { 329 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 330 PetscInt *ai = a->i; 331 PetscScalar *aa; 332 333 PetscFunctionBegin; 334 PetscCall(MatSeqAIJGetArray(A, &aa)); 335 PetscCall(PetscArraycpy(aa + ai[row], v, ai[row + 1] - ai[row])); 336 PetscCall(MatSeqAIJRestoreArray(A, &aa)); 337 PetscFunctionReturn(0); 338 } 339 340 /* 341 MatSeqAIJSetValuesLocalFast - An optimized version of MatSetValuesLocal() for SeqAIJ matrices with several assumptions 342 343 - a single row of values is set with each call 344 - no row or column indices are negative or (in error) larger than the number of rows or columns 345 - the values are always added to the matrix, not set 346 - no new locations are introduced in the nonzero structure of the matrix 347 348 This does NOT assume the global column indices are sorted 349 350 */ 351 352 #include <petsc/private/isimpl.h> 353 PetscErrorCode MatSeqAIJSetValuesLocalFast(Mat A, PetscInt m, const PetscInt im[], PetscInt n, const PetscInt in[], const PetscScalar v[], InsertMode is) 354 { 355 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 356 PetscInt low, high, t, row, nrow, i, col, l; 357 const PetscInt *rp, *ai = a->i, *ailen = a->ilen, *aj = a->j; 358 PetscInt lastcol = -1; 359 MatScalar *ap, value, *aa; 360 const PetscInt *ridx = A->rmap->mapping->indices, *cidx = A->cmap->mapping->indices; 361 362 PetscFunctionBegin; 363 PetscCall(MatSeqAIJGetArray(A, &aa)); 364 row = ridx[im[0]]; 365 rp = aj + ai[row]; 366 ap = aa + ai[row]; 367 nrow = ailen[row]; 368 low = 0; 369 high = nrow; 370 for (l = 0; l < n; l++) { /* loop over added columns */ 371 col = cidx[in[l]]; 372 value = v[l]; 373 374 if (col <= lastcol) low = 0; 375 else high = nrow; 376 lastcol = col; 377 while (high - low > 5) { 378 t = (low + high) / 2; 379 if (rp[t] > col) high = t; 380 else low = t; 381 } 382 for (i = low; i < high; i++) { 383 if (rp[i] == col) { 384 ap[i] += value; 385 low = i + 1; 386 break; 387 } 388 } 389 } 390 PetscCall(MatSeqAIJRestoreArray(A, &aa)); 391 return 0; 392 } 393 394 PetscErrorCode MatSetValues_SeqAIJ(Mat A, PetscInt m, const PetscInt im[], PetscInt n, const PetscInt in[], const PetscScalar v[], InsertMode is) 395 { 396 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 397 PetscInt *rp, k, low, high, t, ii, row, nrow, i, col, l, rmax, N; 398 PetscInt *imax = a->imax, *ai = a->i, *ailen = a->ilen; 399 PetscInt *aj = a->j, nonew = a->nonew, lastcol = -1; 400 MatScalar *ap = NULL, value = 0.0, *aa; 401 PetscBool ignorezeroentries = a->ignorezeroentries; 402 PetscBool roworiented = a->roworiented; 403 404 PetscFunctionBegin; 405 PetscCall(MatSeqAIJGetArray(A, &aa)); 406 for (k = 0; k < m; k++) { /* loop over added rows */ 407 row = im[k]; 408 if (row < 0) continue; 409 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); 410 rp = aj + ai[row]; 411 if (!A->structure_only) ap = aa + ai[row]; 412 rmax = imax[row]; 413 nrow = ailen[row]; 414 low = 0; 415 high = nrow; 416 for (l = 0; l < n; l++) { /* loop over added columns */ 417 if (in[l] < 0) continue; 418 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); 419 col = in[l]; 420 if (v && !A->structure_only) value = roworiented ? v[l + k * n] : v[k + l * m]; 421 if (!A->structure_only && value == 0.0 && ignorezeroentries && is == ADD_VALUES && row != col) continue; 422 423 if (col <= lastcol) low = 0; 424 else high = nrow; 425 lastcol = col; 426 while (high - low > 5) { 427 t = (low + high) / 2; 428 if (rp[t] > col) high = t; 429 else low = t; 430 } 431 for (i = low; i < high; i++) { 432 if (rp[i] > col) break; 433 if (rp[i] == col) { 434 if (!A->structure_only) { 435 if (is == ADD_VALUES) { 436 ap[i] += value; 437 (void)PetscLogFlops(1.0); 438 } else ap[i] = value; 439 } 440 low = i + 1; 441 goto noinsert; 442 } 443 } 444 if (value == 0.0 && ignorezeroentries && row != col) goto noinsert; 445 if (nonew == 1) goto noinsert; 446 PetscCheck(nonew != -1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Inserting a new nonzero at (%" PetscInt_FMT ",%" PetscInt_FMT ") in the matrix", row, col); 447 if (A->structure_only) { 448 MatSeqXAIJReallocateAIJ_structure_only(A, A->rmap->n, 1, nrow, row, col, rmax, ai, aj, rp, imax, nonew, MatScalar); 449 } else { 450 MatSeqXAIJReallocateAIJ(A, A->rmap->n, 1, nrow, row, col, rmax, aa, ai, aj, rp, ap, imax, nonew, MatScalar); 451 } 452 N = nrow++ - 1; 453 a->nz++; 454 high++; 455 /* shift up all the later entries in this row */ 456 PetscCall(PetscArraymove(rp + i + 1, rp + i, N - i + 1)); 457 rp[i] = col; 458 if (!A->structure_only) { 459 PetscCall(PetscArraymove(ap + i + 1, ap + i, N - i + 1)); 460 ap[i] = value; 461 } 462 low = i + 1; 463 A->nonzerostate++; 464 noinsert:; 465 } 466 ailen[row] = nrow; 467 } 468 PetscCall(MatSeqAIJRestoreArray(A, &aa)); 469 PetscFunctionReturn(0); 470 } 471 472 PetscErrorCode MatSetValues_SeqAIJ_SortedFullNoPreallocation(Mat A, PetscInt m, const PetscInt im[], PetscInt n, const PetscInt in[], const PetscScalar v[], InsertMode is) 473 { 474 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 475 PetscInt *rp, k, row; 476 PetscInt *ai = a->i; 477 PetscInt *aj = a->j; 478 MatScalar *aa, *ap; 479 480 PetscFunctionBegin; 481 PetscCheck(!A->was_assembled, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cannot call on assembled matrix."); 482 PetscCheck(m * n + a->nz <= a->maxnz, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of entries in matrix will be larger than maximum nonzeros allocated for %" PetscInt_FMT " in MatSeqAIJSetTotalPreallocation()", a->maxnz); 483 484 PetscCall(MatSeqAIJGetArray(A, &aa)); 485 for (k = 0; k < m; k++) { /* loop over added rows */ 486 row = im[k]; 487 rp = aj + ai[row]; 488 ap = aa + ai[row]; 489 490 PetscCall(PetscMemcpy(rp, in, n * sizeof(PetscInt))); 491 if (!A->structure_only) { 492 if (v) { 493 PetscCall(PetscMemcpy(ap, v, n * sizeof(PetscScalar))); 494 v += n; 495 } else { 496 PetscCall(PetscMemzero(ap, n * sizeof(PetscScalar))); 497 } 498 } 499 a->ilen[row] = n; 500 a->imax[row] = n; 501 a->i[row + 1] = a->i[row] + n; 502 a->nz += n; 503 } 504 PetscCall(MatSeqAIJRestoreArray(A, &aa)); 505 PetscFunctionReturn(0); 506 } 507 508 /*@ 509 MatSeqAIJSetTotalPreallocation - Sets an upper bound on the total number of expected nonzeros in the matrix. 510 511 Input Parameters: 512 + A - the `MATSEQAIJ` matrix 513 - nztotal - bound on the number of nonzeros 514 515 Level: advanced 516 517 Notes: 518 This can be called if you will be provided the matrix row by row (from row zero) with sorted column indices for each row. 519 Simply call `MatSetValues()` after this call to provide the matrix entries in the usual manner. This matrix may be used 520 as always with multiple matrix assemblies. 521 522 .seealso: `MatSetOption()`, `MAT_SORTED_FULL`, `MatSetValues()`, `MatSeqAIJSetPreallocation()` 523 @*/ 524 525 PetscErrorCode MatSeqAIJSetTotalPreallocation(Mat A, PetscInt nztotal) 526 { 527 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 528 529 PetscFunctionBegin; 530 PetscCall(PetscLayoutSetUp(A->rmap)); 531 PetscCall(PetscLayoutSetUp(A->cmap)); 532 a->maxnz = nztotal; 533 if (!a->imax) { PetscCall(PetscMalloc1(A->rmap->n, &a->imax)); } 534 if (!a->ilen) { 535 PetscCall(PetscMalloc1(A->rmap->n, &a->ilen)); 536 } else { 537 PetscCall(PetscMemzero(a->ilen, A->rmap->n * sizeof(PetscInt))); 538 } 539 540 /* allocate the matrix space */ 541 if (A->structure_only) { 542 PetscCall(PetscMalloc1(nztotal, &a->j)); 543 PetscCall(PetscMalloc1(A->rmap->n + 1, &a->i)); 544 } else { 545 PetscCall(PetscMalloc3(nztotal, &a->a, nztotal, &a->j, A->rmap->n + 1, &a->i)); 546 } 547 a->i[0] = 0; 548 if (A->structure_only) { 549 a->singlemalloc = PETSC_FALSE; 550 a->free_a = PETSC_FALSE; 551 } else { 552 a->singlemalloc = PETSC_TRUE; 553 a->free_a = PETSC_TRUE; 554 } 555 a->free_ij = PETSC_TRUE; 556 A->ops->setvalues = MatSetValues_SeqAIJ_SortedFullNoPreallocation; 557 A->preallocated = PETSC_TRUE; 558 PetscFunctionReturn(0); 559 } 560 561 PetscErrorCode MatSetValues_SeqAIJ_SortedFull(Mat A, PetscInt m, const PetscInt im[], PetscInt n, const PetscInt in[], const PetscScalar v[], InsertMode is) 562 { 563 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 564 PetscInt *rp, k, row; 565 PetscInt *ai = a->i, *ailen = a->ilen; 566 PetscInt *aj = a->j; 567 MatScalar *aa, *ap; 568 569 PetscFunctionBegin; 570 PetscCall(MatSeqAIJGetArray(A, &aa)); 571 for (k = 0; k < m; k++) { /* loop over added rows */ 572 row = im[k]; 573 PetscCheck(n <= a->imax[row], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Preallocation for row %" PetscInt_FMT " does not match number of columns provided", n); 574 rp = aj + ai[row]; 575 ap = aa + ai[row]; 576 if (!A->was_assembled) PetscCall(PetscMemcpy(rp, in, n * sizeof(PetscInt))); 577 if (!A->structure_only) { 578 if (v) { 579 PetscCall(PetscMemcpy(ap, v, n * sizeof(PetscScalar))); 580 v += n; 581 } else { 582 PetscCall(PetscMemzero(ap, n * sizeof(PetscScalar))); 583 } 584 } 585 ailen[row] = n; 586 a->nz += n; 587 } 588 PetscCall(MatSeqAIJRestoreArray(A, &aa)); 589 PetscFunctionReturn(0); 590 } 591 592 PetscErrorCode MatGetValues_SeqAIJ(Mat A, PetscInt m, const PetscInt im[], PetscInt n, const PetscInt in[], PetscScalar v[]) 593 { 594 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 595 PetscInt *rp, k, low, high, t, row, nrow, i, col, l, *aj = a->j; 596 PetscInt *ai = a->i, *ailen = a->ilen; 597 MatScalar *ap, *aa; 598 599 PetscFunctionBegin; 600 PetscCall(MatSeqAIJGetArray(A, &aa)); 601 for (k = 0; k < m; k++) { /* loop over rows */ 602 row = im[k]; 603 if (row < 0) { 604 v += n; 605 continue; 606 } /* negative row */ 607 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); 608 rp = aj + ai[row]; 609 ap = aa + ai[row]; 610 nrow = ailen[row]; 611 for (l = 0; l < n; l++) { /* loop over columns */ 612 if (in[l] < 0) { 613 v++; 614 continue; 615 } /* negative column */ 616 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); 617 col = in[l]; 618 high = nrow; 619 low = 0; /* assume unsorted */ 620 while (high - low > 5) { 621 t = (low + high) / 2; 622 if (rp[t] > col) high = t; 623 else low = t; 624 } 625 for (i = low; i < high; i++) { 626 if (rp[i] > col) break; 627 if (rp[i] == col) { 628 *v++ = ap[i]; 629 goto finished; 630 } 631 } 632 *v++ = 0.0; 633 finished:; 634 } 635 } 636 PetscCall(MatSeqAIJRestoreArray(A, &aa)); 637 PetscFunctionReturn(0); 638 } 639 640 PetscErrorCode MatView_SeqAIJ_Binary(Mat mat, PetscViewer viewer) 641 { 642 Mat_SeqAIJ *A = (Mat_SeqAIJ *)mat->data; 643 const PetscScalar *av; 644 PetscInt header[4], M, N, m, nz, i; 645 PetscInt *rowlens; 646 647 PetscFunctionBegin; 648 PetscCall(PetscViewerSetUp(viewer)); 649 650 M = mat->rmap->N; 651 N = mat->cmap->N; 652 m = mat->rmap->n; 653 nz = A->nz; 654 655 /* write matrix header */ 656 header[0] = MAT_FILE_CLASSID; 657 header[1] = M; 658 header[2] = N; 659 header[3] = nz; 660 PetscCall(PetscViewerBinaryWrite(viewer, header, 4, PETSC_INT)); 661 662 /* fill in and store row lengths */ 663 PetscCall(PetscMalloc1(m, &rowlens)); 664 for (i = 0; i < m; i++) rowlens[i] = A->i[i + 1] - A->i[i]; 665 PetscCall(PetscViewerBinaryWrite(viewer, rowlens, m, PETSC_INT)); 666 PetscCall(PetscFree(rowlens)); 667 /* store column indices */ 668 PetscCall(PetscViewerBinaryWrite(viewer, A->j, nz, PETSC_INT)); 669 /* store nonzero values */ 670 PetscCall(MatSeqAIJGetArrayRead(mat, &av)); 671 PetscCall(PetscViewerBinaryWrite(viewer, av, nz, PETSC_SCALAR)); 672 PetscCall(MatSeqAIJRestoreArrayRead(mat, &av)); 673 674 /* write block size option to the viewer's .info file */ 675 PetscCall(MatView_Binary_BlockSizes(mat, viewer)); 676 PetscFunctionReturn(0); 677 } 678 679 static PetscErrorCode MatView_SeqAIJ_ASCII_structonly(Mat A, PetscViewer viewer) 680 { 681 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 682 PetscInt i, k, m = A->rmap->N; 683 684 PetscFunctionBegin; 685 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE)); 686 for (i = 0; i < m; i++) { 687 PetscCall(PetscViewerASCIIPrintf(viewer, "row %" PetscInt_FMT ":", i)); 688 for (k = a->i[i]; k < a->i[i + 1]; k++) PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ") ", a->j[k])); 689 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 690 } 691 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE)); 692 PetscFunctionReturn(0); 693 } 694 695 extern PetscErrorCode MatSeqAIJFactorInfo_Matlab(Mat, PetscViewer); 696 697 PetscErrorCode MatView_SeqAIJ_ASCII(Mat A, PetscViewer viewer) 698 { 699 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 700 const PetscScalar *av; 701 PetscInt i, j, m = A->rmap->n; 702 const char *name; 703 PetscViewerFormat format; 704 705 PetscFunctionBegin; 706 if (A->structure_only) { 707 PetscCall(MatView_SeqAIJ_ASCII_structonly(A, viewer)); 708 PetscFunctionReturn(0); 709 } 710 711 PetscCall(PetscViewerGetFormat(viewer, &format)); 712 if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscFunctionReturn(0); 713 714 /* trigger copy to CPU if needed */ 715 PetscCall(MatSeqAIJGetArrayRead(A, &av)); 716 PetscCall(MatSeqAIJRestoreArrayRead(A, &av)); 717 if (format == PETSC_VIEWER_ASCII_MATLAB) { 718 PetscInt nofinalvalue = 0; 719 if (m && ((a->i[m] == a->i[m - 1]) || (a->j[a->nz - 1] != A->cmap->n - 1))) { 720 /* Need a dummy value to ensure the dimension of the matrix. */ 721 nofinalvalue = 1; 722 } 723 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE)); 724 PetscCall(PetscViewerASCIIPrintf(viewer, "%% Size = %" PetscInt_FMT " %" PetscInt_FMT " \n", m, A->cmap->n)); 725 PetscCall(PetscViewerASCIIPrintf(viewer, "%% Nonzeros = %" PetscInt_FMT " \n", a->nz)); 726 #if defined(PETSC_USE_COMPLEX) 727 PetscCall(PetscViewerASCIIPrintf(viewer, "zzz = zeros(%" PetscInt_FMT ",4);\n", a->nz + nofinalvalue)); 728 #else 729 PetscCall(PetscViewerASCIIPrintf(viewer, "zzz = zeros(%" PetscInt_FMT ",3);\n", a->nz + nofinalvalue)); 730 #endif 731 PetscCall(PetscViewerASCIIPrintf(viewer, "zzz = [\n")); 732 733 for (i = 0; i < m; i++) { 734 for (j = a->i[i]; j < a->i[i + 1]; j++) { 735 #if defined(PETSC_USE_COMPLEX) 736 PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " %" PetscInt_FMT " %18.16e %18.16e\n", i + 1, a->j[j] + 1, (double)PetscRealPart(a->a[j]), (double)PetscImaginaryPart(a->a[j]))); 737 #else 738 PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n", i + 1, a->j[j] + 1, (double)a->a[j])); 739 #endif 740 } 741 } 742 if (nofinalvalue) { 743 #if defined(PETSC_USE_COMPLEX) 744 PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " %" PetscInt_FMT " %18.16e %18.16e\n", m, A->cmap->n, 0., 0.)); 745 #else 746 PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n", m, A->cmap->n, 0.0)); 747 #endif 748 } 749 PetscCall(PetscObjectGetName((PetscObject)A, &name)); 750 PetscCall(PetscViewerASCIIPrintf(viewer, "];\n %s = spconvert(zzz);\n", name)); 751 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE)); 752 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 753 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE)); 754 for (i = 0; i < m; i++) { 755 PetscCall(PetscViewerASCIIPrintf(viewer, "row %" PetscInt_FMT ":", i)); 756 for (j = a->i[i]; j < a->i[i + 1]; j++) { 757 #if defined(PETSC_USE_COMPLEX) 758 if (PetscImaginaryPart(a->a[j]) > 0.0 && PetscRealPart(a->a[j]) != 0.0) { 759 PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g + %g i)", a->j[j], (double)PetscRealPart(a->a[j]), (double)PetscImaginaryPart(a->a[j]))); 760 } else if (PetscImaginaryPart(a->a[j]) < 0.0 && PetscRealPart(a->a[j]) != 0.0) { 761 PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g - %g i)", a->j[j], (double)PetscRealPart(a->a[j]), (double)-PetscImaginaryPart(a->a[j]))); 762 } else if (PetscRealPart(a->a[j]) != 0.0) { 763 PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", a->j[j], (double)PetscRealPart(a->a[j]))); 764 } 765 #else 766 if (a->a[j] != 0.0) PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", a->j[j], (double)a->a[j])); 767 #endif 768 } 769 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 770 } 771 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE)); 772 } else if (format == PETSC_VIEWER_ASCII_SYMMODU) { 773 PetscInt nzd = 0, fshift = 1, *sptr; 774 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE)); 775 PetscCall(PetscMalloc1(m + 1, &sptr)); 776 for (i = 0; i < m; i++) { 777 sptr[i] = nzd + 1; 778 for (j = a->i[i]; j < a->i[i + 1]; j++) { 779 if (a->j[j] >= i) { 780 #if defined(PETSC_USE_COMPLEX) 781 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) nzd++; 782 #else 783 if (a->a[j] != 0.0) nzd++; 784 #endif 785 } 786 } 787 } 788 sptr[m] = nzd + 1; 789 PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT " %" PetscInt_FMT "\n\n", m, nzd)); 790 for (i = 0; i < m + 1; i += 6) { 791 if (i + 4 < m) { 792 PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", sptr[i], sptr[i + 1], sptr[i + 2], sptr[i + 3], sptr[i + 4], sptr[i + 5])); 793 } else if (i + 3 < m) { 794 PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", sptr[i], sptr[i + 1], sptr[i + 2], sptr[i + 3], sptr[i + 4])); 795 } else if (i + 2 < m) { 796 PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", sptr[i], sptr[i + 1], sptr[i + 2], sptr[i + 3])); 797 } else if (i + 1 < m) { 798 PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", sptr[i], sptr[i + 1], sptr[i + 2])); 799 } else if (i < m) { 800 PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT " %" PetscInt_FMT "\n", sptr[i], sptr[i + 1])); 801 } else { 802 PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT "\n", sptr[i])); 803 } 804 } 805 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 806 PetscCall(PetscFree(sptr)); 807 for (i = 0; i < m; i++) { 808 for (j = a->i[i]; j < a->i[i + 1]; j++) { 809 if (a->j[j] >= i) PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT " ", a->j[j] + fshift)); 810 } 811 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 812 } 813 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 814 for (i = 0; i < m; i++) { 815 for (j = a->i[i]; j < a->i[i + 1]; j++) { 816 if (a->j[j] >= i) { 817 #if defined(PETSC_USE_COMPLEX) 818 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) PetscCall(PetscViewerASCIIPrintf(viewer, " %18.16e %18.16e ", (double)PetscRealPart(a->a[j]), (double)PetscImaginaryPart(a->a[j]))); 819 #else 820 if (a->a[j] != 0.0) PetscCall(PetscViewerASCIIPrintf(viewer, " %18.16e ", (double)a->a[j])); 821 #endif 822 } 823 } 824 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 825 } 826 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE)); 827 } else if (format == PETSC_VIEWER_ASCII_DENSE) { 828 PetscInt cnt = 0, jcnt; 829 PetscScalar value; 830 #if defined(PETSC_USE_COMPLEX) 831 PetscBool realonly = PETSC_TRUE; 832 833 for (i = 0; i < a->i[m]; i++) { 834 if (PetscImaginaryPart(a->a[i]) != 0.0) { 835 realonly = PETSC_FALSE; 836 break; 837 } 838 } 839 #endif 840 841 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE)); 842 for (i = 0; i < m; i++) { 843 jcnt = 0; 844 for (j = 0; j < A->cmap->n; j++) { 845 if (jcnt < a->i[i + 1] - a->i[i] && j == a->j[cnt]) { 846 value = a->a[cnt++]; 847 jcnt++; 848 } else { 849 value = 0.0; 850 } 851 #if defined(PETSC_USE_COMPLEX) 852 if (realonly) { 853 PetscCall(PetscViewerASCIIPrintf(viewer, " %7.5e ", (double)PetscRealPart(value))); 854 } else { 855 PetscCall(PetscViewerASCIIPrintf(viewer, " %7.5e+%7.5e i ", (double)PetscRealPart(value), (double)PetscImaginaryPart(value))); 856 } 857 #else 858 PetscCall(PetscViewerASCIIPrintf(viewer, " %7.5e ", (double)value)); 859 #endif 860 } 861 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 862 } 863 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE)); 864 } else if (format == PETSC_VIEWER_ASCII_MATRIXMARKET) { 865 PetscInt fshift = 1; 866 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE)); 867 #if defined(PETSC_USE_COMPLEX) 868 PetscCall(PetscViewerASCIIPrintf(viewer, "%%%%MatrixMarket matrix coordinate complex general\n")); 869 #else 870 PetscCall(PetscViewerASCIIPrintf(viewer, "%%%%MatrixMarket matrix coordinate real general\n")); 871 #endif 872 PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", m, A->cmap->n, a->nz)); 873 for (i = 0; i < m; i++) { 874 for (j = a->i[i]; j < a->i[i + 1]; j++) { 875 #if defined(PETSC_USE_COMPLEX) 876 PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " %" PetscInt_FMT " %g %g\n", i + fshift, a->j[j] + fshift, (double)PetscRealPart(a->a[j]), (double)PetscImaginaryPart(a->a[j]))); 877 #else 878 PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " %" PetscInt_FMT " %g\n", i + fshift, a->j[j] + fshift, (double)a->a[j])); 879 #endif 880 } 881 } 882 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE)); 883 } else { 884 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE)); 885 if (A->factortype) { 886 for (i = 0; i < m; i++) { 887 PetscCall(PetscViewerASCIIPrintf(viewer, "row %" PetscInt_FMT ":", i)); 888 /* L part */ 889 for (j = a->i[i]; j < a->i[i + 1]; j++) { 890 #if defined(PETSC_USE_COMPLEX) 891 if (PetscImaginaryPart(a->a[j]) > 0.0) { 892 PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g + %g i)", a->j[j], (double)PetscRealPart(a->a[j]), (double)PetscImaginaryPart(a->a[j]))); 893 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 894 PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g - %g i)", a->j[j], (double)PetscRealPart(a->a[j]), (double)(-PetscImaginaryPart(a->a[j])))); 895 } else { 896 PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", a->j[j], (double)PetscRealPart(a->a[j]))); 897 } 898 #else 899 PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", a->j[j], (double)a->a[j])); 900 #endif 901 } 902 /* diagonal */ 903 j = a->diag[i]; 904 #if defined(PETSC_USE_COMPLEX) 905 if (PetscImaginaryPart(a->a[j]) > 0.0) { 906 PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g + %g i)", a->j[j], (double)PetscRealPart(1.0 / a->a[j]), (double)PetscImaginaryPart(1.0 / a->a[j]))); 907 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 908 PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g - %g i)", a->j[j], (double)PetscRealPart(1.0 / a->a[j]), (double)(-PetscImaginaryPart(1.0 / a->a[j])))); 909 } else { 910 PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", a->j[j], (double)PetscRealPart(1.0 / a->a[j]))); 911 } 912 #else 913 PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", a->j[j], (double)(1.0 / a->a[j]))); 914 #endif 915 916 /* U part */ 917 for (j = a->diag[i + 1] + 1; j < a->diag[i]; j++) { 918 #if defined(PETSC_USE_COMPLEX) 919 if (PetscImaginaryPart(a->a[j]) > 0.0) { 920 PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g + %g i)", a->j[j], (double)PetscRealPart(a->a[j]), (double)PetscImaginaryPart(a->a[j]))); 921 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 922 PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g - %g i)", a->j[j], (double)PetscRealPart(a->a[j]), (double)(-PetscImaginaryPart(a->a[j])))); 923 } else { 924 PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", a->j[j], (double)PetscRealPart(a->a[j]))); 925 } 926 #else 927 PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", a->j[j], (double)a->a[j])); 928 #endif 929 } 930 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 931 } 932 } else { 933 for (i = 0; i < m; i++) { 934 PetscCall(PetscViewerASCIIPrintf(viewer, "row %" PetscInt_FMT ":", i)); 935 for (j = a->i[i]; j < a->i[i + 1]; j++) { 936 #if defined(PETSC_USE_COMPLEX) 937 if (PetscImaginaryPart(a->a[j]) > 0.0) { 938 PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g + %g i)", a->j[j], (double)PetscRealPart(a->a[j]), (double)PetscImaginaryPart(a->a[j]))); 939 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 940 PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g - %g i)", a->j[j], (double)PetscRealPart(a->a[j]), (double)-PetscImaginaryPart(a->a[j]))); 941 } else { 942 PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", a->j[j], (double)PetscRealPart(a->a[j]))); 943 } 944 #else 945 PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", a->j[j], (double)a->a[j])); 946 #endif 947 } 948 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 949 } 950 } 951 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE)); 952 } 953 PetscCall(PetscViewerFlush(viewer)); 954 PetscFunctionReturn(0); 955 } 956 957 #include <petscdraw.h> 958 PetscErrorCode MatView_SeqAIJ_Draw_Zoom(PetscDraw draw, void *Aa) 959 { 960 Mat A = (Mat)Aa; 961 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 962 PetscInt i, j, m = A->rmap->n; 963 int color; 964 PetscReal xl, yl, xr, yr, x_l, x_r, y_l, y_r; 965 PetscViewer viewer; 966 PetscViewerFormat format; 967 const PetscScalar *aa; 968 969 PetscFunctionBegin; 970 PetscCall(PetscObjectQuery((PetscObject)A, "Zoomviewer", (PetscObject *)&viewer)); 971 PetscCall(PetscViewerGetFormat(viewer, &format)); 972 PetscCall(PetscDrawGetCoordinates(draw, &xl, &yl, &xr, &yr)); 973 974 /* loop over matrix elements drawing boxes */ 975 PetscCall(MatSeqAIJGetArrayRead(A, &aa)); 976 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 977 PetscDrawCollectiveBegin(draw); 978 /* Blue for negative, Cyan for zero and Red for positive */ 979 color = PETSC_DRAW_BLUE; 980 for (i = 0; i < m; i++) { 981 y_l = m - i - 1.0; 982 y_r = y_l + 1.0; 983 for (j = a->i[i]; j < a->i[i + 1]; j++) { 984 x_l = a->j[j]; 985 x_r = x_l + 1.0; 986 if (PetscRealPart(aa[j]) >= 0.) continue; 987 PetscCall(PetscDrawRectangle(draw, x_l, y_l, x_r, y_r, color, color, color, color)); 988 } 989 } 990 color = PETSC_DRAW_CYAN; 991 for (i = 0; i < m; i++) { 992 y_l = m - i - 1.0; 993 y_r = y_l + 1.0; 994 for (j = a->i[i]; j < a->i[i + 1]; j++) { 995 x_l = a->j[j]; 996 x_r = x_l + 1.0; 997 if (aa[j] != 0.) continue; 998 PetscCall(PetscDrawRectangle(draw, x_l, y_l, x_r, y_r, color, color, color, color)); 999 } 1000 } 1001 color = PETSC_DRAW_RED; 1002 for (i = 0; i < m; i++) { 1003 y_l = m - i - 1.0; 1004 y_r = y_l + 1.0; 1005 for (j = a->i[i]; j < a->i[i + 1]; j++) { 1006 x_l = a->j[j]; 1007 x_r = x_l + 1.0; 1008 if (PetscRealPart(aa[j]) <= 0.) continue; 1009 PetscCall(PetscDrawRectangle(draw, x_l, y_l, x_r, y_r, color, color, color, color)); 1010 } 1011 } 1012 PetscDrawCollectiveEnd(draw); 1013 } else { 1014 /* use contour shading to indicate magnitude of values */ 1015 /* first determine max of all nonzero values */ 1016 PetscReal minv = 0.0, maxv = 0.0; 1017 PetscInt nz = a->nz, count = 0; 1018 PetscDraw popup; 1019 1020 for (i = 0; i < nz; i++) { 1021 if (PetscAbsScalar(aa[i]) > maxv) maxv = PetscAbsScalar(aa[i]); 1022 } 1023 if (minv >= maxv) maxv = minv + PETSC_SMALL; 1024 PetscCall(PetscDrawGetPopup(draw, &popup)); 1025 PetscCall(PetscDrawScalePopup(popup, minv, maxv)); 1026 1027 PetscDrawCollectiveBegin(draw); 1028 for (i = 0; i < m; i++) { 1029 y_l = m - i - 1.0; 1030 y_r = y_l + 1.0; 1031 for (j = a->i[i]; j < a->i[i + 1]; j++) { 1032 x_l = a->j[j]; 1033 x_r = x_l + 1.0; 1034 color = PetscDrawRealToColor(PetscAbsScalar(aa[count]), minv, maxv); 1035 PetscCall(PetscDrawRectangle(draw, x_l, y_l, x_r, y_r, color, color, color, color)); 1036 count++; 1037 } 1038 } 1039 PetscDrawCollectiveEnd(draw); 1040 } 1041 PetscCall(MatSeqAIJRestoreArrayRead(A, &aa)); 1042 PetscFunctionReturn(0); 1043 } 1044 1045 #include <petscdraw.h> 1046 PetscErrorCode MatView_SeqAIJ_Draw(Mat A, PetscViewer viewer) 1047 { 1048 PetscDraw draw; 1049 PetscReal xr, yr, xl, yl, h, w; 1050 PetscBool isnull; 1051 1052 PetscFunctionBegin; 1053 PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw)); 1054 PetscCall(PetscDrawIsNull(draw, &isnull)); 1055 if (isnull) PetscFunctionReturn(0); 1056 1057 xr = A->cmap->n; 1058 yr = A->rmap->n; 1059 h = yr / 10.0; 1060 w = xr / 10.0; 1061 xr += w; 1062 yr += h; 1063 xl = -w; 1064 yl = -h; 1065 PetscCall(PetscDrawSetCoordinates(draw, xl, yl, xr, yr)); 1066 PetscCall(PetscObjectCompose((PetscObject)A, "Zoomviewer", (PetscObject)viewer)); 1067 PetscCall(PetscDrawZoom(draw, MatView_SeqAIJ_Draw_Zoom, A)); 1068 PetscCall(PetscObjectCompose((PetscObject)A, "Zoomviewer", NULL)); 1069 PetscCall(PetscDrawSave(draw)); 1070 PetscFunctionReturn(0); 1071 } 1072 1073 PetscErrorCode MatView_SeqAIJ(Mat A, PetscViewer viewer) 1074 { 1075 PetscBool iascii, isbinary, isdraw; 1076 1077 PetscFunctionBegin; 1078 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 1079 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary)); 1080 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw)); 1081 if (iascii) PetscCall(MatView_SeqAIJ_ASCII(A, viewer)); 1082 else if (isbinary) PetscCall(MatView_SeqAIJ_Binary(A, viewer)); 1083 else if (isdraw) PetscCall(MatView_SeqAIJ_Draw(A, viewer)); 1084 PetscCall(MatView_SeqAIJ_Inode(A, viewer)); 1085 PetscFunctionReturn(0); 1086 } 1087 1088 PetscErrorCode MatAssemblyEnd_SeqAIJ(Mat A, MatAssemblyType mode) 1089 { 1090 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 1091 PetscInt fshift = 0, i, *ai = a->i, *aj = a->j, *imax = a->imax; 1092 PetscInt m = A->rmap->n, *ip, N, *ailen = a->ilen, rmax = 0; 1093 MatScalar *aa = a->a, *ap; 1094 PetscReal ratio = 0.6; 1095 1096 PetscFunctionBegin; 1097 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 1098 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 1099 if (A->was_assembled && A->ass_nonzerostate == A->nonzerostate) { 1100 /* we need to respect users asking to use or not the inodes routine in between matrix assemblies */ 1101 PetscCall(MatAssemblyEnd_SeqAIJ_Inode(A, mode)); 1102 PetscFunctionReturn(0); 1103 } 1104 1105 if (m) rmax = ailen[0]; /* determine row with most nonzeros */ 1106 for (i = 1; i < m; i++) { 1107 /* move each row back by the amount of empty slots (fshift) before it*/ 1108 fshift += imax[i - 1] - ailen[i - 1]; 1109 rmax = PetscMax(rmax, ailen[i]); 1110 if (fshift) { 1111 ip = aj + ai[i]; 1112 ap = aa + ai[i]; 1113 N = ailen[i]; 1114 PetscCall(PetscArraymove(ip - fshift, ip, N)); 1115 if (!A->structure_only) PetscCall(PetscArraymove(ap - fshift, ap, N)); 1116 } 1117 ai[i] = ai[i - 1] + ailen[i - 1]; 1118 } 1119 if (m) { 1120 fshift += imax[m - 1] - ailen[m - 1]; 1121 ai[m] = ai[m - 1] + ailen[m - 1]; 1122 } 1123 /* reset ilen and imax for each row */ 1124 a->nonzerorowcnt = 0; 1125 if (A->structure_only) { 1126 PetscCall(PetscFree(a->imax)); 1127 PetscCall(PetscFree(a->ilen)); 1128 } else { /* !A->structure_only */ 1129 for (i = 0; i < m; i++) { 1130 ailen[i] = imax[i] = ai[i + 1] - ai[i]; 1131 a->nonzerorowcnt += ((ai[i + 1] - ai[i]) > 0); 1132 } 1133 } 1134 a->nz = ai[m]; 1135 PetscCheck(!fshift || a->nounused != -1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unused space detected in matrix: %" PetscInt_FMT " X %" PetscInt_FMT ", %" PetscInt_FMT " unneeded", m, A->cmap->n, fshift); 1136 1137 PetscCall(MatMarkDiagonal_SeqAIJ(A)); 1138 PetscCall(PetscInfo(A, "Matrix size: %" PetscInt_FMT " X %" PetscInt_FMT "; storage space: %" PetscInt_FMT " unneeded,%" PetscInt_FMT " used\n", m, A->cmap->n, fshift, a->nz)); 1139 PetscCall(PetscInfo(A, "Number of mallocs during MatSetValues() is %" PetscInt_FMT "\n", a->reallocs)); 1140 PetscCall(PetscInfo(A, "Maximum nonzeros in any row is %" PetscInt_FMT "\n", rmax)); 1141 1142 A->info.mallocs += a->reallocs; 1143 a->reallocs = 0; 1144 A->info.nz_unneeded = (PetscReal)fshift; 1145 a->rmax = rmax; 1146 1147 if (!A->structure_only) PetscCall(MatCheckCompressedRow(A, a->nonzerorowcnt, &a->compressedrow, a->i, m, ratio)); 1148 PetscCall(MatAssemblyEnd_SeqAIJ_Inode(A, mode)); 1149 PetscFunctionReturn(0); 1150 } 1151 1152 PetscErrorCode MatRealPart_SeqAIJ(Mat A) 1153 { 1154 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 1155 PetscInt i, nz = a->nz; 1156 MatScalar *aa; 1157 1158 PetscFunctionBegin; 1159 PetscCall(MatSeqAIJGetArray(A, &aa)); 1160 for (i = 0; i < nz; i++) aa[i] = PetscRealPart(aa[i]); 1161 PetscCall(MatSeqAIJRestoreArray(A, &aa)); 1162 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 1163 PetscFunctionReturn(0); 1164 } 1165 1166 PetscErrorCode MatImaginaryPart_SeqAIJ(Mat A) 1167 { 1168 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 1169 PetscInt i, nz = a->nz; 1170 MatScalar *aa; 1171 1172 PetscFunctionBegin; 1173 PetscCall(MatSeqAIJGetArray(A, &aa)); 1174 for (i = 0; i < nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1175 PetscCall(MatSeqAIJRestoreArray(A, &aa)); 1176 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 1177 PetscFunctionReturn(0); 1178 } 1179 1180 PetscErrorCode MatZeroEntries_SeqAIJ(Mat A) 1181 { 1182 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 1183 MatScalar *aa; 1184 1185 PetscFunctionBegin; 1186 PetscCall(MatSeqAIJGetArrayWrite(A, &aa)); 1187 PetscCall(PetscArrayzero(aa, a->i[A->rmap->n])); 1188 PetscCall(MatSeqAIJRestoreArrayWrite(A, &aa)); 1189 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 1190 PetscFunctionReturn(0); 1191 } 1192 1193 PETSC_INTERN PetscErrorCode MatResetPreallocationCOO_SeqAIJ(Mat A) 1194 { 1195 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 1196 1197 PetscFunctionBegin; 1198 PetscCall(PetscFree(a->perm)); 1199 PetscCall(PetscFree(a->jmap)); 1200 PetscFunctionReturn(0); 1201 } 1202 1203 PetscErrorCode MatDestroy_SeqAIJ(Mat A) 1204 { 1205 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 1206 1207 PetscFunctionBegin; 1208 #if defined(PETSC_USE_LOG) 1209 PetscLogObjectState((PetscObject)A, "Rows=%" PetscInt_FMT ", Cols=%" PetscInt_FMT ", NZ=%" PetscInt_FMT, A->rmap->n, A->cmap->n, a->nz); 1210 #endif 1211 PetscCall(MatSeqXAIJFreeAIJ(A, &a->a, &a->j, &a->i)); 1212 PetscCall(MatResetPreallocationCOO_SeqAIJ(A)); 1213 PetscCall(ISDestroy(&a->row)); 1214 PetscCall(ISDestroy(&a->col)); 1215 PetscCall(PetscFree(a->diag)); 1216 PetscCall(PetscFree(a->ibdiag)); 1217 PetscCall(PetscFree(a->imax)); 1218 PetscCall(PetscFree(a->ilen)); 1219 PetscCall(PetscFree(a->ipre)); 1220 PetscCall(PetscFree3(a->idiag, a->mdiag, a->ssor_work)); 1221 PetscCall(PetscFree(a->solve_work)); 1222 PetscCall(ISDestroy(&a->icol)); 1223 PetscCall(PetscFree(a->saved_values)); 1224 PetscCall(PetscFree2(a->compressedrow.i, a->compressedrow.rindex)); 1225 PetscCall(MatDestroy_SeqAIJ_Inode(A)); 1226 PetscCall(PetscFree(A->data)); 1227 1228 /* MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted may allocate this. 1229 That function is so heavily used (sometimes in an hidden way through multnumeric function pointers) 1230 that is hard to properly add this data to the MatProduct data. We free it here to avoid 1231 users reusing the matrix object with different data to incur in obscure segmentation faults 1232 due to different matrix sizes */ 1233 PetscCall(PetscObjectCompose((PetscObject)A, "__PETSc__ab_dense", NULL)); 1234 1235 PetscCall(PetscObjectChangeTypeName((PetscObject)A, NULL)); 1236 PetscCall(PetscObjectComposeFunction((PetscObject)A, "PetscMatlabEnginePut_C", NULL)); 1237 PetscCall(PetscObjectComposeFunction((PetscObject)A, "PetscMatlabEngineGet_C", NULL)); 1238 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSeqAIJSetColumnIndices_C", NULL)); 1239 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatStoreValues_C", NULL)); 1240 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatRetrieveValues_C", NULL)); 1241 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqaij_seqsbaij_C", NULL)); 1242 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqaij_seqbaij_C", NULL)); 1243 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqaij_seqaijperm_C", NULL)); 1244 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqaij_seqaijsell_C", NULL)); 1245 #if defined(PETSC_HAVE_MKL_SPARSE) 1246 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqaij_seqaijmkl_C", NULL)); 1247 #endif 1248 #if defined(PETSC_HAVE_CUDA) 1249 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqaij_seqaijcusparse_C", NULL)); 1250 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_seqaijcusparse_seqaij_C", NULL)); 1251 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_seqaij_seqaijcusparse_C", NULL)); 1252 #endif 1253 #if defined(PETSC_HAVE_HIP) 1254 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqaij_seqaijhipsparse_C", NULL)); 1255 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_seqaijhipsparse_seqaij_C", NULL)); 1256 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_seqaij_seqaijhipsparse_C", NULL)); 1257 #endif 1258 #if defined(PETSC_HAVE_KOKKOS_KERNELS) 1259 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqaij_seqaijkokkos_C", NULL)); 1260 #endif 1261 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqaij_seqaijcrl_C", NULL)); 1262 #if defined(PETSC_HAVE_ELEMENTAL) 1263 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqaij_elemental_C", NULL)); 1264 #endif 1265 #if defined(PETSC_HAVE_SCALAPACK) 1266 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqaij_scalapack_C", NULL)); 1267 #endif 1268 #if defined(PETSC_HAVE_HYPRE) 1269 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqaij_hypre_C", NULL)); 1270 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_transpose_seqaij_seqaij_C", NULL)); 1271 #endif 1272 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqaij_seqdense_C", NULL)); 1273 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqaij_seqsell_C", NULL)); 1274 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqaij_is_C", NULL)); 1275 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatIsTranspose_C", NULL)); 1276 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatIsHermitianTranspose_C", NULL)); 1277 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSeqAIJSetPreallocation_C", NULL)); 1278 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatResetPreallocation_C", NULL)); 1279 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSeqAIJSetPreallocationCSR_C", NULL)); 1280 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatReorderForNonzeroDiagonal_C", NULL)); 1281 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_is_seqaij_C", NULL)); 1282 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_seqdense_seqaij_C", NULL)); 1283 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_seqaij_seqaij_C", NULL)); 1284 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSeqAIJKron_C", NULL)); 1285 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOO_C", NULL)); 1286 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetValuesCOO_C", NULL)); 1287 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatFactorGetSolverType_C", NULL)); 1288 /* these calls do not belong here: the subclasses Duplicate/Destroy are wrong */ 1289 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqaijsell_seqaij_C", NULL)); 1290 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqaijperm_seqaij_C", NULL)); 1291 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqaij_seqaijviennacl_C", NULL)); 1292 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_seqaijviennacl_seqdense_C", NULL)); 1293 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_seqaijviennacl_seqaij_C", NULL)); 1294 PetscFunctionReturn(0); 1295 } 1296 1297 PetscErrorCode MatSetOption_SeqAIJ(Mat A, MatOption op, PetscBool flg) 1298 { 1299 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 1300 1301 PetscFunctionBegin; 1302 switch (op) { 1303 case MAT_ROW_ORIENTED: 1304 a->roworiented = flg; 1305 break; 1306 case MAT_KEEP_NONZERO_PATTERN: 1307 a->keepnonzeropattern = flg; 1308 break; 1309 case MAT_NEW_NONZERO_LOCATIONS: 1310 a->nonew = (flg ? 0 : 1); 1311 break; 1312 case MAT_NEW_NONZERO_LOCATION_ERR: 1313 a->nonew = (flg ? -1 : 0); 1314 break; 1315 case MAT_NEW_NONZERO_ALLOCATION_ERR: 1316 a->nonew = (flg ? -2 : 0); 1317 break; 1318 case MAT_UNUSED_NONZERO_LOCATION_ERR: 1319 a->nounused = (flg ? -1 : 0); 1320 break; 1321 case MAT_IGNORE_ZERO_ENTRIES: 1322 a->ignorezeroentries = flg; 1323 break; 1324 case MAT_SPD: 1325 case MAT_SYMMETRIC: 1326 case MAT_STRUCTURALLY_SYMMETRIC: 1327 case MAT_HERMITIAN: 1328 case MAT_SYMMETRY_ETERNAL: 1329 case MAT_STRUCTURE_ONLY: 1330 case MAT_STRUCTURAL_SYMMETRY_ETERNAL: 1331 case MAT_SPD_ETERNAL: 1332 /* if the diagonal matrix is square it inherits some of the properties above */ 1333 break; 1334 case MAT_FORCE_DIAGONAL_ENTRIES: 1335 case MAT_IGNORE_OFF_PROC_ENTRIES: 1336 case MAT_USE_HASH_TABLE: 1337 PetscCall(PetscInfo(A, "Option %s ignored\n", MatOptions[op])); 1338 break; 1339 case MAT_USE_INODES: 1340 PetscCall(MatSetOption_SeqAIJ_Inode(A, MAT_USE_INODES, flg)); 1341 break; 1342 case MAT_SUBMAT_SINGLEIS: 1343 A->submat_singleis = flg; 1344 break; 1345 case MAT_SORTED_FULL: 1346 if (flg) A->ops->setvalues = MatSetValues_SeqAIJ_SortedFull; 1347 else A->ops->setvalues = MatSetValues_SeqAIJ; 1348 break; 1349 case MAT_FORM_EXPLICIT_TRANSPOSE: 1350 A->form_explicit_transpose = flg; 1351 break; 1352 default: 1353 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "unknown option %d", op); 1354 } 1355 PetscFunctionReturn(0); 1356 } 1357 1358 PetscErrorCode MatGetDiagonal_SeqAIJ(Mat A, Vec v) 1359 { 1360 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 1361 PetscInt i, j, n, *ai = a->i, *aj = a->j; 1362 PetscScalar *x; 1363 const PetscScalar *aa; 1364 1365 PetscFunctionBegin; 1366 PetscCall(VecGetLocalSize(v, &n)); 1367 PetscCheck(n == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector"); 1368 PetscCall(MatSeqAIJGetArrayRead(A, &aa)); 1369 if (A->factortype == MAT_FACTOR_ILU || A->factortype == MAT_FACTOR_LU) { 1370 PetscInt *diag = a->diag; 1371 PetscCall(VecGetArrayWrite(v, &x)); 1372 for (i = 0; i < n; i++) x[i] = 1.0 / aa[diag[i]]; 1373 PetscCall(VecRestoreArrayWrite(v, &x)); 1374 PetscCall(MatSeqAIJRestoreArrayRead(A, &aa)); 1375 PetscFunctionReturn(0); 1376 } 1377 1378 PetscCall(VecGetArrayWrite(v, &x)); 1379 for (i = 0; i < n; i++) { 1380 x[i] = 0.0; 1381 for (j = ai[i]; j < ai[i + 1]; j++) { 1382 if (aj[j] == i) { 1383 x[i] = aa[j]; 1384 break; 1385 } 1386 } 1387 } 1388 PetscCall(VecRestoreArrayWrite(v, &x)); 1389 PetscCall(MatSeqAIJRestoreArrayRead(A, &aa)); 1390 PetscFunctionReturn(0); 1391 } 1392 1393 #include <../src/mat/impls/aij/seq/ftn-kernels/fmult.h> 1394 PetscErrorCode MatMultTransposeAdd_SeqAIJ(Mat A, Vec xx, Vec zz, Vec yy) 1395 { 1396 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 1397 const MatScalar *aa; 1398 PetscScalar *y; 1399 const PetscScalar *x; 1400 PetscInt m = A->rmap->n; 1401 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 1402 const MatScalar *v; 1403 PetscScalar alpha; 1404 PetscInt n, i, j; 1405 const PetscInt *idx, *ii, *ridx = NULL; 1406 Mat_CompressedRow cprow = a->compressedrow; 1407 PetscBool usecprow = cprow.use; 1408 #endif 1409 1410 PetscFunctionBegin; 1411 if (zz != yy) PetscCall(VecCopy(zz, yy)); 1412 PetscCall(VecGetArrayRead(xx, &x)); 1413 PetscCall(VecGetArray(yy, &y)); 1414 PetscCall(MatSeqAIJGetArrayRead(A, &aa)); 1415 1416 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 1417 fortranmulttransposeaddaij_(&m, x, a->i, a->j, aa, y); 1418 #else 1419 if (usecprow) { 1420 m = cprow.nrows; 1421 ii = cprow.i; 1422 ridx = cprow.rindex; 1423 } else { 1424 ii = a->i; 1425 } 1426 for (i = 0; i < m; i++) { 1427 idx = a->j + ii[i]; 1428 v = aa + ii[i]; 1429 n = ii[i + 1] - ii[i]; 1430 if (usecprow) { 1431 alpha = x[ridx[i]]; 1432 } else { 1433 alpha = x[i]; 1434 } 1435 for (j = 0; j < n; j++) y[idx[j]] += alpha * v[j]; 1436 } 1437 #endif 1438 PetscCall(PetscLogFlops(2.0 * a->nz)); 1439 PetscCall(VecRestoreArrayRead(xx, &x)); 1440 PetscCall(VecRestoreArray(yy, &y)); 1441 PetscCall(MatSeqAIJRestoreArrayRead(A, &aa)); 1442 PetscFunctionReturn(0); 1443 } 1444 1445 PetscErrorCode MatMultTranspose_SeqAIJ(Mat A, Vec xx, Vec yy) 1446 { 1447 PetscFunctionBegin; 1448 PetscCall(VecSet(yy, 0.0)); 1449 PetscCall(MatMultTransposeAdd_SeqAIJ(A, xx, yy, yy)); 1450 PetscFunctionReturn(0); 1451 } 1452 1453 #include <../src/mat/impls/aij/seq/ftn-kernels/fmult.h> 1454 1455 PetscErrorCode MatMult_SeqAIJ(Mat A, Vec xx, Vec yy) 1456 { 1457 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 1458 PetscScalar *y; 1459 const PetscScalar *x; 1460 const MatScalar *aa, *a_a; 1461 PetscInt m = A->rmap->n; 1462 const PetscInt *aj, *ii, *ridx = NULL; 1463 PetscInt n, i; 1464 PetscScalar sum; 1465 PetscBool usecprow = a->compressedrow.use; 1466 1467 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 1468 #pragma disjoint(*x, *y, *aa) 1469 #endif 1470 1471 PetscFunctionBegin; 1472 if (a->inode.use && a->inode.checked) { 1473 PetscCall(MatMult_SeqAIJ_Inode(A, xx, yy)); 1474 PetscFunctionReturn(0); 1475 } 1476 PetscCall(MatSeqAIJGetArrayRead(A, &a_a)); 1477 PetscCall(VecGetArrayRead(xx, &x)); 1478 PetscCall(VecGetArray(yy, &y)); 1479 ii = a->i; 1480 if (usecprow) { /* use compressed row format */ 1481 PetscCall(PetscArrayzero(y, m)); 1482 m = a->compressedrow.nrows; 1483 ii = a->compressedrow.i; 1484 ridx = a->compressedrow.rindex; 1485 for (i = 0; i < m; i++) { 1486 n = ii[i + 1] - ii[i]; 1487 aj = a->j + ii[i]; 1488 aa = a_a + ii[i]; 1489 sum = 0.0; 1490 PetscSparseDensePlusDot(sum, x, aa, aj, n); 1491 /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */ 1492 y[*ridx++] = sum; 1493 } 1494 } else { /* do not use compressed row format */ 1495 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 1496 aj = a->j; 1497 aa = a_a; 1498 fortranmultaij_(&m, x, ii, aj, aa, y); 1499 #else 1500 for (i = 0; i < m; i++) { 1501 n = ii[i + 1] - ii[i]; 1502 aj = a->j + ii[i]; 1503 aa = a_a + ii[i]; 1504 sum = 0.0; 1505 PetscSparseDensePlusDot(sum, x, aa, aj, n); 1506 y[i] = sum; 1507 } 1508 #endif 1509 } 1510 PetscCall(PetscLogFlops(2.0 * a->nz - a->nonzerorowcnt)); 1511 PetscCall(VecRestoreArrayRead(xx, &x)); 1512 PetscCall(VecRestoreArray(yy, &y)); 1513 PetscCall(MatSeqAIJRestoreArrayRead(A, &a_a)); 1514 PetscFunctionReturn(0); 1515 } 1516 1517 PetscErrorCode MatMultMax_SeqAIJ(Mat A, Vec xx, Vec yy) 1518 { 1519 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 1520 PetscScalar *y; 1521 const PetscScalar *x; 1522 const MatScalar *aa, *a_a; 1523 PetscInt m = A->rmap->n; 1524 const PetscInt *aj, *ii, *ridx = NULL; 1525 PetscInt n, i, nonzerorow = 0; 1526 PetscScalar sum; 1527 PetscBool usecprow = a->compressedrow.use; 1528 1529 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 1530 #pragma disjoint(*x, *y, *aa) 1531 #endif 1532 1533 PetscFunctionBegin; 1534 PetscCall(MatSeqAIJGetArrayRead(A, &a_a)); 1535 PetscCall(VecGetArrayRead(xx, &x)); 1536 PetscCall(VecGetArray(yy, &y)); 1537 if (usecprow) { /* use compressed row format */ 1538 m = a->compressedrow.nrows; 1539 ii = a->compressedrow.i; 1540 ridx = a->compressedrow.rindex; 1541 for (i = 0; i < m; i++) { 1542 n = ii[i + 1] - ii[i]; 1543 aj = a->j + ii[i]; 1544 aa = a_a + ii[i]; 1545 sum = 0.0; 1546 nonzerorow += (n > 0); 1547 PetscSparseDenseMaxDot(sum, x, aa, aj, n); 1548 /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */ 1549 y[*ridx++] = sum; 1550 } 1551 } else { /* do not use compressed row format */ 1552 ii = a->i; 1553 for (i = 0; i < m; i++) { 1554 n = ii[i + 1] - ii[i]; 1555 aj = a->j + ii[i]; 1556 aa = a_a + ii[i]; 1557 sum = 0.0; 1558 nonzerorow += (n > 0); 1559 PetscSparseDenseMaxDot(sum, x, aa, aj, n); 1560 y[i] = sum; 1561 } 1562 } 1563 PetscCall(PetscLogFlops(2.0 * a->nz - nonzerorow)); 1564 PetscCall(VecRestoreArrayRead(xx, &x)); 1565 PetscCall(VecRestoreArray(yy, &y)); 1566 PetscCall(MatSeqAIJRestoreArrayRead(A, &a_a)); 1567 PetscFunctionReturn(0); 1568 } 1569 1570 PetscErrorCode MatMultAddMax_SeqAIJ(Mat A, Vec xx, Vec yy, Vec zz) 1571 { 1572 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 1573 PetscScalar *y, *z; 1574 const PetscScalar *x; 1575 const MatScalar *aa, *a_a; 1576 PetscInt m = A->rmap->n, *aj, *ii; 1577 PetscInt n, i, *ridx = NULL; 1578 PetscScalar sum; 1579 PetscBool usecprow = a->compressedrow.use; 1580 1581 PetscFunctionBegin; 1582 PetscCall(MatSeqAIJGetArrayRead(A, &a_a)); 1583 PetscCall(VecGetArrayRead(xx, &x)); 1584 PetscCall(VecGetArrayPair(yy, zz, &y, &z)); 1585 if (usecprow) { /* use compressed row format */ 1586 if (zz != yy) PetscCall(PetscArraycpy(z, y, m)); 1587 m = a->compressedrow.nrows; 1588 ii = a->compressedrow.i; 1589 ridx = a->compressedrow.rindex; 1590 for (i = 0; i < m; i++) { 1591 n = ii[i + 1] - ii[i]; 1592 aj = a->j + ii[i]; 1593 aa = a_a + ii[i]; 1594 sum = y[*ridx]; 1595 PetscSparseDenseMaxDot(sum, x, aa, aj, n); 1596 z[*ridx++] = sum; 1597 } 1598 } else { /* do not use compressed row format */ 1599 ii = a->i; 1600 for (i = 0; i < m; i++) { 1601 n = ii[i + 1] - ii[i]; 1602 aj = a->j + ii[i]; 1603 aa = a_a + ii[i]; 1604 sum = y[i]; 1605 PetscSparseDenseMaxDot(sum, x, aa, aj, n); 1606 z[i] = sum; 1607 } 1608 } 1609 PetscCall(PetscLogFlops(2.0 * a->nz)); 1610 PetscCall(VecRestoreArrayRead(xx, &x)); 1611 PetscCall(VecRestoreArrayPair(yy, zz, &y, &z)); 1612 PetscCall(MatSeqAIJRestoreArrayRead(A, &a_a)); 1613 PetscFunctionReturn(0); 1614 } 1615 1616 #include <../src/mat/impls/aij/seq/ftn-kernels/fmultadd.h> 1617 PetscErrorCode MatMultAdd_SeqAIJ(Mat A, Vec xx, Vec yy, Vec zz) 1618 { 1619 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 1620 PetscScalar *y, *z; 1621 const PetscScalar *x; 1622 const MatScalar *aa, *a_a; 1623 const PetscInt *aj, *ii, *ridx = NULL; 1624 PetscInt m = A->rmap->n, n, i; 1625 PetscScalar sum; 1626 PetscBool usecprow = a->compressedrow.use; 1627 1628 PetscFunctionBegin; 1629 if (a->inode.use && a->inode.checked) { 1630 PetscCall(MatMultAdd_SeqAIJ_Inode(A, xx, yy, zz)); 1631 PetscFunctionReturn(0); 1632 } 1633 PetscCall(MatSeqAIJGetArrayRead(A, &a_a)); 1634 PetscCall(VecGetArrayRead(xx, &x)); 1635 PetscCall(VecGetArrayPair(yy, zz, &y, &z)); 1636 if (usecprow) { /* use compressed row format */ 1637 if (zz != yy) PetscCall(PetscArraycpy(z, y, m)); 1638 m = a->compressedrow.nrows; 1639 ii = a->compressedrow.i; 1640 ridx = a->compressedrow.rindex; 1641 for (i = 0; i < m; i++) { 1642 n = ii[i + 1] - ii[i]; 1643 aj = a->j + ii[i]; 1644 aa = a_a + ii[i]; 1645 sum = y[*ridx]; 1646 PetscSparseDensePlusDot(sum, x, aa, aj, n); 1647 z[*ridx++] = sum; 1648 } 1649 } else { /* do not use compressed row format */ 1650 ii = a->i; 1651 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 1652 aj = a->j; 1653 aa = a_a; 1654 fortranmultaddaij_(&m, x, ii, aj, aa, y, z); 1655 #else 1656 for (i = 0; i < m; i++) { 1657 n = ii[i + 1] - ii[i]; 1658 aj = a->j + ii[i]; 1659 aa = a_a + ii[i]; 1660 sum = y[i]; 1661 PetscSparseDensePlusDot(sum, x, aa, aj, n); 1662 z[i] = sum; 1663 } 1664 #endif 1665 } 1666 PetscCall(PetscLogFlops(2.0 * a->nz)); 1667 PetscCall(VecRestoreArrayRead(xx, &x)); 1668 PetscCall(VecRestoreArrayPair(yy, zz, &y, &z)); 1669 PetscCall(MatSeqAIJRestoreArrayRead(A, &a_a)); 1670 PetscFunctionReturn(0); 1671 } 1672 1673 /* 1674 Adds diagonal pointers to sparse matrix structure. 1675 */ 1676 PetscErrorCode MatMarkDiagonal_SeqAIJ(Mat A) 1677 { 1678 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 1679 PetscInt i, j, m = A->rmap->n; 1680 PetscBool alreadySet = PETSC_TRUE; 1681 1682 PetscFunctionBegin; 1683 if (!a->diag) { 1684 PetscCall(PetscMalloc1(m, &a->diag)); 1685 alreadySet = PETSC_FALSE; 1686 } 1687 for (i = 0; i < A->rmap->n; i++) { 1688 /* If A's diagonal is already correctly set, this fast track enables cheap and repeated MatMarkDiagonal_SeqAIJ() calls */ 1689 if (alreadySet) { 1690 PetscInt pos = a->diag[i]; 1691 if (pos >= a->i[i] && pos < a->i[i + 1] && a->j[pos] == i) continue; 1692 } 1693 1694 a->diag[i] = a->i[i + 1]; 1695 for (j = a->i[i]; j < a->i[i + 1]; j++) { 1696 if (a->j[j] == i) { 1697 a->diag[i] = j; 1698 break; 1699 } 1700 } 1701 } 1702 PetscFunctionReturn(0); 1703 } 1704 1705 PetscErrorCode MatShift_SeqAIJ(Mat A, PetscScalar v) 1706 { 1707 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 1708 const PetscInt *diag = (const PetscInt *)a->diag; 1709 const PetscInt *ii = (const PetscInt *)a->i; 1710 PetscInt i, *mdiag = NULL; 1711 PetscInt cnt = 0; /* how many diagonals are missing */ 1712 1713 PetscFunctionBegin; 1714 if (!A->preallocated || !a->nz) { 1715 PetscCall(MatSeqAIJSetPreallocation(A, 1, NULL)); 1716 PetscCall(MatShift_Basic(A, v)); 1717 PetscFunctionReturn(0); 1718 } 1719 1720 if (a->diagonaldense) { 1721 cnt = 0; 1722 } else { 1723 PetscCall(PetscCalloc1(A->rmap->n, &mdiag)); 1724 for (i = 0; i < A->rmap->n; i++) { 1725 if (i < A->cmap->n && diag[i] >= ii[i + 1]) { /* 'out of range' rows never have diagonals */ 1726 cnt++; 1727 mdiag[i] = 1; 1728 } 1729 } 1730 } 1731 if (!cnt) { 1732 PetscCall(MatShift_Basic(A, v)); 1733 } else { 1734 PetscScalar *olda = a->a; /* preserve pointers to current matrix nonzeros structure and values */ 1735 PetscInt *oldj = a->j, *oldi = a->i; 1736 PetscBool singlemalloc = a->singlemalloc, free_a = a->free_a, free_ij = a->free_ij; 1737 1738 a->a = NULL; 1739 a->j = NULL; 1740 a->i = NULL; 1741 /* increase the values in imax for each row where a diagonal is being inserted then reallocate the matrix data structures */ 1742 for (i = 0; i < PetscMin(A->rmap->n, A->cmap->n); i++) a->imax[i] += mdiag[i]; 1743 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(A, 0, a->imax)); 1744 1745 /* copy old values into new matrix data structure */ 1746 for (i = 0; i < A->rmap->n; i++) { 1747 PetscCall(MatSetValues(A, 1, &i, a->imax[i] - mdiag[i], &oldj[oldi[i]], &olda[oldi[i]], ADD_VALUES)); 1748 if (i < A->cmap->n) PetscCall(MatSetValue(A, i, i, v, ADD_VALUES)); 1749 } 1750 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 1751 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 1752 if (singlemalloc) { 1753 PetscCall(PetscFree3(olda, oldj, oldi)); 1754 } else { 1755 if (free_a) PetscCall(PetscFree(olda)); 1756 if (free_ij) PetscCall(PetscFree(oldj)); 1757 if (free_ij) PetscCall(PetscFree(oldi)); 1758 } 1759 } 1760 PetscCall(PetscFree(mdiag)); 1761 a->diagonaldense = PETSC_TRUE; 1762 PetscFunctionReturn(0); 1763 } 1764 1765 /* 1766 Checks for missing diagonals 1767 */ 1768 PetscErrorCode MatMissingDiagonal_SeqAIJ(Mat A, PetscBool *missing, PetscInt *d) 1769 { 1770 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 1771 PetscInt *diag, *ii = a->i, i; 1772 1773 PetscFunctionBegin; 1774 *missing = PETSC_FALSE; 1775 if (A->rmap->n > 0 && !ii) { 1776 *missing = PETSC_TRUE; 1777 if (d) *d = 0; 1778 PetscCall(PetscInfo(A, "Matrix has no entries therefore is missing diagonal\n")); 1779 } else { 1780 PetscInt n; 1781 n = PetscMin(A->rmap->n, A->cmap->n); 1782 diag = a->diag; 1783 for (i = 0; i < n; i++) { 1784 if (diag[i] >= ii[i + 1]) { 1785 *missing = PETSC_TRUE; 1786 if (d) *d = i; 1787 PetscCall(PetscInfo(A, "Matrix is missing diagonal number %" PetscInt_FMT "\n", i)); 1788 break; 1789 } 1790 } 1791 } 1792 PetscFunctionReturn(0); 1793 } 1794 1795 #include <petscblaslapack.h> 1796 #include <petsc/private/kernels/blockinvert.h> 1797 1798 /* 1799 Note that values is allocated externally by the PC and then passed into this routine 1800 */ 1801 PetscErrorCode MatInvertVariableBlockDiagonal_SeqAIJ(Mat A, PetscInt nblocks, const PetscInt *bsizes, PetscScalar *diag) 1802 { 1803 PetscInt n = A->rmap->n, i, ncnt = 0, *indx, j, bsizemax = 0, *v_pivots; 1804 PetscBool allowzeropivot, zeropivotdetected = PETSC_FALSE; 1805 const PetscReal shift = 0.0; 1806 PetscInt ipvt[5]; 1807 PetscScalar work[25], *v_work; 1808 1809 PetscFunctionBegin; 1810 allowzeropivot = PetscNot(A->erroriffailure); 1811 for (i = 0; i < nblocks; i++) ncnt += bsizes[i]; 1812 PetscCheck(ncnt == n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Total blocksizes %" PetscInt_FMT " doesn't match number matrix rows %" PetscInt_FMT, ncnt, n); 1813 for (i = 0; i < nblocks; i++) bsizemax = PetscMax(bsizemax, bsizes[i]); 1814 PetscCall(PetscMalloc1(bsizemax, &indx)); 1815 if (bsizemax > 7) PetscCall(PetscMalloc2(bsizemax, &v_work, bsizemax, &v_pivots)); 1816 ncnt = 0; 1817 for (i = 0; i < nblocks; i++) { 1818 for (j = 0; j < bsizes[i]; j++) indx[j] = ncnt + j; 1819 PetscCall(MatGetValues(A, bsizes[i], indx, bsizes[i], indx, diag)); 1820 switch (bsizes[i]) { 1821 case 1: 1822 *diag = 1.0 / (*diag); 1823 break; 1824 case 2: 1825 PetscCall(PetscKernel_A_gets_inverse_A_2(diag, shift, allowzeropivot, &zeropivotdetected)); 1826 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1827 PetscCall(PetscKernel_A_gets_transpose_A_2(diag)); 1828 break; 1829 case 3: 1830 PetscCall(PetscKernel_A_gets_inverse_A_3(diag, shift, allowzeropivot, &zeropivotdetected)); 1831 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1832 PetscCall(PetscKernel_A_gets_transpose_A_3(diag)); 1833 break; 1834 case 4: 1835 PetscCall(PetscKernel_A_gets_inverse_A_4(diag, shift, allowzeropivot, &zeropivotdetected)); 1836 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1837 PetscCall(PetscKernel_A_gets_transpose_A_4(diag)); 1838 break; 1839 case 5: 1840 PetscCall(PetscKernel_A_gets_inverse_A_5(diag, ipvt, work, shift, allowzeropivot, &zeropivotdetected)); 1841 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1842 PetscCall(PetscKernel_A_gets_transpose_A_5(diag)); 1843 break; 1844 case 6: 1845 PetscCall(PetscKernel_A_gets_inverse_A_6(diag, shift, allowzeropivot, &zeropivotdetected)); 1846 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1847 PetscCall(PetscKernel_A_gets_transpose_A_6(diag)); 1848 break; 1849 case 7: 1850 PetscCall(PetscKernel_A_gets_inverse_A_7(diag, shift, allowzeropivot, &zeropivotdetected)); 1851 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1852 PetscCall(PetscKernel_A_gets_transpose_A_7(diag)); 1853 break; 1854 default: 1855 PetscCall(PetscKernel_A_gets_inverse_A(bsizes[i], diag, v_pivots, v_work, allowzeropivot, &zeropivotdetected)); 1856 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1857 PetscCall(PetscKernel_A_gets_transpose_A_N(diag, bsizes[i])); 1858 } 1859 ncnt += bsizes[i]; 1860 diag += bsizes[i] * bsizes[i]; 1861 } 1862 if (bsizemax > 7) PetscCall(PetscFree2(v_work, v_pivots)); 1863 PetscCall(PetscFree(indx)); 1864 PetscFunctionReturn(0); 1865 } 1866 1867 /* 1868 Negative shift indicates do not generate an error if there is a zero diagonal, just invert it anyways 1869 */ 1870 PetscErrorCode MatInvertDiagonal_SeqAIJ(Mat A, PetscScalar omega, PetscScalar fshift) 1871 { 1872 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 1873 PetscInt i, *diag, m = A->rmap->n; 1874 const MatScalar *v; 1875 PetscScalar *idiag, *mdiag; 1876 1877 PetscFunctionBegin; 1878 if (a->idiagvalid) PetscFunctionReturn(0); 1879 PetscCall(MatMarkDiagonal_SeqAIJ(A)); 1880 diag = a->diag; 1881 if (!a->idiag) { PetscCall(PetscMalloc3(m, &a->idiag, m, &a->mdiag, m, &a->ssor_work)); } 1882 1883 mdiag = a->mdiag; 1884 idiag = a->idiag; 1885 PetscCall(MatSeqAIJGetArrayRead(A, &v)); 1886 if (omega == 1.0 && PetscRealPart(fshift) <= 0.0) { 1887 for (i = 0; i < m; i++) { 1888 mdiag[i] = v[diag[i]]; 1889 if (!PetscAbsScalar(mdiag[i])) { /* zero diagonal */ 1890 if (PetscRealPart(fshift)) { 1891 PetscCall(PetscInfo(A, "Zero diagonal on row %" PetscInt_FMT "\n", i)); 1892 A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1893 A->factorerror_zeropivot_value = 0.0; 1894 A->factorerror_zeropivot_row = i; 1895 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Zero diagonal on row %" PetscInt_FMT, i); 1896 } 1897 idiag[i] = 1.0 / v[diag[i]]; 1898 } 1899 PetscCall(PetscLogFlops(m)); 1900 } else { 1901 for (i = 0; i < m; i++) { 1902 mdiag[i] = v[diag[i]]; 1903 idiag[i] = omega / (fshift + v[diag[i]]); 1904 } 1905 PetscCall(PetscLogFlops(2.0 * m)); 1906 } 1907 a->idiagvalid = PETSC_TRUE; 1908 PetscCall(MatSeqAIJRestoreArrayRead(A, &v)); 1909 PetscFunctionReturn(0); 1910 } 1911 1912 #include <../src/mat/impls/aij/seq/ftn-kernels/frelax.h> 1913 PetscErrorCode MatSOR_SeqAIJ(Mat A, Vec bb, PetscReal omega, MatSORType flag, PetscReal fshift, PetscInt its, PetscInt lits, Vec xx) 1914 { 1915 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 1916 PetscScalar *x, d, sum, *t, scale; 1917 const MatScalar *v, *idiag = NULL, *mdiag, *aa; 1918 const PetscScalar *b, *bs, *xb, *ts; 1919 PetscInt n, m = A->rmap->n, i; 1920 const PetscInt *idx, *diag; 1921 1922 PetscFunctionBegin; 1923 if (a->inode.use && a->inode.checked && omega == 1.0 && fshift == 0.0) { 1924 PetscCall(MatSOR_SeqAIJ_Inode(A, bb, omega, flag, fshift, its, lits, xx)); 1925 PetscFunctionReturn(0); 1926 } 1927 its = its * lits; 1928 1929 if (fshift != a->fshift || omega != a->omega) a->idiagvalid = PETSC_FALSE; /* must recompute idiag[] */ 1930 if (!a->idiagvalid) PetscCall(MatInvertDiagonal_SeqAIJ(A, omega, fshift)); 1931 a->fshift = fshift; 1932 a->omega = omega; 1933 1934 diag = a->diag; 1935 t = a->ssor_work; 1936 idiag = a->idiag; 1937 mdiag = a->mdiag; 1938 1939 PetscCall(MatSeqAIJGetArrayRead(A, &aa)); 1940 PetscCall(VecGetArray(xx, &x)); 1941 PetscCall(VecGetArrayRead(bb, &b)); 1942 /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */ 1943 if (flag == SOR_APPLY_UPPER) { 1944 /* apply (U + D/omega) to the vector */ 1945 bs = b; 1946 for (i = 0; i < m; i++) { 1947 d = fshift + mdiag[i]; 1948 n = a->i[i + 1] - diag[i] - 1; 1949 idx = a->j + diag[i] + 1; 1950 v = aa + diag[i] + 1; 1951 sum = b[i] * d / omega; 1952 PetscSparseDensePlusDot(sum, bs, v, idx, n); 1953 x[i] = sum; 1954 } 1955 PetscCall(VecRestoreArray(xx, &x)); 1956 PetscCall(VecRestoreArrayRead(bb, &b)); 1957 PetscCall(MatSeqAIJRestoreArrayRead(A, &aa)); 1958 PetscCall(PetscLogFlops(a->nz)); 1959 PetscFunctionReturn(0); 1960 } 1961 1962 PetscCheck(flag != SOR_APPLY_LOWER, PETSC_COMM_SELF, PETSC_ERR_SUP, "SOR_APPLY_LOWER is not implemented"); 1963 if (flag & SOR_EISENSTAT) { 1964 /* Let A = L + U + D; where L is lower triangular, 1965 U is upper triangular, E = D/omega; This routine applies 1966 1967 (L + E)^{-1} A (U + E)^{-1} 1968 1969 to a vector efficiently using Eisenstat's trick. 1970 */ 1971 scale = (2.0 / omega) - 1.0; 1972 1973 /* x = (E + U)^{-1} b */ 1974 for (i = m - 1; i >= 0; i--) { 1975 n = a->i[i + 1] - diag[i] - 1; 1976 idx = a->j + diag[i] + 1; 1977 v = aa + diag[i] + 1; 1978 sum = b[i]; 1979 PetscSparseDenseMinusDot(sum, x, v, idx, n); 1980 x[i] = sum * idiag[i]; 1981 } 1982 1983 /* t = b - (2*E - D)x */ 1984 v = aa; 1985 for (i = 0; i < m; i++) t[i] = b[i] - scale * (v[*diag++]) * x[i]; 1986 1987 /* t = (E + L)^{-1}t */ 1988 ts = t; 1989 diag = a->diag; 1990 for (i = 0; i < m; i++) { 1991 n = diag[i] - a->i[i]; 1992 idx = a->j + a->i[i]; 1993 v = aa + a->i[i]; 1994 sum = t[i]; 1995 PetscSparseDenseMinusDot(sum, ts, v, idx, n); 1996 t[i] = sum * idiag[i]; 1997 /* x = x + t */ 1998 x[i] += t[i]; 1999 } 2000 2001 PetscCall(PetscLogFlops(6.0 * m - 1 + 2.0 * a->nz)); 2002 PetscCall(VecRestoreArray(xx, &x)); 2003 PetscCall(VecRestoreArrayRead(bb, &b)); 2004 PetscFunctionReturn(0); 2005 } 2006 if (flag & SOR_ZERO_INITIAL_GUESS) { 2007 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 2008 for (i = 0; i < m; i++) { 2009 n = diag[i] - a->i[i]; 2010 idx = a->j + a->i[i]; 2011 v = aa + a->i[i]; 2012 sum = b[i]; 2013 PetscSparseDenseMinusDot(sum, x, v, idx, n); 2014 t[i] = sum; 2015 x[i] = sum * idiag[i]; 2016 } 2017 xb = t; 2018 PetscCall(PetscLogFlops(a->nz)); 2019 } else xb = b; 2020 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 2021 for (i = m - 1; i >= 0; i--) { 2022 n = a->i[i + 1] - diag[i] - 1; 2023 idx = a->j + diag[i] + 1; 2024 v = aa + diag[i] + 1; 2025 sum = xb[i]; 2026 PetscSparseDenseMinusDot(sum, x, v, idx, n); 2027 if (xb == b) { 2028 x[i] = sum * idiag[i]; 2029 } else { 2030 x[i] = (1 - omega) * x[i] + sum * idiag[i]; /* omega in idiag */ 2031 } 2032 } 2033 PetscCall(PetscLogFlops(a->nz)); /* assumes 1/2 in upper */ 2034 } 2035 its--; 2036 } 2037 while (its--) { 2038 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 2039 for (i = 0; i < m; i++) { 2040 /* lower */ 2041 n = diag[i] - a->i[i]; 2042 idx = a->j + a->i[i]; 2043 v = aa + a->i[i]; 2044 sum = b[i]; 2045 PetscSparseDenseMinusDot(sum, x, v, idx, n); 2046 t[i] = sum; /* save application of the lower-triangular part */ 2047 /* upper */ 2048 n = a->i[i + 1] - diag[i] - 1; 2049 idx = a->j + diag[i] + 1; 2050 v = aa + diag[i] + 1; 2051 PetscSparseDenseMinusDot(sum, x, v, idx, n); 2052 x[i] = (1. - omega) * x[i] + sum * idiag[i]; /* omega in idiag */ 2053 } 2054 xb = t; 2055 PetscCall(PetscLogFlops(2.0 * a->nz)); 2056 } else xb = b; 2057 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 2058 for (i = m - 1; i >= 0; i--) { 2059 sum = xb[i]; 2060 if (xb == b) { 2061 /* whole matrix (no checkpointing available) */ 2062 n = a->i[i + 1] - a->i[i]; 2063 idx = a->j + a->i[i]; 2064 v = aa + a->i[i]; 2065 PetscSparseDenseMinusDot(sum, x, v, idx, n); 2066 x[i] = (1. - omega) * x[i] + (sum + mdiag[i] * x[i]) * idiag[i]; 2067 } else { /* lower-triangular part has been saved, so only apply upper-triangular */ 2068 n = a->i[i + 1] - diag[i] - 1; 2069 idx = a->j + diag[i] + 1; 2070 v = aa + diag[i] + 1; 2071 PetscSparseDenseMinusDot(sum, x, v, idx, n); 2072 x[i] = (1. - omega) * x[i] + sum * idiag[i]; /* omega in idiag */ 2073 } 2074 } 2075 if (xb == b) { 2076 PetscCall(PetscLogFlops(2.0 * a->nz)); 2077 } else { 2078 PetscCall(PetscLogFlops(a->nz)); /* assumes 1/2 in upper */ 2079 } 2080 } 2081 } 2082 PetscCall(MatSeqAIJRestoreArrayRead(A, &aa)); 2083 PetscCall(VecRestoreArray(xx, &x)); 2084 PetscCall(VecRestoreArrayRead(bb, &b)); 2085 PetscFunctionReturn(0); 2086 } 2087 2088 PetscErrorCode MatGetInfo_SeqAIJ(Mat A, MatInfoType flag, MatInfo *info) 2089 { 2090 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 2091 2092 PetscFunctionBegin; 2093 info->block_size = 1.0; 2094 info->nz_allocated = a->maxnz; 2095 info->nz_used = a->nz; 2096 info->nz_unneeded = (a->maxnz - a->nz); 2097 info->assemblies = A->num_ass; 2098 info->mallocs = A->info.mallocs; 2099 info->memory = 0; /* REVIEW ME */ 2100 if (A->factortype) { 2101 info->fill_ratio_given = A->info.fill_ratio_given; 2102 info->fill_ratio_needed = A->info.fill_ratio_needed; 2103 info->factor_mallocs = A->info.factor_mallocs; 2104 } else { 2105 info->fill_ratio_given = 0; 2106 info->fill_ratio_needed = 0; 2107 info->factor_mallocs = 0; 2108 } 2109 PetscFunctionReturn(0); 2110 } 2111 2112 PetscErrorCode MatZeroRows_SeqAIJ(Mat A, PetscInt N, const PetscInt rows[], PetscScalar diag, Vec x, Vec b) 2113 { 2114 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 2115 PetscInt i, m = A->rmap->n - 1; 2116 const PetscScalar *xx; 2117 PetscScalar *bb, *aa; 2118 PetscInt d = 0; 2119 2120 PetscFunctionBegin; 2121 if (x && b) { 2122 PetscCall(VecGetArrayRead(x, &xx)); 2123 PetscCall(VecGetArray(b, &bb)); 2124 for (i = 0; i < N; i++) { 2125 PetscCheck(rows[i] >= 0 && rows[i] <= m, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "row %" PetscInt_FMT " out of range", rows[i]); 2126 if (rows[i] >= A->cmap->n) continue; 2127 bb[rows[i]] = diag * xx[rows[i]]; 2128 } 2129 PetscCall(VecRestoreArrayRead(x, &xx)); 2130 PetscCall(VecRestoreArray(b, &bb)); 2131 } 2132 2133 PetscCall(MatSeqAIJGetArray(A, &aa)); 2134 if (a->keepnonzeropattern) { 2135 for (i = 0; i < N; i++) { 2136 PetscCheck(rows[i] >= 0 && rows[i] <= m, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "row %" PetscInt_FMT " out of range", rows[i]); 2137 PetscCall(PetscArrayzero(&aa[a->i[rows[i]]], a->ilen[rows[i]])); 2138 } 2139 if (diag != 0.0) { 2140 for (i = 0; i < N; i++) { 2141 d = rows[i]; 2142 if (rows[i] >= A->cmap->n) continue; 2143 PetscCheck(a->diag[d] < a->i[d + 1], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Matrix is missing diagonal entry in the zeroed row %" PetscInt_FMT, d); 2144 } 2145 for (i = 0; i < N; i++) { 2146 if (rows[i] >= A->cmap->n) continue; 2147 aa[a->diag[rows[i]]] = diag; 2148 } 2149 } 2150 } else { 2151 if (diag != 0.0) { 2152 for (i = 0; i < N; i++) { 2153 PetscCheck(rows[i] >= 0 && rows[i] <= m, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "row %" PetscInt_FMT " out of range", rows[i]); 2154 if (a->ilen[rows[i]] > 0) { 2155 if (rows[i] >= A->cmap->n) { 2156 a->ilen[rows[i]] = 0; 2157 } else { 2158 a->ilen[rows[i]] = 1; 2159 aa[a->i[rows[i]]] = diag; 2160 a->j[a->i[rows[i]]] = rows[i]; 2161 } 2162 } else if (rows[i] < A->cmap->n) { /* in case row was completely empty */ 2163 PetscCall(MatSetValues_SeqAIJ(A, 1, &rows[i], 1, &rows[i], &diag, INSERT_VALUES)); 2164 } 2165 } 2166 } else { 2167 for (i = 0; i < N; i++) { 2168 PetscCheck(rows[i] >= 0 && rows[i] <= m, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "row %" PetscInt_FMT " out of range", rows[i]); 2169 a->ilen[rows[i]] = 0; 2170 } 2171 } 2172 A->nonzerostate++; 2173 } 2174 PetscCall(MatSeqAIJRestoreArray(A, &aa)); 2175 PetscUseTypeMethod(A, assemblyend, MAT_FINAL_ASSEMBLY); 2176 PetscFunctionReturn(0); 2177 } 2178 2179 PetscErrorCode MatZeroRowsColumns_SeqAIJ(Mat A, PetscInt N, const PetscInt rows[], PetscScalar diag, Vec x, Vec b) 2180 { 2181 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 2182 PetscInt i, j, m = A->rmap->n - 1, d = 0; 2183 PetscBool missing, *zeroed, vecs = PETSC_FALSE; 2184 const PetscScalar *xx; 2185 PetscScalar *bb, *aa; 2186 2187 PetscFunctionBegin; 2188 if (!N) PetscFunctionReturn(0); 2189 PetscCall(MatSeqAIJGetArray(A, &aa)); 2190 if (x && b) { 2191 PetscCall(VecGetArrayRead(x, &xx)); 2192 PetscCall(VecGetArray(b, &bb)); 2193 vecs = PETSC_TRUE; 2194 } 2195 PetscCall(PetscCalloc1(A->rmap->n, &zeroed)); 2196 for (i = 0; i < N; i++) { 2197 PetscCheck(rows[i] >= 0 && rows[i] <= m, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "row %" PetscInt_FMT " out of range", rows[i]); 2198 PetscCall(PetscArrayzero(&aa[a->i[rows[i]]], a->ilen[rows[i]])); 2199 2200 zeroed[rows[i]] = PETSC_TRUE; 2201 } 2202 for (i = 0; i < A->rmap->n; i++) { 2203 if (!zeroed[i]) { 2204 for (j = a->i[i]; j < a->i[i + 1]; j++) { 2205 if (a->j[j] < A->rmap->n && zeroed[a->j[j]]) { 2206 if (vecs) bb[i] -= aa[j] * xx[a->j[j]]; 2207 aa[j] = 0.0; 2208 } 2209 } 2210 } else if (vecs && i < A->cmap->N) bb[i] = diag * xx[i]; 2211 } 2212 if (x && b) { 2213 PetscCall(VecRestoreArrayRead(x, &xx)); 2214 PetscCall(VecRestoreArray(b, &bb)); 2215 } 2216 PetscCall(PetscFree(zeroed)); 2217 if (diag != 0.0) { 2218 PetscCall(MatMissingDiagonal_SeqAIJ(A, &missing, &d)); 2219 if (missing) { 2220 for (i = 0; i < N; i++) { 2221 if (rows[i] >= A->cmap->N) continue; 2222 PetscCheck(!a->nonew || rows[i] < d, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Matrix is missing diagonal entry in row %" PetscInt_FMT " (%" PetscInt_FMT ")", d, rows[i]); 2223 PetscCall(MatSetValues_SeqAIJ(A, 1, &rows[i], 1, &rows[i], &diag, INSERT_VALUES)); 2224 } 2225 } else { 2226 for (i = 0; i < N; i++) aa[a->diag[rows[i]]] = diag; 2227 } 2228 } 2229 PetscCall(MatSeqAIJRestoreArray(A, &aa)); 2230 PetscUseTypeMethod(A, assemblyend, MAT_FINAL_ASSEMBLY); 2231 PetscFunctionReturn(0); 2232 } 2233 2234 PetscErrorCode MatGetRow_SeqAIJ(Mat A, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v) 2235 { 2236 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 2237 const PetscScalar *aa; 2238 PetscInt *itmp; 2239 2240 PetscFunctionBegin; 2241 PetscCall(MatSeqAIJGetArrayRead(A, &aa)); 2242 *nz = a->i[row + 1] - a->i[row]; 2243 if (v) *v = (PetscScalar *)(aa + a->i[row]); 2244 if (idx) { 2245 itmp = a->j + a->i[row]; 2246 if (*nz) *idx = itmp; 2247 else *idx = NULL; 2248 } 2249 PetscCall(MatSeqAIJRestoreArrayRead(A, &aa)); 2250 PetscFunctionReturn(0); 2251 } 2252 2253 PetscErrorCode MatRestoreRow_SeqAIJ(Mat A, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v) 2254 { 2255 PetscFunctionBegin; 2256 if (nz) *nz = 0; 2257 if (idx) *idx = NULL; 2258 if (v) *v = NULL; 2259 PetscFunctionReturn(0); 2260 } 2261 2262 PetscErrorCode MatNorm_SeqAIJ(Mat A, NormType type, PetscReal *nrm) 2263 { 2264 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 2265 const MatScalar *v; 2266 PetscReal sum = 0.0; 2267 PetscInt i, j; 2268 2269 PetscFunctionBegin; 2270 PetscCall(MatSeqAIJGetArrayRead(A, &v)); 2271 if (type == NORM_FROBENIUS) { 2272 #if defined(PETSC_USE_REAL___FP16) 2273 PetscBLASInt one = 1, nz = a->nz; 2274 PetscCallBLAS("BLASnrm2", *nrm = BLASnrm2_(&nz, v, &one)); 2275 #else 2276 for (i = 0; i < a->nz; i++) { 2277 sum += PetscRealPart(PetscConj(*v) * (*v)); 2278 v++; 2279 } 2280 *nrm = PetscSqrtReal(sum); 2281 #endif 2282 PetscCall(PetscLogFlops(2.0 * a->nz)); 2283 } else if (type == NORM_1) { 2284 PetscReal *tmp; 2285 PetscInt *jj = a->j; 2286 PetscCall(PetscCalloc1(A->cmap->n + 1, &tmp)); 2287 *nrm = 0.0; 2288 for (j = 0; j < a->nz; j++) { 2289 tmp[*jj++] += PetscAbsScalar(*v); 2290 v++; 2291 } 2292 for (j = 0; j < A->cmap->n; j++) { 2293 if (tmp[j] > *nrm) *nrm = tmp[j]; 2294 } 2295 PetscCall(PetscFree(tmp)); 2296 PetscCall(PetscLogFlops(PetscMax(a->nz - 1, 0))); 2297 } else if (type == NORM_INFINITY) { 2298 *nrm = 0.0; 2299 for (j = 0; j < A->rmap->n; j++) { 2300 const PetscScalar *v2 = v + a->i[j]; 2301 sum = 0.0; 2302 for (i = 0; i < a->i[j + 1] - a->i[j]; i++) { 2303 sum += PetscAbsScalar(*v2); 2304 v2++; 2305 } 2306 if (sum > *nrm) *nrm = sum; 2307 } 2308 PetscCall(PetscLogFlops(PetscMax(a->nz - 1, 0))); 2309 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for two norm"); 2310 PetscCall(MatSeqAIJRestoreArrayRead(A, &v)); 2311 PetscFunctionReturn(0); 2312 } 2313 2314 PetscErrorCode MatIsTranspose_SeqAIJ(Mat A, Mat B, PetscReal tol, PetscBool *f) 2315 { 2316 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)A->data, *bij = (Mat_SeqAIJ *)B->data; 2317 PetscInt *adx, *bdx, *aii, *bii, *aptr, *bptr; 2318 const MatScalar *va, *vb; 2319 PetscInt ma, na, mb, nb, i; 2320 2321 PetscFunctionBegin; 2322 PetscCall(MatGetSize(A, &ma, &na)); 2323 PetscCall(MatGetSize(B, &mb, &nb)); 2324 if (ma != nb || na != mb) { 2325 *f = PETSC_FALSE; 2326 PetscFunctionReturn(0); 2327 } 2328 PetscCall(MatSeqAIJGetArrayRead(A, &va)); 2329 PetscCall(MatSeqAIJGetArrayRead(B, &vb)); 2330 aii = aij->i; 2331 bii = bij->i; 2332 adx = aij->j; 2333 bdx = bij->j; 2334 PetscCall(PetscMalloc1(ma, &aptr)); 2335 PetscCall(PetscMalloc1(mb, &bptr)); 2336 for (i = 0; i < ma; i++) aptr[i] = aii[i]; 2337 for (i = 0; i < mb; i++) bptr[i] = bii[i]; 2338 2339 *f = PETSC_TRUE; 2340 for (i = 0; i < ma; i++) { 2341 while (aptr[i] < aii[i + 1]) { 2342 PetscInt idc, idr; 2343 PetscScalar vc, vr; 2344 /* column/row index/value */ 2345 idc = adx[aptr[i]]; 2346 idr = bdx[bptr[idc]]; 2347 vc = va[aptr[i]]; 2348 vr = vb[bptr[idc]]; 2349 if (i != idr || PetscAbsScalar(vc - vr) > tol) { 2350 *f = PETSC_FALSE; 2351 goto done; 2352 } else { 2353 aptr[i]++; 2354 if (B || i != idc) bptr[idc]++; 2355 } 2356 } 2357 } 2358 done: 2359 PetscCall(PetscFree(aptr)); 2360 PetscCall(PetscFree(bptr)); 2361 PetscCall(MatSeqAIJRestoreArrayRead(A, &va)); 2362 PetscCall(MatSeqAIJRestoreArrayRead(B, &vb)); 2363 PetscFunctionReturn(0); 2364 } 2365 2366 PetscErrorCode MatIsHermitianTranspose_SeqAIJ(Mat A, Mat B, PetscReal tol, PetscBool *f) 2367 { 2368 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)A->data, *bij = (Mat_SeqAIJ *)B->data; 2369 PetscInt *adx, *bdx, *aii, *bii, *aptr, *bptr; 2370 MatScalar *va, *vb; 2371 PetscInt ma, na, mb, nb, i; 2372 2373 PetscFunctionBegin; 2374 PetscCall(MatGetSize(A, &ma, &na)); 2375 PetscCall(MatGetSize(B, &mb, &nb)); 2376 if (ma != nb || na != mb) { 2377 *f = PETSC_FALSE; 2378 PetscFunctionReturn(0); 2379 } 2380 aii = aij->i; 2381 bii = bij->i; 2382 adx = aij->j; 2383 bdx = bij->j; 2384 va = aij->a; 2385 vb = bij->a; 2386 PetscCall(PetscMalloc1(ma, &aptr)); 2387 PetscCall(PetscMalloc1(mb, &bptr)); 2388 for (i = 0; i < ma; i++) aptr[i] = aii[i]; 2389 for (i = 0; i < mb; i++) bptr[i] = bii[i]; 2390 2391 *f = PETSC_TRUE; 2392 for (i = 0; i < ma; i++) { 2393 while (aptr[i] < aii[i + 1]) { 2394 PetscInt idc, idr; 2395 PetscScalar vc, vr; 2396 /* column/row index/value */ 2397 idc = adx[aptr[i]]; 2398 idr = bdx[bptr[idc]]; 2399 vc = va[aptr[i]]; 2400 vr = vb[bptr[idc]]; 2401 if (i != idr || PetscAbsScalar(vc - PetscConj(vr)) > tol) { 2402 *f = PETSC_FALSE; 2403 goto done; 2404 } else { 2405 aptr[i]++; 2406 if (B || i != idc) bptr[idc]++; 2407 } 2408 } 2409 } 2410 done: 2411 PetscCall(PetscFree(aptr)); 2412 PetscCall(PetscFree(bptr)); 2413 PetscFunctionReturn(0); 2414 } 2415 2416 PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A, PetscReal tol, PetscBool *f) 2417 { 2418 PetscFunctionBegin; 2419 PetscCall(MatIsTranspose_SeqAIJ(A, A, tol, f)); 2420 PetscFunctionReturn(0); 2421 } 2422 2423 PetscErrorCode MatIsHermitian_SeqAIJ(Mat A, PetscReal tol, PetscBool *f) 2424 { 2425 PetscFunctionBegin; 2426 PetscCall(MatIsHermitianTranspose_SeqAIJ(A, A, tol, f)); 2427 PetscFunctionReturn(0); 2428 } 2429 2430 PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A, Vec ll, Vec rr) 2431 { 2432 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 2433 const PetscScalar *l, *r; 2434 PetscScalar x; 2435 MatScalar *v; 2436 PetscInt i, j, m = A->rmap->n, n = A->cmap->n, M, nz = a->nz; 2437 const PetscInt *jj; 2438 2439 PetscFunctionBegin; 2440 if (ll) { 2441 /* The local size is used so that VecMPI can be passed to this routine 2442 by MatDiagonalScale_MPIAIJ */ 2443 PetscCall(VecGetLocalSize(ll, &m)); 2444 PetscCheck(m == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Left scaling vector wrong length"); 2445 PetscCall(VecGetArrayRead(ll, &l)); 2446 PetscCall(MatSeqAIJGetArray(A, &v)); 2447 for (i = 0; i < m; i++) { 2448 x = l[i]; 2449 M = a->i[i + 1] - a->i[i]; 2450 for (j = 0; j < M; j++) (*v++) *= x; 2451 } 2452 PetscCall(VecRestoreArrayRead(ll, &l)); 2453 PetscCall(PetscLogFlops(nz)); 2454 PetscCall(MatSeqAIJRestoreArray(A, &v)); 2455 } 2456 if (rr) { 2457 PetscCall(VecGetLocalSize(rr, &n)); 2458 PetscCheck(n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Right scaling vector wrong length"); 2459 PetscCall(VecGetArrayRead(rr, &r)); 2460 PetscCall(MatSeqAIJGetArray(A, &v)); 2461 jj = a->j; 2462 for (i = 0; i < nz; i++) (*v++) *= r[*jj++]; 2463 PetscCall(MatSeqAIJRestoreArray(A, &v)); 2464 PetscCall(VecRestoreArrayRead(rr, &r)); 2465 PetscCall(PetscLogFlops(nz)); 2466 } 2467 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 2468 PetscFunctionReturn(0); 2469 } 2470 2471 PetscErrorCode MatCreateSubMatrix_SeqAIJ(Mat A, IS isrow, IS iscol, PetscInt csize, MatReuse scall, Mat *B) 2472 { 2473 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *c; 2474 PetscInt *smap, i, k, kstart, kend, oldcols = A->cmap->n, *lens; 2475 PetscInt row, mat_i, *mat_j, tcol, first, step, *mat_ilen, sum, lensi; 2476 const PetscInt *irow, *icol; 2477 const PetscScalar *aa; 2478 PetscInt nrows, ncols; 2479 PetscInt *starts, *j_new, *i_new, *aj = a->j, *ai = a->i, ii, *ailen = a->ilen; 2480 MatScalar *a_new, *mat_a; 2481 Mat C; 2482 PetscBool stride; 2483 2484 PetscFunctionBegin; 2485 PetscCall(ISGetIndices(isrow, &irow)); 2486 PetscCall(ISGetLocalSize(isrow, &nrows)); 2487 PetscCall(ISGetLocalSize(iscol, &ncols)); 2488 2489 PetscCall(PetscObjectTypeCompare((PetscObject)iscol, ISSTRIDE, &stride)); 2490 if (stride) { 2491 PetscCall(ISStrideGetInfo(iscol, &first, &step)); 2492 } else { 2493 first = 0; 2494 step = 0; 2495 } 2496 if (stride && step == 1) { 2497 /* special case of contiguous rows */ 2498 PetscCall(PetscMalloc2(nrows, &lens, nrows, &starts)); 2499 /* loop over new rows determining lens and starting points */ 2500 for (i = 0; i < nrows; i++) { 2501 kstart = ai[irow[i]]; 2502 kend = kstart + ailen[irow[i]]; 2503 starts[i] = kstart; 2504 for (k = kstart; k < kend; k++) { 2505 if (aj[k] >= first) { 2506 starts[i] = k; 2507 break; 2508 } 2509 } 2510 sum = 0; 2511 while (k < kend) { 2512 if (aj[k++] >= first + ncols) break; 2513 sum++; 2514 } 2515 lens[i] = sum; 2516 } 2517 /* create submatrix */ 2518 if (scall == MAT_REUSE_MATRIX) { 2519 PetscInt n_cols, n_rows; 2520 PetscCall(MatGetSize(*B, &n_rows, &n_cols)); 2521 PetscCheck(n_rows == nrows && n_cols == ncols, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Reused submatrix wrong size"); 2522 PetscCall(MatZeroEntries(*B)); 2523 C = *B; 2524 } else { 2525 PetscInt rbs, cbs; 2526 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &C)); 2527 PetscCall(MatSetSizes(C, nrows, ncols, PETSC_DETERMINE, PETSC_DETERMINE)); 2528 PetscCall(ISGetBlockSize(isrow, &rbs)); 2529 PetscCall(ISGetBlockSize(iscol, &cbs)); 2530 PetscCall(MatSetBlockSizes(C, rbs, cbs)); 2531 PetscCall(MatSetType(C, ((PetscObject)A)->type_name)); 2532 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(C, 0, lens)); 2533 } 2534 c = (Mat_SeqAIJ *)C->data; 2535 2536 /* loop over rows inserting into submatrix */ 2537 a_new = c->a; 2538 j_new = c->j; 2539 i_new = c->i; 2540 PetscCall(MatSeqAIJGetArrayRead(A, &aa)); 2541 for (i = 0; i < nrows; i++) { 2542 ii = starts[i]; 2543 lensi = lens[i]; 2544 for (k = 0; k < lensi; k++) *j_new++ = aj[ii + k] - first; 2545 PetscCall(PetscArraycpy(a_new, aa + starts[i], lensi)); 2546 a_new += lensi; 2547 i_new[i + 1] = i_new[i] + lensi; 2548 c->ilen[i] = lensi; 2549 } 2550 PetscCall(MatSeqAIJRestoreArrayRead(A, &aa)); 2551 PetscCall(PetscFree2(lens, starts)); 2552 } else { 2553 PetscCall(ISGetIndices(iscol, &icol)); 2554 PetscCall(PetscCalloc1(oldcols, &smap)); 2555 PetscCall(PetscMalloc1(1 + nrows, &lens)); 2556 for (i = 0; i < ncols; i++) { 2557 PetscCheck(icol[i] < oldcols, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Requesting column beyond largest column icol[%" PetscInt_FMT "] %" PetscInt_FMT " >= A->cmap->n %" PetscInt_FMT, i, icol[i], oldcols); 2558 smap[icol[i]] = i + 1; 2559 } 2560 2561 /* determine lens of each row */ 2562 for (i = 0; i < nrows; i++) { 2563 kstart = ai[irow[i]]; 2564 kend = kstart + a->ilen[irow[i]]; 2565 lens[i] = 0; 2566 for (k = kstart; k < kend; k++) { 2567 if (smap[aj[k]]) lens[i]++; 2568 } 2569 } 2570 /* Create and fill new matrix */ 2571 if (scall == MAT_REUSE_MATRIX) { 2572 PetscBool equal; 2573 2574 c = (Mat_SeqAIJ *)((*B)->data); 2575 PetscCheck((*B)->rmap->n == nrows && (*B)->cmap->n == ncols, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Cannot reuse matrix. wrong size"); 2576 PetscCall(PetscArraycmp(c->ilen, lens, (*B)->rmap->n, &equal)); 2577 PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Cannot reuse matrix. wrong no of nonzeros"); 2578 PetscCall(PetscArrayzero(c->ilen, (*B)->rmap->n)); 2579 C = *B; 2580 } else { 2581 PetscInt rbs, cbs; 2582 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &C)); 2583 PetscCall(MatSetSizes(C, nrows, ncols, PETSC_DETERMINE, PETSC_DETERMINE)); 2584 PetscCall(ISGetBlockSize(isrow, &rbs)); 2585 PetscCall(ISGetBlockSize(iscol, &cbs)); 2586 PetscCall(MatSetBlockSizes(C, rbs, cbs)); 2587 PetscCall(MatSetType(C, ((PetscObject)A)->type_name)); 2588 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(C, 0, lens)); 2589 } 2590 PetscCall(MatSeqAIJGetArrayRead(A, &aa)); 2591 c = (Mat_SeqAIJ *)(C->data); 2592 for (i = 0; i < nrows; i++) { 2593 row = irow[i]; 2594 kstart = ai[row]; 2595 kend = kstart + a->ilen[row]; 2596 mat_i = c->i[i]; 2597 mat_j = c->j + mat_i; 2598 mat_a = c->a + mat_i; 2599 mat_ilen = c->ilen + i; 2600 for (k = kstart; k < kend; k++) { 2601 if ((tcol = smap[a->j[k]])) { 2602 *mat_j++ = tcol - 1; 2603 *mat_a++ = aa[k]; 2604 (*mat_ilen)++; 2605 } 2606 } 2607 } 2608 PetscCall(MatSeqAIJRestoreArrayRead(A, &aa)); 2609 /* Free work space */ 2610 PetscCall(ISRestoreIndices(iscol, &icol)); 2611 PetscCall(PetscFree(smap)); 2612 PetscCall(PetscFree(lens)); 2613 /* sort */ 2614 for (i = 0; i < nrows; i++) { 2615 PetscInt ilen; 2616 2617 mat_i = c->i[i]; 2618 mat_j = c->j + mat_i; 2619 mat_a = c->a + mat_i; 2620 ilen = c->ilen[i]; 2621 PetscCall(PetscSortIntWithScalarArray(ilen, mat_j, mat_a)); 2622 } 2623 } 2624 #if defined(PETSC_HAVE_DEVICE) 2625 PetscCall(MatBindToCPU(C, A->boundtocpu)); 2626 #endif 2627 PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY)); 2628 PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY)); 2629 2630 PetscCall(ISRestoreIndices(isrow, &irow)); 2631 *B = C; 2632 PetscFunctionReturn(0); 2633 } 2634 2635 PetscErrorCode MatGetMultiProcBlock_SeqAIJ(Mat mat, MPI_Comm subComm, MatReuse scall, Mat *subMat) 2636 { 2637 Mat B; 2638 2639 PetscFunctionBegin; 2640 if (scall == MAT_INITIAL_MATRIX) { 2641 PetscCall(MatCreate(subComm, &B)); 2642 PetscCall(MatSetSizes(B, mat->rmap->n, mat->cmap->n, mat->rmap->n, mat->cmap->n)); 2643 PetscCall(MatSetBlockSizesFromMats(B, mat, mat)); 2644 PetscCall(MatSetType(B, MATSEQAIJ)); 2645 PetscCall(MatDuplicateNoCreate_SeqAIJ(B, mat, MAT_COPY_VALUES, PETSC_TRUE)); 2646 *subMat = B; 2647 } else { 2648 PetscCall(MatCopy_SeqAIJ(mat, *subMat, SAME_NONZERO_PATTERN)); 2649 } 2650 PetscFunctionReturn(0); 2651 } 2652 2653 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA, IS row, IS col, const MatFactorInfo *info) 2654 { 2655 Mat_SeqAIJ *a = (Mat_SeqAIJ *)inA->data; 2656 Mat outA; 2657 PetscBool row_identity, col_identity; 2658 2659 PetscFunctionBegin; 2660 PetscCheck(info->levels == 0, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only levels=0 supported for in-place ilu"); 2661 2662 PetscCall(ISIdentity(row, &row_identity)); 2663 PetscCall(ISIdentity(col, &col_identity)); 2664 2665 outA = inA; 2666 outA->factortype = MAT_FACTOR_LU; 2667 PetscCall(PetscFree(inA->solvertype)); 2668 PetscCall(PetscStrallocpy(MATSOLVERPETSC, &inA->solvertype)); 2669 2670 PetscCall(PetscObjectReference((PetscObject)row)); 2671 PetscCall(ISDestroy(&a->row)); 2672 2673 a->row = row; 2674 2675 PetscCall(PetscObjectReference((PetscObject)col)); 2676 PetscCall(ISDestroy(&a->col)); 2677 2678 a->col = col; 2679 2680 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 2681 PetscCall(ISDestroy(&a->icol)); 2682 PetscCall(ISInvertPermutation(col, PETSC_DECIDE, &a->icol)); 2683 2684 if (!a->solve_work) { /* this matrix may have been factored before */ 2685 PetscCall(PetscMalloc1(inA->rmap->n + 1, &a->solve_work)); 2686 } 2687 2688 PetscCall(MatMarkDiagonal_SeqAIJ(inA)); 2689 if (row_identity && col_identity) { 2690 PetscCall(MatLUFactorNumeric_SeqAIJ_inplace(outA, inA, info)); 2691 } else { 2692 PetscCall(MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(outA, inA, info)); 2693 } 2694 PetscFunctionReturn(0); 2695 } 2696 2697 PetscErrorCode MatScale_SeqAIJ(Mat inA, PetscScalar alpha) 2698 { 2699 Mat_SeqAIJ *a = (Mat_SeqAIJ *)inA->data; 2700 PetscScalar *v; 2701 PetscBLASInt one = 1, bnz; 2702 2703 PetscFunctionBegin; 2704 PetscCall(MatSeqAIJGetArray(inA, &v)); 2705 PetscCall(PetscBLASIntCast(a->nz, &bnz)); 2706 PetscCallBLAS("BLASscal", BLASscal_(&bnz, &alpha, v, &one)); 2707 PetscCall(PetscLogFlops(a->nz)); 2708 PetscCall(MatSeqAIJRestoreArray(inA, &v)); 2709 PetscCall(MatSeqAIJInvalidateDiagonal(inA)); 2710 PetscFunctionReturn(0); 2711 } 2712 2713 PetscErrorCode MatDestroySubMatrix_Private(Mat_SubSppt *submatj) 2714 { 2715 PetscInt i; 2716 2717 PetscFunctionBegin; 2718 if (!submatj->id) { /* delete data that are linked only to submats[id=0] */ 2719 PetscCall(PetscFree4(submatj->sbuf1, submatj->ptr, submatj->tmp, submatj->ctr)); 2720 2721 for (i = 0; i < submatj->nrqr; ++i) PetscCall(PetscFree(submatj->sbuf2[i])); 2722 PetscCall(PetscFree3(submatj->sbuf2, submatj->req_size, submatj->req_source1)); 2723 2724 if (submatj->rbuf1) { 2725 PetscCall(PetscFree(submatj->rbuf1[0])); 2726 PetscCall(PetscFree(submatj->rbuf1)); 2727 } 2728 2729 for (i = 0; i < submatj->nrqs; ++i) PetscCall(PetscFree(submatj->rbuf3[i])); 2730 PetscCall(PetscFree3(submatj->req_source2, submatj->rbuf2, submatj->rbuf3)); 2731 PetscCall(PetscFree(submatj->pa)); 2732 } 2733 2734 #if defined(PETSC_USE_CTABLE) 2735 PetscCall(PetscHMapIDestroy(&submatj->rmap)); 2736 if (submatj->cmap_loc) PetscCall(PetscFree(submatj->cmap_loc)); 2737 PetscCall(PetscFree(submatj->rmap_loc)); 2738 #else 2739 PetscCall(PetscFree(submatj->rmap)); 2740 #endif 2741 2742 if (!submatj->allcolumns) { 2743 #if defined(PETSC_USE_CTABLE) 2744 PetscCall(PetscHMapIDestroy((PetscHMapI *)&submatj->cmap)); 2745 #else 2746 PetscCall(PetscFree(submatj->cmap)); 2747 #endif 2748 } 2749 PetscCall(PetscFree(submatj->row2proc)); 2750 2751 PetscCall(PetscFree(submatj)); 2752 PetscFunctionReturn(0); 2753 } 2754 2755 PetscErrorCode MatDestroySubMatrix_SeqAIJ(Mat C) 2756 { 2757 Mat_SeqAIJ *c = (Mat_SeqAIJ *)C->data; 2758 Mat_SubSppt *submatj = c->submatis1; 2759 2760 PetscFunctionBegin; 2761 PetscCall((*submatj->destroy)(C)); 2762 PetscCall(MatDestroySubMatrix_Private(submatj)); 2763 PetscFunctionReturn(0); 2764 } 2765 2766 /* Note this has code duplication with MatDestroySubMatrices_SeqBAIJ() */ 2767 PetscErrorCode MatDestroySubMatrices_SeqAIJ(PetscInt n, Mat *mat[]) 2768 { 2769 PetscInt i; 2770 Mat C; 2771 Mat_SeqAIJ *c; 2772 Mat_SubSppt *submatj; 2773 2774 PetscFunctionBegin; 2775 for (i = 0; i < n; i++) { 2776 C = (*mat)[i]; 2777 c = (Mat_SeqAIJ *)C->data; 2778 submatj = c->submatis1; 2779 if (submatj) { 2780 if (--((PetscObject)C)->refct <= 0) { 2781 PetscCall(PetscFree(C->factorprefix)); 2782 PetscCall((*submatj->destroy)(C)); 2783 PetscCall(MatDestroySubMatrix_Private(submatj)); 2784 PetscCall(PetscFree(C->defaultvectype)); 2785 PetscCall(PetscFree(C->defaultrandtype)); 2786 PetscCall(PetscLayoutDestroy(&C->rmap)); 2787 PetscCall(PetscLayoutDestroy(&C->cmap)); 2788 PetscCall(PetscHeaderDestroy(&C)); 2789 } 2790 } else { 2791 PetscCall(MatDestroy(&C)); 2792 } 2793 } 2794 2795 /* Destroy Dummy submatrices created for reuse */ 2796 PetscCall(MatDestroySubMatrices_Dummy(n, mat)); 2797 2798 PetscCall(PetscFree(*mat)); 2799 PetscFunctionReturn(0); 2800 } 2801 2802 PetscErrorCode MatCreateSubMatrices_SeqAIJ(Mat A, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *B[]) 2803 { 2804 PetscInt i; 2805 2806 PetscFunctionBegin; 2807 if (scall == MAT_INITIAL_MATRIX) PetscCall(PetscCalloc1(n + 1, B)); 2808 2809 for (i = 0; i < n; i++) PetscCall(MatCreateSubMatrix_SeqAIJ(A, irow[i], icol[i], PETSC_DECIDE, scall, &(*B)[i])); 2810 PetscFunctionReturn(0); 2811 } 2812 2813 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A, PetscInt is_max, IS is[], PetscInt ov) 2814 { 2815 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 2816 PetscInt row, i, j, k, l, m, n, *nidx, isz, val; 2817 const PetscInt *idx; 2818 PetscInt start, end, *ai, *aj; 2819 PetscBT table; 2820 2821 PetscFunctionBegin; 2822 m = A->rmap->n; 2823 ai = a->i; 2824 aj = a->j; 2825 2826 PetscCheck(ov >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "illegal negative overlap value used"); 2827 2828 PetscCall(PetscMalloc1(m + 1, &nidx)); 2829 PetscCall(PetscBTCreate(m, &table)); 2830 2831 for (i = 0; i < is_max; i++) { 2832 /* Initialize the two local arrays */ 2833 isz = 0; 2834 PetscCall(PetscBTMemzero(m, table)); 2835 2836 /* Extract the indices, assume there can be duplicate entries */ 2837 PetscCall(ISGetIndices(is[i], &idx)); 2838 PetscCall(ISGetLocalSize(is[i], &n)); 2839 2840 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 2841 for (j = 0; j < n; ++j) { 2842 if (!PetscBTLookupSet(table, idx[j])) nidx[isz++] = idx[j]; 2843 } 2844 PetscCall(ISRestoreIndices(is[i], &idx)); 2845 PetscCall(ISDestroy(&is[i])); 2846 2847 k = 0; 2848 for (j = 0; j < ov; j++) { /* for each overlap */ 2849 n = isz; 2850 for (; k < n; k++) { /* do only those rows in nidx[k], which are not done yet */ 2851 row = nidx[k]; 2852 start = ai[row]; 2853 end = ai[row + 1]; 2854 for (l = start; l < end; l++) { 2855 val = aj[l]; 2856 if (!PetscBTLookupSet(table, val)) nidx[isz++] = val; 2857 } 2858 } 2859 } 2860 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, isz, nidx, PETSC_COPY_VALUES, (is + i))); 2861 } 2862 PetscCall(PetscBTDestroy(&table)); 2863 PetscCall(PetscFree(nidx)); 2864 PetscFunctionReturn(0); 2865 } 2866 2867 /* -------------------------------------------------------------- */ 2868 PetscErrorCode MatPermute_SeqAIJ(Mat A, IS rowp, IS colp, Mat *B) 2869 { 2870 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 2871 PetscInt i, nz = 0, m = A->rmap->n, n = A->cmap->n; 2872 const PetscInt *row, *col; 2873 PetscInt *cnew, j, *lens; 2874 IS icolp, irowp; 2875 PetscInt *cwork = NULL; 2876 PetscScalar *vwork = NULL; 2877 2878 PetscFunctionBegin; 2879 PetscCall(ISInvertPermutation(rowp, PETSC_DECIDE, &irowp)); 2880 PetscCall(ISGetIndices(irowp, &row)); 2881 PetscCall(ISInvertPermutation(colp, PETSC_DECIDE, &icolp)); 2882 PetscCall(ISGetIndices(icolp, &col)); 2883 2884 /* determine lengths of permuted rows */ 2885 PetscCall(PetscMalloc1(m + 1, &lens)); 2886 for (i = 0; i < m; i++) lens[row[i]] = a->i[i + 1] - a->i[i]; 2887 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), B)); 2888 PetscCall(MatSetSizes(*B, m, n, m, n)); 2889 PetscCall(MatSetBlockSizesFromMats(*B, A, A)); 2890 PetscCall(MatSetType(*B, ((PetscObject)A)->type_name)); 2891 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(*B, 0, lens)); 2892 PetscCall(PetscFree(lens)); 2893 2894 PetscCall(PetscMalloc1(n, &cnew)); 2895 for (i = 0; i < m; i++) { 2896 PetscCall(MatGetRow_SeqAIJ(A, i, &nz, &cwork, &vwork)); 2897 for (j = 0; j < nz; j++) cnew[j] = col[cwork[j]]; 2898 PetscCall(MatSetValues_SeqAIJ(*B, 1, &row[i], nz, cnew, vwork, INSERT_VALUES)); 2899 PetscCall(MatRestoreRow_SeqAIJ(A, i, &nz, &cwork, &vwork)); 2900 } 2901 PetscCall(PetscFree(cnew)); 2902 2903 (*B)->assembled = PETSC_FALSE; 2904 2905 #if defined(PETSC_HAVE_DEVICE) 2906 PetscCall(MatBindToCPU(*B, A->boundtocpu)); 2907 #endif 2908 PetscCall(MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY)); 2909 PetscCall(MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY)); 2910 PetscCall(ISRestoreIndices(irowp, &row)); 2911 PetscCall(ISRestoreIndices(icolp, &col)); 2912 PetscCall(ISDestroy(&irowp)); 2913 PetscCall(ISDestroy(&icolp)); 2914 if (rowp == colp) PetscCall(MatPropagateSymmetryOptions(A, *B)); 2915 PetscFunctionReturn(0); 2916 } 2917 2918 PetscErrorCode MatCopy_SeqAIJ(Mat A, Mat B, MatStructure str) 2919 { 2920 PetscFunctionBegin; 2921 /* If the two matrices have the same copy implementation, use fast copy. */ 2922 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 2923 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 2924 Mat_SeqAIJ *b = (Mat_SeqAIJ *)B->data; 2925 const PetscScalar *aa; 2926 2927 PetscCall(MatSeqAIJGetArrayRead(A, &aa)); 2928 PetscCheck(a->i[A->rmap->n] == b->i[B->rmap->n], PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of nonzeros in two matrices are different %" PetscInt_FMT " != %" PetscInt_FMT, a->i[A->rmap->n], b->i[B->rmap->n]); 2929 PetscCall(PetscArraycpy(b->a, aa, a->i[A->rmap->n])); 2930 PetscCall(PetscObjectStateIncrease((PetscObject)B)); 2931 PetscCall(MatSeqAIJRestoreArrayRead(A, &aa)); 2932 } else { 2933 PetscCall(MatCopy_Basic(A, B, str)); 2934 } 2935 PetscFunctionReturn(0); 2936 } 2937 2938 PetscErrorCode MatSetUp_SeqAIJ(Mat A) 2939 { 2940 PetscFunctionBegin; 2941 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(A, PETSC_DEFAULT, NULL)); 2942 PetscFunctionReturn(0); 2943 } 2944 2945 PETSC_INTERN PetscErrorCode MatSeqAIJGetArray_SeqAIJ(Mat A, PetscScalar *array[]) 2946 { 2947 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 2948 2949 PetscFunctionBegin; 2950 *array = a->a; 2951 PetscFunctionReturn(0); 2952 } 2953 2954 PETSC_INTERN PetscErrorCode MatSeqAIJRestoreArray_SeqAIJ(Mat A, PetscScalar *array[]) 2955 { 2956 PetscFunctionBegin; 2957 *array = NULL; 2958 PetscFunctionReturn(0); 2959 } 2960 2961 /* 2962 Computes the number of nonzeros per row needed for preallocation when X and Y 2963 have different nonzero structure. 2964 */ 2965 PetscErrorCode MatAXPYGetPreallocation_SeqX_private(PetscInt m, const PetscInt *xi, const PetscInt *xj, const PetscInt *yi, const PetscInt *yj, PetscInt *nnz) 2966 { 2967 PetscInt i, j, k, nzx, nzy; 2968 2969 PetscFunctionBegin; 2970 /* Set the number of nonzeros in the new matrix */ 2971 for (i = 0; i < m; i++) { 2972 const PetscInt *xjj = xj + xi[i], *yjj = yj + yi[i]; 2973 nzx = xi[i + 1] - xi[i]; 2974 nzy = yi[i + 1] - yi[i]; 2975 nnz[i] = 0; 2976 for (j = 0, k = 0; j < nzx; j++) { /* Point in X */ 2977 for (; k < nzy && yjj[k] < xjj[j]; k++) nnz[i]++; /* Catch up to X */ 2978 if (k < nzy && yjj[k] == xjj[j]) k++; /* Skip duplicate */ 2979 nnz[i]++; 2980 } 2981 for (; k < nzy; k++) nnz[i]++; 2982 } 2983 PetscFunctionReturn(0); 2984 } 2985 2986 PetscErrorCode MatAXPYGetPreallocation_SeqAIJ(Mat Y, Mat X, PetscInt *nnz) 2987 { 2988 PetscInt m = Y->rmap->N; 2989 Mat_SeqAIJ *x = (Mat_SeqAIJ *)X->data; 2990 Mat_SeqAIJ *y = (Mat_SeqAIJ *)Y->data; 2991 2992 PetscFunctionBegin; 2993 /* Set the number of nonzeros in the new matrix */ 2994 PetscCall(MatAXPYGetPreallocation_SeqX_private(m, x->i, x->j, y->i, y->j, nnz)); 2995 PetscFunctionReturn(0); 2996 } 2997 2998 PetscErrorCode MatAXPY_SeqAIJ(Mat Y, PetscScalar a, Mat X, MatStructure str) 2999 { 3000 Mat_SeqAIJ *x = (Mat_SeqAIJ *)X->data, *y = (Mat_SeqAIJ *)Y->data; 3001 3002 PetscFunctionBegin; 3003 if (str == UNKNOWN_NONZERO_PATTERN || (PetscDefined(USE_DEBUG) && str == SAME_NONZERO_PATTERN)) { 3004 PetscBool e = x->nz == y->nz ? PETSC_TRUE : PETSC_FALSE; 3005 if (e) { 3006 PetscCall(PetscArraycmp(x->i, y->i, Y->rmap->n + 1, &e)); 3007 if (e) { 3008 PetscCall(PetscArraycmp(x->j, y->j, y->nz, &e)); 3009 if (e) str = SAME_NONZERO_PATTERN; 3010 } 3011 } 3012 if (!e) PetscCheck(str != SAME_NONZERO_PATTERN, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "MatStructure is not SAME_NONZERO_PATTERN"); 3013 } 3014 if (str == SAME_NONZERO_PATTERN) { 3015 const PetscScalar *xa; 3016 PetscScalar *ya, alpha = a; 3017 PetscBLASInt one = 1, bnz; 3018 3019 PetscCall(PetscBLASIntCast(x->nz, &bnz)); 3020 PetscCall(MatSeqAIJGetArray(Y, &ya)); 3021 PetscCall(MatSeqAIJGetArrayRead(X, &xa)); 3022 PetscCallBLAS("BLASaxpy", BLASaxpy_(&bnz, &alpha, xa, &one, ya, &one)); 3023 PetscCall(MatSeqAIJRestoreArrayRead(X, &xa)); 3024 PetscCall(MatSeqAIJRestoreArray(Y, &ya)); 3025 PetscCall(PetscLogFlops(2.0 * bnz)); 3026 PetscCall(MatSeqAIJInvalidateDiagonal(Y)); 3027 PetscCall(PetscObjectStateIncrease((PetscObject)Y)); 3028 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 3029 PetscCall(MatAXPY_Basic(Y, a, X, str)); 3030 } else { 3031 Mat B; 3032 PetscInt *nnz; 3033 PetscCall(PetscMalloc1(Y->rmap->N, &nnz)); 3034 PetscCall(MatCreate(PetscObjectComm((PetscObject)Y), &B)); 3035 PetscCall(PetscObjectSetName((PetscObject)B, ((PetscObject)Y)->name)); 3036 PetscCall(MatSetLayouts(B, Y->rmap, Y->cmap)); 3037 PetscCall(MatSetType(B, ((PetscObject)Y)->type_name)); 3038 PetscCall(MatAXPYGetPreallocation_SeqAIJ(Y, X, nnz)); 3039 PetscCall(MatSeqAIJSetPreallocation(B, 0, nnz)); 3040 PetscCall(MatAXPY_BasicWithPreallocation(B, Y, a, X, str)); 3041 PetscCall(MatHeaderMerge(Y, &B)); 3042 PetscCall(MatSeqAIJCheckInode(Y)); 3043 PetscCall(PetscFree(nnz)); 3044 } 3045 PetscFunctionReturn(0); 3046 } 3047 3048 PETSC_INTERN PetscErrorCode MatConjugate_SeqAIJ(Mat mat) 3049 { 3050 #if defined(PETSC_USE_COMPLEX) 3051 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 3052 PetscInt i, nz; 3053 PetscScalar *a; 3054 3055 PetscFunctionBegin; 3056 nz = aij->nz; 3057 PetscCall(MatSeqAIJGetArray(mat, &a)); 3058 for (i = 0; i < nz; i++) a[i] = PetscConj(a[i]); 3059 PetscCall(MatSeqAIJRestoreArray(mat, &a)); 3060 #else 3061 PetscFunctionBegin; 3062 #endif 3063 PetscFunctionReturn(0); 3064 } 3065 3066 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A, Vec v, PetscInt idx[]) 3067 { 3068 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 3069 PetscInt i, j, m = A->rmap->n, *ai, *aj, ncols, n; 3070 PetscReal atmp; 3071 PetscScalar *x; 3072 const MatScalar *aa, *av; 3073 3074 PetscFunctionBegin; 3075 PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix"); 3076 PetscCall(MatSeqAIJGetArrayRead(A, &av)); 3077 aa = av; 3078 ai = a->i; 3079 aj = a->j; 3080 3081 PetscCall(VecSet(v, 0.0)); 3082 PetscCall(VecGetArrayWrite(v, &x)); 3083 PetscCall(VecGetLocalSize(v, &n)); 3084 PetscCheck(n == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector"); 3085 for (i = 0; i < m; i++) { 3086 ncols = ai[1] - ai[0]; 3087 ai++; 3088 for (j = 0; j < ncols; j++) { 3089 atmp = PetscAbsScalar(*aa); 3090 if (PetscAbsScalar(x[i]) < atmp) { 3091 x[i] = atmp; 3092 if (idx) idx[i] = *aj; 3093 } 3094 aa++; 3095 aj++; 3096 } 3097 } 3098 PetscCall(VecRestoreArrayWrite(v, &x)); 3099 PetscCall(MatSeqAIJRestoreArrayRead(A, &av)); 3100 PetscFunctionReturn(0); 3101 } 3102 3103 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A, Vec v, PetscInt idx[]) 3104 { 3105 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 3106 PetscInt i, j, m = A->rmap->n, *ai, *aj, ncols, n; 3107 PetscScalar *x; 3108 const MatScalar *aa, *av; 3109 3110 PetscFunctionBegin; 3111 PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix"); 3112 PetscCall(MatSeqAIJGetArrayRead(A, &av)); 3113 aa = av; 3114 ai = a->i; 3115 aj = a->j; 3116 3117 PetscCall(VecSet(v, 0.0)); 3118 PetscCall(VecGetArrayWrite(v, &x)); 3119 PetscCall(VecGetLocalSize(v, &n)); 3120 PetscCheck(n == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector"); 3121 for (i = 0; i < m; i++) { 3122 ncols = ai[1] - ai[0]; 3123 ai++; 3124 if (ncols == A->cmap->n) { /* row is dense */ 3125 x[i] = *aa; 3126 if (idx) idx[i] = 0; 3127 } else { /* row is sparse so already KNOW maximum is 0.0 or higher */ 3128 x[i] = 0.0; 3129 if (idx) { 3130 for (j = 0; j < ncols; j++) { /* find first implicit 0.0 in the row */ 3131 if (aj[j] > j) { 3132 idx[i] = j; 3133 break; 3134 } 3135 } 3136 /* in case first implicit 0.0 in the row occurs at ncols-th column */ 3137 if (j == ncols && j < A->cmap->n) idx[i] = j; 3138 } 3139 } 3140 for (j = 0; j < ncols; j++) { 3141 if (PetscRealPart(x[i]) < PetscRealPart(*aa)) { 3142 x[i] = *aa; 3143 if (idx) idx[i] = *aj; 3144 } 3145 aa++; 3146 aj++; 3147 } 3148 } 3149 PetscCall(VecRestoreArrayWrite(v, &x)); 3150 PetscCall(MatSeqAIJRestoreArrayRead(A, &av)); 3151 PetscFunctionReturn(0); 3152 } 3153 3154 PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A, Vec v, PetscInt idx[]) 3155 { 3156 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 3157 PetscInt i, j, m = A->rmap->n, *ai, *aj, ncols, n; 3158 PetscScalar *x; 3159 const MatScalar *aa, *av; 3160 3161 PetscFunctionBegin; 3162 PetscCall(MatSeqAIJGetArrayRead(A, &av)); 3163 aa = av; 3164 ai = a->i; 3165 aj = a->j; 3166 3167 PetscCall(VecSet(v, 0.0)); 3168 PetscCall(VecGetArrayWrite(v, &x)); 3169 PetscCall(VecGetLocalSize(v, &n)); 3170 PetscCheck(n == m, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector, %" PetscInt_FMT " vs. %" PetscInt_FMT " rows", m, n); 3171 for (i = 0; i < m; i++) { 3172 ncols = ai[1] - ai[0]; 3173 ai++; 3174 if (ncols == A->cmap->n) { /* row is dense */ 3175 x[i] = *aa; 3176 if (idx) idx[i] = 0; 3177 } else { /* row is sparse so already KNOW minimum is 0.0 or higher */ 3178 x[i] = 0.0; 3179 if (idx) { /* find first implicit 0.0 in the row */ 3180 for (j = 0; j < ncols; j++) { 3181 if (aj[j] > j) { 3182 idx[i] = j; 3183 break; 3184 } 3185 } 3186 /* in case first implicit 0.0 in the row occurs at ncols-th column */ 3187 if (j == ncols && j < A->cmap->n) idx[i] = j; 3188 } 3189 } 3190 for (j = 0; j < ncols; j++) { 3191 if (PetscAbsScalar(x[i]) > PetscAbsScalar(*aa)) { 3192 x[i] = *aa; 3193 if (idx) idx[i] = *aj; 3194 } 3195 aa++; 3196 aj++; 3197 } 3198 } 3199 PetscCall(VecRestoreArrayWrite(v, &x)); 3200 PetscCall(MatSeqAIJRestoreArrayRead(A, &av)); 3201 PetscFunctionReturn(0); 3202 } 3203 3204 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A, Vec v, PetscInt idx[]) 3205 { 3206 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 3207 PetscInt i, j, m = A->rmap->n, ncols, n; 3208 const PetscInt *ai, *aj; 3209 PetscScalar *x; 3210 const MatScalar *aa, *av; 3211 3212 PetscFunctionBegin; 3213 PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix"); 3214 PetscCall(MatSeqAIJGetArrayRead(A, &av)); 3215 aa = av; 3216 ai = a->i; 3217 aj = a->j; 3218 3219 PetscCall(VecSet(v, 0.0)); 3220 PetscCall(VecGetArrayWrite(v, &x)); 3221 PetscCall(VecGetLocalSize(v, &n)); 3222 PetscCheck(n == m, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector"); 3223 for (i = 0; i < m; i++) { 3224 ncols = ai[1] - ai[0]; 3225 ai++; 3226 if (ncols == A->cmap->n) { /* row is dense */ 3227 x[i] = *aa; 3228 if (idx) idx[i] = 0; 3229 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 3230 x[i] = 0.0; 3231 if (idx) { /* find first implicit 0.0 in the row */ 3232 for (j = 0; j < ncols; j++) { 3233 if (aj[j] > j) { 3234 idx[i] = j; 3235 break; 3236 } 3237 } 3238 /* in case first implicit 0.0 in the row occurs at ncols-th column */ 3239 if (j == ncols && j < A->cmap->n) idx[i] = j; 3240 } 3241 } 3242 for (j = 0; j < ncols; j++) { 3243 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) { 3244 x[i] = *aa; 3245 if (idx) idx[i] = *aj; 3246 } 3247 aa++; 3248 aj++; 3249 } 3250 } 3251 PetscCall(VecRestoreArrayWrite(v, &x)); 3252 PetscCall(MatSeqAIJRestoreArrayRead(A, &av)); 3253 PetscFunctionReturn(0); 3254 } 3255 3256 PetscErrorCode MatInvertBlockDiagonal_SeqAIJ(Mat A, const PetscScalar **values) 3257 { 3258 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 3259 PetscInt i, bs = PetscAbs(A->rmap->bs), mbs = A->rmap->n / bs, ipvt[5], bs2 = bs * bs, *v_pivots, ij[7], *IJ, j; 3260 MatScalar *diag, work[25], *v_work; 3261 const PetscReal shift = 0.0; 3262 PetscBool allowzeropivot, zeropivotdetected = PETSC_FALSE; 3263 3264 PetscFunctionBegin; 3265 allowzeropivot = PetscNot(A->erroriffailure); 3266 if (a->ibdiagvalid) { 3267 if (values) *values = a->ibdiag; 3268 PetscFunctionReturn(0); 3269 } 3270 PetscCall(MatMarkDiagonal_SeqAIJ(A)); 3271 if (!a->ibdiag) { PetscCall(PetscMalloc1(bs2 * mbs, &a->ibdiag)); } 3272 diag = a->ibdiag; 3273 if (values) *values = a->ibdiag; 3274 /* factor and invert each block */ 3275 switch (bs) { 3276 case 1: 3277 for (i = 0; i < mbs; i++) { 3278 PetscCall(MatGetValues(A, 1, &i, 1, &i, diag + i)); 3279 if (PetscAbsScalar(diag[i] + shift) < PETSC_MACHINE_EPSILON) { 3280 if (allowzeropivot) { 3281 A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3282 A->factorerror_zeropivot_value = PetscAbsScalar(diag[i]); 3283 A->factorerror_zeropivot_row = i; 3284 PetscCall(PetscInfo(A, "Zero pivot, row %" PetscInt_FMT " pivot %g tolerance %g\n", i, (double)PetscAbsScalar(diag[i]), (double)PETSC_MACHINE_EPSILON)); 3285 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_MAT_LU_ZRPVT, "Zero pivot, row %" PetscInt_FMT " pivot %g tolerance %g", i, (double)PetscAbsScalar(diag[i]), (double)PETSC_MACHINE_EPSILON); 3286 } 3287 diag[i] = (PetscScalar)1.0 / (diag[i] + shift); 3288 } 3289 break; 3290 case 2: 3291 for (i = 0; i < mbs; i++) { 3292 ij[0] = 2 * i; 3293 ij[1] = 2 * i + 1; 3294 PetscCall(MatGetValues(A, 2, ij, 2, ij, diag)); 3295 PetscCall(PetscKernel_A_gets_inverse_A_2(diag, shift, allowzeropivot, &zeropivotdetected)); 3296 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3297 PetscCall(PetscKernel_A_gets_transpose_A_2(diag)); 3298 diag += 4; 3299 } 3300 break; 3301 case 3: 3302 for (i = 0; i < mbs; i++) { 3303 ij[0] = 3 * i; 3304 ij[1] = 3 * i + 1; 3305 ij[2] = 3 * i + 2; 3306 PetscCall(MatGetValues(A, 3, ij, 3, ij, diag)); 3307 PetscCall(PetscKernel_A_gets_inverse_A_3(diag, shift, allowzeropivot, &zeropivotdetected)); 3308 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3309 PetscCall(PetscKernel_A_gets_transpose_A_3(diag)); 3310 diag += 9; 3311 } 3312 break; 3313 case 4: 3314 for (i = 0; i < mbs; i++) { 3315 ij[0] = 4 * i; 3316 ij[1] = 4 * i + 1; 3317 ij[2] = 4 * i + 2; 3318 ij[3] = 4 * i + 3; 3319 PetscCall(MatGetValues(A, 4, ij, 4, ij, diag)); 3320 PetscCall(PetscKernel_A_gets_inverse_A_4(diag, shift, allowzeropivot, &zeropivotdetected)); 3321 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3322 PetscCall(PetscKernel_A_gets_transpose_A_4(diag)); 3323 diag += 16; 3324 } 3325 break; 3326 case 5: 3327 for (i = 0; i < mbs; i++) { 3328 ij[0] = 5 * i; 3329 ij[1] = 5 * i + 1; 3330 ij[2] = 5 * i + 2; 3331 ij[3] = 5 * i + 3; 3332 ij[4] = 5 * i + 4; 3333 PetscCall(MatGetValues(A, 5, ij, 5, ij, diag)); 3334 PetscCall(PetscKernel_A_gets_inverse_A_5(diag, ipvt, work, shift, allowzeropivot, &zeropivotdetected)); 3335 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3336 PetscCall(PetscKernel_A_gets_transpose_A_5(diag)); 3337 diag += 25; 3338 } 3339 break; 3340 case 6: 3341 for (i = 0; i < mbs; i++) { 3342 ij[0] = 6 * i; 3343 ij[1] = 6 * i + 1; 3344 ij[2] = 6 * i + 2; 3345 ij[3] = 6 * i + 3; 3346 ij[4] = 6 * i + 4; 3347 ij[5] = 6 * i + 5; 3348 PetscCall(MatGetValues(A, 6, ij, 6, ij, diag)); 3349 PetscCall(PetscKernel_A_gets_inverse_A_6(diag, shift, allowzeropivot, &zeropivotdetected)); 3350 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3351 PetscCall(PetscKernel_A_gets_transpose_A_6(diag)); 3352 diag += 36; 3353 } 3354 break; 3355 case 7: 3356 for (i = 0; i < mbs; i++) { 3357 ij[0] = 7 * i; 3358 ij[1] = 7 * i + 1; 3359 ij[2] = 7 * i + 2; 3360 ij[3] = 7 * i + 3; 3361 ij[4] = 7 * i + 4; 3362 ij[5] = 7 * i + 5; 3363 ij[5] = 7 * i + 6; 3364 PetscCall(MatGetValues(A, 7, ij, 7, ij, diag)); 3365 PetscCall(PetscKernel_A_gets_inverse_A_7(diag, shift, allowzeropivot, &zeropivotdetected)); 3366 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3367 PetscCall(PetscKernel_A_gets_transpose_A_7(diag)); 3368 diag += 49; 3369 } 3370 break; 3371 default: 3372 PetscCall(PetscMalloc3(bs, &v_work, bs, &v_pivots, bs, &IJ)); 3373 for (i = 0; i < mbs; i++) { 3374 for (j = 0; j < bs; j++) IJ[j] = bs * i + j; 3375 PetscCall(MatGetValues(A, bs, IJ, bs, IJ, diag)); 3376 PetscCall(PetscKernel_A_gets_inverse_A(bs, diag, v_pivots, v_work, allowzeropivot, &zeropivotdetected)); 3377 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3378 PetscCall(PetscKernel_A_gets_transpose_A_N(diag, bs)); 3379 diag += bs2; 3380 } 3381 PetscCall(PetscFree3(v_work, v_pivots, IJ)); 3382 } 3383 a->ibdiagvalid = PETSC_TRUE; 3384 PetscFunctionReturn(0); 3385 } 3386 3387 static PetscErrorCode MatSetRandom_SeqAIJ(Mat x, PetscRandom rctx) 3388 { 3389 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)x->data; 3390 PetscScalar a, *aa; 3391 PetscInt m, n, i, j, col; 3392 3393 PetscFunctionBegin; 3394 if (!x->assembled) { 3395 PetscCall(MatGetSize(x, &m, &n)); 3396 for (i = 0; i < m; i++) { 3397 for (j = 0; j < aij->imax[i]; j++) { 3398 PetscCall(PetscRandomGetValue(rctx, &a)); 3399 col = (PetscInt)(n * PetscRealPart(a)); 3400 PetscCall(MatSetValues(x, 1, &i, 1, &col, &a, ADD_VALUES)); 3401 } 3402 } 3403 } else { 3404 PetscCall(MatSeqAIJGetArrayWrite(x, &aa)); 3405 for (i = 0; i < aij->nz; i++) PetscCall(PetscRandomGetValue(rctx, aa + i)); 3406 PetscCall(MatSeqAIJRestoreArrayWrite(x, &aa)); 3407 } 3408 PetscCall(MatAssemblyBegin(x, MAT_FINAL_ASSEMBLY)); 3409 PetscCall(MatAssemblyEnd(x, MAT_FINAL_ASSEMBLY)); 3410 PetscFunctionReturn(0); 3411 } 3412 3413 /* Like MatSetRandom_SeqAIJ, but do not set values on columns in range of [low, high) */ 3414 PetscErrorCode MatSetRandomSkipColumnRange_SeqAIJ_Private(Mat x, PetscInt low, PetscInt high, PetscRandom rctx) 3415 { 3416 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)x->data; 3417 PetscScalar a; 3418 PetscInt m, n, i, j, col, nskip; 3419 3420 PetscFunctionBegin; 3421 nskip = high - low; 3422 PetscCall(MatGetSize(x, &m, &n)); 3423 n -= nskip; /* shrink number of columns where nonzeros can be set */ 3424 for (i = 0; i < m; i++) { 3425 for (j = 0; j < aij->imax[i]; j++) { 3426 PetscCall(PetscRandomGetValue(rctx, &a)); 3427 col = (PetscInt)(n * PetscRealPart(a)); 3428 if (col >= low) col += nskip; /* shift col rightward to skip the hole */ 3429 PetscCall(MatSetValues(x, 1, &i, 1, &col, &a, ADD_VALUES)); 3430 } 3431 } 3432 PetscCall(MatAssemblyBegin(x, MAT_FINAL_ASSEMBLY)); 3433 PetscCall(MatAssemblyEnd(x, MAT_FINAL_ASSEMBLY)); 3434 PetscFunctionReturn(0); 3435 } 3436 3437 /* -------------------------------------------------------------------*/ 3438 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 3439 MatGetRow_SeqAIJ, 3440 MatRestoreRow_SeqAIJ, 3441 MatMult_SeqAIJ, 3442 /* 4*/ MatMultAdd_SeqAIJ, 3443 MatMultTranspose_SeqAIJ, 3444 MatMultTransposeAdd_SeqAIJ, 3445 NULL, 3446 NULL, 3447 NULL, 3448 /* 10*/ NULL, 3449 MatLUFactor_SeqAIJ, 3450 NULL, 3451 MatSOR_SeqAIJ, 3452 MatTranspose_SeqAIJ, 3453 /*1 5*/ MatGetInfo_SeqAIJ, 3454 MatEqual_SeqAIJ, 3455 MatGetDiagonal_SeqAIJ, 3456 MatDiagonalScale_SeqAIJ, 3457 MatNorm_SeqAIJ, 3458 /* 20*/ NULL, 3459 MatAssemblyEnd_SeqAIJ, 3460 MatSetOption_SeqAIJ, 3461 MatZeroEntries_SeqAIJ, 3462 /* 24*/ MatZeroRows_SeqAIJ, 3463 NULL, 3464 NULL, 3465 NULL, 3466 NULL, 3467 /* 29*/ MatSetUp_SeqAIJ, 3468 NULL, 3469 NULL, 3470 NULL, 3471 NULL, 3472 /* 34*/ MatDuplicate_SeqAIJ, 3473 NULL, 3474 NULL, 3475 MatILUFactor_SeqAIJ, 3476 NULL, 3477 /* 39*/ MatAXPY_SeqAIJ, 3478 MatCreateSubMatrices_SeqAIJ, 3479 MatIncreaseOverlap_SeqAIJ, 3480 MatGetValues_SeqAIJ, 3481 MatCopy_SeqAIJ, 3482 /* 44*/ MatGetRowMax_SeqAIJ, 3483 MatScale_SeqAIJ, 3484 MatShift_SeqAIJ, 3485 MatDiagonalSet_SeqAIJ, 3486 MatZeroRowsColumns_SeqAIJ, 3487 /* 49*/ MatSetRandom_SeqAIJ, 3488 MatGetRowIJ_SeqAIJ, 3489 MatRestoreRowIJ_SeqAIJ, 3490 MatGetColumnIJ_SeqAIJ, 3491 MatRestoreColumnIJ_SeqAIJ, 3492 /* 54*/ MatFDColoringCreate_SeqXAIJ, 3493 NULL, 3494 NULL, 3495 MatPermute_SeqAIJ, 3496 NULL, 3497 /* 59*/ NULL, 3498 MatDestroy_SeqAIJ, 3499 MatView_SeqAIJ, 3500 NULL, 3501 NULL, 3502 /* 64*/ NULL, 3503 MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ, 3504 NULL, 3505 NULL, 3506 NULL, 3507 /* 69*/ MatGetRowMaxAbs_SeqAIJ, 3508 MatGetRowMinAbs_SeqAIJ, 3509 NULL, 3510 NULL, 3511 NULL, 3512 /* 74*/ NULL, 3513 MatFDColoringApply_AIJ, 3514 NULL, 3515 NULL, 3516 NULL, 3517 /* 79*/ MatFindZeroDiagonals_SeqAIJ, 3518 NULL, 3519 NULL, 3520 NULL, 3521 MatLoad_SeqAIJ, 3522 /* 84*/ MatIsSymmetric_SeqAIJ, 3523 MatIsHermitian_SeqAIJ, 3524 NULL, 3525 NULL, 3526 NULL, 3527 /* 89*/ NULL, 3528 NULL, 3529 MatMatMultNumeric_SeqAIJ_SeqAIJ, 3530 NULL, 3531 NULL, 3532 /* 94*/ MatPtAPNumeric_SeqAIJ_SeqAIJ_SparseAxpy, 3533 NULL, 3534 NULL, 3535 MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ, 3536 NULL, 3537 /* 99*/ MatProductSetFromOptions_SeqAIJ, 3538 NULL, 3539 NULL, 3540 MatConjugate_SeqAIJ, 3541 NULL, 3542 /*104*/ MatSetValuesRow_SeqAIJ, 3543 MatRealPart_SeqAIJ, 3544 MatImaginaryPart_SeqAIJ, 3545 NULL, 3546 NULL, 3547 /*109*/ MatMatSolve_SeqAIJ, 3548 NULL, 3549 MatGetRowMin_SeqAIJ, 3550 NULL, 3551 MatMissingDiagonal_SeqAIJ, 3552 /*114*/ NULL, 3553 NULL, 3554 NULL, 3555 NULL, 3556 NULL, 3557 /*119*/ NULL, 3558 NULL, 3559 NULL, 3560 NULL, 3561 MatGetMultiProcBlock_SeqAIJ, 3562 /*124*/ MatFindNonzeroRows_SeqAIJ, 3563 MatGetColumnReductions_SeqAIJ, 3564 MatInvertBlockDiagonal_SeqAIJ, 3565 MatInvertVariableBlockDiagonal_SeqAIJ, 3566 NULL, 3567 /*129*/ NULL, 3568 NULL, 3569 NULL, 3570 MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ, 3571 MatTransposeColoringCreate_SeqAIJ, 3572 /*134*/ MatTransColoringApplySpToDen_SeqAIJ, 3573 MatTransColoringApplyDenToSp_SeqAIJ, 3574 NULL, 3575 NULL, 3576 MatRARtNumeric_SeqAIJ_SeqAIJ, 3577 /*139*/ NULL, 3578 NULL, 3579 NULL, 3580 MatFDColoringSetUp_SeqXAIJ, 3581 MatFindOffBlockDiagonalEntries_SeqAIJ, 3582 MatCreateMPIMatConcatenateSeqMat_SeqAIJ, 3583 /*145*/ MatDestroySubMatrices_SeqAIJ, 3584 NULL, 3585 NULL, 3586 MatCreateGraph_Simple_AIJ, 3587 NULL, 3588 /*150*/ MatTransposeSymbolic_SeqAIJ, 3589 MatEliminateZeros_SeqAIJ}; 3590 3591 PetscErrorCode MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat, PetscInt *indices) 3592 { 3593 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 3594 PetscInt i, nz, n; 3595 3596 PetscFunctionBegin; 3597 nz = aij->maxnz; 3598 n = mat->rmap->n; 3599 for (i = 0; i < nz; i++) aij->j[i] = indices[i]; 3600 aij->nz = nz; 3601 for (i = 0; i < n; i++) aij->ilen[i] = aij->imax[i]; 3602 PetscFunctionReturn(0); 3603 } 3604 3605 /* 3606 * Given a sparse matrix with global column indices, compact it by using a local column space. 3607 * The result matrix helps saving memory in other algorithms, such as MatPtAPSymbolic_MPIAIJ_MPIAIJ_scalable() 3608 */ 3609 PetscErrorCode MatSeqAIJCompactOutExtraColumns_SeqAIJ(Mat mat, ISLocalToGlobalMapping *mapping) 3610 { 3611 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 3612 PetscHMapI gid1_lid1; 3613 PetscHashIter tpos; 3614 PetscInt gid, lid, i, ec, nz = aij->nz; 3615 PetscInt *garray, *jj = aij->j; 3616 3617 PetscFunctionBegin; 3618 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 3619 PetscValidPointer(mapping, 2); 3620 /* use a table */ 3621 PetscCall(PetscHMapICreateWithSize(mat->rmap->n, &gid1_lid1)); 3622 ec = 0; 3623 for (i = 0; i < nz; i++) { 3624 PetscInt data, gid1 = jj[i] + 1; 3625 PetscCall(PetscHMapIGetWithDefault(gid1_lid1, gid1, 0, &data)); 3626 if (!data) { 3627 /* one based table */ 3628 PetscCall(PetscHMapISet(gid1_lid1, gid1, ++ec)); 3629 } 3630 } 3631 /* form array of columns we need */ 3632 PetscCall(PetscMalloc1(ec, &garray)); 3633 PetscHashIterBegin(gid1_lid1, tpos); 3634 while (!PetscHashIterAtEnd(gid1_lid1, tpos)) { 3635 PetscHashIterGetKey(gid1_lid1, tpos, gid); 3636 PetscHashIterGetVal(gid1_lid1, tpos, lid); 3637 PetscHashIterNext(gid1_lid1, tpos); 3638 gid--; 3639 lid--; 3640 garray[lid] = gid; 3641 } 3642 PetscCall(PetscSortInt(ec, garray)); /* sort, and rebuild */ 3643 PetscCall(PetscHMapIClear(gid1_lid1)); 3644 for (i = 0; i < ec; i++) PetscCall(PetscHMapISet(gid1_lid1, garray[i] + 1, i + 1)); 3645 /* compact out the extra columns in B */ 3646 for (i = 0; i < nz; i++) { 3647 PetscInt gid1 = jj[i] + 1; 3648 PetscCall(PetscHMapIGetWithDefault(gid1_lid1, gid1, 0, &lid)); 3649 lid--; 3650 jj[i] = lid; 3651 } 3652 PetscCall(PetscLayoutDestroy(&mat->cmap)); 3653 PetscCall(PetscHMapIDestroy(&gid1_lid1)); 3654 PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)mat), ec, ec, 1, &mat->cmap)); 3655 PetscCall(ISLocalToGlobalMappingCreate(PETSC_COMM_SELF, mat->cmap->bs, mat->cmap->n, garray, PETSC_OWN_POINTER, mapping)); 3656 PetscCall(ISLocalToGlobalMappingSetType(*mapping, ISLOCALTOGLOBALMAPPINGHASH)); 3657 PetscFunctionReturn(0); 3658 } 3659 3660 /*@ 3661 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 3662 in the matrix. 3663 3664 Input Parameters: 3665 + mat - the `MATSEQAIJ` matrix 3666 - indices - the column indices 3667 3668 Level: advanced 3669 3670 Notes: 3671 This can be called if you have precomputed the nonzero structure of the 3672 matrix and want to provide it to the matrix object to improve the performance 3673 of the `MatSetValues()` operation. 3674 3675 You MUST have set the correct numbers of nonzeros per row in the call to 3676 `MatCreateSeqAIJ()`, and the columns indices MUST be sorted. 3677 3678 MUST be called before any calls to `MatSetValues()` 3679 3680 The indices should start with zero, not one. 3681 3682 @*/ 3683 PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat, PetscInt *indices) 3684 { 3685 PetscFunctionBegin; 3686 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 3687 PetscValidIntPointer(indices, 2); 3688 PetscUseMethod(mat, "MatSeqAIJSetColumnIndices_C", (Mat, PetscInt *), (mat, indices)); 3689 PetscFunctionReturn(0); 3690 } 3691 3692 /* ----------------------------------------------------------------------------------------*/ 3693 3694 PetscErrorCode MatStoreValues_SeqAIJ(Mat mat) 3695 { 3696 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 3697 size_t nz = aij->i[mat->rmap->n]; 3698 3699 PetscFunctionBegin; 3700 PetscCheck(aij->nonew, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3701 3702 /* allocate space for values if not already there */ 3703 if (!aij->saved_values) { PetscCall(PetscMalloc1(nz + 1, &aij->saved_values)); } 3704 3705 /* copy values over */ 3706 PetscCall(PetscArraycpy(aij->saved_values, aij->a, nz)); 3707 PetscFunctionReturn(0); 3708 } 3709 3710 /*@ 3711 MatStoreValues - Stashes a copy of the matrix values; this allows, for 3712 example, reuse of the linear part of a Jacobian, while recomputing the 3713 nonlinear portion. 3714 3715 Logically Collect 3716 3717 Input Parameters: 3718 . mat - the matrix (currently only `MATAIJ` matrices support this option) 3719 3720 Level: advanced 3721 3722 Common Usage, with `SNESSolve()`: 3723 $ Create Jacobian matrix 3724 $ Set linear terms into matrix 3725 $ Apply boundary conditions to matrix, at this time matrix must have 3726 $ final nonzero structure (i.e. setting the nonlinear terms and applying 3727 $ boundary conditions again will not change the nonzero structure 3728 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3729 $ ierr = MatStoreValues(mat); 3730 $ Call SNESSetJacobian() with matrix 3731 $ In your Jacobian routine 3732 $ ierr = MatRetrieveValues(mat); 3733 $ Set nonlinear terms in matrix 3734 3735 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 3736 $ // build linear portion of Jacobian 3737 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3738 $ ierr = MatStoreValues(mat); 3739 $ loop over nonlinear iterations 3740 $ ierr = MatRetrieveValues(mat); 3741 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 3742 $ // call MatAssemblyBegin/End() on matrix 3743 $ Solve linear system with Jacobian 3744 $ endloop 3745 3746 Notes: 3747 Matrix must already be assemblied before calling this routine 3748 Must set the matrix option `MatSetOption`(mat,`MAT_NEW_NONZERO_LOCATIONS`,`PETSC_FALSE`); before 3749 calling this routine. 3750 3751 When this is called multiple times it overwrites the previous set of stored values 3752 and does not allocated additional space. 3753 3754 .seealso: `MatRetrieveValues()` 3755 @*/ 3756 PetscErrorCode MatStoreValues(Mat mat) 3757 { 3758 PetscFunctionBegin; 3759 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 3760 PetscCheck(mat->assembled, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix"); 3761 PetscCheck(!mat->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix"); 3762 PetscUseMethod(mat, "MatStoreValues_C", (Mat), (mat)); 3763 PetscFunctionReturn(0); 3764 } 3765 3766 PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat) 3767 { 3768 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 3769 PetscInt nz = aij->i[mat->rmap->n]; 3770 3771 PetscFunctionBegin; 3772 PetscCheck(aij->nonew, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3773 PetscCheck(aij->saved_values, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call MatStoreValues(A);first"); 3774 /* copy values over */ 3775 PetscCall(PetscArraycpy(aij->a, aij->saved_values, nz)); 3776 PetscFunctionReturn(0); 3777 } 3778 3779 /*@ 3780 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 3781 example, reuse of the linear part of a Jacobian, while recomputing the 3782 nonlinear portion. 3783 3784 Logically Collect 3785 3786 Input Parameters: 3787 . mat - the matrix (currently only `MATAIJ` matrices support this option) 3788 3789 Level: advanced 3790 3791 .seealso: `MatStoreValues()` 3792 @*/ 3793 PetscErrorCode MatRetrieveValues(Mat mat) 3794 { 3795 PetscFunctionBegin; 3796 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 3797 PetscCheck(mat->assembled, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix"); 3798 PetscCheck(!mat->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix"); 3799 PetscUseMethod(mat, "MatRetrieveValues_C", (Mat), (mat)); 3800 PetscFunctionReturn(0); 3801 } 3802 3803 /* --------------------------------------------------------------------------------*/ 3804 /*@C 3805 MatCreateSeqAIJ - Creates a sparse matrix in `MATSEQAIJ` (compressed row) format 3806 (the default parallel PETSc format). For good matrix assembly performance 3807 the user should preallocate the matrix storage by setting the parameter nz 3808 (or the array nnz). By setting these parameters accurately, performance 3809 during matrix assembly can be increased by more than a factor of 50. 3810 3811 Collective 3812 3813 Input Parameters: 3814 + comm - MPI communicator, set to `PETSC_COMM_SELF` 3815 . m - number of rows 3816 . n - number of columns 3817 . nz - number of nonzeros per row (same for all rows) 3818 - nnz - array containing the number of nonzeros in the various rows 3819 (possibly different for each row) or NULL 3820 3821 Output Parameter: 3822 . A - the matrix 3823 3824 It is recommended that one use the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`, 3825 MatXXXXSetPreallocation() paradigm instead of this routine directly. 3826 [MatXXXXSetPreallocation() is, for example, `MatSeqAIJSetPreallocation()`] 3827 3828 Notes: 3829 If nnz is given then nz is ignored 3830 3831 The AIJ format, also called 3832 compressed row storage, is fully compatible with standard Fortran 77 3833 storage. That is, the stored row and column indices can begin at 3834 either one (as in Fortran) or zero. See the users' manual for details. 3835 3836 Specify the preallocated storage with either nz or nnz (not both). 3837 Set nz = `PETSC_DEFAULT` and nnz = NULL for PETSc to control dynamic memory 3838 allocation. For large problems you MUST preallocate memory or you 3839 will get TERRIBLE performance, see the users' manual chapter on matrices. 3840 3841 By default, this format uses inodes (identical nodes) when possible, to 3842 improve numerical efficiency of matrix-vector products and solves. We 3843 search for consecutive rows with the same nonzero structure, thereby 3844 reusing matrix information to achieve increased efficiency. 3845 3846 Options Database Keys: 3847 + -mat_no_inode - Do not use inodes 3848 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3849 3850 Level: intermediate 3851 3852 .seealso: [Sparse Matrix Creation](sec_matsparse), `MatCreate()`, `MatCreateAIJ()`, `MatSetValues()`, `MatSeqAIJSetColumnIndices()`, `MatCreateSeqAIJWithArrays()` 3853 @*/ 3854 PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt nz, const PetscInt nnz[], Mat *A) 3855 { 3856 PetscFunctionBegin; 3857 PetscCall(MatCreate(comm, A)); 3858 PetscCall(MatSetSizes(*A, m, n, m, n)); 3859 PetscCall(MatSetType(*A, MATSEQAIJ)); 3860 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(*A, nz, nnz)); 3861 PetscFunctionReturn(0); 3862 } 3863 3864 /*@C 3865 MatSeqAIJSetPreallocation - For good matrix assembly performance 3866 the user should preallocate the matrix storage by setting the parameter nz 3867 (or the array nnz). By setting these parameters accurately, performance 3868 during matrix assembly can be increased by more than a factor of 50. 3869 3870 Collective 3871 3872 Input Parameters: 3873 + B - The matrix 3874 . nz - number of nonzeros per row (same for all rows) 3875 - nnz - array containing the number of nonzeros in the various rows 3876 (possibly different for each row) or NULL 3877 3878 Notes: 3879 If nnz is given then nz is ignored 3880 3881 The `MATSEQAIJ` format also called 3882 compressed row storage, is fully compatible with standard Fortran 77 3883 storage. That is, the stored row and column indices can begin at 3884 either one (as in Fortran) or zero. See the users' manual for details. 3885 3886 Specify the preallocated storage with either nz or nnz (not both). 3887 Set nz = `PETSC_DEFAULT` and nnz = NULL for PETSc to control dynamic memory 3888 allocation. For large problems you MUST preallocate memory or you 3889 will get TERRIBLE performance, see the users' manual chapter on matrices. 3890 3891 You can call `MatGetInfo()` to get information on how effective the preallocation was; 3892 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 3893 You can also run with the option -info and look for messages with the string 3894 malloc in them to see if additional memory allocation was needed. 3895 3896 Developer Notes: 3897 Use nz of `MAT_SKIP_ALLOCATION` to not allocate any space for the matrix 3898 entries or columns indices 3899 3900 By default, this format uses inodes (identical nodes) when possible, to 3901 improve numerical efficiency of matrix-vector products and solves. We 3902 search for consecutive rows with the same nonzero structure, thereby 3903 reusing matrix information to achieve increased efficiency. 3904 3905 Options Database Keys: 3906 + -mat_no_inode - Do not use inodes 3907 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3908 3909 Level: intermediate 3910 3911 .seealso: `MatCreate()`, `MatCreateAIJ()`, `MatSetValues()`, `MatSeqAIJSetColumnIndices()`, `MatCreateSeqAIJWithArrays()`, `MatGetInfo()`, 3912 `MatSeqAIJSetTotalPreallocation()` 3913 @*/ 3914 PetscErrorCode MatSeqAIJSetPreallocation(Mat B, PetscInt nz, const PetscInt nnz[]) 3915 { 3916 PetscFunctionBegin; 3917 PetscValidHeaderSpecific(B, MAT_CLASSID, 1); 3918 PetscValidType(B, 1); 3919 PetscTryMethod(B, "MatSeqAIJSetPreallocation_C", (Mat, PetscInt, const PetscInt[]), (B, nz, nnz)); 3920 PetscFunctionReturn(0); 3921 } 3922 3923 PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B, PetscInt nz, const PetscInt *nnz) 3924 { 3925 Mat_SeqAIJ *b; 3926 PetscBool skipallocation = PETSC_FALSE, realalloc = PETSC_FALSE; 3927 PetscInt i; 3928 3929 PetscFunctionBegin; 3930 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 3931 if (nz == MAT_SKIP_ALLOCATION) { 3932 skipallocation = PETSC_TRUE; 3933 nz = 0; 3934 } 3935 PetscCall(PetscLayoutSetUp(B->rmap)); 3936 PetscCall(PetscLayoutSetUp(B->cmap)); 3937 3938 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3939 PetscCheck(nz >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "nz cannot be less than 0: value %" PetscInt_FMT, nz); 3940 if (PetscUnlikelyDebug(nnz)) { 3941 for (i = 0; i < B->rmap->n; i++) { 3942 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]); 3943 PetscCheck(nnz[i] <= B->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "nnz cannot be greater than row length: local row %" PetscInt_FMT " value %" PetscInt_FMT " rowlength %" PetscInt_FMT, i, nnz[i], B->cmap->n); 3944 } 3945 } 3946 3947 B->preallocated = PETSC_TRUE; 3948 3949 b = (Mat_SeqAIJ *)B->data; 3950 3951 if (!skipallocation) { 3952 if (!b->imax) { PetscCall(PetscMalloc1(B->rmap->n, &b->imax)); } 3953 if (!b->ilen) { 3954 /* b->ilen will count nonzeros in each row so far. */ 3955 PetscCall(PetscCalloc1(B->rmap->n, &b->ilen)); 3956 } else { 3957 PetscCall(PetscMemzero(b->ilen, B->rmap->n * sizeof(PetscInt))); 3958 } 3959 if (!b->ipre) { PetscCall(PetscMalloc1(B->rmap->n, &b->ipre)); } 3960 if (!nnz) { 3961 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 3962 else if (nz < 0) nz = 1; 3963 nz = PetscMin(nz, B->cmap->n); 3964 for (i = 0; i < B->rmap->n; i++) b->imax[i] = nz; 3965 nz = nz * B->rmap->n; 3966 } else { 3967 PetscInt64 nz64 = 0; 3968 for (i = 0; i < B->rmap->n; i++) { 3969 b->imax[i] = nnz[i]; 3970 nz64 += nnz[i]; 3971 } 3972 PetscCall(PetscIntCast(nz64, &nz)); 3973 } 3974 3975 /* allocate the matrix space */ 3976 /* FIXME: should B's old memory be unlogged? */ 3977 PetscCall(MatSeqXAIJFreeAIJ(B, &b->a, &b->j, &b->i)); 3978 if (B->structure_only) { 3979 PetscCall(PetscMalloc1(nz, &b->j)); 3980 PetscCall(PetscMalloc1(B->rmap->n + 1, &b->i)); 3981 } else { 3982 PetscCall(PetscMalloc3(nz, &b->a, nz, &b->j, B->rmap->n + 1, &b->i)); 3983 } 3984 b->i[0] = 0; 3985 for (i = 1; i < B->rmap->n + 1; i++) b->i[i] = b->i[i - 1] + b->imax[i - 1]; 3986 if (B->structure_only) { 3987 b->singlemalloc = PETSC_FALSE; 3988 b->free_a = PETSC_FALSE; 3989 } else { 3990 b->singlemalloc = PETSC_TRUE; 3991 b->free_a = PETSC_TRUE; 3992 } 3993 b->free_ij = PETSC_TRUE; 3994 } else { 3995 b->free_a = PETSC_FALSE; 3996 b->free_ij = PETSC_FALSE; 3997 } 3998 3999 if (b->ipre && nnz != b->ipre && b->imax) { 4000 /* reserve user-requested sparsity */ 4001 PetscCall(PetscArraycpy(b->ipre, b->imax, B->rmap->n)); 4002 } 4003 4004 b->nz = 0; 4005 b->maxnz = nz; 4006 B->info.nz_unneeded = (double)b->maxnz; 4007 if (realalloc) PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE)); 4008 B->was_assembled = PETSC_FALSE; 4009 B->assembled = PETSC_FALSE; 4010 /* We simply deem preallocation has changed nonzero state. Updating the state 4011 will give clients (like AIJKokkos) a chance to know something has happened. 4012 */ 4013 B->nonzerostate++; 4014 PetscFunctionReturn(0); 4015 } 4016 4017 PetscErrorCode MatResetPreallocation_SeqAIJ(Mat A) 4018 { 4019 Mat_SeqAIJ *a; 4020 PetscInt i; 4021 4022 PetscFunctionBegin; 4023 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 4024 4025 /* Check local size. If zero, then return */ 4026 if (!A->rmap->n) PetscFunctionReturn(0); 4027 4028 a = (Mat_SeqAIJ *)A->data; 4029 /* if no saved info, we error out */ 4030 PetscCheck(a->ipre, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "No saved preallocation info "); 4031 4032 PetscCheck(a->i && a->j && a->a && a->imax && a->ilen, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Memory info is incomplete, and can not reset preallocation "); 4033 4034 PetscCall(PetscArraycpy(a->imax, a->ipre, A->rmap->n)); 4035 PetscCall(PetscArrayzero(a->ilen, A->rmap->n)); 4036 a->i[0] = 0; 4037 for (i = 1; i < A->rmap->n + 1; i++) a->i[i] = a->i[i - 1] + a->imax[i - 1]; 4038 A->preallocated = PETSC_TRUE; 4039 a->nz = 0; 4040 a->maxnz = a->i[A->rmap->n]; 4041 A->info.nz_unneeded = (double)a->maxnz; 4042 A->was_assembled = PETSC_FALSE; 4043 A->assembled = PETSC_FALSE; 4044 PetscFunctionReturn(0); 4045 } 4046 4047 /*@ 4048 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in `MATSEQAIJ` format. 4049 4050 Input Parameters: 4051 + B - the matrix 4052 . i - the indices into j for the start of each row (starts with zero) 4053 . j - the column indices for each row (starts with zero) these must be sorted for each row 4054 - v - optional values in the matrix 4055 4056 Level: developer 4057 4058 Notes: 4059 The i,j,v values are COPIED with this routine; to avoid the copy use `MatCreateSeqAIJWithArrays()` 4060 4061 This routine may be called multiple times with different nonzero patterns (or the same nonzero pattern). The nonzero 4062 structure will be the union of all the previous nonzero structures. 4063 4064 Developer Notes: 4065 An optimization could be added to the implementation where it checks if the i, and j are identical to the current i and j and 4066 then just copies the v values directly with `PetscMemcpy()`. 4067 4068 This routine could also take a `PetscCopyMode` argument to allow sharing the values instead of always copying them. 4069 4070 .seealso: `MatCreate()`, `MatCreateSeqAIJ()`, `MatSetValues()`, `MatSeqAIJSetPreallocation()`, `MatCreateSeqAIJ()`, `MATSEQAIJ`, `MatResetPreallocation()` 4071 @*/ 4072 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B, const PetscInt i[], const PetscInt j[], const PetscScalar v[]) 4073 { 4074 PetscFunctionBegin; 4075 PetscValidHeaderSpecific(B, MAT_CLASSID, 1); 4076 PetscValidType(B, 1); 4077 PetscTryMethod(B, "MatSeqAIJSetPreallocationCSR_C", (Mat, const PetscInt[], const PetscInt[], const PetscScalar[]), (B, i, j, v)); 4078 PetscFunctionReturn(0); 4079 } 4080 4081 PetscErrorCode MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B, const PetscInt Ii[], const PetscInt J[], const PetscScalar v[]) 4082 { 4083 PetscInt i; 4084 PetscInt m, n; 4085 PetscInt nz; 4086 PetscInt *nnz; 4087 4088 PetscFunctionBegin; 4089 PetscCheck(Ii[0] == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %" PetscInt_FMT, Ii[0]); 4090 4091 PetscCall(PetscLayoutSetUp(B->rmap)); 4092 PetscCall(PetscLayoutSetUp(B->cmap)); 4093 4094 PetscCall(MatGetSize(B, &m, &n)); 4095 PetscCall(PetscMalloc1(m + 1, &nnz)); 4096 for (i = 0; i < m; i++) { 4097 nz = Ii[i + 1] - Ii[i]; 4098 PetscCheck(nz >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local row %" PetscInt_FMT " has a negative number of columns %" PetscInt_FMT, i, nz); 4099 nnz[i] = nz; 4100 } 4101 PetscCall(MatSeqAIJSetPreallocation(B, 0, nnz)); 4102 PetscCall(PetscFree(nnz)); 4103 4104 for (i = 0; i < m; i++) PetscCall(MatSetValues_SeqAIJ(B, 1, &i, Ii[i + 1] - Ii[i], J + Ii[i], v ? v + Ii[i] : NULL, INSERT_VALUES)); 4105 4106 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 4107 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 4108 4109 PetscCall(MatSetOption(B, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE)); 4110 PetscFunctionReturn(0); 4111 } 4112 4113 /*@ 4114 MatSeqAIJKron - Computes C, the Kronecker product of A and B. 4115 4116 Input Parameters: 4117 + A - left-hand side matrix 4118 . B - right-hand side matrix 4119 - reuse - either `MAT_INITIAL_MATRIX` or `MAT_REUSE_MATRIX` 4120 4121 Output Parameter: 4122 . C - Kronecker product of A and B 4123 4124 Level: intermediate 4125 4126 Note: 4127 `MAT_REUSE_MATRIX` can only be used when the nonzero structure of the product matrix has not changed from that last call to `MatSeqAIJKron()`. 4128 4129 .seealso: `MatCreateSeqAIJ()`, `MATSEQAIJ`, `MATKAIJ`, `MatReuse` 4130 @*/ 4131 PetscErrorCode MatSeqAIJKron(Mat A, Mat B, MatReuse reuse, Mat *C) 4132 { 4133 PetscFunctionBegin; 4134 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 4135 PetscValidType(A, 1); 4136 PetscValidHeaderSpecific(B, MAT_CLASSID, 2); 4137 PetscValidType(B, 2); 4138 PetscValidPointer(C, 4); 4139 if (reuse == MAT_REUSE_MATRIX) { 4140 PetscValidHeaderSpecific(*C, MAT_CLASSID, 4); 4141 PetscValidType(*C, 4); 4142 } 4143 PetscTryMethod(A, "MatSeqAIJKron_C", (Mat, Mat, MatReuse, Mat *), (A, B, reuse, C)); 4144 PetscFunctionReturn(0); 4145 } 4146 4147 PetscErrorCode MatSeqAIJKron_SeqAIJ(Mat A, Mat B, MatReuse reuse, Mat *C) 4148 { 4149 Mat newmat; 4150 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 4151 Mat_SeqAIJ *b = (Mat_SeqAIJ *)B->data; 4152 PetscScalar *v; 4153 const PetscScalar *aa, *ba; 4154 PetscInt *i, *j, m, n, p, q, nnz = 0, am = A->rmap->n, bm = B->rmap->n, an = A->cmap->n, bn = B->cmap->n; 4155 PetscBool flg; 4156 4157 PetscFunctionBegin; 4158 PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix"); 4159 PetscCheck(A->assembled, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix"); 4160 PetscCheck(!B->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix"); 4161 PetscCheck(B->assembled, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix"); 4162 PetscCall(PetscObjectTypeCompare((PetscObject)B, MATSEQAIJ, &flg)); 4163 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "MatType %s", ((PetscObject)B)->type_name); 4164 PetscCheck(reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX, PETSC_COMM_SELF, PETSC_ERR_SUP, "MatReuse %d", (int)reuse); 4165 if (reuse == MAT_INITIAL_MATRIX) { 4166 PetscCall(PetscMalloc2(am * bm + 1, &i, a->i[am] * b->i[bm], &j)); 4167 PetscCall(MatCreate(PETSC_COMM_SELF, &newmat)); 4168 PetscCall(MatSetSizes(newmat, am * bm, an * bn, am * bm, an * bn)); 4169 PetscCall(MatSetType(newmat, MATAIJ)); 4170 i[0] = 0; 4171 for (m = 0; m < am; ++m) { 4172 for (p = 0; p < bm; ++p) { 4173 i[m * bm + p + 1] = i[m * bm + p] + (a->i[m + 1] - a->i[m]) * (b->i[p + 1] - b->i[p]); 4174 for (n = a->i[m]; n < a->i[m + 1]; ++n) { 4175 for (q = b->i[p]; q < b->i[p + 1]; ++q) j[nnz++] = a->j[n] * bn + b->j[q]; 4176 } 4177 } 4178 } 4179 PetscCall(MatSeqAIJSetPreallocationCSR(newmat, i, j, NULL)); 4180 *C = newmat; 4181 PetscCall(PetscFree2(i, j)); 4182 nnz = 0; 4183 } 4184 PetscCall(MatSeqAIJGetArray(*C, &v)); 4185 PetscCall(MatSeqAIJGetArrayRead(A, &aa)); 4186 PetscCall(MatSeqAIJGetArrayRead(B, &ba)); 4187 for (m = 0; m < am; ++m) { 4188 for (p = 0; p < bm; ++p) { 4189 for (n = a->i[m]; n < a->i[m + 1]; ++n) { 4190 for (q = b->i[p]; q < b->i[p + 1]; ++q) v[nnz++] = aa[n] * ba[q]; 4191 } 4192 } 4193 } 4194 PetscCall(MatSeqAIJRestoreArray(*C, &v)); 4195 PetscCall(MatSeqAIJRestoreArrayRead(A, &aa)); 4196 PetscCall(MatSeqAIJRestoreArrayRead(B, &ba)); 4197 PetscFunctionReturn(0); 4198 } 4199 4200 #include <../src/mat/impls/dense/seq/dense.h> 4201 #include <petsc/private/kernels/petscaxpy.h> 4202 4203 /* 4204 Computes (B'*A')' since computing B*A directly is untenable 4205 4206 n p p 4207 [ ] [ ] [ ] 4208 m [ A ] * n [ B ] = m [ C ] 4209 [ ] [ ] [ ] 4210 4211 */ 4212 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A, Mat B, Mat C) 4213 { 4214 Mat_SeqDense *sub_a = (Mat_SeqDense *)A->data; 4215 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ *)B->data; 4216 Mat_SeqDense *sub_c = (Mat_SeqDense *)C->data; 4217 PetscInt i, j, n, m, q, p; 4218 const PetscInt *ii, *idx; 4219 const PetscScalar *b, *a, *a_q; 4220 PetscScalar *c, *c_q; 4221 PetscInt clda = sub_c->lda; 4222 PetscInt alda = sub_a->lda; 4223 4224 PetscFunctionBegin; 4225 m = A->rmap->n; 4226 n = A->cmap->n; 4227 p = B->cmap->n; 4228 a = sub_a->v; 4229 b = sub_b->a; 4230 c = sub_c->v; 4231 if (clda == m) { 4232 PetscCall(PetscArrayzero(c, m * p)); 4233 } else { 4234 for (j = 0; j < p; j++) 4235 for (i = 0; i < m; i++) c[j * clda + i] = 0.0; 4236 } 4237 ii = sub_b->i; 4238 idx = sub_b->j; 4239 for (i = 0; i < n; i++) { 4240 q = ii[i + 1] - ii[i]; 4241 while (q-- > 0) { 4242 c_q = c + clda * (*idx); 4243 a_q = a + alda * i; 4244 PetscKernelAXPY(c_q, *b, a_q, m); 4245 idx++; 4246 b++; 4247 } 4248 } 4249 PetscFunctionReturn(0); 4250 } 4251 4252 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A, Mat B, PetscReal fill, Mat C) 4253 { 4254 PetscInt m = A->rmap->n, n = B->cmap->n; 4255 PetscBool cisdense; 4256 4257 PetscFunctionBegin; 4258 PetscCheck(A->cmap->n == B->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "A->cmap->n %" PetscInt_FMT " != B->rmap->n %" PetscInt_FMT, A->cmap->n, B->rmap->n); 4259 PetscCall(MatSetSizes(C, m, n, m, n)); 4260 PetscCall(MatSetBlockSizesFromMats(C, A, B)); 4261 PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSECUDA, MATSEQDENSEHIP, "")); 4262 if (!cisdense) PetscCall(MatSetType(C, MATDENSE)); 4263 PetscCall(MatSetUp(C)); 4264 4265 C->ops->matmultnumeric = MatMatMultNumeric_SeqDense_SeqAIJ; 4266 PetscFunctionReturn(0); 4267 } 4268 4269 /* ----------------------------------------------------------------*/ 4270 /*MC 4271 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 4272 based on compressed sparse row format. 4273 4274 Options Database Keys: 4275 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 4276 4277 Level: beginner 4278 4279 Notes: 4280 `MatSetValues()` may be called for this matrix type with a NULL argument for the numerical values, 4281 in this case the values associated with the rows and columns one passes in are set to zero 4282 in the matrix 4283 4284 `MatSetOptions`(,`MAT_STRUCTURE_ONLY`,`PETSC_TRUE`) may be called for this matrix type. In this no 4285 space is allocated for the nonzero entries and any entries passed with `MatSetValues()` are ignored 4286 4287 Developer Note: 4288 It would be nice if all matrix formats supported passing NULL in for the numerical values 4289 4290 .seealso: `MatCreateSeqAIJ()`, `MatSetFromOptions()`, `MatSetType()`, `MatCreate()`, `MatType`, `MATSELL`, `MATSEQSELL`, `MATMPISELL` 4291 M*/ 4292 4293 /*MC 4294 MATAIJ - MATAIJ = "aij" - A matrix type to be used for sparse matrices. 4295 4296 This matrix type is identical to `MATSEQAIJ` when constructed with a single process communicator, 4297 and `MATMPIAIJ` otherwise. As a result, for single process communicators, 4298 `MatSeqAIJSetPreallocation()` is supported, and similarly `MatMPIAIJSetPreallocation()` is supported 4299 for communicators controlling multiple processes. It is recommended that you call both of 4300 the above preallocation routines for simplicity. 4301 4302 Options Database Keys: 4303 . -mat_type aij - sets the matrix type to "aij" during a call to `MatSetFromOptions()` 4304 4305 Note: 4306 Subclasses include `MATAIJCUSPARSE`, `MATAIJPERM`, `MATAIJSELL`, `MATAIJMKL`, `MATAIJCRL`, and also automatically switches over to use inodes when 4307 enough exist. 4308 4309 Level: beginner 4310 4311 .seealso: `MatCreateAIJ()`, `MatCreateSeqAIJ()`, `MATSEQAIJ`, `MATMPIAIJ`, `MATSELL`, `MATSEQSELL`, `MATMPISELL` 4312 M*/ 4313 4314 /*MC 4315 MATAIJCRL - MATAIJCRL = "aijcrl" - A matrix type to be used for sparse matrices. 4316 4317 This matrix type is identical to `MATSEQAIJCRL` when constructed with a single process communicator, 4318 and `MATMPIAIJCRL` otherwise. As a result, for single process communicators, 4319 `MatSeqAIJSetPreallocation()` is supported, and similarly `MatMPIAIJSetPreallocation()` is supported 4320 for communicators controlling multiple processes. It is recommended that you call both of 4321 the above preallocation routines for simplicity. 4322 4323 Options Database Keys: 4324 . -mat_type aijcrl - sets the matrix type to "aijcrl" during a call to `MatSetFromOptions()` 4325 4326 Level: beginner 4327 4328 .seealso: `MatCreateMPIAIJCRL`, `MATSEQAIJCRL`, `MATMPIAIJCRL`, `MATSEQAIJCRL`, `MATMPIAIJCRL` 4329 M*/ 4330 4331 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCRL(Mat, MatType, MatReuse, Mat *); 4332 #if defined(PETSC_HAVE_ELEMENTAL) 4333 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_Elemental(Mat, MatType, MatReuse, Mat *); 4334 #endif 4335 #if defined(PETSC_HAVE_SCALAPACK) 4336 PETSC_INTERN PetscErrorCode MatConvert_AIJ_ScaLAPACK(Mat, MatType, MatReuse, Mat *); 4337 #endif 4338 #if defined(PETSC_HAVE_HYPRE) 4339 PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A, MatType, MatReuse, Mat *); 4340 #endif 4341 4342 PETSC_EXTERN PetscErrorCode MatConvert_SeqAIJ_SeqSELL(Mat, MatType, MatReuse, Mat *); 4343 PETSC_INTERN PetscErrorCode MatConvert_XAIJ_IS(Mat, MatType, MatReuse, Mat *); 4344 PETSC_INTERN PetscErrorCode MatProductSetFromOptions_IS_XAIJ(Mat); 4345 4346 /*@C 4347 MatSeqAIJGetArray - gives read/write access to the array where the data for a `MATSEQAIJ` matrix is stored 4348 4349 Not Collective 4350 4351 Input Parameter: 4352 . mat - a `MATSEQAIJ` matrix 4353 4354 Output Parameter: 4355 . array - pointer to the data 4356 4357 Level: intermediate 4358 4359 .seealso: `MatSeqAIJRestoreArray()`, `MatSeqAIJGetArrayF90()` 4360 @*/ 4361 PetscErrorCode MatSeqAIJGetArray(Mat A, PetscScalar **array) 4362 { 4363 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)A->data; 4364 4365 PetscFunctionBegin; 4366 if (aij->ops->getarray) { 4367 PetscCall((*aij->ops->getarray)(A, array)); 4368 } else { 4369 *array = aij->a; 4370 } 4371 PetscFunctionReturn(0); 4372 } 4373 4374 /*@C 4375 MatSeqAIJRestoreArray - returns access to the array where the data for a `MATSEQAIJ` matrix is stored obtained by `MatSeqAIJGetArray()` 4376 4377 Not Collective 4378 4379 Input Parameters: 4380 + mat - a `MATSEQAIJ` matrix 4381 - array - pointer to the data 4382 4383 Level: intermediate 4384 4385 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArrayF90()` 4386 @*/ 4387 PetscErrorCode MatSeqAIJRestoreArray(Mat A, PetscScalar **array) 4388 { 4389 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)A->data; 4390 4391 PetscFunctionBegin; 4392 if (aij->ops->restorearray) { 4393 PetscCall((*aij->ops->restorearray)(A, array)); 4394 } else { 4395 *array = NULL; 4396 } 4397 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 4398 PetscCall(PetscObjectStateIncrease((PetscObject)A)); 4399 PetscFunctionReturn(0); 4400 } 4401 4402 /*@C 4403 MatSeqAIJGetArrayRead - gives read-only access to the array where the data for a `MATSEQAIJ` matrix is stored 4404 4405 Not Collective 4406 4407 Input Parameter: 4408 . mat - a `MATSEQAIJ` matrix 4409 4410 Output Parameter: 4411 . array - pointer to the data 4412 4413 Level: intermediate 4414 4415 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArrayRead()` 4416 @*/ 4417 PetscErrorCode MatSeqAIJGetArrayRead(Mat A, const PetscScalar **array) 4418 { 4419 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)A->data; 4420 4421 PetscFunctionBegin; 4422 if (aij->ops->getarrayread) { 4423 PetscCall((*aij->ops->getarrayread)(A, array)); 4424 } else { 4425 *array = aij->a; 4426 } 4427 PetscFunctionReturn(0); 4428 } 4429 4430 /*@C 4431 MatSeqAIJRestoreArrayRead - restore the read-only access array obtained from `MatSeqAIJGetArrayRead()` 4432 4433 Not Collective 4434 4435 Input Parameter: 4436 . mat - a `MATSEQAIJ` matrix 4437 4438 Output Parameter: 4439 . array - pointer to the data 4440 4441 Level: intermediate 4442 4443 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJGetArrayRead()` 4444 @*/ 4445 PetscErrorCode MatSeqAIJRestoreArrayRead(Mat A, const PetscScalar **array) 4446 { 4447 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)A->data; 4448 4449 PetscFunctionBegin; 4450 if (aij->ops->restorearrayread) { 4451 PetscCall((*aij->ops->restorearrayread)(A, array)); 4452 } else { 4453 *array = NULL; 4454 } 4455 PetscFunctionReturn(0); 4456 } 4457 4458 /*@C 4459 MatSeqAIJGetArrayWrite - gives write-only access to the array where the data for a `MATSEQAIJ` matrix is stored 4460 4461 Not Collective 4462 4463 Input Parameter: 4464 . mat - a `MATSEQAIJ` matrix 4465 4466 Output Parameter: 4467 . array - pointer to the data 4468 4469 Level: intermediate 4470 4471 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArrayRead()` 4472 @*/ 4473 PetscErrorCode MatSeqAIJGetArrayWrite(Mat A, PetscScalar **array) 4474 { 4475 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)A->data; 4476 4477 PetscFunctionBegin; 4478 if (aij->ops->getarraywrite) { 4479 PetscCall((*aij->ops->getarraywrite)(A, array)); 4480 } else { 4481 *array = aij->a; 4482 } 4483 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 4484 PetscCall(PetscObjectStateIncrease((PetscObject)A)); 4485 PetscFunctionReturn(0); 4486 } 4487 4488 /*@C 4489 MatSeqAIJRestoreArrayWrite - restore the read-only access array obtained from MatSeqAIJGetArrayRead 4490 4491 Not Collective 4492 4493 Input Parameter: 4494 . mat - a MATSEQAIJ matrix 4495 4496 Output Parameter: 4497 . array - pointer to the data 4498 4499 Level: intermediate 4500 4501 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJGetArrayRead()` 4502 @*/ 4503 PetscErrorCode MatSeqAIJRestoreArrayWrite(Mat A, PetscScalar **array) 4504 { 4505 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)A->data; 4506 4507 PetscFunctionBegin; 4508 if (aij->ops->restorearraywrite) { 4509 PetscCall((*aij->ops->restorearraywrite)(A, array)); 4510 } else { 4511 *array = NULL; 4512 } 4513 PetscFunctionReturn(0); 4514 } 4515 4516 /*@C 4517 MatSeqAIJGetCSRAndMemType - Get the CSR arrays and the memory type of the `MATSEQAIJ` matrix 4518 4519 Not Collective 4520 4521 Input Parameter: 4522 . mat - a matrix of type `MATSEQAIJ` or its subclasses 4523 4524 Output Parameters: 4525 + i - row map array of the matrix 4526 . j - column index array of the matrix 4527 . a - data array of the matrix 4528 - memtype - memory type of the arrays 4529 4530 Notes: 4531 Any of the output parameters can be NULL, in which case the corresponding value is not returned. 4532 If mat is a device matrix, the arrays are on the device. Otherwise, they are on the host. 4533 4534 One can call this routine on a preallocated but not assembled matrix to just get the memory of the CSR underneath the matrix. 4535 If the matrix is assembled, the data array 'a' is guaranteed to have the latest values of the matrix. 4536 4537 Level: Developer 4538 4539 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJGetArrayRead()` 4540 @*/ 4541 PetscErrorCode MatSeqAIJGetCSRAndMemType(Mat mat, const PetscInt **i, const PetscInt **j, PetscScalar **a, PetscMemType *mtype) 4542 { 4543 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 4544 4545 PetscFunctionBegin; 4546 PetscCheck(mat->preallocated, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "matrix is not preallocated"); 4547 if (aij->ops->getcsrandmemtype) { 4548 PetscCall((*aij->ops->getcsrandmemtype)(mat, i, j, a, mtype)); 4549 } else { 4550 if (i) *i = aij->i; 4551 if (j) *j = aij->j; 4552 if (a) *a = aij->a; 4553 if (mtype) *mtype = PETSC_MEMTYPE_HOST; 4554 } 4555 PetscFunctionReturn(0); 4556 } 4557 4558 /*@C 4559 MatSeqAIJGetMaxRowNonzeros - returns the maximum number of nonzeros in any row 4560 4561 Not Collective 4562 4563 Input Parameter: 4564 . mat - a `MATSEQAIJ` matrix 4565 4566 Output Parameter: 4567 . nz - the maximum number of nonzeros in any row 4568 4569 Level: intermediate 4570 4571 .seealso: `MatSeqAIJRestoreArray()`, `MatSeqAIJGetArrayF90()` 4572 @*/ 4573 PetscErrorCode MatSeqAIJGetMaxRowNonzeros(Mat A, PetscInt *nz) 4574 { 4575 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)A->data; 4576 4577 PetscFunctionBegin; 4578 *nz = aij->rmax; 4579 PetscFunctionReturn(0); 4580 } 4581 4582 PetscErrorCode MatSetPreallocationCOO_SeqAIJ(Mat mat, PetscCount coo_n, PetscInt coo_i[], PetscInt coo_j[]) 4583 { 4584 MPI_Comm comm; 4585 PetscInt *i, *j; 4586 PetscInt M, N, row; 4587 PetscCount k, p, q, nneg, nnz, start, end; /* Index the coo array, so use PetscCount as their type */ 4588 PetscInt *Ai; /* Change to PetscCount once we use it for row pointers */ 4589 PetscInt *Aj; 4590 PetscScalar *Aa; 4591 Mat_SeqAIJ *seqaij = (Mat_SeqAIJ *)(mat->data); 4592 MatType rtype; 4593 PetscCount *perm, *jmap; 4594 4595 PetscFunctionBegin; 4596 PetscCall(MatResetPreallocationCOO_SeqAIJ(mat)); 4597 PetscCall(PetscObjectGetComm((PetscObject)mat, &comm)); 4598 PetscCall(MatGetSize(mat, &M, &N)); 4599 i = coo_i; 4600 j = coo_j; 4601 PetscCall(PetscMalloc1(coo_n, &perm)); 4602 for (k = 0; k < coo_n; k++) { /* Ignore entries with negative row or col indices */ 4603 if (j[k] < 0) i[k] = -1; 4604 perm[k] = k; 4605 } 4606 4607 /* Sort by row */ 4608 PetscCall(PetscSortIntWithIntCountArrayPair(coo_n, i, j, perm)); 4609 for (k = 0; k < coo_n; k++) { 4610 if (i[k] >= 0) break; 4611 } /* Advance k to the first row with a non-negative index */ 4612 nneg = k; 4613 PetscCall(PetscMalloc1(coo_n - nneg + 1, &jmap)); /* +1 to make a CSR-like data structure. jmap[i] originally is the number of repeats for i-th nonzero */ 4614 nnz = 0; /* Total number of unique nonzeros to be counted */ 4615 jmap++; /* Inc jmap by 1 for convinience */ 4616 4617 PetscCall(PetscCalloc1(M + 1, &Ai)); /* CSR of A */ 4618 PetscCall(PetscMalloc1(coo_n - nneg, &Aj)); /* We have at most coo_n-nneg unique nonzeros */ 4619 4620 /* In each row, sort by column, then unique column indices to get row length */ 4621 Ai++; /* Inc by 1 for convinience */ 4622 q = 0; /* q-th unique nonzero, with q starting from 0 */ 4623 while (k < coo_n) { 4624 row = i[k]; 4625 start = k; /* [start,end) indices for this row */ 4626 while (k < coo_n && i[k] == row) k++; 4627 end = k; 4628 PetscCall(PetscSortIntWithCountArray(end - start, j + start, perm + start)); 4629 /* Find number of unique col entries in this row */ 4630 Aj[q] = j[start]; /* Log the first nonzero in this row */ 4631 jmap[q] = 1; /* Number of repeats of this nozero entry */ 4632 Ai[row] = 1; 4633 nnz++; 4634 4635 for (p = start + 1; p < end; p++) { /* Scan remaining nonzero in this row */ 4636 if (j[p] != j[p - 1]) { /* Meet a new nonzero */ 4637 q++; 4638 jmap[q] = 1; 4639 Aj[q] = j[p]; 4640 Ai[row]++; 4641 nnz++; 4642 } else { 4643 jmap[q]++; 4644 } 4645 } 4646 q++; /* Move to next row and thus next unique nonzero */ 4647 } 4648 4649 Ai--; /* Back to the beginning of Ai[] */ 4650 for (k = 0; k < M; k++) Ai[k + 1] += Ai[k]; 4651 jmap--; /* Back to the beginning of jmap[] */ 4652 jmap[0] = 0; 4653 for (k = 0; k < nnz; k++) jmap[k + 1] += jmap[k]; 4654 if (nnz < coo_n - nneg) { /* Realloc with actual number of unique nonzeros */ 4655 PetscCount *jmap_new; 4656 PetscInt *Aj_new; 4657 4658 PetscCall(PetscMalloc1(nnz + 1, &jmap_new)); 4659 PetscCall(PetscArraycpy(jmap_new, jmap, nnz + 1)); 4660 PetscCall(PetscFree(jmap)); 4661 jmap = jmap_new; 4662 4663 PetscCall(PetscMalloc1(nnz, &Aj_new)); 4664 PetscCall(PetscArraycpy(Aj_new, Aj, nnz)); 4665 PetscCall(PetscFree(Aj)); 4666 Aj = Aj_new; 4667 } 4668 4669 if (nneg) { /* Discard heading entries with negative indices in perm[], as we'll access it from index 0 in MatSetValuesCOO */ 4670 PetscCount *perm_new; 4671 4672 PetscCall(PetscMalloc1(coo_n - nneg, &perm_new)); 4673 PetscCall(PetscArraycpy(perm_new, perm + nneg, coo_n - nneg)); 4674 PetscCall(PetscFree(perm)); 4675 perm = perm_new; 4676 } 4677 4678 PetscCall(MatGetRootType_Private(mat, &rtype)); 4679 PetscCall(PetscCalloc1(nnz, &Aa)); /* Zero the matrix */ 4680 PetscCall(MatSetSeqAIJWithArrays_private(PETSC_COMM_SELF, M, N, Ai, Aj, Aa, rtype, mat)); 4681 4682 seqaij->singlemalloc = PETSC_FALSE; /* Ai, Aj and Aa are not allocated in one big malloc */ 4683 seqaij->free_a = seqaij->free_ij = PETSC_TRUE; /* Let newmat own Ai, Aj and Aa */ 4684 /* Record COO fields */ 4685 seqaij->coo_n = coo_n; 4686 seqaij->Atot = coo_n - nneg; /* Annz is seqaij->nz, so no need to record that again */ 4687 seqaij->jmap = jmap; /* of length nnz+1 */ 4688 seqaij->perm = perm; 4689 PetscFunctionReturn(0); 4690 } 4691 4692 static PetscErrorCode MatSetValuesCOO_SeqAIJ(Mat A, const PetscScalar v[], InsertMode imode) 4693 { 4694 Mat_SeqAIJ *aseq = (Mat_SeqAIJ *)A->data; 4695 PetscCount i, j, Annz = aseq->nz; 4696 PetscCount *perm = aseq->perm, *jmap = aseq->jmap; 4697 PetscScalar *Aa; 4698 4699 PetscFunctionBegin; 4700 PetscCall(MatSeqAIJGetArray(A, &Aa)); 4701 for (i = 0; i < Annz; i++) { 4702 PetscScalar sum = 0.0; 4703 for (j = jmap[i]; j < jmap[i + 1]; j++) sum += v[perm[j]]; 4704 Aa[i] = (imode == INSERT_VALUES ? 0.0 : Aa[i]) + sum; 4705 } 4706 PetscCall(MatSeqAIJRestoreArray(A, &Aa)); 4707 PetscFunctionReturn(0); 4708 } 4709 4710 #if defined(PETSC_HAVE_CUDA) 4711 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCUSPARSE(Mat, MatType, MatReuse, Mat *); 4712 #endif 4713 #if defined(PETSC_HAVE_HIP) 4714 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJHIPSPARSE(Mat, MatType, MatReuse, Mat *); 4715 #endif 4716 #if defined(PETSC_HAVE_KOKKOS_KERNELS) 4717 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJKokkos(Mat, MatType, MatReuse, Mat *); 4718 #endif 4719 4720 PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJ(Mat B) 4721 { 4722 Mat_SeqAIJ *b; 4723 PetscMPIInt size; 4724 4725 PetscFunctionBegin; 4726 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B), &size)); 4727 PetscCheck(size <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Comm must be of size 1"); 4728 4729 PetscCall(PetscNew(&b)); 4730 4731 B->data = (void *)b; 4732 4733 PetscCall(PetscMemcpy(B->ops, &MatOps_Values, sizeof(struct _MatOps))); 4734 if (B->sortedfull) B->ops->setvalues = MatSetValues_SeqAIJ_SortedFull; 4735 4736 b->row = NULL; 4737 b->col = NULL; 4738 b->icol = NULL; 4739 b->reallocs = 0; 4740 b->ignorezeroentries = PETSC_FALSE; 4741 b->roworiented = PETSC_TRUE; 4742 b->nonew = 0; 4743 b->diag = NULL; 4744 b->solve_work = NULL; 4745 B->spptr = NULL; 4746 b->saved_values = NULL; 4747 b->idiag = NULL; 4748 b->mdiag = NULL; 4749 b->ssor_work = NULL; 4750 b->omega = 1.0; 4751 b->fshift = 0.0; 4752 b->idiagvalid = PETSC_FALSE; 4753 b->ibdiagvalid = PETSC_FALSE; 4754 b->keepnonzeropattern = PETSC_FALSE; 4755 4756 PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATSEQAIJ)); 4757 #if defined(PETSC_HAVE_MATLAB) 4758 PetscCall(PetscObjectComposeFunction((PetscObject)B, "PetscMatlabEnginePut_C", MatlabEnginePut_SeqAIJ)); 4759 PetscCall(PetscObjectComposeFunction((PetscObject)B, "PetscMatlabEngineGet_C", MatlabEngineGet_SeqAIJ)); 4760 #endif 4761 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqAIJSetColumnIndices_C", MatSeqAIJSetColumnIndices_SeqAIJ)); 4762 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatStoreValues_C", MatStoreValues_SeqAIJ)); 4763 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatRetrieveValues_C", MatRetrieveValues_SeqAIJ)); 4764 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqsbaij_C", MatConvert_SeqAIJ_SeqSBAIJ)); 4765 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqbaij_C", MatConvert_SeqAIJ_SeqBAIJ)); 4766 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqaijperm_C", MatConvert_SeqAIJ_SeqAIJPERM)); 4767 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqaijsell_C", MatConvert_SeqAIJ_SeqAIJSELL)); 4768 #if defined(PETSC_HAVE_MKL_SPARSE) 4769 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqaijmkl_C", MatConvert_SeqAIJ_SeqAIJMKL)); 4770 #endif 4771 #if defined(PETSC_HAVE_CUDA) 4772 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqaijcusparse_C", MatConvert_SeqAIJ_SeqAIJCUSPARSE)); 4773 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaijcusparse_seqaij_C", MatProductSetFromOptions_SeqAIJ)); 4774 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaij_seqaijcusparse_C", MatProductSetFromOptions_SeqAIJ)); 4775 #endif 4776 #if defined(PETSC_HAVE_HIP) 4777 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqaijhipsparse_C", MatConvert_SeqAIJ_SeqAIJHIPSPARSE)); 4778 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaijhipsparse_seqaij_C", MatProductSetFromOptions_SeqAIJ)); 4779 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaij_seqaijhipsparse_C", MatProductSetFromOptions_SeqAIJ)); 4780 #endif 4781 #if defined(PETSC_HAVE_KOKKOS_KERNELS) 4782 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqaijkokkos_C", MatConvert_SeqAIJ_SeqAIJKokkos)); 4783 #endif 4784 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqaijcrl_C", MatConvert_SeqAIJ_SeqAIJCRL)); 4785 #if defined(PETSC_HAVE_ELEMENTAL) 4786 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_elemental_C", MatConvert_SeqAIJ_Elemental)); 4787 #endif 4788 #if defined(PETSC_HAVE_SCALAPACK) 4789 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_scalapack_C", MatConvert_AIJ_ScaLAPACK)); 4790 #endif 4791 #if defined(PETSC_HAVE_HYPRE) 4792 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_hypre_C", MatConvert_AIJ_HYPRE)); 4793 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_transpose_seqaij_seqaij_C", MatProductSetFromOptions_Transpose_AIJ_AIJ)); 4794 #endif 4795 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqdense_C", MatConvert_SeqAIJ_SeqDense)); 4796 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqsell_C", MatConvert_SeqAIJ_SeqSELL)); 4797 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_is_C", MatConvert_XAIJ_IS)); 4798 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatIsTranspose_C", MatIsTranspose_SeqAIJ)); 4799 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatIsHermitianTranspose_C", MatIsTranspose_SeqAIJ)); 4800 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqAIJSetPreallocation_C", MatSeqAIJSetPreallocation_SeqAIJ)); 4801 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatResetPreallocation_C", MatResetPreallocation_SeqAIJ)); 4802 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqAIJSetPreallocationCSR_C", MatSeqAIJSetPreallocationCSR_SeqAIJ)); 4803 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatReorderForNonzeroDiagonal_C", MatReorderForNonzeroDiagonal_SeqAIJ)); 4804 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_is_seqaij_C", MatProductSetFromOptions_IS_XAIJ)); 4805 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdense_seqaij_C", MatProductSetFromOptions_SeqDense_SeqAIJ)); 4806 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaij_seqaij_C", MatProductSetFromOptions_SeqAIJ)); 4807 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqAIJKron_C", MatSeqAIJKron_SeqAIJ)); 4808 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSetPreallocationCOO_C", MatSetPreallocationCOO_SeqAIJ)); 4809 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSetValuesCOO_C", MatSetValuesCOO_SeqAIJ)); 4810 PetscCall(MatCreate_SeqAIJ_Inode(B)); 4811 PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATSEQAIJ)); 4812 PetscCall(MatSeqAIJSetTypeFromOptions(B)); /* this allows changing the matrix subtype to say MATSEQAIJPERM */ 4813 PetscFunctionReturn(0); 4814 } 4815 4816 /* 4817 Given a matrix generated with MatGetFactor() duplicates all the information in A into C 4818 */ 4819 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C, Mat A, MatDuplicateOption cpvalues, PetscBool mallocmatspace) 4820 { 4821 Mat_SeqAIJ *c = (Mat_SeqAIJ *)C->data, *a = (Mat_SeqAIJ *)A->data; 4822 PetscInt m = A->rmap->n, i; 4823 4824 PetscFunctionBegin; 4825 PetscCheck(A->assembled || cpvalues == MAT_DO_NOT_COPY_VALUES, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cannot duplicate unassembled matrix"); 4826 4827 C->factortype = A->factortype; 4828 c->row = NULL; 4829 c->col = NULL; 4830 c->icol = NULL; 4831 c->reallocs = 0; 4832 4833 C->assembled = A->assembled; 4834 C->preallocated = A->preallocated; 4835 4836 if (A->preallocated) { 4837 PetscCall(PetscLayoutReference(A->rmap, &C->rmap)); 4838 PetscCall(PetscLayoutReference(A->cmap, &C->cmap)); 4839 4840 PetscCall(PetscMalloc1(m, &c->imax)); 4841 PetscCall(PetscMemcpy(c->imax, a->imax, m * sizeof(PetscInt))); 4842 PetscCall(PetscMalloc1(m, &c->ilen)); 4843 PetscCall(PetscMemcpy(c->ilen, a->ilen, m * sizeof(PetscInt))); 4844 4845 /* allocate the matrix space */ 4846 if (mallocmatspace) { 4847 PetscCall(PetscMalloc3(a->i[m], &c->a, a->i[m], &c->j, m + 1, &c->i)); 4848 4849 c->singlemalloc = PETSC_TRUE; 4850 4851 PetscCall(PetscArraycpy(c->i, a->i, m + 1)); 4852 if (m > 0) { 4853 PetscCall(PetscArraycpy(c->j, a->j, a->i[m])); 4854 if (cpvalues == MAT_COPY_VALUES) { 4855 const PetscScalar *aa; 4856 4857 PetscCall(MatSeqAIJGetArrayRead(A, &aa)); 4858 PetscCall(PetscArraycpy(c->a, aa, a->i[m])); 4859 PetscCall(MatSeqAIJGetArrayRead(A, &aa)); 4860 } else { 4861 PetscCall(PetscArrayzero(c->a, a->i[m])); 4862 } 4863 } 4864 } 4865 4866 c->ignorezeroentries = a->ignorezeroentries; 4867 c->roworiented = a->roworiented; 4868 c->nonew = a->nonew; 4869 if (a->diag) { 4870 PetscCall(PetscMalloc1(m + 1, &c->diag)); 4871 PetscCall(PetscMemcpy(c->diag, a->diag, m * sizeof(PetscInt))); 4872 } else c->diag = NULL; 4873 4874 c->solve_work = NULL; 4875 c->saved_values = NULL; 4876 c->idiag = NULL; 4877 c->ssor_work = NULL; 4878 c->keepnonzeropattern = a->keepnonzeropattern; 4879 c->free_a = PETSC_TRUE; 4880 c->free_ij = PETSC_TRUE; 4881 4882 c->rmax = a->rmax; 4883 c->nz = a->nz; 4884 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 4885 4886 c->compressedrow.use = a->compressedrow.use; 4887 c->compressedrow.nrows = a->compressedrow.nrows; 4888 if (a->compressedrow.use) { 4889 i = a->compressedrow.nrows; 4890 PetscCall(PetscMalloc2(i + 1, &c->compressedrow.i, i, &c->compressedrow.rindex)); 4891 PetscCall(PetscArraycpy(c->compressedrow.i, a->compressedrow.i, i + 1)); 4892 PetscCall(PetscArraycpy(c->compressedrow.rindex, a->compressedrow.rindex, i)); 4893 } else { 4894 c->compressedrow.use = PETSC_FALSE; 4895 c->compressedrow.i = NULL; 4896 c->compressedrow.rindex = NULL; 4897 } 4898 c->nonzerorowcnt = a->nonzerorowcnt; 4899 C->nonzerostate = A->nonzerostate; 4900 4901 PetscCall(MatDuplicate_SeqAIJ_Inode(A, cpvalues, &C)); 4902 } 4903 PetscCall(PetscFunctionListDuplicate(((PetscObject)A)->qlist, &((PetscObject)C)->qlist)); 4904 PetscFunctionReturn(0); 4905 } 4906 4907 PetscErrorCode MatDuplicate_SeqAIJ(Mat A, MatDuplicateOption cpvalues, Mat *B) 4908 { 4909 PetscFunctionBegin; 4910 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), B)); 4911 PetscCall(MatSetSizes(*B, A->rmap->n, A->cmap->n, A->rmap->n, A->cmap->n)); 4912 if (!(A->rmap->n % A->rmap->bs) && !(A->cmap->n % A->cmap->bs)) PetscCall(MatSetBlockSizesFromMats(*B, A, A)); 4913 PetscCall(MatSetType(*B, ((PetscObject)A)->type_name)); 4914 PetscCall(MatDuplicateNoCreate_SeqAIJ(*B, A, cpvalues, PETSC_TRUE)); 4915 PetscFunctionReturn(0); 4916 } 4917 4918 PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer) 4919 { 4920 PetscBool isbinary, ishdf5; 4921 4922 PetscFunctionBegin; 4923 PetscValidHeaderSpecific(newMat, MAT_CLASSID, 1); 4924 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 4925 /* force binary viewer to load .info file if it has not yet done so */ 4926 PetscCall(PetscViewerSetUp(viewer)); 4927 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary)); 4928 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5)); 4929 if (isbinary) { 4930 PetscCall(MatLoad_SeqAIJ_Binary(newMat, viewer)); 4931 } else if (ishdf5) { 4932 #if defined(PETSC_HAVE_HDF5) 4933 PetscCall(MatLoad_AIJ_HDF5(newMat, viewer)); 4934 #else 4935 SETERRQ(PetscObjectComm((PetscObject)newMat), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5"); 4936 #endif 4937 } else { 4938 SETERRQ(PetscObjectComm((PetscObject)newMat), PETSC_ERR_SUP, "Viewer type %s not yet supported for reading %s matrices", ((PetscObject)viewer)->type_name, ((PetscObject)newMat)->type_name); 4939 } 4940 PetscFunctionReturn(0); 4941 } 4942 4943 PetscErrorCode MatLoad_SeqAIJ_Binary(Mat mat, PetscViewer viewer) 4944 { 4945 Mat_SeqAIJ *a = (Mat_SeqAIJ *)mat->data; 4946 PetscInt header[4], *rowlens, M, N, nz, sum, rows, cols, i; 4947 4948 PetscFunctionBegin; 4949 PetscCall(PetscViewerSetUp(viewer)); 4950 4951 /* read in matrix header */ 4952 PetscCall(PetscViewerBinaryRead(viewer, header, 4, NULL, PETSC_INT)); 4953 PetscCheck(header[0] == MAT_FILE_CLASSID, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Not a matrix object in file"); 4954 M = header[1]; 4955 N = header[2]; 4956 nz = header[3]; 4957 PetscCheck(M >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Matrix row size (%" PetscInt_FMT ") in file is negative", M); 4958 PetscCheck(N >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Matrix column size (%" PetscInt_FMT ") in file is negative", N); 4959 PetscCheck(nz >= 0, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Matrix stored in special format on disk, cannot load as SeqAIJ"); 4960 4961 /* set block sizes from the viewer's .info file */ 4962 PetscCall(MatLoad_Binary_BlockSizes(mat, viewer)); 4963 /* set local and global sizes if not set already */ 4964 if (mat->rmap->n < 0) mat->rmap->n = M; 4965 if (mat->cmap->n < 0) mat->cmap->n = N; 4966 if (mat->rmap->N < 0) mat->rmap->N = M; 4967 if (mat->cmap->N < 0) mat->cmap->N = N; 4968 PetscCall(PetscLayoutSetUp(mat->rmap)); 4969 PetscCall(PetscLayoutSetUp(mat->cmap)); 4970 4971 /* check if the matrix sizes are correct */ 4972 PetscCall(MatGetSize(mat, &rows, &cols)); 4973 PetscCheck(M == rows && N == cols, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Matrix in file of different sizes (%" PetscInt_FMT ", %" PetscInt_FMT ") than the input matrix (%" PetscInt_FMT ", %" PetscInt_FMT ")", M, N, rows, cols); 4974 4975 /* read in row lengths */ 4976 PetscCall(PetscMalloc1(M, &rowlens)); 4977 PetscCall(PetscViewerBinaryRead(viewer, rowlens, M, NULL, PETSC_INT)); 4978 /* check if sum(rowlens) is same as nz */ 4979 sum = 0; 4980 for (i = 0; i < M; i++) sum += rowlens[i]; 4981 PetscCheck(sum == nz, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Inconsistent matrix data in file: nonzeros = %" PetscInt_FMT ", sum-row-lengths = %" PetscInt_FMT, nz, sum); 4982 /* preallocate and check sizes */ 4983 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(mat, 0, rowlens)); 4984 PetscCall(MatGetSize(mat, &rows, &cols)); 4985 PetscCheck(M == rows && N == cols, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Matrix in file of different length (%" PetscInt_FMT ", %" PetscInt_FMT ") than the input matrix (%" PetscInt_FMT ", %" PetscInt_FMT ")", M, N, rows, cols); 4986 /* store row lengths */ 4987 PetscCall(PetscArraycpy(a->ilen, rowlens, M)); 4988 PetscCall(PetscFree(rowlens)); 4989 4990 /* fill in "i" row pointers */ 4991 a->i[0] = 0; 4992 for (i = 0; i < M; i++) a->i[i + 1] = a->i[i] + a->ilen[i]; 4993 /* read in "j" column indices */ 4994 PetscCall(PetscViewerBinaryRead(viewer, a->j, nz, NULL, PETSC_INT)); 4995 /* read in "a" nonzero values */ 4996 PetscCall(PetscViewerBinaryRead(viewer, a->a, nz, NULL, PETSC_SCALAR)); 4997 4998 PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY)); 4999 PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY)); 5000 PetscFunctionReturn(0); 5001 } 5002 5003 PetscErrorCode MatEqual_SeqAIJ(Mat A, Mat B, PetscBool *flg) 5004 { 5005 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data; 5006 const PetscScalar *aa, *ba; 5007 #if defined(PETSC_USE_COMPLEX) 5008 PetscInt k; 5009 #endif 5010 5011 PetscFunctionBegin; 5012 /* If the matrix dimensions are not equal,or no of nonzeros */ 5013 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) || (a->nz != b->nz)) { 5014 *flg = PETSC_FALSE; 5015 PetscFunctionReturn(0); 5016 } 5017 5018 /* if the a->i are the same */ 5019 PetscCall(PetscArraycmp(a->i, b->i, A->rmap->n + 1, flg)); 5020 if (!*flg) PetscFunctionReturn(0); 5021 5022 /* if a->j are the same */ 5023 PetscCall(PetscArraycmp(a->j, b->j, a->nz, flg)); 5024 if (!*flg) PetscFunctionReturn(0); 5025 5026 PetscCall(MatSeqAIJGetArrayRead(A, &aa)); 5027 PetscCall(MatSeqAIJGetArrayRead(B, &ba)); 5028 /* if a->a are the same */ 5029 #if defined(PETSC_USE_COMPLEX) 5030 for (k = 0; k < a->nz; k++) { 5031 if (PetscRealPart(aa[k]) != PetscRealPart(ba[k]) || PetscImaginaryPart(aa[k]) != PetscImaginaryPart(ba[k])) { 5032 *flg = PETSC_FALSE; 5033 PetscFunctionReturn(0); 5034 } 5035 } 5036 #else 5037 PetscCall(PetscArraycmp(aa, ba, a->nz, flg)); 5038 #endif 5039 PetscCall(MatSeqAIJRestoreArrayRead(A, &aa)); 5040 PetscCall(MatSeqAIJRestoreArrayRead(B, &ba)); 5041 PetscFunctionReturn(0); 5042 } 5043 5044 /*@ 5045 MatCreateSeqAIJWithArrays - Creates an sequential `MATSEQAIJ` matrix using matrix elements (in CSR format) 5046 provided by the user. 5047 5048 Collective 5049 5050 Input Parameters: 5051 + comm - must be an MPI communicator of size 1 5052 . m - number of rows 5053 . n - number of columns 5054 . i - row indices; that is i[0] = 0, i[row] = i[row-1] + number of elements in that row of the matrix 5055 . j - column indices 5056 - a - matrix values 5057 5058 Output Parameter: 5059 . mat - the matrix 5060 5061 Level: intermediate 5062 5063 Notes: 5064 The i, j, and a arrays are not copied by this routine, the user must free these arrays 5065 once the matrix is destroyed and not before 5066 5067 You cannot set new nonzero locations into this matrix, that will generate an error. 5068 5069 The i and j indices are 0 based 5070 5071 The format which is used for the sparse matrix input, is equivalent to a 5072 row-major ordering.. i.e for the following matrix, the input data expected is 5073 as shown 5074 5075 $ 1 0 0 5076 $ 2 0 3 5077 $ 4 5 6 5078 $ 5079 $ i = {0,1,3,6} [size = nrow+1 = 3+1] 5080 $ j = {0,0,2,0,1,2} [size = 6]; values must be sorted for each row 5081 $ v = {1,2,3,4,5,6} [size = 6] 5082 5083 .seealso: `MatCreate()`, `MatCreateAIJ()`, `MatCreateSeqAIJ()`, `MatCreateMPIAIJWithArrays()`, `MatMPIAIJSetPreallocationCSR()` 5084 @*/ 5085 PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt i[], PetscInt j[], PetscScalar a[], Mat *mat) 5086 { 5087 PetscInt ii; 5088 Mat_SeqAIJ *aij; 5089 PetscInt jj; 5090 5091 PetscFunctionBegin; 5092 PetscCheck(m <= 0 || i[0] == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "i (row indices) must start with 0"); 5093 PetscCall(MatCreate(comm, mat)); 5094 PetscCall(MatSetSizes(*mat, m, n, m, n)); 5095 /* PetscCall(MatSetBlockSizes(*mat,,)); */ 5096 PetscCall(MatSetType(*mat, MATSEQAIJ)); 5097 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(*mat, MAT_SKIP_ALLOCATION, NULL)); 5098 aij = (Mat_SeqAIJ *)(*mat)->data; 5099 PetscCall(PetscMalloc1(m, &aij->imax)); 5100 PetscCall(PetscMalloc1(m, &aij->ilen)); 5101 5102 aij->i = i; 5103 aij->j = j; 5104 aij->a = a; 5105 aij->singlemalloc = PETSC_FALSE; 5106 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 5107 aij->free_a = PETSC_FALSE; 5108 aij->free_ij = PETSC_FALSE; 5109 5110 for (ii = 0, aij->nonzerorowcnt = 0, aij->rmax = 0; ii < m; ii++) { 5111 aij->ilen[ii] = aij->imax[ii] = i[ii + 1] - i[ii]; 5112 if (PetscDefined(USE_DEBUG)) { 5113 PetscCheck(i[ii + 1] - i[ii] >= 0, 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]); 5114 for (jj = i[ii] + 1; jj < i[ii + 1]; jj++) { 5115 PetscCheck(j[jj] >= j[jj - 1], PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column entry number %" PetscInt_FMT " (actual column %" PetscInt_FMT ") in row %" PetscInt_FMT " is not sorted", jj - i[ii], j[jj], ii); 5116 PetscCheck(j[jj] != j[jj - 1], PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column entry number %" PetscInt_FMT " (actual column %" PetscInt_FMT ") in row %" PetscInt_FMT " is identical to previous entry", jj - i[ii], j[jj], ii); 5117 } 5118 } 5119 } 5120 if (PetscDefined(USE_DEBUG)) { 5121 for (ii = 0; ii < aij->i[m]; ii++) { 5122 PetscCheck(j[ii] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative column index at location = %" PetscInt_FMT " index = %" PetscInt_FMT, ii, j[ii]); 5123 PetscCheck(j[ii] <= n - 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column index to large at location = %" PetscInt_FMT " index = %" PetscInt_FMT, ii, j[ii]); 5124 } 5125 } 5126 5127 PetscCall(MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY)); 5128 PetscCall(MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY)); 5129 PetscFunctionReturn(0); 5130 } 5131 5132 /*@ 5133 MatCreateSeqAIJFromTriple - Creates an sequential `MATSEQAIJ` matrix using matrix elements (in COO format) 5134 provided by the user. 5135 5136 Collective 5137 5138 Input Parameters: 5139 + comm - must be an MPI communicator of size 1 5140 . m - number of rows 5141 . n - number of columns 5142 . i - row indices 5143 . j - column indices 5144 . a - matrix values 5145 . nz - number of nonzeros 5146 - idx - if the i and j indices start with 1 use `PETSC_TRUE` otherwise use `PETSC_FALSE` 5147 5148 Output Parameter: 5149 . mat - the matrix 5150 5151 Level: intermediate 5152 5153 Example: 5154 For the following matrix, the input data expected is as shown (using 0 based indexing) 5155 .vb 5156 1 0 0 5157 2 0 3 5158 4 5 6 5159 5160 i = {0,1,1,2,2,2} 5161 j = {0,0,2,0,1,2} 5162 v = {1,2,3,4,5,6} 5163 .ve 5164 Notes: 5165 Instead of using this function, users should also consider `MatSetPreallocationCOO()` and `MatSetValuesCOO()`, which allow repeated or remote entries, 5166 and are particularly useful in iterative applications. 5167 5168 .seealso: `MatCreate()`, `MatCreateAIJ()`, `MatCreateSeqAIJ()`, `MatCreateSeqAIJWithArrays()`, `MatMPIAIJSetPreallocationCSR()`, `MatSetValuesCOO()`, `MatSetPreallocationCOO()` 5169 @*/ 5170 PetscErrorCode MatCreateSeqAIJFromTriple(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt i[], PetscInt j[], PetscScalar a[], Mat *mat, PetscInt nz, PetscBool idx) 5171 { 5172 PetscInt ii, *nnz, one = 1, row, col; 5173 5174 PetscFunctionBegin; 5175 PetscCall(PetscCalloc1(m, &nnz)); 5176 for (ii = 0; ii < nz; ii++) nnz[i[ii] - !!idx] += 1; 5177 PetscCall(MatCreate(comm, mat)); 5178 PetscCall(MatSetSizes(*mat, m, n, m, n)); 5179 PetscCall(MatSetType(*mat, MATSEQAIJ)); 5180 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(*mat, 0, nnz)); 5181 for (ii = 0; ii < nz; ii++) { 5182 if (idx) { 5183 row = i[ii] - 1; 5184 col = j[ii] - 1; 5185 } else { 5186 row = i[ii]; 5187 col = j[ii]; 5188 } 5189 PetscCall(MatSetValues(*mat, one, &row, one, &col, &a[ii], ADD_VALUES)); 5190 } 5191 PetscCall(MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY)); 5192 PetscCall(MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY)); 5193 PetscCall(PetscFree(nnz)); 5194 PetscFunctionReturn(0); 5195 } 5196 5197 PetscErrorCode MatSeqAIJInvalidateDiagonal(Mat A) 5198 { 5199 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 5200 5201 PetscFunctionBegin; 5202 a->idiagvalid = PETSC_FALSE; 5203 a->ibdiagvalid = PETSC_FALSE; 5204 5205 PetscCall(MatSeqAIJInvalidateDiagonal_Inode(A)); 5206 PetscFunctionReturn(0); 5207 } 5208 5209 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqAIJ(MPI_Comm comm, Mat inmat, PetscInt n, MatReuse scall, Mat *outmat) 5210 { 5211 PetscFunctionBegin; 5212 PetscCall(MatCreateMPIMatConcatenateSeqMat_MPIAIJ(comm, inmat, n, scall, outmat)); 5213 PetscFunctionReturn(0); 5214 } 5215 5216 /* 5217 Permute A into C's *local* index space using rowemb,colemb. 5218 The embedding are supposed to be injections and the above implies that the range of rowemb is a subset 5219 of [0,m), colemb is in [0,n). 5220 If pattern == DIFFERENT_NONZERO_PATTERN, C is preallocated according to A. 5221 */ 5222 PetscErrorCode MatSetSeqMat_SeqAIJ(Mat C, IS rowemb, IS colemb, MatStructure pattern, Mat B) 5223 { 5224 /* If making this function public, change the error returned in this function away from _PLIB. */ 5225 Mat_SeqAIJ *Baij; 5226 PetscBool seqaij; 5227 PetscInt m, n, *nz, i, j, count; 5228 PetscScalar v; 5229 const PetscInt *rowindices, *colindices; 5230 5231 PetscFunctionBegin; 5232 if (!B) PetscFunctionReturn(0); 5233 /* Check to make sure the target matrix (and embeddings) are compatible with C and each other. */ 5234 PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATSEQAIJ, &seqaij)); 5235 PetscCheck(seqaij, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Input matrix is of wrong type"); 5236 if (rowemb) { 5237 PetscCall(ISGetLocalSize(rowemb, &m)); 5238 PetscCheck(m == B->rmap->n, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Row IS of size %" PetscInt_FMT " is incompatible with matrix row size %" PetscInt_FMT, m, B->rmap->n); 5239 } else { 5240 PetscCheck(C->rmap->n == B->rmap->n, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Input matrix is row-incompatible with the target matrix"); 5241 } 5242 if (colemb) { 5243 PetscCall(ISGetLocalSize(colemb, &n)); 5244 PetscCheck(n == B->cmap->n, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Diag col IS of size %" PetscInt_FMT " is incompatible with input matrix col size %" PetscInt_FMT, n, B->cmap->n); 5245 } else { 5246 PetscCheck(C->cmap->n == B->cmap->n, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Input matrix is col-incompatible with the target matrix"); 5247 } 5248 5249 Baij = (Mat_SeqAIJ *)(B->data); 5250 if (pattern == DIFFERENT_NONZERO_PATTERN) { 5251 PetscCall(PetscMalloc1(B->rmap->n, &nz)); 5252 for (i = 0; i < B->rmap->n; i++) nz[i] = Baij->i[i + 1] - Baij->i[i]; 5253 PetscCall(MatSeqAIJSetPreallocation(C, 0, nz)); 5254 PetscCall(PetscFree(nz)); 5255 } 5256 if (pattern == SUBSET_NONZERO_PATTERN) PetscCall(MatZeroEntries(C)); 5257 count = 0; 5258 rowindices = NULL; 5259 colindices = NULL; 5260 if (rowemb) PetscCall(ISGetIndices(rowemb, &rowindices)); 5261 if (colemb) PetscCall(ISGetIndices(colemb, &colindices)); 5262 for (i = 0; i < B->rmap->n; i++) { 5263 PetscInt row; 5264 row = i; 5265 if (rowindices) row = rowindices[i]; 5266 for (j = Baij->i[i]; j < Baij->i[i + 1]; j++) { 5267 PetscInt col; 5268 col = Baij->j[count]; 5269 if (colindices) col = colindices[col]; 5270 v = Baij->a[count]; 5271 PetscCall(MatSetValues(C, 1, &row, 1, &col, &v, INSERT_VALUES)); 5272 ++count; 5273 } 5274 } 5275 /* FIXME: set C's nonzerostate correctly. */ 5276 /* Assembly for C is necessary. */ 5277 C->preallocated = PETSC_TRUE; 5278 C->assembled = PETSC_TRUE; 5279 C->was_assembled = PETSC_FALSE; 5280 PetscFunctionReturn(0); 5281 } 5282 5283 PetscErrorCode MatEliminateZeros_SeqAIJ(Mat A) 5284 { 5285 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 5286 MatScalar *aa = a->a; 5287 PetscInt m = A->rmap->n, fshift = 0, fshift_prev = 0, i, k; 5288 PetscInt *ailen = a->ilen, *imax = a->imax, *ai = a->i, *aj = a->j, rmax = 0; 5289 5290 PetscFunctionBegin; 5291 PetscCheck(A->assembled, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cannot eliminate zeros for unassembled matrix"); 5292 if (m) rmax = ailen[0]; /* determine row with most nonzeros */ 5293 for (i = 1; i <= m; i++) { 5294 /* move each nonzero entry back by the amount of zero slots (fshift) before it*/ 5295 for (k = ai[i - 1]; k < ai[i]; k++) { 5296 if (aa[k] == 0 && aj[k] != i - 1) fshift++; 5297 else { 5298 if (aa[k] == 0 && aj[k] == i - 1) PetscCall(PetscInfo(A, "Keep the diagonal zero at row %" PetscInt_FMT "\n", i - 1)); 5299 aa[k - fshift] = aa[k]; 5300 aj[k - fshift] = aj[k]; 5301 } 5302 } 5303 ai[i - 1] -= fshift_prev; // safe to update ai[i-1] now since it will not be used in the next iteration 5304 fshift_prev = fshift; 5305 /* reset ilen and imax for each row */ 5306 ailen[i - 1] = imax[i - 1] = ai[i] - fshift - ai[i - 1]; 5307 a->nonzerorowcnt += ((ai[i] - fshift - ai[i - 1]) > 0); 5308 rmax = PetscMax(rmax, ailen[i - 1]); 5309 } 5310 if (m) { 5311 ai[m] -= fshift; 5312 a->nz = ai[m]; 5313 } 5314 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)); 5315 A->nonzerostate -= fshift; 5316 A->info.nz_unneeded += (PetscReal)fshift; 5317 a->rmax = rmax; 5318 if (a->inode.use && a->inode.checked) PetscCall(MatSeqAIJCheckInode(A)); 5319 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 5320 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 5321 PetscFunctionReturn(0); 5322 } 5323 5324 PetscFunctionList MatSeqAIJList = NULL; 5325 5326 /*@C 5327 MatSeqAIJSetType - Converts a `MATSEQAIJ` matrix to a subtype 5328 5329 Collective 5330 5331 Input Parameters: 5332 + mat - the matrix object 5333 - matype - matrix type 5334 5335 Options Database Key: 5336 . -mat_seqai_type <method> - for example seqaijcrl 5337 5338 Level: intermediate 5339 5340 .seealso: `PCSetType()`, `VecSetType()`, `MatCreate()`, `MatType`, `Mat` 5341 @*/ 5342 PetscErrorCode MatSeqAIJSetType(Mat mat, MatType matype) 5343 { 5344 PetscBool sametype; 5345 PetscErrorCode (*r)(Mat, MatType, MatReuse, Mat *); 5346 5347 PetscFunctionBegin; 5348 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 5349 PetscCall(PetscObjectTypeCompare((PetscObject)mat, matype, &sametype)); 5350 if (sametype) PetscFunctionReturn(0); 5351 5352 PetscCall(PetscFunctionListFind(MatSeqAIJList, matype, &r)); 5353 PetscCheck(r, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown Mat type given: %s", matype); 5354 PetscCall((*r)(mat, matype, MAT_INPLACE_MATRIX, &mat)); 5355 PetscFunctionReturn(0); 5356 } 5357 5358 /*@C 5359 MatSeqAIJRegister - - Adds a new sub-matrix type for sequential `MATSEQAIJ` matrices 5360 5361 Not Collective 5362 5363 Input Parameters: 5364 + name - name of a new user-defined matrix type, for example `MATSEQAIJCRL` 5365 - function - routine to convert to subtype 5366 5367 Notes: 5368 `MatSeqAIJRegister()` may be called multiple times to add several user-defined solvers. 5369 5370 Then, your matrix can be chosen with the procedural interface at runtime via the option 5371 $ -mat_seqaij_type my_mat 5372 5373 Level: advanced 5374 5375 .seealso: `MatSeqAIJRegisterAll()` 5376 @*/ 5377 PetscErrorCode MatSeqAIJRegister(const char sname[], PetscErrorCode (*function)(Mat, MatType, MatReuse, Mat *)) 5378 { 5379 PetscFunctionBegin; 5380 PetscCall(MatInitializePackage()); 5381 PetscCall(PetscFunctionListAdd(&MatSeqAIJList, sname, function)); 5382 PetscFunctionReturn(0); 5383 } 5384 5385 PetscBool MatSeqAIJRegisterAllCalled = PETSC_FALSE; 5386 5387 /*@C 5388 MatSeqAIJRegisterAll - Registers all of the matrix subtypes of `MATSSEQAIJ` 5389 5390 Not Collective 5391 5392 Level: advanced 5393 5394 .seealso: `MatRegisterAll()`, `MatSeqAIJRegister()` 5395 @*/ 5396 PetscErrorCode MatSeqAIJRegisterAll(void) 5397 { 5398 PetscFunctionBegin; 5399 if (MatSeqAIJRegisterAllCalled) PetscFunctionReturn(0); 5400 MatSeqAIJRegisterAllCalled = PETSC_TRUE; 5401 5402 PetscCall(MatSeqAIJRegister(MATSEQAIJCRL, MatConvert_SeqAIJ_SeqAIJCRL)); 5403 PetscCall(MatSeqAIJRegister(MATSEQAIJPERM, MatConvert_SeqAIJ_SeqAIJPERM)); 5404 PetscCall(MatSeqAIJRegister(MATSEQAIJSELL, MatConvert_SeqAIJ_SeqAIJSELL)); 5405 #if defined(PETSC_HAVE_MKL_SPARSE) 5406 PetscCall(MatSeqAIJRegister(MATSEQAIJMKL, MatConvert_SeqAIJ_SeqAIJMKL)); 5407 #endif 5408 #if defined(PETSC_HAVE_CUDA) 5409 PetscCall(MatSeqAIJRegister(MATSEQAIJCUSPARSE, MatConvert_SeqAIJ_SeqAIJCUSPARSE)); 5410 #endif 5411 #if defined(PETSC_HAVE_HIP) 5412 PetscCall(MatSeqAIJRegister(MATSEQAIJHIPSPARSE, MatConvert_SeqAIJ_SeqAIJHIPSPARSE)); 5413 #endif 5414 #if defined(PETSC_HAVE_KOKKOS_KERNELS) 5415 PetscCall(MatSeqAIJRegister(MATSEQAIJKOKKOS, MatConvert_SeqAIJ_SeqAIJKokkos)); 5416 #endif 5417 #if defined(PETSC_HAVE_VIENNACL) && defined(PETSC_HAVE_VIENNACL_NO_CUDA) 5418 PetscCall(MatSeqAIJRegister(MATMPIAIJVIENNACL, MatConvert_SeqAIJ_SeqAIJViennaCL)); 5419 #endif 5420 PetscFunctionReturn(0); 5421 } 5422 5423 /* 5424 Special version for direct calls from Fortran 5425 */ 5426 #include <petsc/private/fortranimpl.h> 5427 #if defined(PETSC_HAVE_FORTRAN_CAPS) 5428 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 5429 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 5430 #define matsetvaluesseqaij_ matsetvaluesseqaij 5431 #endif 5432 5433 /* Change these macros so can be used in void function */ 5434 5435 /* Change these macros so can be used in void function */ 5436 /* Identical to PetscCallVoid, except it assigns to *_ierr */ 5437 #undef PetscCall 5438 #define PetscCall(...) \ 5439 do { \ 5440 PetscErrorCode ierr_msv_mpiaij = __VA_ARGS__; \ 5441 if (PetscUnlikely(ierr_msv_mpiaij)) { \ 5442 *_ierr = PetscError(PETSC_COMM_SELF, __LINE__, PETSC_FUNCTION_NAME, __FILE__, ierr_msv_mpiaij, PETSC_ERROR_REPEAT, " "); \ 5443 return; \ 5444 } \ 5445 } while (0) 5446 5447 #undef SETERRQ 5448 #define SETERRQ(comm, ierr, ...) \ 5449 do { \ 5450 *_ierr = PetscError(comm, __LINE__, PETSC_FUNCTION_NAME, __FILE__, ierr, PETSC_ERROR_INITIAL, __VA_ARGS__); \ 5451 return; \ 5452 } while (0) 5453 5454 PETSC_EXTERN void matsetvaluesseqaij_(Mat *AA, PetscInt *mm, const PetscInt im[], PetscInt *nn, const PetscInt in[], const PetscScalar v[], InsertMode *isis, PetscErrorCode *_ierr) 5455 { 5456 Mat A = *AA; 5457 PetscInt m = *mm, n = *nn; 5458 InsertMode is = *isis; 5459 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 5460 PetscInt *rp, k, low, high, t, ii, row, nrow, i, col, l, rmax, N; 5461 PetscInt *imax, *ai, *ailen; 5462 PetscInt *aj, nonew = a->nonew, lastcol = -1; 5463 MatScalar *ap, value, *aa; 5464 PetscBool ignorezeroentries = a->ignorezeroentries; 5465 PetscBool roworiented = a->roworiented; 5466 5467 PetscFunctionBegin; 5468 MatCheckPreallocated(A, 1); 5469 imax = a->imax; 5470 ai = a->i; 5471 ailen = a->ilen; 5472 aj = a->j; 5473 aa = a->a; 5474 5475 for (k = 0; k < m; k++) { /* loop over added rows */ 5476 row = im[k]; 5477 if (row < 0) continue; 5478 PetscCheck(row < A->rmap->n, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Row too large"); 5479 rp = aj + ai[row]; 5480 ap = aa + ai[row]; 5481 rmax = imax[row]; 5482 nrow = ailen[row]; 5483 low = 0; 5484 high = nrow; 5485 for (l = 0; l < n; l++) { /* loop over added columns */ 5486 if (in[l] < 0) continue; 5487 PetscCheck(in[l] < A->cmap->n, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Column too large"); 5488 col = in[l]; 5489 if (roworiented) value = v[l + k * n]; 5490 else value = v[k + l * m]; 5491 5492 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 5493 5494 if (col <= lastcol) low = 0; 5495 else high = nrow; 5496 lastcol = col; 5497 while (high - low > 5) { 5498 t = (low + high) / 2; 5499 if (rp[t] > col) high = t; 5500 else low = t; 5501 } 5502 for (i = low; i < high; i++) { 5503 if (rp[i] > col) break; 5504 if (rp[i] == col) { 5505 if (is == ADD_VALUES) ap[i] += value; 5506 else ap[i] = value; 5507 goto noinsert; 5508 } 5509 } 5510 if (value == 0.0 && ignorezeroentries) goto noinsert; 5511 if (nonew == 1) goto noinsert; 5512 PetscCheck(nonew != -1, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Inserting a new nonzero in the matrix"); 5513 MatSeqXAIJReallocateAIJ(A, A->rmap->n, 1, nrow, row, col, rmax, aa, ai, aj, rp, ap, imax, nonew, MatScalar); 5514 N = nrow++ - 1; 5515 a->nz++; 5516 high++; 5517 /* shift up all the later entries in this row */ 5518 for (ii = N; ii >= i; ii--) { 5519 rp[ii + 1] = rp[ii]; 5520 ap[ii + 1] = ap[ii]; 5521 } 5522 rp[i] = col; 5523 ap[i] = value; 5524 A->nonzerostate++; 5525 noinsert:; 5526 low = i + 1; 5527 } 5528 ailen[row] = nrow; 5529 } 5530 PetscFunctionReturnVoid(); 5531 } 5532 /* Undefining these here since they were redefined from their original definition above! No 5533 * other PETSc functions should be defined past this point, as it is impossible to recover the 5534 * original definitions */ 5535 #undef PetscCall 5536 #undef SETERRQ 5537