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, ll, m, n, *nidx, isz, val; 2817 const PetscInt *idx; 2818 PetscInt start, end, *ai, *aj, bs = (A->rmap->bs > 0 && A->rmap->bs == A->cmap->bs) ? A->rmap->bs : 1; 2819 PetscBT table; 2820 2821 PetscFunctionBegin; 2822 m = A->rmap->n / bs; 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 if (bs > 1) { 2841 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 2842 for (j = 0; j < n; ++j) { 2843 if (!PetscBTLookupSet(table, idx[j] / bs)) nidx[isz++] = idx[j] / bs; 2844 } 2845 PetscCall(ISRestoreIndices(is[i], &idx)); 2846 PetscCall(ISDestroy(&is[i])); 2847 2848 k = 0; 2849 for (j = 0; j < ov; j++) { /* for each overlap */ 2850 n = isz; 2851 for (; k < n; k++) { /* do only those rows in nidx[k], which are not done yet */ 2852 for (ll = 0; ll < bs; ll++) { 2853 row = bs * nidx[k] + ll; 2854 start = ai[row]; 2855 end = ai[row + 1]; 2856 for (l = start; l < end; l++) { 2857 val = aj[l] / bs; 2858 if (!PetscBTLookupSet(table, val)) nidx[isz++] = val; 2859 } 2860 } 2861 } 2862 } 2863 PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, isz, nidx, PETSC_COPY_VALUES, (is + i))); 2864 } else { 2865 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 2866 for (j = 0; j < n; ++j) { 2867 if (!PetscBTLookupSet(table, idx[j])) nidx[isz++] = idx[j]; 2868 } 2869 PetscCall(ISRestoreIndices(is[i], &idx)); 2870 PetscCall(ISDestroy(&is[i])); 2871 2872 k = 0; 2873 for (j = 0; j < ov; j++) { /* for each overlap */ 2874 n = isz; 2875 for (; k < n; k++) { /* do only those rows in nidx[k], which are not done yet */ 2876 row = nidx[k]; 2877 start = ai[row]; 2878 end = ai[row + 1]; 2879 for (l = start; l < end; l++) { 2880 val = aj[l]; 2881 if (!PetscBTLookupSet(table, val)) nidx[isz++] = val; 2882 } 2883 } 2884 } 2885 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, isz, nidx, PETSC_COPY_VALUES, (is + i))); 2886 } 2887 } 2888 PetscCall(PetscBTDestroy(&table)); 2889 PetscCall(PetscFree(nidx)); 2890 PetscFunctionReturn(0); 2891 } 2892 2893 /* -------------------------------------------------------------- */ 2894 PetscErrorCode MatPermute_SeqAIJ(Mat A, IS rowp, IS colp, Mat *B) 2895 { 2896 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 2897 PetscInt i, nz = 0, m = A->rmap->n, n = A->cmap->n; 2898 const PetscInt *row, *col; 2899 PetscInt *cnew, j, *lens; 2900 IS icolp, irowp; 2901 PetscInt *cwork = NULL; 2902 PetscScalar *vwork = NULL; 2903 2904 PetscFunctionBegin; 2905 PetscCall(ISInvertPermutation(rowp, PETSC_DECIDE, &irowp)); 2906 PetscCall(ISGetIndices(irowp, &row)); 2907 PetscCall(ISInvertPermutation(colp, PETSC_DECIDE, &icolp)); 2908 PetscCall(ISGetIndices(icolp, &col)); 2909 2910 /* determine lengths of permuted rows */ 2911 PetscCall(PetscMalloc1(m + 1, &lens)); 2912 for (i = 0; i < m; i++) lens[row[i]] = a->i[i + 1] - a->i[i]; 2913 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), B)); 2914 PetscCall(MatSetSizes(*B, m, n, m, n)); 2915 PetscCall(MatSetBlockSizesFromMats(*B, A, A)); 2916 PetscCall(MatSetType(*B, ((PetscObject)A)->type_name)); 2917 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(*B, 0, lens)); 2918 PetscCall(PetscFree(lens)); 2919 2920 PetscCall(PetscMalloc1(n, &cnew)); 2921 for (i = 0; i < m; i++) { 2922 PetscCall(MatGetRow_SeqAIJ(A, i, &nz, &cwork, &vwork)); 2923 for (j = 0; j < nz; j++) cnew[j] = col[cwork[j]]; 2924 PetscCall(MatSetValues_SeqAIJ(*B, 1, &row[i], nz, cnew, vwork, INSERT_VALUES)); 2925 PetscCall(MatRestoreRow_SeqAIJ(A, i, &nz, &cwork, &vwork)); 2926 } 2927 PetscCall(PetscFree(cnew)); 2928 2929 (*B)->assembled = PETSC_FALSE; 2930 2931 #if defined(PETSC_HAVE_DEVICE) 2932 PetscCall(MatBindToCPU(*B, A->boundtocpu)); 2933 #endif 2934 PetscCall(MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY)); 2935 PetscCall(MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY)); 2936 PetscCall(ISRestoreIndices(irowp, &row)); 2937 PetscCall(ISRestoreIndices(icolp, &col)); 2938 PetscCall(ISDestroy(&irowp)); 2939 PetscCall(ISDestroy(&icolp)); 2940 if (rowp == colp) PetscCall(MatPropagateSymmetryOptions(A, *B)); 2941 PetscFunctionReturn(0); 2942 } 2943 2944 PetscErrorCode MatCopy_SeqAIJ(Mat A, Mat B, MatStructure str) 2945 { 2946 PetscFunctionBegin; 2947 /* If the two matrices have the same copy implementation, use fast copy. */ 2948 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 2949 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 2950 Mat_SeqAIJ *b = (Mat_SeqAIJ *)B->data; 2951 const PetscScalar *aa; 2952 2953 PetscCall(MatSeqAIJGetArrayRead(A, &aa)); 2954 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]); 2955 PetscCall(PetscArraycpy(b->a, aa, a->i[A->rmap->n])); 2956 PetscCall(PetscObjectStateIncrease((PetscObject)B)); 2957 PetscCall(MatSeqAIJRestoreArrayRead(A, &aa)); 2958 } else { 2959 PetscCall(MatCopy_Basic(A, B, str)); 2960 } 2961 PetscFunctionReturn(0); 2962 } 2963 2964 PetscErrorCode MatSetUp_SeqAIJ(Mat A) 2965 { 2966 PetscFunctionBegin; 2967 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(A, PETSC_DEFAULT, NULL)); 2968 PetscFunctionReturn(0); 2969 } 2970 2971 PETSC_INTERN PetscErrorCode MatSeqAIJGetArray_SeqAIJ(Mat A, PetscScalar *array[]) 2972 { 2973 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 2974 2975 PetscFunctionBegin; 2976 *array = a->a; 2977 PetscFunctionReturn(0); 2978 } 2979 2980 PETSC_INTERN PetscErrorCode MatSeqAIJRestoreArray_SeqAIJ(Mat A, PetscScalar *array[]) 2981 { 2982 PetscFunctionBegin; 2983 *array = NULL; 2984 PetscFunctionReturn(0); 2985 } 2986 2987 /* 2988 Computes the number of nonzeros per row needed for preallocation when X and Y 2989 have different nonzero structure. 2990 */ 2991 PetscErrorCode MatAXPYGetPreallocation_SeqX_private(PetscInt m, const PetscInt *xi, const PetscInt *xj, const PetscInt *yi, const PetscInt *yj, PetscInt *nnz) 2992 { 2993 PetscInt i, j, k, nzx, nzy; 2994 2995 PetscFunctionBegin; 2996 /* Set the number of nonzeros in the new matrix */ 2997 for (i = 0; i < m; i++) { 2998 const PetscInt *xjj = xj + xi[i], *yjj = yj + yi[i]; 2999 nzx = xi[i + 1] - xi[i]; 3000 nzy = yi[i + 1] - yi[i]; 3001 nnz[i] = 0; 3002 for (j = 0, k = 0; j < nzx; j++) { /* Point in X */ 3003 for (; k < nzy && yjj[k] < xjj[j]; k++) nnz[i]++; /* Catch up to X */ 3004 if (k < nzy && yjj[k] == xjj[j]) k++; /* Skip duplicate */ 3005 nnz[i]++; 3006 } 3007 for (; k < nzy; k++) nnz[i]++; 3008 } 3009 PetscFunctionReturn(0); 3010 } 3011 3012 PetscErrorCode MatAXPYGetPreallocation_SeqAIJ(Mat Y, Mat X, PetscInt *nnz) 3013 { 3014 PetscInt m = Y->rmap->N; 3015 Mat_SeqAIJ *x = (Mat_SeqAIJ *)X->data; 3016 Mat_SeqAIJ *y = (Mat_SeqAIJ *)Y->data; 3017 3018 PetscFunctionBegin; 3019 /* Set the number of nonzeros in the new matrix */ 3020 PetscCall(MatAXPYGetPreallocation_SeqX_private(m, x->i, x->j, y->i, y->j, nnz)); 3021 PetscFunctionReturn(0); 3022 } 3023 3024 PetscErrorCode MatAXPY_SeqAIJ(Mat Y, PetscScalar a, Mat X, MatStructure str) 3025 { 3026 Mat_SeqAIJ *x = (Mat_SeqAIJ *)X->data, *y = (Mat_SeqAIJ *)Y->data; 3027 3028 PetscFunctionBegin; 3029 if (str == UNKNOWN_NONZERO_PATTERN || (PetscDefined(USE_DEBUG) && str == SAME_NONZERO_PATTERN)) { 3030 PetscBool e = x->nz == y->nz ? PETSC_TRUE : PETSC_FALSE; 3031 if (e) { 3032 PetscCall(PetscArraycmp(x->i, y->i, Y->rmap->n + 1, &e)); 3033 if (e) { 3034 PetscCall(PetscArraycmp(x->j, y->j, y->nz, &e)); 3035 if (e) str = SAME_NONZERO_PATTERN; 3036 } 3037 } 3038 if (!e) PetscCheck(str != SAME_NONZERO_PATTERN, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "MatStructure is not SAME_NONZERO_PATTERN"); 3039 } 3040 if (str == SAME_NONZERO_PATTERN) { 3041 const PetscScalar *xa; 3042 PetscScalar *ya, alpha = a; 3043 PetscBLASInt one = 1, bnz; 3044 3045 PetscCall(PetscBLASIntCast(x->nz, &bnz)); 3046 PetscCall(MatSeqAIJGetArray(Y, &ya)); 3047 PetscCall(MatSeqAIJGetArrayRead(X, &xa)); 3048 PetscCallBLAS("BLASaxpy", BLASaxpy_(&bnz, &alpha, xa, &one, ya, &one)); 3049 PetscCall(MatSeqAIJRestoreArrayRead(X, &xa)); 3050 PetscCall(MatSeqAIJRestoreArray(Y, &ya)); 3051 PetscCall(PetscLogFlops(2.0 * bnz)); 3052 PetscCall(MatSeqAIJInvalidateDiagonal(Y)); 3053 PetscCall(PetscObjectStateIncrease((PetscObject)Y)); 3054 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 3055 PetscCall(MatAXPY_Basic(Y, a, X, str)); 3056 } else { 3057 Mat B; 3058 PetscInt *nnz; 3059 PetscCall(PetscMalloc1(Y->rmap->N, &nnz)); 3060 PetscCall(MatCreate(PetscObjectComm((PetscObject)Y), &B)); 3061 PetscCall(PetscObjectSetName((PetscObject)B, ((PetscObject)Y)->name)); 3062 PetscCall(MatSetLayouts(B, Y->rmap, Y->cmap)); 3063 PetscCall(MatSetType(B, ((PetscObject)Y)->type_name)); 3064 PetscCall(MatAXPYGetPreallocation_SeqAIJ(Y, X, nnz)); 3065 PetscCall(MatSeqAIJSetPreallocation(B, 0, nnz)); 3066 PetscCall(MatAXPY_BasicWithPreallocation(B, Y, a, X, str)); 3067 PetscCall(MatHeaderMerge(Y, &B)); 3068 PetscCall(MatSeqAIJCheckInode(Y)); 3069 PetscCall(PetscFree(nnz)); 3070 } 3071 PetscFunctionReturn(0); 3072 } 3073 3074 PETSC_INTERN PetscErrorCode MatConjugate_SeqAIJ(Mat mat) 3075 { 3076 #if defined(PETSC_USE_COMPLEX) 3077 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 3078 PetscInt i, nz; 3079 PetscScalar *a; 3080 3081 PetscFunctionBegin; 3082 nz = aij->nz; 3083 PetscCall(MatSeqAIJGetArray(mat, &a)); 3084 for (i = 0; i < nz; i++) a[i] = PetscConj(a[i]); 3085 PetscCall(MatSeqAIJRestoreArray(mat, &a)); 3086 #else 3087 PetscFunctionBegin; 3088 #endif 3089 PetscFunctionReturn(0); 3090 } 3091 3092 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A, Vec v, PetscInt idx[]) 3093 { 3094 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 3095 PetscInt i, j, m = A->rmap->n, *ai, *aj, ncols, n; 3096 PetscReal atmp; 3097 PetscScalar *x; 3098 const MatScalar *aa, *av; 3099 3100 PetscFunctionBegin; 3101 PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix"); 3102 PetscCall(MatSeqAIJGetArrayRead(A, &av)); 3103 aa = av; 3104 ai = a->i; 3105 aj = a->j; 3106 3107 PetscCall(VecSet(v, 0.0)); 3108 PetscCall(VecGetArrayWrite(v, &x)); 3109 PetscCall(VecGetLocalSize(v, &n)); 3110 PetscCheck(n == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector"); 3111 for (i = 0; i < m; i++) { 3112 ncols = ai[1] - ai[0]; 3113 ai++; 3114 for (j = 0; j < ncols; j++) { 3115 atmp = PetscAbsScalar(*aa); 3116 if (PetscAbsScalar(x[i]) < atmp) { 3117 x[i] = atmp; 3118 if (idx) idx[i] = *aj; 3119 } 3120 aa++; 3121 aj++; 3122 } 3123 } 3124 PetscCall(VecRestoreArrayWrite(v, &x)); 3125 PetscCall(MatSeqAIJRestoreArrayRead(A, &av)); 3126 PetscFunctionReturn(0); 3127 } 3128 3129 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A, Vec v, PetscInt idx[]) 3130 { 3131 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 3132 PetscInt i, j, m = A->rmap->n, *ai, *aj, ncols, n; 3133 PetscScalar *x; 3134 const MatScalar *aa, *av; 3135 3136 PetscFunctionBegin; 3137 PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix"); 3138 PetscCall(MatSeqAIJGetArrayRead(A, &av)); 3139 aa = av; 3140 ai = a->i; 3141 aj = a->j; 3142 3143 PetscCall(VecSet(v, 0.0)); 3144 PetscCall(VecGetArrayWrite(v, &x)); 3145 PetscCall(VecGetLocalSize(v, &n)); 3146 PetscCheck(n == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector"); 3147 for (i = 0; i < m; i++) { 3148 ncols = ai[1] - ai[0]; 3149 ai++; 3150 if (ncols == A->cmap->n) { /* row is dense */ 3151 x[i] = *aa; 3152 if (idx) idx[i] = 0; 3153 } else { /* row is sparse so already KNOW maximum is 0.0 or higher */ 3154 x[i] = 0.0; 3155 if (idx) { 3156 for (j = 0; j < ncols; j++) { /* find first implicit 0.0 in the row */ 3157 if (aj[j] > j) { 3158 idx[i] = j; 3159 break; 3160 } 3161 } 3162 /* in case first implicit 0.0 in the row occurs at ncols-th column */ 3163 if (j == ncols && j < A->cmap->n) idx[i] = j; 3164 } 3165 } 3166 for (j = 0; j < ncols; j++) { 3167 if (PetscRealPart(x[i]) < PetscRealPart(*aa)) { 3168 x[i] = *aa; 3169 if (idx) idx[i] = *aj; 3170 } 3171 aa++; 3172 aj++; 3173 } 3174 } 3175 PetscCall(VecRestoreArrayWrite(v, &x)); 3176 PetscCall(MatSeqAIJRestoreArrayRead(A, &av)); 3177 PetscFunctionReturn(0); 3178 } 3179 3180 PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A, Vec v, PetscInt idx[]) 3181 { 3182 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 3183 PetscInt i, j, m = A->rmap->n, *ai, *aj, ncols, n; 3184 PetscScalar *x; 3185 const MatScalar *aa, *av; 3186 3187 PetscFunctionBegin; 3188 PetscCall(MatSeqAIJGetArrayRead(A, &av)); 3189 aa = av; 3190 ai = a->i; 3191 aj = a->j; 3192 3193 PetscCall(VecSet(v, 0.0)); 3194 PetscCall(VecGetArrayWrite(v, &x)); 3195 PetscCall(VecGetLocalSize(v, &n)); 3196 PetscCheck(n == m, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector, %" PetscInt_FMT " vs. %" PetscInt_FMT " rows", m, n); 3197 for (i = 0; i < m; i++) { 3198 ncols = ai[1] - ai[0]; 3199 ai++; 3200 if (ncols == A->cmap->n) { /* row is dense */ 3201 x[i] = *aa; 3202 if (idx) idx[i] = 0; 3203 } else { /* row is sparse so already KNOW minimum is 0.0 or higher */ 3204 x[i] = 0.0; 3205 if (idx) { /* find first implicit 0.0 in the row */ 3206 for (j = 0; j < ncols; j++) { 3207 if (aj[j] > j) { 3208 idx[i] = j; 3209 break; 3210 } 3211 } 3212 /* in case first implicit 0.0 in the row occurs at ncols-th column */ 3213 if (j == ncols && j < A->cmap->n) idx[i] = j; 3214 } 3215 } 3216 for (j = 0; j < ncols; j++) { 3217 if (PetscAbsScalar(x[i]) > PetscAbsScalar(*aa)) { 3218 x[i] = *aa; 3219 if (idx) idx[i] = *aj; 3220 } 3221 aa++; 3222 aj++; 3223 } 3224 } 3225 PetscCall(VecRestoreArrayWrite(v, &x)); 3226 PetscCall(MatSeqAIJRestoreArrayRead(A, &av)); 3227 PetscFunctionReturn(0); 3228 } 3229 3230 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A, Vec v, PetscInt idx[]) 3231 { 3232 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 3233 PetscInt i, j, m = A->rmap->n, ncols, n; 3234 const PetscInt *ai, *aj; 3235 PetscScalar *x; 3236 const MatScalar *aa, *av; 3237 3238 PetscFunctionBegin; 3239 PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix"); 3240 PetscCall(MatSeqAIJGetArrayRead(A, &av)); 3241 aa = av; 3242 ai = a->i; 3243 aj = a->j; 3244 3245 PetscCall(VecSet(v, 0.0)); 3246 PetscCall(VecGetArrayWrite(v, &x)); 3247 PetscCall(VecGetLocalSize(v, &n)); 3248 PetscCheck(n == m, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector"); 3249 for (i = 0; i < m; i++) { 3250 ncols = ai[1] - ai[0]; 3251 ai++; 3252 if (ncols == A->cmap->n) { /* row is dense */ 3253 x[i] = *aa; 3254 if (idx) idx[i] = 0; 3255 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 3256 x[i] = 0.0; 3257 if (idx) { /* find first implicit 0.0 in the row */ 3258 for (j = 0; j < ncols; j++) { 3259 if (aj[j] > j) { 3260 idx[i] = j; 3261 break; 3262 } 3263 } 3264 /* in case first implicit 0.0 in the row occurs at ncols-th column */ 3265 if (j == ncols && j < A->cmap->n) idx[i] = j; 3266 } 3267 } 3268 for (j = 0; j < ncols; j++) { 3269 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) { 3270 x[i] = *aa; 3271 if (idx) idx[i] = *aj; 3272 } 3273 aa++; 3274 aj++; 3275 } 3276 } 3277 PetscCall(VecRestoreArrayWrite(v, &x)); 3278 PetscCall(MatSeqAIJRestoreArrayRead(A, &av)); 3279 PetscFunctionReturn(0); 3280 } 3281 3282 PetscErrorCode MatInvertBlockDiagonal_SeqAIJ(Mat A, const PetscScalar **values) 3283 { 3284 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 3285 PetscInt i, bs = PetscAbs(A->rmap->bs), mbs = A->rmap->n / bs, ipvt[5], bs2 = bs * bs, *v_pivots, ij[7], *IJ, j; 3286 MatScalar *diag, work[25], *v_work; 3287 const PetscReal shift = 0.0; 3288 PetscBool allowzeropivot, zeropivotdetected = PETSC_FALSE; 3289 3290 PetscFunctionBegin; 3291 allowzeropivot = PetscNot(A->erroriffailure); 3292 if (a->ibdiagvalid) { 3293 if (values) *values = a->ibdiag; 3294 PetscFunctionReturn(0); 3295 } 3296 PetscCall(MatMarkDiagonal_SeqAIJ(A)); 3297 if (!a->ibdiag) { PetscCall(PetscMalloc1(bs2 * mbs, &a->ibdiag)); } 3298 diag = a->ibdiag; 3299 if (values) *values = a->ibdiag; 3300 /* factor and invert each block */ 3301 switch (bs) { 3302 case 1: 3303 for (i = 0; i < mbs; i++) { 3304 PetscCall(MatGetValues(A, 1, &i, 1, &i, diag + i)); 3305 if (PetscAbsScalar(diag[i] + shift) < PETSC_MACHINE_EPSILON) { 3306 if (allowzeropivot) { 3307 A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3308 A->factorerror_zeropivot_value = PetscAbsScalar(diag[i]); 3309 A->factorerror_zeropivot_row = i; 3310 PetscCall(PetscInfo(A, "Zero pivot, row %" PetscInt_FMT " pivot %g tolerance %g\n", i, (double)PetscAbsScalar(diag[i]), (double)PETSC_MACHINE_EPSILON)); 3311 } 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); 3312 } 3313 diag[i] = (PetscScalar)1.0 / (diag[i] + shift); 3314 } 3315 break; 3316 case 2: 3317 for (i = 0; i < mbs; i++) { 3318 ij[0] = 2 * i; 3319 ij[1] = 2 * i + 1; 3320 PetscCall(MatGetValues(A, 2, ij, 2, ij, diag)); 3321 PetscCall(PetscKernel_A_gets_inverse_A_2(diag, shift, allowzeropivot, &zeropivotdetected)); 3322 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3323 PetscCall(PetscKernel_A_gets_transpose_A_2(diag)); 3324 diag += 4; 3325 } 3326 break; 3327 case 3: 3328 for (i = 0; i < mbs; i++) { 3329 ij[0] = 3 * i; 3330 ij[1] = 3 * i + 1; 3331 ij[2] = 3 * i + 2; 3332 PetscCall(MatGetValues(A, 3, ij, 3, ij, diag)); 3333 PetscCall(PetscKernel_A_gets_inverse_A_3(diag, shift, allowzeropivot, &zeropivotdetected)); 3334 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3335 PetscCall(PetscKernel_A_gets_transpose_A_3(diag)); 3336 diag += 9; 3337 } 3338 break; 3339 case 4: 3340 for (i = 0; i < mbs; i++) { 3341 ij[0] = 4 * i; 3342 ij[1] = 4 * i + 1; 3343 ij[2] = 4 * i + 2; 3344 ij[3] = 4 * i + 3; 3345 PetscCall(MatGetValues(A, 4, ij, 4, ij, diag)); 3346 PetscCall(PetscKernel_A_gets_inverse_A_4(diag, shift, allowzeropivot, &zeropivotdetected)); 3347 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3348 PetscCall(PetscKernel_A_gets_transpose_A_4(diag)); 3349 diag += 16; 3350 } 3351 break; 3352 case 5: 3353 for (i = 0; i < mbs; i++) { 3354 ij[0] = 5 * i; 3355 ij[1] = 5 * i + 1; 3356 ij[2] = 5 * i + 2; 3357 ij[3] = 5 * i + 3; 3358 ij[4] = 5 * i + 4; 3359 PetscCall(MatGetValues(A, 5, ij, 5, ij, diag)); 3360 PetscCall(PetscKernel_A_gets_inverse_A_5(diag, ipvt, work, shift, allowzeropivot, &zeropivotdetected)); 3361 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3362 PetscCall(PetscKernel_A_gets_transpose_A_5(diag)); 3363 diag += 25; 3364 } 3365 break; 3366 case 6: 3367 for (i = 0; i < mbs; i++) { 3368 ij[0] = 6 * i; 3369 ij[1] = 6 * i + 1; 3370 ij[2] = 6 * i + 2; 3371 ij[3] = 6 * i + 3; 3372 ij[4] = 6 * i + 4; 3373 ij[5] = 6 * i + 5; 3374 PetscCall(MatGetValues(A, 6, ij, 6, ij, diag)); 3375 PetscCall(PetscKernel_A_gets_inverse_A_6(diag, shift, allowzeropivot, &zeropivotdetected)); 3376 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3377 PetscCall(PetscKernel_A_gets_transpose_A_6(diag)); 3378 diag += 36; 3379 } 3380 break; 3381 case 7: 3382 for (i = 0; i < mbs; i++) { 3383 ij[0] = 7 * i; 3384 ij[1] = 7 * i + 1; 3385 ij[2] = 7 * i + 2; 3386 ij[3] = 7 * i + 3; 3387 ij[4] = 7 * i + 4; 3388 ij[5] = 7 * i + 5; 3389 ij[5] = 7 * i + 6; 3390 PetscCall(MatGetValues(A, 7, ij, 7, ij, diag)); 3391 PetscCall(PetscKernel_A_gets_inverse_A_7(diag, shift, allowzeropivot, &zeropivotdetected)); 3392 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3393 PetscCall(PetscKernel_A_gets_transpose_A_7(diag)); 3394 diag += 49; 3395 } 3396 break; 3397 default: 3398 PetscCall(PetscMalloc3(bs, &v_work, bs, &v_pivots, bs, &IJ)); 3399 for (i = 0; i < mbs; i++) { 3400 for (j = 0; j < bs; j++) IJ[j] = bs * i + j; 3401 PetscCall(MatGetValues(A, bs, IJ, bs, IJ, diag)); 3402 PetscCall(PetscKernel_A_gets_inverse_A(bs, diag, v_pivots, v_work, allowzeropivot, &zeropivotdetected)); 3403 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3404 PetscCall(PetscKernel_A_gets_transpose_A_N(diag, bs)); 3405 diag += bs2; 3406 } 3407 PetscCall(PetscFree3(v_work, v_pivots, IJ)); 3408 } 3409 a->ibdiagvalid = PETSC_TRUE; 3410 PetscFunctionReturn(0); 3411 } 3412 3413 static PetscErrorCode MatSetRandom_SeqAIJ(Mat x, PetscRandom rctx) 3414 { 3415 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)x->data; 3416 PetscScalar a, *aa; 3417 PetscInt m, n, i, j, col; 3418 3419 PetscFunctionBegin; 3420 if (!x->assembled) { 3421 PetscCall(MatGetSize(x, &m, &n)); 3422 for (i = 0; i < m; i++) { 3423 for (j = 0; j < aij->imax[i]; j++) { 3424 PetscCall(PetscRandomGetValue(rctx, &a)); 3425 col = (PetscInt)(n * PetscRealPart(a)); 3426 PetscCall(MatSetValues(x, 1, &i, 1, &col, &a, ADD_VALUES)); 3427 } 3428 } 3429 } else { 3430 PetscCall(MatSeqAIJGetArrayWrite(x, &aa)); 3431 for (i = 0; i < aij->nz; i++) PetscCall(PetscRandomGetValue(rctx, aa + i)); 3432 PetscCall(MatSeqAIJRestoreArrayWrite(x, &aa)); 3433 } 3434 PetscCall(MatAssemblyBegin(x, MAT_FINAL_ASSEMBLY)); 3435 PetscCall(MatAssemblyEnd(x, MAT_FINAL_ASSEMBLY)); 3436 PetscFunctionReturn(0); 3437 } 3438 3439 /* Like MatSetRandom_SeqAIJ, but do not set values on columns in range of [low, high) */ 3440 PetscErrorCode MatSetRandomSkipColumnRange_SeqAIJ_Private(Mat x, PetscInt low, PetscInt high, PetscRandom rctx) 3441 { 3442 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)x->data; 3443 PetscScalar a; 3444 PetscInt m, n, i, j, col, nskip; 3445 3446 PetscFunctionBegin; 3447 nskip = high - low; 3448 PetscCall(MatGetSize(x, &m, &n)); 3449 n -= nskip; /* shrink number of columns where nonzeros can be set */ 3450 for (i = 0; i < m; i++) { 3451 for (j = 0; j < aij->imax[i]; j++) { 3452 PetscCall(PetscRandomGetValue(rctx, &a)); 3453 col = (PetscInt)(n * PetscRealPart(a)); 3454 if (col >= low) col += nskip; /* shift col rightward to skip the hole */ 3455 PetscCall(MatSetValues(x, 1, &i, 1, &col, &a, ADD_VALUES)); 3456 } 3457 } 3458 PetscCall(MatAssemblyBegin(x, MAT_FINAL_ASSEMBLY)); 3459 PetscCall(MatAssemblyEnd(x, MAT_FINAL_ASSEMBLY)); 3460 PetscFunctionReturn(0); 3461 } 3462 3463 /* -------------------------------------------------------------------*/ 3464 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 3465 MatGetRow_SeqAIJ, 3466 MatRestoreRow_SeqAIJ, 3467 MatMult_SeqAIJ, 3468 /* 4*/ MatMultAdd_SeqAIJ, 3469 MatMultTranspose_SeqAIJ, 3470 MatMultTransposeAdd_SeqAIJ, 3471 NULL, 3472 NULL, 3473 NULL, 3474 /* 10*/ NULL, 3475 MatLUFactor_SeqAIJ, 3476 NULL, 3477 MatSOR_SeqAIJ, 3478 MatTranspose_SeqAIJ, 3479 /*1 5*/ MatGetInfo_SeqAIJ, 3480 MatEqual_SeqAIJ, 3481 MatGetDiagonal_SeqAIJ, 3482 MatDiagonalScale_SeqAIJ, 3483 MatNorm_SeqAIJ, 3484 /* 20*/ NULL, 3485 MatAssemblyEnd_SeqAIJ, 3486 MatSetOption_SeqAIJ, 3487 MatZeroEntries_SeqAIJ, 3488 /* 24*/ MatZeroRows_SeqAIJ, 3489 NULL, 3490 NULL, 3491 NULL, 3492 NULL, 3493 /* 29*/ MatSetUp_SeqAIJ, 3494 NULL, 3495 NULL, 3496 NULL, 3497 NULL, 3498 /* 34*/ MatDuplicate_SeqAIJ, 3499 NULL, 3500 NULL, 3501 MatILUFactor_SeqAIJ, 3502 NULL, 3503 /* 39*/ MatAXPY_SeqAIJ, 3504 MatCreateSubMatrices_SeqAIJ, 3505 MatIncreaseOverlap_SeqAIJ, 3506 MatGetValues_SeqAIJ, 3507 MatCopy_SeqAIJ, 3508 /* 44*/ MatGetRowMax_SeqAIJ, 3509 MatScale_SeqAIJ, 3510 MatShift_SeqAIJ, 3511 MatDiagonalSet_SeqAIJ, 3512 MatZeroRowsColumns_SeqAIJ, 3513 /* 49*/ MatSetRandom_SeqAIJ, 3514 MatGetRowIJ_SeqAIJ, 3515 MatRestoreRowIJ_SeqAIJ, 3516 MatGetColumnIJ_SeqAIJ, 3517 MatRestoreColumnIJ_SeqAIJ, 3518 /* 54*/ MatFDColoringCreate_SeqXAIJ, 3519 NULL, 3520 NULL, 3521 MatPermute_SeqAIJ, 3522 NULL, 3523 /* 59*/ NULL, 3524 MatDestroy_SeqAIJ, 3525 MatView_SeqAIJ, 3526 NULL, 3527 NULL, 3528 /* 64*/ NULL, 3529 MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ, 3530 NULL, 3531 NULL, 3532 NULL, 3533 /* 69*/ MatGetRowMaxAbs_SeqAIJ, 3534 MatGetRowMinAbs_SeqAIJ, 3535 NULL, 3536 NULL, 3537 NULL, 3538 /* 74*/ NULL, 3539 MatFDColoringApply_AIJ, 3540 NULL, 3541 NULL, 3542 NULL, 3543 /* 79*/ MatFindZeroDiagonals_SeqAIJ, 3544 NULL, 3545 NULL, 3546 NULL, 3547 MatLoad_SeqAIJ, 3548 /* 84*/ MatIsSymmetric_SeqAIJ, 3549 MatIsHermitian_SeqAIJ, 3550 NULL, 3551 NULL, 3552 NULL, 3553 /* 89*/ NULL, 3554 NULL, 3555 MatMatMultNumeric_SeqAIJ_SeqAIJ, 3556 NULL, 3557 NULL, 3558 /* 94*/ MatPtAPNumeric_SeqAIJ_SeqAIJ_SparseAxpy, 3559 NULL, 3560 NULL, 3561 MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ, 3562 NULL, 3563 /* 99*/ MatProductSetFromOptions_SeqAIJ, 3564 NULL, 3565 NULL, 3566 MatConjugate_SeqAIJ, 3567 NULL, 3568 /*104*/ MatSetValuesRow_SeqAIJ, 3569 MatRealPart_SeqAIJ, 3570 MatImaginaryPart_SeqAIJ, 3571 NULL, 3572 NULL, 3573 /*109*/ MatMatSolve_SeqAIJ, 3574 NULL, 3575 MatGetRowMin_SeqAIJ, 3576 NULL, 3577 MatMissingDiagonal_SeqAIJ, 3578 /*114*/ NULL, 3579 NULL, 3580 NULL, 3581 NULL, 3582 NULL, 3583 /*119*/ NULL, 3584 NULL, 3585 NULL, 3586 NULL, 3587 MatGetMultiProcBlock_SeqAIJ, 3588 /*124*/ MatFindNonzeroRows_SeqAIJ, 3589 MatGetColumnReductions_SeqAIJ, 3590 MatInvertBlockDiagonal_SeqAIJ, 3591 MatInvertVariableBlockDiagonal_SeqAIJ, 3592 NULL, 3593 /*129*/ NULL, 3594 NULL, 3595 NULL, 3596 MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ, 3597 MatTransposeColoringCreate_SeqAIJ, 3598 /*134*/ MatTransColoringApplySpToDen_SeqAIJ, 3599 MatTransColoringApplyDenToSp_SeqAIJ, 3600 NULL, 3601 NULL, 3602 MatRARtNumeric_SeqAIJ_SeqAIJ, 3603 /*139*/ NULL, 3604 NULL, 3605 NULL, 3606 MatFDColoringSetUp_SeqXAIJ, 3607 MatFindOffBlockDiagonalEntries_SeqAIJ, 3608 MatCreateMPIMatConcatenateSeqMat_SeqAIJ, 3609 /*145*/ MatDestroySubMatrices_SeqAIJ, 3610 NULL, 3611 NULL, 3612 MatCreateGraph_Simple_AIJ, 3613 NULL, 3614 /*150*/ MatTransposeSymbolic_SeqAIJ, 3615 MatEliminateZeros_SeqAIJ}; 3616 3617 PetscErrorCode MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat, PetscInt *indices) 3618 { 3619 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 3620 PetscInt i, nz, n; 3621 3622 PetscFunctionBegin; 3623 nz = aij->maxnz; 3624 n = mat->rmap->n; 3625 for (i = 0; i < nz; i++) aij->j[i] = indices[i]; 3626 aij->nz = nz; 3627 for (i = 0; i < n; i++) aij->ilen[i] = aij->imax[i]; 3628 PetscFunctionReturn(0); 3629 } 3630 3631 /* 3632 * Given a sparse matrix with global column indices, compact it by using a local column space. 3633 * The result matrix helps saving memory in other algorithms, such as MatPtAPSymbolic_MPIAIJ_MPIAIJ_scalable() 3634 */ 3635 PetscErrorCode MatSeqAIJCompactOutExtraColumns_SeqAIJ(Mat mat, ISLocalToGlobalMapping *mapping) 3636 { 3637 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 3638 PetscHMapI gid1_lid1; 3639 PetscHashIter tpos; 3640 PetscInt gid, lid, i, ec, nz = aij->nz; 3641 PetscInt *garray, *jj = aij->j; 3642 3643 PetscFunctionBegin; 3644 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 3645 PetscValidPointer(mapping, 2); 3646 /* use a table */ 3647 PetscCall(PetscHMapICreateWithSize(mat->rmap->n, &gid1_lid1)); 3648 ec = 0; 3649 for (i = 0; i < nz; i++) { 3650 PetscInt data, gid1 = jj[i] + 1; 3651 PetscCall(PetscHMapIGetWithDefault(gid1_lid1, gid1, 0, &data)); 3652 if (!data) { 3653 /* one based table */ 3654 PetscCall(PetscHMapISet(gid1_lid1, gid1, ++ec)); 3655 } 3656 } 3657 /* form array of columns we need */ 3658 PetscCall(PetscMalloc1(ec, &garray)); 3659 PetscHashIterBegin(gid1_lid1, tpos); 3660 while (!PetscHashIterAtEnd(gid1_lid1, tpos)) { 3661 PetscHashIterGetKey(gid1_lid1, tpos, gid); 3662 PetscHashIterGetVal(gid1_lid1, tpos, lid); 3663 PetscHashIterNext(gid1_lid1, tpos); 3664 gid--; 3665 lid--; 3666 garray[lid] = gid; 3667 } 3668 PetscCall(PetscSortInt(ec, garray)); /* sort, and rebuild */ 3669 PetscCall(PetscHMapIClear(gid1_lid1)); 3670 for (i = 0; i < ec; i++) PetscCall(PetscHMapISet(gid1_lid1, garray[i] + 1, i + 1)); 3671 /* compact out the extra columns in B */ 3672 for (i = 0; i < nz; i++) { 3673 PetscInt gid1 = jj[i] + 1; 3674 PetscCall(PetscHMapIGetWithDefault(gid1_lid1, gid1, 0, &lid)); 3675 lid--; 3676 jj[i] = lid; 3677 } 3678 PetscCall(PetscLayoutDestroy(&mat->cmap)); 3679 PetscCall(PetscHMapIDestroy(&gid1_lid1)); 3680 PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)mat), ec, ec, 1, &mat->cmap)); 3681 PetscCall(ISLocalToGlobalMappingCreate(PETSC_COMM_SELF, mat->cmap->bs, mat->cmap->n, garray, PETSC_OWN_POINTER, mapping)); 3682 PetscCall(ISLocalToGlobalMappingSetType(*mapping, ISLOCALTOGLOBALMAPPINGHASH)); 3683 PetscFunctionReturn(0); 3684 } 3685 3686 /*@ 3687 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 3688 in the matrix. 3689 3690 Input Parameters: 3691 + mat - the `MATSEQAIJ` matrix 3692 - indices - the column indices 3693 3694 Level: advanced 3695 3696 Notes: 3697 This can be called if you have precomputed the nonzero structure of the 3698 matrix and want to provide it to the matrix object to improve the performance 3699 of the `MatSetValues()` operation. 3700 3701 You MUST have set the correct numbers of nonzeros per row in the call to 3702 `MatCreateSeqAIJ()`, and the columns indices MUST be sorted. 3703 3704 MUST be called before any calls to `MatSetValues()` 3705 3706 The indices should start with zero, not one. 3707 3708 @*/ 3709 PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat, PetscInt *indices) 3710 { 3711 PetscFunctionBegin; 3712 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 3713 PetscValidIntPointer(indices, 2); 3714 PetscUseMethod(mat, "MatSeqAIJSetColumnIndices_C", (Mat, PetscInt *), (mat, indices)); 3715 PetscFunctionReturn(0); 3716 } 3717 3718 /* ----------------------------------------------------------------------------------------*/ 3719 3720 PetscErrorCode MatStoreValues_SeqAIJ(Mat mat) 3721 { 3722 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 3723 size_t nz = aij->i[mat->rmap->n]; 3724 3725 PetscFunctionBegin; 3726 PetscCheck(aij->nonew, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3727 3728 /* allocate space for values if not already there */ 3729 if (!aij->saved_values) { PetscCall(PetscMalloc1(nz + 1, &aij->saved_values)); } 3730 3731 /* copy values over */ 3732 PetscCall(PetscArraycpy(aij->saved_values, aij->a, nz)); 3733 PetscFunctionReturn(0); 3734 } 3735 3736 /*@ 3737 MatStoreValues - Stashes a copy of the matrix values; this allows, for 3738 example, reuse of the linear part of a Jacobian, while recomputing the 3739 nonlinear portion. 3740 3741 Logically Collect 3742 3743 Input Parameters: 3744 . mat - the matrix (currently only `MATAIJ` matrices support this option) 3745 3746 Level: advanced 3747 3748 Common Usage, with `SNESSolve()`: 3749 $ Create Jacobian matrix 3750 $ Set linear terms into matrix 3751 $ Apply boundary conditions to matrix, at this time matrix must have 3752 $ final nonzero structure (i.e. setting the nonlinear terms and applying 3753 $ boundary conditions again will not change the nonzero structure 3754 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3755 $ ierr = MatStoreValues(mat); 3756 $ Call SNESSetJacobian() with matrix 3757 $ In your Jacobian routine 3758 $ ierr = MatRetrieveValues(mat); 3759 $ Set nonlinear terms in matrix 3760 3761 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 3762 $ // build linear portion of Jacobian 3763 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3764 $ ierr = MatStoreValues(mat); 3765 $ loop over nonlinear iterations 3766 $ ierr = MatRetrieveValues(mat); 3767 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 3768 $ // call MatAssemblyBegin/End() on matrix 3769 $ Solve linear system with Jacobian 3770 $ endloop 3771 3772 Notes: 3773 Matrix must already be assemblied before calling this routine 3774 Must set the matrix option `MatSetOption`(mat,`MAT_NEW_NONZERO_LOCATIONS`,`PETSC_FALSE`); before 3775 calling this routine. 3776 3777 When this is called multiple times it overwrites the previous set of stored values 3778 and does not allocated additional space. 3779 3780 .seealso: `MatRetrieveValues()` 3781 @*/ 3782 PetscErrorCode MatStoreValues(Mat mat) 3783 { 3784 PetscFunctionBegin; 3785 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 3786 PetscCheck(mat->assembled, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix"); 3787 PetscCheck(!mat->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix"); 3788 PetscUseMethod(mat, "MatStoreValues_C", (Mat), (mat)); 3789 PetscFunctionReturn(0); 3790 } 3791 3792 PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat) 3793 { 3794 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 3795 PetscInt nz = aij->i[mat->rmap->n]; 3796 3797 PetscFunctionBegin; 3798 PetscCheck(aij->nonew, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3799 PetscCheck(aij->saved_values, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call MatStoreValues(A);first"); 3800 /* copy values over */ 3801 PetscCall(PetscArraycpy(aij->a, aij->saved_values, nz)); 3802 PetscFunctionReturn(0); 3803 } 3804 3805 /*@ 3806 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 3807 example, reuse of the linear part of a Jacobian, while recomputing the 3808 nonlinear portion. 3809 3810 Logically Collect 3811 3812 Input Parameters: 3813 . mat - the matrix (currently only `MATAIJ` matrices support this option) 3814 3815 Level: advanced 3816 3817 .seealso: `MatStoreValues()` 3818 @*/ 3819 PetscErrorCode MatRetrieveValues(Mat mat) 3820 { 3821 PetscFunctionBegin; 3822 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 3823 PetscCheck(mat->assembled, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix"); 3824 PetscCheck(!mat->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix"); 3825 PetscUseMethod(mat, "MatRetrieveValues_C", (Mat), (mat)); 3826 PetscFunctionReturn(0); 3827 } 3828 3829 /* --------------------------------------------------------------------------------*/ 3830 /*@C 3831 MatCreateSeqAIJ - Creates a sparse matrix in `MATSEQAIJ` (compressed row) format 3832 (the default parallel PETSc format). For good matrix assembly performance 3833 the user should preallocate the matrix storage by setting the parameter nz 3834 (or the array nnz). By setting these parameters accurately, performance 3835 during matrix assembly can be increased by more than a factor of 50. 3836 3837 Collective 3838 3839 Input Parameters: 3840 + comm - MPI communicator, set to `PETSC_COMM_SELF` 3841 . m - number of rows 3842 . n - number of columns 3843 . nz - number of nonzeros per row (same for all rows) 3844 - nnz - array containing the number of nonzeros in the various rows 3845 (possibly different for each row) or NULL 3846 3847 Output Parameter: 3848 . A - the matrix 3849 3850 It is recommended that one use the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`, 3851 MatXXXXSetPreallocation() paradigm instead of this routine directly. 3852 [MatXXXXSetPreallocation() is, for example, `MatSeqAIJSetPreallocation()`] 3853 3854 Notes: 3855 If nnz is given then nz is ignored 3856 3857 The AIJ format, also called 3858 compressed row storage, is fully compatible with standard Fortran 77 3859 storage. That is, the stored row and column indices can begin at 3860 either one (as in Fortran) or zero. See the users' manual for details. 3861 3862 Specify the preallocated storage with either nz or nnz (not both). 3863 Set nz = `PETSC_DEFAULT` and nnz = NULL for PETSc to control dynamic memory 3864 allocation. For large problems you MUST preallocate memory or you 3865 will get TERRIBLE performance, see the users' manual chapter on matrices. 3866 3867 By default, this format uses inodes (identical nodes) when possible, to 3868 improve numerical efficiency of matrix-vector products and solves. We 3869 search for consecutive rows with the same nonzero structure, thereby 3870 reusing matrix information to achieve increased efficiency. 3871 3872 Options Database Keys: 3873 + -mat_no_inode - Do not use inodes 3874 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3875 3876 Level: intermediate 3877 3878 .seealso: [Sparse Matrix Creation](sec_matsparse), `MatCreate()`, `MatCreateAIJ()`, `MatSetValues()`, `MatSeqAIJSetColumnIndices()`, `MatCreateSeqAIJWithArrays()` 3879 @*/ 3880 PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt nz, const PetscInt nnz[], Mat *A) 3881 { 3882 PetscFunctionBegin; 3883 PetscCall(MatCreate(comm, A)); 3884 PetscCall(MatSetSizes(*A, m, n, m, n)); 3885 PetscCall(MatSetType(*A, MATSEQAIJ)); 3886 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(*A, nz, nnz)); 3887 PetscFunctionReturn(0); 3888 } 3889 3890 /*@C 3891 MatSeqAIJSetPreallocation - For good matrix assembly performance 3892 the user should preallocate the matrix storage by setting the parameter nz 3893 (or the array nnz). By setting these parameters accurately, performance 3894 during matrix assembly can be increased by more than a factor of 50. 3895 3896 Collective 3897 3898 Input Parameters: 3899 + B - The matrix 3900 . nz - number of nonzeros per row (same for all rows) 3901 - nnz - array containing the number of nonzeros in the various rows 3902 (possibly different for each row) or NULL 3903 3904 Notes: 3905 If nnz is given then nz is ignored 3906 3907 The `MATSEQAIJ` format also called 3908 compressed row storage, is fully compatible with standard Fortran 77 3909 storage. That is, the stored row and column indices can begin at 3910 either one (as in Fortran) or zero. See the users' manual for details. 3911 3912 Specify the preallocated storage with either nz or nnz (not both). 3913 Set nz = `PETSC_DEFAULT` and nnz = NULL for PETSc to control dynamic memory 3914 allocation. For large problems you MUST preallocate memory or you 3915 will get TERRIBLE performance, see the users' manual chapter on matrices. 3916 3917 You can call `MatGetInfo()` to get information on how effective the preallocation was; 3918 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 3919 You can also run with the option -info and look for messages with the string 3920 malloc in them to see if additional memory allocation was needed. 3921 3922 Developer Notes: 3923 Use nz of `MAT_SKIP_ALLOCATION` to not allocate any space for the matrix 3924 entries or columns indices 3925 3926 By default, this format uses inodes (identical nodes) when possible, to 3927 improve numerical efficiency of matrix-vector products and solves. We 3928 search for consecutive rows with the same nonzero structure, thereby 3929 reusing matrix information to achieve increased efficiency. 3930 3931 Options Database Keys: 3932 + -mat_no_inode - Do not use inodes 3933 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3934 3935 Level: intermediate 3936 3937 .seealso: `MatCreate()`, `MatCreateAIJ()`, `MatSetValues()`, `MatSeqAIJSetColumnIndices()`, `MatCreateSeqAIJWithArrays()`, `MatGetInfo()`, 3938 `MatSeqAIJSetTotalPreallocation()` 3939 @*/ 3940 PetscErrorCode MatSeqAIJSetPreallocation(Mat B, PetscInt nz, const PetscInt nnz[]) 3941 { 3942 PetscFunctionBegin; 3943 PetscValidHeaderSpecific(B, MAT_CLASSID, 1); 3944 PetscValidType(B, 1); 3945 PetscTryMethod(B, "MatSeqAIJSetPreallocation_C", (Mat, PetscInt, const PetscInt[]), (B, nz, nnz)); 3946 PetscFunctionReturn(0); 3947 } 3948 3949 PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B, PetscInt nz, const PetscInt *nnz) 3950 { 3951 Mat_SeqAIJ *b; 3952 PetscBool skipallocation = PETSC_FALSE, realalloc = PETSC_FALSE; 3953 PetscInt i; 3954 3955 PetscFunctionBegin; 3956 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 3957 if (nz == MAT_SKIP_ALLOCATION) { 3958 skipallocation = PETSC_TRUE; 3959 nz = 0; 3960 } 3961 PetscCall(PetscLayoutSetUp(B->rmap)); 3962 PetscCall(PetscLayoutSetUp(B->cmap)); 3963 3964 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3965 PetscCheck(nz >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "nz cannot be less than 0: value %" PetscInt_FMT, nz); 3966 if (PetscUnlikelyDebug(nnz)) { 3967 for (i = 0; i < B->rmap->n; i++) { 3968 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]); 3969 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); 3970 } 3971 } 3972 3973 B->preallocated = PETSC_TRUE; 3974 3975 b = (Mat_SeqAIJ *)B->data; 3976 3977 if (!skipallocation) { 3978 if (!b->imax) { PetscCall(PetscMalloc1(B->rmap->n, &b->imax)); } 3979 if (!b->ilen) { 3980 /* b->ilen will count nonzeros in each row so far. */ 3981 PetscCall(PetscCalloc1(B->rmap->n, &b->ilen)); 3982 } else { 3983 PetscCall(PetscMemzero(b->ilen, B->rmap->n * sizeof(PetscInt))); 3984 } 3985 if (!b->ipre) { PetscCall(PetscMalloc1(B->rmap->n, &b->ipre)); } 3986 if (!nnz) { 3987 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 3988 else if (nz < 0) nz = 1; 3989 nz = PetscMin(nz, B->cmap->n); 3990 for (i = 0; i < B->rmap->n; i++) b->imax[i] = nz; 3991 nz = nz * B->rmap->n; 3992 } else { 3993 PetscInt64 nz64 = 0; 3994 for (i = 0; i < B->rmap->n; i++) { 3995 b->imax[i] = nnz[i]; 3996 nz64 += nnz[i]; 3997 } 3998 PetscCall(PetscIntCast(nz64, &nz)); 3999 } 4000 4001 /* allocate the matrix space */ 4002 /* FIXME: should B's old memory be unlogged? */ 4003 PetscCall(MatSeqXAIJFreeAIJ(B, &b->a, &b->j, &b->i)); 4004 if (B->structure_only) { 4005 PetscCall(PetscMalloc1(nz, &b->j)); 4006 PetscCall(PetscMalloc1(B->rmap->n + 1, &b->i)); 4007 } else { 4008 PetscCall(PetscMalloc3(nz, &b->a, nz, &b->j, B->rmap->n + 1, &b->i)); 4009 } 4010 b->i[0] = 0; 4011 for (i = 1; i < B->rmap->n + 1; i++) b->i[i] = b->i[i - 1] + b->imax[i - 1]; 4012 if (B->structure_only) { 4013 b->singlemalloc = PETSC_FALSE; 4014 b->free_a = PETSC_FALSE; 4015 } else { 4016 b->singlemalloc = PETSC_TRUE; 4017 b->free_a = PETSC_TRUE; 4018 } 4019 b->free_ij = PETSC_TRUE; 4020 } else { 4021 b->free_a = PETSC_FALSE; 4022 b->free_ij = PETSC_FALSE; 4023 } 4024 4025 if (b->ipre && nnz != b->ipre && b->imax) { 4026 /* reserve user-requested sparsity */ 4027 PetscCall(PetscArraycpy(b->ipre, b->imax, B->rmap->n)); 4028 } 4029 4030 b->nz = 0; 4031 b->maxnz = nz; 4032 B->info.nz_unneeded = (double)b->maxnz; 4033 if (realalloc) PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE)); 4034 B->was_assembled = PETSC_FALSE; 4035 B->assembled = PETSC_FALSE; 4036 /* We simply deem preallocation has changed nonzero state. Updating the state 4037 will give clients (like AIJKokkos) a chance to know something has happened. 4038 */ 4039 B->nonzerostate++; 4040 PetscFunctionReturn(0); 4041 } 4042 4043 PetscErrorCode MatResetPreallocation_SeqAIJ(Mat A) 4044 { 4045 Mat_SeqAIJ *a; 4046 PetscInt i; 4047 4048 PetscFunctionBegin; 4049 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 4050 4051 /* Check local size. If zero, then return */ 4052 if (!A->rmap->n) PetscFunctionReturn(0); 4053 4054 a = (Mat_SeqAIJ *)A->data; 4055 /* if no saved info, we error out */ 4056 PetscCheck(a->ipre, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "No saved preallocation info "); 4057 4058 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 "); 4059 4060 PetscCall(PetscArraycpy(a->imax, a->ipre, A->rmap->n)); 4061 PetscCall(PetscArrayzero(a->ilen, A->rmap->n)); 4062 a->i[0] = 0; 4063 for (i = 1; i < A->rmap->n + 1; i++) a->i[i] = a->i[i - 1] + a->imax[i - 1]; 4064 A->preallocated = PETSC_TRUE; 4065 a->nz = 0; 4066 a->maxnz = a->i[A->rmap->n]; 4067 A->info.nz_unneeded = (double)a->maxnz; 4068 A->was_assembled = PETSC_FALSE; 4069 A->assembled = PETSC_FALSE; 4070 PetscFunctionReturn(0); 4071 } 4072 4073 /*@ 4074 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in `MATSEQAIJ` format. 4075 4076 Input Parameters: 4077 + B - the matrix 4078 . i - the indices into j for the start of each row (starts with zero) 4079 . j - the column indices for each row (starts with zero) these must be sorted for each row 4080 - v - optional values in the matrix 4081 4082 Level: developer 4083 4084 Notes: 4085 The i,j,v values are COPIED with this routine; to avoid the copy use `MatCreateSeqAIJWithArrays()` 4086 4087 This routine may be called multiple times with different nonzero patterns (or the same nonzero pattern). The nonzero 4088 structure will be the union of all the previous nonzero structures. 4089 4090 Developer Notes: 4091 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 4092 then just copies the v values directly with `PetscMemcpy()`. 4093 4094 This routine could also take a `PetscCopyMode` argument to allow sharing the values instead of always copying them. 4095 4096 .seealso: `MatCreate()`, `MatCreateSeqAIJ()`, `MatSetValues()`, `MatSeqAIJSetPreallocation()`, `MatCreateSeqAIJ()`, `MATSEQAIJ`, `MatResetPreallocation()` 4097 @*/ 4098 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B, const PetscInt i[], const PetscInt j[], const PetscScalar v[]) 4099 { 4100 PetscFunctionBegin; 4101 PetscValidHeaderSpecific(B, MAT_CLASSID, 1); 4102 PetscValidType(B, 1); 4103 PetscTryMethod(B, "MatSeqAIJSetPreallocationCSR_C", (Mat, const PetscInt[], const PetscInt[], const PetscScalar[]), (B, i, j, v)); 4104 PetscFunctionReturn(0); 4105 } 4106 4107 PetscErrorCode MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B, const PetscInt Ii[], const PetscInt J[], const PetscScalar v[]) 4108 { 4109 PetscInt i; 4110 PetscInt m, n; 4111 PetscInt nz; 4112 PetscInt *nnz; 4113 4114 PetscFunctionBegin; 4115 PetscCheck(Ii[0] == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %" PetscInt_FMT, Ii[0]); 4116 4117 PetscCall(PetscLayoutSetUp(B->rmap)); 4118 PetscCall(PetscLayoutSetUp(B->cmap)); 4119 4120 PetscCall(MatGetSize(B, &m, &n)); 4121 PetscCall(PetscMalloc1(m + 1, &nnz)); 4122 for (i = 0; i < m; i++) { 4123 nz = Ii[i + 1] - Ii[i]; 4124 PetscCheck(nz >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local row %" PetscInt_FMT " has a negative number of columns %" PetscInt_FMT, i, nz); 4125 nnz[i] = nz; 4126 } 4127 PetscCall(MatSeqAIJSetPreallocation(B, 0, nnz)); 4128 PetscCall(PetscFree(nnz)); 4129 4130 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)); 4131 4132 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 4133 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 4134 4135 PetscCall(MatSetOption(B, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE)); 4136 PetscFunctionReturn(0); 4137 } 4138 4139 /*@ 4140 MatSeqAIJKron - Computes C, the Kronecker product of A and B. 4141 4142 Input Parameters: 4143 + A - left-hand side matrix 4144 . B - right-hand side matrix 4145 - reuse - either `MAT_INITIAL_MATRIX` or `MAT_REUSE_MATRIX` 4146 4147 Output Parameter: 4148 . C - Kronecker product of A and B 4149 4150 Level: intermediate 4151 4152 Note: 4153 `MAT_REUSE_MATRIX` can only be used when the nonzero structure of the product matrix has not changed from that last call to `MatSeqAIJKron()`. 4154 4155 .seealso: `MatCreateSeqAIJ()`, `MATSEQAIJ`, `MATKAIJ`, `MatReuse` 4156 @*/ 4157 PetscErrorCode MatSeqAIJKron(Mat A, Mat B, MatReuse reuse, Mat *C) 4158 { 4159 PetscFunctionBegin; 4160 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 4161 PetscValidType(A, 1); 4162 PetscValidHeaderSpecific(B, MAT_CLASSID, 2); 4163 PetscValidType(B, 2); 4164 PetscValidPointer(C, 4); 4165 if (reuse == MAT_REUSE_MATRIX) { 4166 PetscValidHeaderSpecific(*C, MAT_CLASSID, 4); 4167 PetscValidType(*C, 4); 4168 } 4169 PetscTryMethod(A, "MatSeqAIJKron_C", (Mat, Mat, MatReuse, Mat *), (A, B, reuse, C)); 4170 PetscFunctionReturn(0); 4171 } 4172 4173 PetscErrorCode MatSeqAIJKron_SeqAIJ(Mat A, Mat B, MatReuse reuse, Mat *C) 4174 { 4175 Mat newmat; 4176 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 4177 Mat_SeqAIJ *b = (Mat_SeqAIJ *)B->data; 4178 PetscScalar *v; 4179 const PetscScalar *aa, *ba; 4180 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; 4181 PetscBool flg; 4182 4183 PetscFunctionBegin; 4184 PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix"); 4185 PetscCheck(A->assembled, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix"); 4186 PetscCheck(!B->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix"); 4187 PetscCheck(B->assembled, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix"); 4188 PetscCall(PetscObjectTypeCompare((PetscObject)B, MATSEQAIJ, &flg)); 4189 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "MatType %s", ((PetscObject)B)->type_name); 4190 PetscCheck(reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX, PETSC_COMM_SELF, PETSC_ERR_SUP, "MatReuse %d", (int)reuse); 4191 if (reuse == MAT_INITIAL_MATRIX) { 4192 PetscCall(PetscMalloc2(am * bm + 1, &i, a->i[am] * b->i[bm], &j)); 4193 PetscCall(MatCreate(PETSC_COMM_SELF, &newmat)); 4194 PetscCall(MatSetSizes(newmat, am * bm, an * bn, am * bm, an * bn)); 4195 PetscCall(MatSetType(newmat, MATAIJ)); 4196 i[0] = 0; 4197 for (m = 0; m < am; ++m) { 4198 for (p = 0; p < bm; ++p) { 4199 i[m * bm + p + 1] = i[m * bm + p] + (a->i[m + 1] - a->i[m]) * (b->i[p + 1] - b->i[p]); 4200 for (n = a->i[m]; n < a->i[m + 1]; ++n) { 4201 for (q = b->i[p]; q < b->i[p + 1]; ++q) j[nnz++] = a->j[n] * bn + b->j[q]; 4202 } 4203 } 4204 } 4205 PetscCall(MatSeqAIJSetPreallocationCSR(newmat, i, j, NULL)); 4206 *C = newmat; 4207 PetscCall(PetscFree2(i, j)); 4208 nnz = 0; 4209 } 4210 PetscCall(MatSeqAIJGetArray(*C, &v)); 4211 PetscCall(MatSeqAIJGetArrayRead(A, &aa)); 4212 PetscCall(MatSeqAIJGetArrayRead(B, &ba)); 4213 for (m = 0; m < am; ++m) { 4214 for (p = 0; p < bm; ++p) { 4215 for (n = a->i[m]; n < a->i[m + 1]; ++n) { 4216 for (q = b->i[p]; q < b->i[p + 1]; ++q) v[nnz++] = aa[n] * ba[q]; 4217 } 4218 } 4219 } 4220 PetscCall(MatSeqAIJRestoreArray(*C, &v)); 4221 PetscCall(MatSeqAIJRestoreArrayRead(A, &aa)); 4222 PetscCall(MatSeqAIJRestoreArrayRead(B, &ba)); 4223 PetscFunctionReturn(0); 4224 } 4225 4226 #include <../src/mat/impls/dense/seq/dense.h> 4227 #include <petsc/private/kernels/petscaxpy.h> 4228 4229 /* 4230 Computes (B'*A')' since computing B*A directly is untenable 4231 4232 n p p 4233 [ ] [ ] [ ] 4234 m [ A ] * n [ B ] = m [ C ] 4235 [ ] [ ] [ ] 4236 4237 */ 4238 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A, Mat B, Mat C) 4239 { 4240 Mat_SeqDense *sub_a = (Mat_SeqDense *)A->data; 4241 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ *)B->data; 4242 Mat_SeqDense *sub_c = (Mat_SeqDense *)C->data; 4243 PetscInt i, j, n, m, q, p; 4244 const PetscInt *ii, *idx; 4245 const PetscScalar *b, *a, *a_q; 4246 PetscScalar *c, *c_q; 4247 PetscInt clda = sub_c->lda; 4248 PetscInt alda = sub_a->lda; 4249 4250 PetscFunctionBegin; 4251 m = A->rmap->n; 4252 n = A->cmap->n; 4253 p = B->cmap->n; 4254 a = sub_a->v; 4255 b = sub_b->a; 4256 c = sub_c->v; 4257 if (clda == m) { 4258 PetscCall(PetscArrayzero(c, m * p)); 4259 } else { 4260 for (j = 0; j < p; j++) 4261 for (i = 0; i < m; i++) c[j * clda + i] = 0.0; 4262 } 4263 ii = sub_b->i; 4264 idx = sub_b->j; 4265 for (i = 0; i < n; i++) { 4266 q = ii[i + 1] - ii[i]; 4267 while (q-- > 0) { 4268 c_q = c + clda * (*idx); 4269 a_q = a + alda * i; 4270 PetscKernelAXPY(c_q, *b, a_q, m); 4271 idx++; 4272 b++; 4273 } 4274 } 4275 PetscFunctionReturn(0); 4276 } 4277 4278 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A, Mat B, PetscReal fill, Mat C) 4279 { 4280 PetscInt m = A->rmap->n, n = B->cmap->n; 4281 PetscBool cisdense; 4282 4283 PetscFunctionBegin; 4284 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); 4285 PetscCall(MatSetSizes(C, m, n, m, n)); 4286 PetscCall(MatSetBlockSizesFromMats(C, A, B)); 4287 PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSECUDA, MATSEQDENSEHIP, "")); 4288 if (!cisdense) PetscCall(MatSetType(C, MATDENSE)); 4289 PetscCall(MatSetUp(C)); 4290 4291 C->ops->matmultnumeric = MatMatMultNumeric_SeqDense_SeqAIJ; 4292 PetscFunctionReturn(0); 4293 } 4294 4295 /* ----------------------------------------------------------------*/ 4296 /*MC 4297 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 4298 based on compressed sparse row format. 4299 4300 Options Database Keys: 4301 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 4302 4303 Level: beginner 4304 4305 Notes: 4306 `MatSetValues()` may be called for this matrix type with a NULL argument for the numerical values, 4307 in this case the values associated with the rows and columns one passes in are set to zero 4308 in the matrix 4309 4310 `MatSetOptions`(,`MAT_STRUCTURE_ONLY`,`PETSC_TRUE`) may be called for this matrix type. In this no 4311 space is allocated for the nonzero entries and any entries passed with `MatSetValues()` are ignored 4312 4313 Developer Note: 4314 It would be nice if all matrix formats supported passing NULL in for the numerical values 4315 4316 .seealso: `MatCreateSeqAIJ()`, `MatSetFromOptions()`, `MatSetType()`, `MatCreate()`, `MatType`, `MATSELL`, `MATSEQSELL`, `MATMPISELL` 4317 M*/ 4318 4319 /*MC 4320 MATAIJ - MATAIJ = "aij" - A matrix type to be used for sparse matrices. 4321 4322 This matrix type is identical to `MATSEQAIJ` when constructed with a single process communicator, 4323 and `MATMPIAIJ` otherwise. As a result, for single process communicators, 4324 `MatSeqAIJSetPreallocation()` is supported, and similarly `MatMPIAIJSetPreallocation()` is supported 4325 for communicators controlling multiple processes. It is recommended that you call both of 4326 the above preallocation routines for simplicity. 4327 4328 Options Database Keys: 4329 . -mat_type aij - sets the matrix type to "aij" during a call to `MatSetFromOptions()` 4330 4331 Note: 4332 Subclasses include `MATAIJCUSPARSE`, `MATAIJPERM`, `MATAIJSELL`, `MATAIJMKL`, `MATAIJCRL`, and also automatically switches over to use inodes when 4333 enough exist. 4334 4335 Level: beginner 4336 4337 .seealso: `MatCreateAIJ()`, `MatCreateSeqAIJ()`, `MATSEQAIJ`, `MATMPIAIJ`, `MATSELL`, `MATSEQSELL`, `MATMPISELL` 4338 M*/ 4339 4340 /*MC 4341 MATAIJCRL - MATAIJCRL = "aijcrl" - A matrix type to be used for sparse matrices. 4342 4343 This matrix type is identical to `MATSEQAIJCRL` when constructed with a single process communicator, 4344 and `MATMPIAIJCRL` otherwise. As a result, for single process communicators, 4345 `MatSeqAIJSetPreallocation()` is supported, and similarly `MatMPIAIJSetPreallocation()` is supported 4346 for communicators controlling multiple processes. It is recommended that you call both of 4347 the above preallocation routines for simplicity. 4348 4349 Options Database Keys: 4350 . -mat_type aijcrl - sets the matrix type to "aijcrl" during a call to `MatSetFromOptions()` 4351 4352 Level: beginner 4353 4354 .seealso: `MatCreateMPIAIJCRL`, `MATSEQAIJCRL`, `MATMPIAIJCRL`, `MATSEQAIJCRL`, `MATMPIAIJCRL` 4355 M*/ 4356 4357 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCRL(Mat, MatType, MatReuse, Mat *); 4358 #if defined(PETSC_HAVE_ELEMENTAL) 4359 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_Elemental(Mat, MatType, MatReuse, Mat *); 4360 #endif 4361 #if defined(PETSC_HAVE_SCALAPACK) 4362 PETSC_INTERN PetscErrorCode MatConvert_AIJ_ScaLAPACK(Mat, MatType, MatReuse, Mat *); 4363 #endif 4364 #if defined(PETSC_HAVE_HYPRE) 4365 PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A, MatType, MatReuse, Mat *); 4366 #endif 4367 4368 PETSC_EXTERN PetscErrorCode MatConvert_SeqAIJ_SeqSELL(Mat, MatType, MatReuse, Mat *); 4369 PETSC_INTERN PetscErrorCode MatConvert_XAIJ_IS(Mat, MatType, MatReuse, Mat *); 4370 PETSC_INTERN PetscErrorCode MatProductSetFromOptions_IS_XAIJ(Mat); 4371 4372 /*@C 4373 MatSeqAIJGetArray - gives read/write access to the array where the data for a `MATSEQAIJ` matrix is stored 4374 4375 Not Collective 4376 4377 Input Parameter: 4378 . mat - a `MATSEQAIJ` matrix 4379 4380 Output Parameter: 4381 . array - pointer to the data 4382 4383 Level: intermediate 4384 4385 .seealso: `MatSeqAIJRestoreArray()`, `MatSeqAIJGetArrayF90()` 4386 @*/ 4387 PetscErrorCode MatSeqAIJGetArray(Mat A, PetscScalar **array) 4388 { 4389 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)A->data; 4390 4391 PetscFunctionBegin; 4392 if (aij->ops->getarray) { 4393 PetscCall((*aij->ops->getarray)(A, array)); 4394 } else { 4395 *array = aij->a; 4396 } 4397 PetscFunctionReturn(0); 4398 } 4399 4400 /*@C 4401 MatSeqAIJRestoreArray - returns access to the array where the data for a `MATSEQAIJ` matrix is stored obtained by `MatSeqAIJGetArray()` 4402 4403 Not Collective 4404 4405 Input Parameters: 4406 + mat - a `MATSEQAIJ` matrix 4407 - array - pointer to the data 4408 4409 Level: intermediate 4410 4411 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArrayF90()` 4412 @*/ 4413 PetscErrorCode MatSeqAIJRestoreArray(Mat A, PetscScalar **array) 4414 { 4415 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)A->data; 4416 4417 PetscFunctionBegin; 4418 if (aij->ops->restorearray) { 4419 PetscCall((*aij->ops->restorearray)(A, array)); 4420 } else { 4421 *array = NULL; 4422 } 4423 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 4424 PetscCall(PetscObjectStateIncrease((PetscObject)A)); 4425 PetscFunctionReturn(0); 4426 } 4427 4428 /*@C 4429 MatSeqAIJGetArrayRead - gives read-only access to the array where the data for a `MATSEQAIJ` matrix is stored 4430 4431 Not Collective 4432 4433 Input Parameter: 4434 . mat - a `MATSEQAIJ` matrix 4435 4436 Output Parameter: 4437 . array - pointer to the data 4438 4439 Level: intermediate 4440 4441 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArrayRead()` 4442 @*/ 4443 PetscErrorCode MatSeqAIJGetArrayRead(Mat A, const PetscScalar **array) 4444 { 4445 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)A->data; 4446 4447 PetscFunctionBegin; 4448 if (aij->ops->getarrayread) { 4449 PetscCall((*aij->ops->getarrayread)(A, array)); 4450 } else { 4451 *array = aij->a; 4452 } 4453 PetscFunctionReturn(0); 4454 } 4455 4456 /*@C 4457 MatSeqAIJRestoreArrayRead - restore the read-only access array obtained from `MatSeqAIJGetArrayRead()` 4458 4459 Not Collective 4460 4461 Input Parameter: 4462 . mat - a `MATSEQAIJ` matrix 4463 4464 Output Parameter: 4465 . array - pointer to the data 4466 4467 Level: intermediate 4468 4469 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJGetArrayRead()` 4470 @*/ 4471 PetscErrorCode MatSeqAIJRestoreArrayRead(Mat A, const PetscScalar **array) 4472 { 4473 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)A->data; 4474 4475 PetscFunctionBegin; 4476 if (aij->ops->restorearrayread) { 4477 PetscCall((*aij->ops->restorearrayread)(A, array)); 4478 } else { 4479 *array = NULL; 4480 } 4481 PetscFunctionReturn(0); 4482 } 4483 4484 /*@C 4485 MatSeqAIJGetArrayWrite - gives write-only access to the array where the data for a `MATSEQAIJ` matrix is stored 4486 4487 Not Collective 4488 4489 Input Parameter: 4490 . mat - a `MATSEQAIJ` matrix 4491 4492 Output Parameter: 4493 . array - pointer to the data 4494 4495 Level: intermediate 4496 4497 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArrayRead()` 4498 @*/ 4499 PetscErrorCode MatSeqAIJGetArrayWrite(Mat A, PetscScalar **array) 4500 { 4501 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)A->data; 4502 4503 PetscFunctionBegin; 4504 if (aij->ops->getarraywrite) { 4505 PetscCall((*aij->ops->getarraywrite)(A, array)); 4506 } else { 4507 *array = aij->a; 4508 } 4509 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 4510 PetscCall(PetscObjectStateIncrease((PetscObject)A)); 4511 PetscFunctionReturn(0); 4512 } 4513 4514 /*@C 4515 MatSeqAIJRestoreArrayWrite - restore the read-only access array obtained from MatSeqAIJGetArrayRead 4516 4517 Not Collective 4518 4519 Input Parameter: 4520 . mat - a MATSEQAIJ matrix 4521 4522 Output Parameter: 4523 . array - pointer to the data 4524 4525 Level: intermediate 4526 4527 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJGetArrayRead()` 4528 @*/ 4529 PetscErrorCode MatSeqAIJRestoreArrayWrite(Mat A, PetscScalar **array) 4530 { 4531 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)A->data; 4532 4533 PetscFunctionBegin; 4534 if (aij->ops->restorearraywrite) { 4535 PetscCall((*aij->ops->restorearraywrite)(A, array)); 4536 } else { 4537 *array = NULL; 4538 } 4539 PetscFunctionReturn(0); 4540 } 4541 4542 /*@C 4543 MatSeqAIJGetCSRAndMemType - Get the CSR arrays and the memory type of the `MATSEQAIJ` matrix 4544 4545 Not Collective 4546 4547 Input Parameter: 4548 . mat - a matrix of type `MATSEQAIJ` or its subclasses 4549 4550 Output Parameters: 4551 + i - row map array of the matrix 4552 . j - column index array of the matrix 4553 . a - data array of the matrix 4554 - memtype - memory type of the arrays 4555 4556 Notes: 4557 Any of the output parameters can be NULL, in which case the corresponding value is not returned. 4558 If mat is a device matrix, the arrays are on the device. Otherwise, they are on the host. 4559 4560 One can call this routine on a preallocated but not assembled matrix to just get the memory of the CSR underneath the matrix. 4561 If the matrix is assembled, the data array 'a' is guaranteed to have the latest values of the matrix. 4562 4563 Level: Developer 4564 4565 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJGetArrayRead()` 4566 @*/ 4567 PetscErrorCode MatSeqAIJGetCSRAndMemType(Mat mat, const PetscInt **i, const PetscInt **j, PetscScalar **a, PetscMemType *mtype) 4568 { 4569 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 4570 4571 PetscFunctionBegin; 4572 PetscCheck(mat->preallocated, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "matrix is not preallocated"); 4573 if (aij->ops->getcsrandmemtype) { 4574 PetscCall((*aij->ops->getcsrandmemtype)(mat, i, j, a, mtype)); 4575 } else { 4576 if (i) *i = aij->i; 4577 if (j) *j = aij->j; 4578 if (a) *a = aij->a; 4579 if (mtype) *mtype = PETSC_MEMTYPE_HOST; 4580 } 4581 PetscFunctionReturn(0); 4582 } 4583 4584 /*@C 4585 MatSeqAIJGetMaxRowNonzeros - returns the maximum number of nonzeros in any row 4586 4587 Not Collective 4588 4589 Input Parameter: 4590 . mat - a `MATSEQAIJ` matrix 4591 4592 Output Parameter: 4593 . nz - the maximum number of nonzeros in any row 4594 4595 Level: intermediate 4596 4597 .seealso: `MatSeqAIJRestoreArray()`, `MatSeqAIJGetArrayF90()` 4598 @*/ 4599 PetscErrorCode MatSeqAIJGetMaxRowNonzeros(Mat A, PetscInt *nz) 4600 { 4601 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)A->data; 4602 4603 PetscFunctionBegin; 4604 *nz = aij->rmax; 4605 PetscFunctionReturn(0); 4606 } 4607 4608 PetscErrorCode MatSetPreallocationCOO_SeqAIJ(Mat mat, PetscCount coo_n, PetscInt coo_i[], PetscInt coo_j[]) 4609 { 4610 MPI_Comm comm; 4611 PetscInt *i, *j; 4612 PetscInt M, N, row; 4613 PetscCount k, p, q, nneg, nnz, start, end; /* Index the coo array, so use PetscCount as their type */ 4614 PetscInt *Ai; /* Change to PetscCount once we use it for row pointers */ 4615 PetscInt *Aj; 4616 PetscScalar *Aa; 4617 Mat_SeqAIJ *seqaij = (Mat_SeqAIJ *)(mat->data); 4618 MatType rtype; 4619 PetscCount *perm, *jmap; 4620 4621 PetscFunctionBegin; 4622 PetscCall(MatResetPreallocationCOO_SeqAIJ(mat)); 4623 PetscCall(PetscObjectGetComm((PetscObject)mat, &comm)); 4624 PetscCall(MatGetSize(mat, &M, &N)); 4625 i = coo_i; 4626 j = coo_j; 4627 PetscCall(PetscMalloc1(coo_n, &perm)); 4628 for (k = 0; k < coo_n; k++) { /* Ignore entries with negative row or col indices */ 4629 if (j[k] < 0) i[k] = -1; 4630 perm[k] = k; 4631 } 4632 4633 /* Sort by row */ 4634 PetscCall(PetscSortIntWithIntCountArrayPair(coo_n, i, j, perm)); 4635 for (k = 0; k < coo_n; k++) { 4636 if (i[k] >= 0) break; 4637 } /* Advance k to the first row with a non-negative index */ 4638 nneg = k; 4639 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 */ 4640 nnz = 0; /* Total number of unique nonzeros to be counted */ 4641 jmap++; /* Inc jmap by 1 for convinience */ 4642 4643 PetscCall(PetscCalloc1(M + 1, &Ai)); /* CSR of A */ 4644 PetscCall(PetscMalloc1(coo_n - nneg, &Aj)); /* We have at most coo_n-nneg unique nonzeros */ 4645 4646 /* In each row, sort by column, then unique column indices to get row length */ 4647 Ai++; /* Inc by 1 for convinience */ 4648 q = 0; /* q-th unique nonzero, with q starting from 0 */ 4649 while (k < coo_n) { 4650 row = i[k]; 4651 start = k; /* [start,end) indices for this row */ 4652 while (k < coo_n && i[k] == row) k++; 4653 end = k; 4654 PetscCall(PetscSortIntWithCountArray(end - start, j + start, perm + start)); 4655 /* Find number of unique col entries in this row */ 4656 Aj[q] = j[start]; /* Log the first nonzero in this row */ 4657 jmap[q] = 1; /* Number of repeats of this nozero entry */ 4658 Ai[row] = 1; 4659 nnz++; 4660 4661 for (p = start + 1; p < end; p++) { /* Scan remaining nonzero in this row */ 4662 if (j[p] != j[p - 1]) { /* Meet a new nonzero */ 4663 q++; 4664 jmap[q] = 1; 4665 Aj[q] = j[p]; 4666 Ai[row]++; 4667 nnz++; 4668 } else { 4669 jmap[q]++; 4670 } 4671 } 4672 q++; /* Move to next row and thus next unique nonzero */ 4673 } 4674 4675 Ai--; /* Back to the beginning of Ai[] */ 4676 for (k = 0; k < M; k++) Ai[k + 1] += Ai[k]; 4677 jmap--; /* Back to the beginning of jmap[] */ 4678 jmap[0] = 0; 4679 for (k = 0; k < nnz; k++) jmap[k + 1] += jmap[k]; 4680 if (nnz < coo_n - nneg) { /* Realloc with actual number of unique nonzeros */ 4681 PetscCount *jmap_new; 4682 PetscInt *Aj_new; 4683 4684 PetscCall(PetscMalloc1(nnz + 1, &jmap_new)); 4685 PetscCall(PetscArraycpy(jmap_new, jmap, nnz + 1)); 4686 PetscCall(PetscFree(jmap)); 4687 jmap = jmap_new; 4688 4689 PetscCall(PetscMalloc1(nnz, &Aj_new)); 4690 PetscCall(PetscArraycpy(Aj_new, Aj, nnz)); 4691 PetscCall(PetscFree(Aj)); 4692 Aj = Aj_new; 4693 } 4694 4695 if (nneg) { /* Discard heading entries with negative indices in perm[], as we'll access it from index 0 in MatSetValuesCOO */ 4696 PetscCount *perm_new; 4697 4698 PetscCall(PetscMalloc1(coo_n - nneg, &perm_new)); 4699 PetscCall(PetscArraycpy(perm_new, perm + nneg, coo_n - nneg)); 4700 PetscCall(PetscFree(perm)); 4701 perm = perm_new; 4702 } 4703 4704 PetscCall(MatGetRootType_Private(mat, &rtype)); 4705 PetscCall(PetscCalloc1(nnz, &Aa)); /* Zero the matrix */ 4706 PetscCall(MatSetSeqAIJWithArrays_private(PETSC_COMM_SELF, M, N, Ai, Aj, Aa, rtype, mat)); 4707 4708 seqaij->singlemalloc = PETSC_FALSE; /* Ai, Aj and Aa are not allocated in one big malloc */ 4709 seqaij->free_a = seqaij->free_ij = PETSC_TRUE; /* Let newmat own Ai, Aj and Aa */ 4710 /* Record COO fields */ 4711 seqaij->coo_n = coo_n; 4712 seqaij->Atot = coo_n - nneg; /* Annz is seqaij->nz, so no need to record that again */ 4713 seqaij->jmap = jmap; /* of length nnz+1 */ 4714 seqaij->perm = perm; 4715 PetscFunctionReturn(0); 4716 } 4717 4718 static PetscErrorCode MatSetValuesCOO_SeqAIJ(Mat A, const PetscScalar v[], InsertMode imode) 4719 { 4720 Mat_SeqAIJ *aseq = (Mat_SeqAIJ *)A->data; 4721 PetscCount i, j, Annz = aseq->nz; 4722 PetscCount *perm = aseq->perm, *jmap = aseq->jmap; 4723 PetscScalar *Aa; 4724 4725 PetscFunctionBegin; 4726 PetscCall(MatSeqAIJGetArray(A, &Aa)); 4727 for (i = 0; i < Annz; i++) { 4728 PetscScalar sum = 0.0; 4729 for (j = jmap[i]; j < jmap[i + 1]; j++) sum += v[perm[j]]; 4730 Aa[i] = (imode == INSERT_VALUES ? 0.0 : Aa[i]) + sum; 4731 } 4732 PetscCall(MatSeqAIJRestoreArray(A, &Aa)); 4733 PetscFunctionReturn(0); 4734 } 4735 4736 #if defined(PETSC_HAVE_CUDA) 4737 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCUSPARSE(Mat, MatType, MatReuse, Mat *); 4738 #endif 4739 #if defined(PETSC_HAVE_HIP) 4740 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJHIPSPARSE(Mat, MatType, MatReuse, Mat *); 4741 #endif 4742 #if defined(PETSC_HAVE_KOKKOS_KERNELS) 4743 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJKokkos(Mat, MatType, MatReuse, Mat *); 4744 #endif 4745 4746 PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJ(Mat B) 4747 { 4748 Mat_SeqAIJ *b; 4749 PetscMPIInt size; 4750 4751 PetscFunctionBegin; 4752 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B), &size)); 4753 PetscCheck(size <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Comm must be of size 1"); 4754 4755 PetscCall(PetscNew(&b)); 4756 4757 B->data = (void *)b; 4758 4759 PetscCall(PetscMemcpy(B->ops, &MatOps_Values, sizeof(struct _MatOps))); 4760 if (B->sortedfull) B->ops->setvalues = MatSetValues_SeqAIJ_SortedFull; 4761 4762 b->row = NULL; 4763 b->col = NULL; 4764 b->icol = NULL; 4765 b->reallocs = 0; 4766 b->ignorezeroentries = PETSC_FALSE; 4767 b->roworiented = PETSC_TRUE; 4768 b->nonew = 0; 4769 b->diag = NULL; 4770 b->solve_work = NULL; 4771 B->spptr = NULL; 4772 b->saved_values = NULL; 4773 b->idiag = NULL; 4774 b->mdiag = NULL; 4775 b->ssor_work = NULL; 4776 b->omega = 1.0; 4777 b->fshift = 0.0; 4778 b->idiagvalid = PETSC_FALSE; 4779 b->ibdiagvalid = PETSC_FALSE; 4780 b->keepnonzeropattern = PETSC_FALSE; 4781 4782 PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATSEQAIJ)); 4783 #if defined(PETSC_HAVE_MATLAB) 4784 PetscCall(PetscObjectComposeFunction((PetscObject)B, "PetscMatlabEnginePut_C", MatlabEnginePut_SeqAIJ)); 4785 PetscCall(PetscObjectComposeFunction((PetscObject)B, "PetscMatlabEngineGet_C", MatlabEngineGet_SeqAIJ)); 4786 #endif 4787 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqAIJSetColumnIndices_C", MatSeqAIJSetColumnIndices_SeqAIJ)); 4788 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatStoreValues_C", MatStoreValues_SeqAIJ)); 4789 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatRetrieveValues_C", MatRetrieveValues_SeqAIJ)); 4790 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqsbaij_C", MatConvert_SeqAIJ_SeqSBAIJ)); 4791 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqbaij_C", MatConvert_SeqAIJ_SeqBAIJ)); 4792 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqaijperm_C", MatConvert_SeqAIJ_SeqAIJPERM)); 4793 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqaijsell_C", MatConvert_SeqAIJ_SeqAIJSELL)); 4794 #if defined(PETSC_HAVE_MKL_SPARSE) 4795 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqaijmkl_C", MatConvert_SeqAIJ_SeqAIJMKL)); 4796 #endif 4797 #if defined(PETSC_HAVE_CUDA) 4798 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqaijcusparse_C", MatConvert_SeqAIJ_SeqAIJCUSPARSE)); 4799 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaijcusparse_seqaij_C", MatProductSetFromOptions_SeqAIJ)); 4800 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaij_seqaijcusparse_C", MatProductSetFromOptions_SeqAIJ)); 4801 #endif 4802 #if defined(PETSC_HAVE_HIP) 4803 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqaijhipsparse_C", MatConvert_SeqAIJ_SeqAIJHIPSPARSE)); 4804 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaijhipsparse_seqaij_C", MatProductSetFromOptions_SeqAIJ)); 4805 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaij_seqaijhipsparse_C", MatProductSetFromOptions_SeqAIJ)); 4806 #endif 4807 #if defined(PETSC_HAVE_KOKKOS_KERNELS) 4808 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqaijkokkos_C", MatConvert_SeqAIJ_SeqAIJKokkos)); 4809 #endif 4810 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqaijcrl_C", MatConvert_SeqAIJ_SeqAIJCRL)); 4811 #if defined(PETSC_HAVE_ELEMENTAL) 4812 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_elemental_C", MatConvert_SeqAIJ_Elemental)); 4813 #endif 4814 #if defined(PETSC_HAVE_SCALAPACK) 4815 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_scalapack_C", MatConvert_AIJ_ScaLAPACK)); 4816 #endif 4817 #if defined(PETSC_HAVE_HYPRE) 4818 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_hypre_C", MatConvert_AIJ_HYPRE)); 4819 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_transpose_seqaij_seqaij_C", MatProductSetFromOptions_Transpose_AIJ_AIJ)); 4820 #endif 4821 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqdense_C", MatConvert_SeqAIJ_SeqDense)); 4822 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqsell_C", MatConvert_SeqAIJ_SeqSELL)); 4823 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_is_C", MatConvert_XAIJ_IS)); 4824 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatIsTranspose_C", MatIsTranspose_SeqAIJ)); 4825 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatIsHermitianTranspose_C", MatIsTranspose_SeqAIJ)); 4826 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqAIJSetPreallocation_C", MatSeqAIJSetPreallocation_SeqAIJ)); 4827 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatResetPreallocation_C", MatResetPreallocation_SeqAIJ)); 4828 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqAIJSetPreallocationCSR_C", MatSeqAIJSetPreallocationCSR_SeqAIJ)); 4829 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatReorderForNonzeroDiagonal_C", MatReorderForNonzeroDiagonal_SeqAIJ)); 4830 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_is_seqaij_C", MatProductSetFromOptions_IS_XAIJ)); 4831 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdense_seqaij_C", MatProductSetFromOptions_SeqDense_SeqAIJ)); 4832 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaij_seqaij_C", MatProductSetFromOptions_SeqAIJ)); 4833 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqAIJKron_C", MatSeqAIJKron_SeqAIJ)); 4834 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSetPreallocationCOO_C", MatSetPreallocationCOO_SeqAIJ)); 4835 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSetValuesCOO_C", MatSetValuesCOO_SeqAIJ)); 4836 PetscCall(MatCreate_SeqAIJ_Inode(B)); 4837 PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATSEQAIJ)); 4838 PetscCall(MatSeqAIJSetTypeFromOptions(B)); /* this allows changing the matrix subtype to say MATSEQAIJPERM */ 4839 PetscFunctionReturn(0); 4840 } 4841 4842 /* 4843 Given a matrix generated with MatGetFactor() duplicates all the information in A into C 4844 */ 4845 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C, Mat A, MatDuplicateOption cpvalues, PetscBool mallocmatspace) 4846 { 4847 Mat_SeqAIJ *c = (Mat_SeqAIJ *)C->data, *a = (Mat_SeqAIJ *)A->data; 4848 PetscInt m = A->rmap->n, i; 4849 4850 PetscFunctionBegin; 4851 PetscCheck(A->assembled || cpvalues == MAT_DO_NOT_COPY_VALUES, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cannot duplicate unassembled matrix"); 4852 4853 C->factortype = A->factortype; 4854 c->row = NULL; 4855 c->col = NULL; 4856 c->icol = NULL; 4857 c->reallocs = 0; 4858 4859 C->assembled = A->assembled; 4860 C->preallocated = A->preallocated; 4861 4862 if (A->preallocated) { 4863 PetscCall(PetscLayoutReference(A->rmap, &C->rmap)); 4864 PetscCall(PetscLayoutReference(A->cmap, &C->cmap)); 4865 4866 PetscCall(PetscMalloc1(m, &c->imax)); 4867 PetscCall(PetscMemcpy(c->imax, a->imax, m * sizeof(PetscInt))); 4868 PetscCall(PetscMalloc1(m, &c->ilen)); 4869 PetscCall(PetscMemcpy(c->ilen, a->ilen, m * sizeof(PetscInt))); 4870 4871 /* allocate the matrix space */ 4872 if (mallocmatspace) { 4873 PetscCall(PetscMalloc3(a->i[m], &c->a, a->i[m], &c->j, m + 1, &c->i)); 4874 4875 c->singlemalloc = PETSC_TRUE; 4876 4877 PetscCall(PetscArraycpy(c->i, a->i, m + 1)); 4878 if (m > 0) { 4879 PetscCall(PetscArraycpy(c->j, a->j, a->i[m])); 4880 if (cpvalues == MAT_COPY_VALUES) { 4881 const PetscScalar *aa; 4882 4883 PetscCall(MatSeqAIJGetArrayRead(A, &aa)); 4884 PetscCall(PetscArraycpy(c->a, aa, a->i[m])); 4885 PetscCall(MatSeqAIJGetArrayRead(A, &aa)); 4886 } else { 4887 PetscCall(PetscArrayzero(c->a, a->i[m])); 4888 } 4889 } 4890 } 4891 4892 c->ignorezeroentries = a->ignorezeroentries; 4893 c->roworiented = a->roworiented; 4894 c->nonew = a->nonew; 4895 if (a->diag) { 4896 PetscCall(PetscMalloc1(m + 1, &c->diag)); 4897 PetscCall(PetscMemcpy(c->diag, a->diag, m * sizeof(PetscInt))); 4898 } else c->diag = NULL; 4899 4900 c->solve_work = NULL; 4901 c->saved_values = NULL; 4902 c->idiag = NULL; 4903 c->ssor_work = NULL; 4904 c->keepnonzeropattern = a->keepnonzeropattern; 4905 c->free_a = PETSC_TRUE; 4906 c->free_ij = PETSC_TRUE; 4907 4908 c->rmax = a->rmax; 4909 c->nz = a->nz; 4910 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 4911 4912 c->compressedrow.use = a->compressedrow.use; 4913 c->compressedrow.nrows = a->compressedrow.nrows; 4914 if (a->compressedrow.use) { 4915 i = a->compressedrow.nrows; 4916 PetscCall(PetscMalloc2(i + 1, &c->compressedrow.i, i, &c->compressedrow.rindex)); 4917 PetscCall(PetscArraycpy(c->compressedrow.i, a->compressedrow.i, i + 1)); 4918 PetscCall(PetscArraycpy(c->compressedrow.rindex, a->compressedrow.rindex, i)); 4919 } else { 4920 c->compressedrow.use = PETSC_FALSE; 4921 c->compressedrow.i = NULL; 4922 c->compressedrow.rindex = NULL; 4923 } 4924 c->nonzerorowcnt = a->nonzerorowcnt; 4925 C->nonzerostate = A->nonzerostate; 4926 4927 PetscCall(MatDuplicate_SeqAIJ_Inode(A, cpvalues, &C)); 4928 } 4929 PetscCall(PetscFunctionListDuplicate(((PetscObject)A)->qlist, &((PetscObject)C)->qlist)); 4930 PetscFunctionReturn(0); 4931 } 4932 4933 PetscErrorCode MatDuplicate_SeqAIJ(Mat A, MatDuplicateOption cpvalues, Mat *B) 4934 { 4935 PetscFunctionBegin; 4936 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), B)); 4937 PetscCall(MatSetSizes(*B, A->rmap->n, A->cmap->n, A->rmap->n, A->cmap->n)); 4938 if (!(A->rmap->n % A->rmap->bs) && !(A->cmap->n % A->cmap->bs)) PetscCall(MatSetBlockSizesFromMats(*B, A, A)); 4939 PetscCall(MatSetType(*B, ((PetscObject)A)->type_name)); 4940 PetscCall(MatDuplicateNoCreate_SeqAIJ(*B, A, cpvalues, PETSC_TRUE)); 4941 PetscFunctionReturn(0); 4942 } 4943 4944 PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer) 4945 { 4946 PetscBool isbinary, ishdf5; 4947 4948 PetscFunctionBegin; 4949 PetscValidHeaderSpecific(newMat, MAT_CLASSID, 1); 4950 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 4951 /* force binary viewer to load .info file if it has not yet done so */ 4952 PetscCall(PetscViewerSetUp(viewer)); 4953 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary)); 4954 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5)); 4955 if (isbinary) { 4956 PetscCall(MatLoad_SeqAIJ_Binary(newMat, viewer)); 4957 } else if (ishdf5) { 4958 #if defined(PETSC_HAVE_HDF5) 4959 PetscCall(MatLoad_AIJ_HDF5(newMat, viewer)); 4960 #else 4961 SETERRQ(PetscObjectComm((PetscObject)newMat), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5"); 4962 #endif 4963 } else { 4964 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); 4965 } 4966 PetscFunctionReturn(0); 4967 } 4968 4969 PetscErrorCode MatLoad_SeqAIJ_Binary(Mat mat, PetscViewer viewer) 4970 { 4971 Mat_SeqAIJ *a = (Mat_SeqAIJ *)mat->data; 4972 PetscInt header[4], *rowlens, M, N, nz, sum, rows, cols, i; 4973 4974 PetscFunctionBegin; 4975 PetscCall(PetscViewerSetUp(viewer)); 4976 4977 /* read in matrix header */ 4978 PetscCall(PetscViewerBinaryRead(viewer, header, 4, NULL, PETSC_INT)); 4979 PetscCheck(header[0] == MAT_FILE_CLASSID, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Not a matrix object in file"); 4980 M = header[1]; 4981 N = header[2]; 4982 nz = header[3]; 4983 PetscCheck(M >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Matrix row size (%" PetscInt_FMT ") in file is negative", M); 4984 PetscCheck(N >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Matrix column size (%" PetscInt_FMT ") in file is negative", N); 4985 PetscCheck(nz >= 0, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Matrix stored in special format on disk, cannot load as SeqAIJ"); 4986 4987 /* set block sizes from the viewer's .info file */ 4988 PetscCall(MatLoad_Binary_BlockSizes(mat, viewer)); 4989 /* set local and global sizes if not set already */ 4990 if (mat->rmap->n < 0) mat->rmap->n = M; 4991 if (mat->cmap->n < 0) mat->cmap->n = N; 4992 if (mat->rmap->N < 0) mat->rmap->N = M; 4993 if (mat->cmap->N < 0) mat->cmap->N = N; 4994 PetscCall(PetscLayoutSetUp(mat->rmap)); 4995 PetscCall(PetscLayoutSetUp(mat->cmap)); 4996 4997 /* check if the matrix sizes are correct */ 4998 PetscCall(MatGetSize(mat, &rows, &cols)); 4999 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); 5000 5001 /* read in row lengths */ 5002 PetscCall(PetscMalloc1(M, &rowlens)); 5003 PetscCall(PetscViewerBinaryRead(viewer, rowlens, M, NULL, PETSC_INT)); 5004 /* check if sum(rowlens) is same as nz */ 5005 sum = 0; 5006 for (i = 0; i < M; i++) sum += rowlens[i]; 5007 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); 5008 /* preallocate and check sizes */ 5009 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(mat, 0, rowlens)); 5010 PetscCall(MatGetSize(mat, &rows, &cols)); 5011 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); 5012 /* store row lengths */ 5013 PetscCall(PetscArraycpy(a->ilen, rowlens, M)); 5014 PetscCall(PetscFree(rowlens)); 5015 5016 /* fill in "i" row pointers */ 5017 a->i[0] = 0; 5018 for (i = 0; i < M; i++) a->i[i + 1] = a->i[i] + a->ilen[i]; 5019 /* read in "j" column indices */ 5020 PetscCall(PetscViewerBinaryRead(viewer, a->j, nz, NULL, PETSC_INT)); 5021 /* read in "a" nonzero values */ 5022 PetscCall(PetscViewerBinaryRead(viewer, a->a, nz, NULL, PETSC_SCALAR)); 5023 5024 PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY)); 5025 PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY)); 5026 PetscFunctionReturn(0); 5027 } 5028 5029 PetscErrorCode MatEqual_SeqAIJ(Mat A, Mat B, PetscBool *flg) 5030 { 5031 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data; 5032 const PetscScalar *aa, *ba; 5033 #if defined(PETSC_USE_COMPLEX) 5034 PetscInt k; 5035 #endif 5036 5037 PetscFunctionBegin; 5038 /* If the matrix dimensions are not equal,or no of nonzeros */ 5039 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) || (a->nz != b->nz)) { 5040 *flg = PETSC_FALSE; 5041 PetscFunctionReturn(0); 5042 } 5043 5044 /* if the a->i are the same */ 5045 PetscCall(PetscArraycmp(a->i, b->i, A->rmap->n + 1, flg)); 5046 if (!*flg) PetscFunctionReturn(0); 5047 5048 /* if a->j are the same */ 5049 PetscCall(PetscArraycmp(a->j, b->j, a->nz, flg)); 5050 if (!*flg) PetscFunctionReturn(0); 5051 5052 PetscCall(MatSeqAIJGetArrayRead(A, &aa)); 5053 PetscCall(MatSeqAIJGetArrayRead(B, &ba)); 5054 /* if a->a are the same */ 5055 #if defined(PETSC_USE_COMPLEX) 5056 for (k = 0; k < a->nz; k++) { 5057 if (PetscRealPart(aa[k]) != PetscRealPart(ba[k]) || PetscImaginaryPart(aa[k]) != PetscImaginaryPart(ba[k])) { 5058 *flg = PETSC_FALSE; 5059 PetscFunctionReturn(0); 5060 } 5061 } 5062 #else 5063 PetscCall(PetscArraycmp(aa, ba, a->nz, flg)); 5064 #endif 5065 PetscCall(MatSeqAIJRestoreArrayRead(A, &aa)); 5066 PetscCall(MatSeqAIJRestoreArrayRead(B, &ba)); 5067 PetscFunctionReturn(0); 5068 } 5069 5070 /*@ 5071 MatCreateSeqAIJWithArrays - Creates an sequential `MATSEQAIJ` matrix using matrix elements (in CSR format) 5072 provided by the user. 5073 5074 Collective 5075 5076 Input Parameters: 5077 + comm - must be an MPI communicator of size 1 5078 . m - number of rows 5079 . n - number of columns 5080 . i - row indices; that is i[0] = 0, i[row] = i[row-1] + number of elements in that row of the matrix 5081 . j - column indices 5082 - a - matrix values 5083 5084 Output Parameter: 5085 . mat - the matrix 5086 5087 Level: intermediate 5088 5089 Notes: 5090 The i, j, and a arrays are not copied by this routine, the user must free these arrays 5091 once the matrix is destroyed and not before 5092 5093 You cannot set new nonzero locations into this matrix, that will generate an error. 5094 5095 The i and j indices are 0 based 5096 5097 The format which is used for the sparse matrix input, is equivalent to a 5098 row-major ordering.. i.e for the following matrix, the input data expected is 5099 as shown 5100 5101 $ 1 0 0 5102 $ 2 0 3 5103 $ 4 5 6 5104 $ 5105 $ i = {0,1,3,6} [size = nrow+1 = 3+1] 5106 $ j = {0,0,2,0,1,2} [size = 6]; values must be sorted for each row 5107 $ v = {1,2,3,4,5,6} [size = 6] 5108 5109 .seealso: `MatCreate()`, `MatCreateAIJ()`, `MatCreateSeqAIJ()`, `MatCreateMPIAIJWithArrays()`, `MatMPIAIJSetPreallocationCSR()` 5110 @*/ 5111 PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt i[], PetscInt j[], PetscScalar a[], Mat *mat) 5112 { 5113 PetscInt ii; 5114 Mat_SeqAIJ *aij; 5115 PetscInt jj; 5116 5117 PetscFunctionBegin; 5118 PetscCheck(m <= 0 || i[0] == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "i (row indices) must start with 0"); 5119 PetscCall(MatCreate(comm, mat)); 5120 PetscCall(MatSetSizes(*mat, m, n, m, n)); 5121 /* PetscCall(MatSetBlockSizes(*mat,,)); */ 5122 PetscCall(MatSetType(*mat, MATSEQAIJ)); 5123 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(*mat, MAT_SKIP_ALLOCATION, NULL)); 5124 aij = (Mat_SeqAIJ *)(*mat)->data; 5125 PetscCall(PetscMalloc1(m, &aij->imax)); 5126 PetscCall(PetscMalloc1(m, &aij->ilen)); 5127 5128 aij->i = i; 5129 aij->j = j; 5130 aij->a = a; 5131 aij->singlemalloc = PETSC_FALSE; 5132 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 5133 aij->free_a = PETSC_FALSE; 5134 aij->free_ij = PETSC_FALSE; 5135 5136 for (ii = 0, aij->nonzerorowcnt = 0, aij->rmax = 0; ii < m; ii++) { 5137 aij->ilen[ii] = aij->imax[ii] = i[ii + 1] - i[ii]; 5138 if (PetscDefined(USE_DEBUG)) { 5139 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]); 5140 for (jj = i[ii] + 1; jj < i[ii + 1]; jj++) { 5141 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); 5142 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); 5143 } 5144 } 5145 } 5146 if (PetscDefined(USE_DEBUG)) { 5147 for (ii = 0; ii < aij->i[m]; ii++) { 5148 PetscCheck(j[ii] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative column index at location = %" PetscInt_FMT " index = %" PetscInt_FMT, ii, j[ii]); 5149 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]); 5150 } 5151 } 5152 5153 PetscCall(MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY)); 5154 PetscCall(MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY)); 5155 PetscFunctionReturn(0); 5156 } 5157 5158 /*@ 5159 MatCreateSeqAIJFromTriple - Creates an sequential `MATSEQAIJ` matrix using matrix elements (in COO format) 5160 provided by the user. 5161 5162 Collective 5163 5164 Input Parameters: 5165 + comm - must be an MPI communicator of size 1 5166 . m - number of rows 5167 . n - number of columns 5168 . i - row indices 5169 . j - column indices 5170 . a - matrix values 5171 . nz - number of nonzeros 5172 - idx - if the i and j indices start with 1 use `PETSC_TRUE` otherwise use `PETSC_FALSE` 5173 5174 Output Parameter: 5175 . mat - the matrix 5176 5177 Level: intermediate 5178 5179 Example: 5180 For the following matrix, the input data expected is as shown (using 0 based indexing) 5181 .vb 5182 1 0 0 5183 2 0 3 5184 4 5 6 5185 5186 i = {0,1,1,2,2,2} 5187 j = {0,0,2,0,1,2} 5188 v = {1,2,3,4,5,6} 5189 .ve 5190 Notes: 5191 Instead of using this function, users should also consider `MatSetPreallocationCOO()` and `MatSetValuesCOO()`, which allow repeated or remote entries, 5192 and are particularly useful in iterative applications. 5193 5194 .seealso: `MatCreate()`, `MatCreateAIJ()`, `MatCreateSeqAIJ()`, `MatCreateSeqAIJWithArrays()`, `MatMPIAIJSetPreallocationCSR()`, `MatSetValuesCOO()`, `MatSetPreallocationCOO()` 5195 @*/ 5196 PetscErrorCode MatCreateSeqAIJFromTriple(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt i[], PetscInt j[], PetscScalar a[], Mat *mat, PetscInt nz, PetscBool idx) 5197 { 5198 PetscInt ii, *nnz, one = 1, row, col; 5199 5200 PetscFunctionBegin; 5201 PetscCall(PetscCalloc1(m, &nnz)); 5202 for (ii = 0; ii < nz; ii++) nnz[i[ii] - !!idx] += 1; 5203 PetscCall(MatCreate(comm, mat)); 5204 PetscCall(MatSetSizes(*mat, m, n, m, n)); 5205 PetscCall(MatSetType(*mat, MATSEQAIJ)); 5206 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(*mat, 0, nnz)); 5207 for (ii = 0; ii < nz; ii++) { 5208 if (idx) { 5209 row = i[ii] - 1; 5210 col = j[ii] - 1; 5211 } else { 5212 row = i[ii]; 5213 col = j[ii]; 5214 } 5215 PetscCall(MatSetValues(*mat, one, &row, one, &col, &a[ii], ADD_VALUES)); 5216 } 5217 PetscCall(MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY)); 5218 PetscCall(MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY)); 5219 PetscCall(PetscFree(nnz)); 5220 PetscFunctionReturn(0); 5221 } 5222 5223 PetscErrorCode MatSeqAIJInvalidateDiagonal(Mat A) 5224 { 5225 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 5226 5227 PetscFunctionBegin; 5228 a->idiagvalid = PETSC_FALSE; 5229 a->ibdiagvalid = PETSC_FALSE; 5230 5231 PetscCall(MatSeqAIJInvalidateDiagonal_Inode(A)); 5232 PetscFunctionReturn(0); 5233 } 5234 5235 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqAIJ(MPI_Comm comm, Mat inmat, PetscInt n, MatReuse scall, Mat *outmat) 5236 { 5237 PetscFunctionBegin; 5238 PetscCall(MatCreateMPIMatConcatenateSeqMat_MPIAIJ(comm, inmat, n, scall, outmat)); 5239 PetscFunctionReturn(0); 5240 } 5241 5242 /* 5243 Permute A into C's *local* index space using rowemb,colemb. 5244 The embedding are supposed to be injections and the above implies that the range of rowemb is a subset 5245 of [0,m), colemb is in [0,n). 5246 If pattern == DIFFERENT_NONZERO_PATTERN, C is preallocated according to A. 5247 */ 5248 PetscErrorCode MatSetSeqMat_SeqAIJ(Mat C, IS rowemb, IS colemb, MatStructure pattern, Mat B) 5249 { 5250 /* If making this function public, change the error returned in this function away from _PLIB. */ 5251 Mat_SeqAIJ *Baij; 5252 PetscBool seqaij; 5253 PetscInt m, n, *nz, i, j, count; 5254 PetscScalar v; 5255 const PetscInt *rowindices, *colindices; 5256 5257 PetscFunctionBegin; 5258 if (!B) PetscFunctionReturn(0); 5259 /* Check to make sure the target matrix (and embeddings) are compatible with C and each other. */ 5260 PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATSEQAIJ, &seqaij)); 5261 PetscCheck(seqaij, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Input matrix is of wrong type"); 5262 if (rowemb) { 5263 PetscCall(ISGetLocalSize(rowemb, &m)); 5264 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); 5265 } else { 5266 PetscCheck(C->rmap->n == B->rmap->n, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Input matrix is row-incompatible with the target matrix"); 5267 } 5268 if (colemb) { 5269 PetscCall(ISGetLocalSize(colemb, &n)); 5270 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); 5271 } else { 5272 PetscCheck(C->cmap->n == B->cmap->n, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Input matrix is col-incompatible with the target matrix"); 5273 } 5274 5275 Baij = (Mat_SeqAIJ *)(B->data); 5276 if (pattern == DIFFERENT_NONZERO_PATTERN) { 5277 PetscCall(PetscMalloc1(B->rmap->n, &nz)); 5278 for (i = 0; i < B->rmap->n; i++) nz[i] = Baij->i[i + 1] - Baij->i[i]; 5279 PetscCall(MatSeqAIJSetPreallocation(C, 0, nz)); 5280 PetscCall(PetscFree(nz)); 5281 } 5282 if (pattern == SUBSET_NONZERO_PATTERN) PetscCall(MatZeroEntries(C)); 5283 count = 0; 5284 rowindices = NULL; 5285 colindices = NULL; 5286 if (rowemb) PetscCall(ISGetIndices(rowemb, &rowindices)); 5287 if (colemb) PetscCall(ISGetIndices(colemb, &colindices)); 5288 for (i = 0; i < B->rmap->n; i++) { 5289 PetscInt row; 5290 row = i; 5291 if (rowindices) row = rowindices[i]; 5292 for (j = Baij->i[i]; j < Baij->i[i + 1]; j++) { 5293 PetscInt col; 5294 col = Baij->j[count]; 5295 if (colindices) col = colindices[col]; 5296 v = Baij->a[count]; 5297 PetscCall(MatSetValues(C, 1, &row, 1, &col, &v, INSERT_VALUES)); 5298 ++count; 5299 } 5300 } 5301 /* FIXME: set C's nonzerostate correctly. */ 5302 /* Assembly for C is necessary. */ 5303 C->preallocated = PETSC_TRUE; 5304 C->assembled = PETSC_TRUE; 5305 C->was_assembled = PETSC_FALSE; 5306 PetscFunctionReturn(0); 5307 } 5308 5309 PetscErrorCode MatEliminateZeros_SeqAIJ(Mat A) 5310 { 5311 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 5312 MatScalar *aa = a->a; 5313 PetscInt m = A->rmap->n, fshift = 0, fshift_prev = 0, i, k; 5314 PetscInt *ailen = a->ilen, *imax = a->imax, *ai = a->i, *aj = a->j, rmax = 0; 5315 5316 PetscFunctionBegin; 5317 PetscCheck(A->assembled, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cannot eliminate zeros for unassembled matrix"); 5318 if (m) rmax = ailen[0]; /* determine row with most nonzeros */ 5319 for (i = 1; i <= m; i++) { 5320 /* move each nonzero entry back by the amount of zero slots (fshift) before it*/ 5321 for (k = ai[i - 1]; k < ai[i]; k++) { 5322 if (aa[k] == 0 && aj[k] != i - 1) fshift++; 5323 else { 5324 if (aa[k] == 0 && aj[k] == i - 1) PetscCall(PetscInfo(A, "Keep the diagonal zero at row %" PetscInt_FMT "\n", i - 1)); 5325 aa[k - fshift] = aa[k]; 5326 aj[k - fshift] = aj[k]; 5327 } 5328 } 5329 ai[i - 1] -= fshift_prev; // safe to update ai[i-1] now since it will not be used in the next iteration 5330 fshift_prev = fshift; 5331 /* reset ilen and imax for each row */ 5332 ailen[i - 1] = imax[i - 1] = ai[i] - fshift - ai[i - 1]; 5333 a->nonzerorowcnt += ((ai[i] - fshift - ai[i - 1]) > 0); 5334 rmax = PetscMax(rmax, ailen[i - 1]); 5335 } 5336 if (m) { 5337 ai[m] -= fshift; 5338 a->nz = ai[m]; 5339 } 5340 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)); 5341 A->nonzerostate -= fshift; 5342 A->info.nz_unneeded += (PetscReal)fshift; 5343 a->rmax = rmax; 5344 if (a->inode.use && a->inode.checked) PetscCall(MatSeqAIJCheckInode(A)); 5345 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 5346 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 5347 PetscFunctionReturn(0); 5348 } 5349 5350 PetscFunctionList MatSeqAIJList = NULL; 5351 5352 /*@C 5353 MatSeqAIJSetType - Converts a `MATSEQAIJ` matrix to a subtype 5354 5355 Collective 5356 5357 Input Parameters: 5358 + mat - the matrix object 5359 - matype - matrix type 5360 5361 Options Database Key: 5362 . -mat_seqai_type <method> - for example seqaijcrl 5363 5364 Level: intermediate 5365 5366 .seealso: `PCSetType()`, `VecSetType()`, `MatCreate()`, `MatType`, `Mat` 5367 @*/ 5368 PetscErrorCode MatSeqAIJSetType(Mat mat, MatType matype) 5369 { 5370 PetscBool sametype; 5371 PetscErrorCode (*r)(Mat, MatType, MatReuse, Mat *); 5372 5373 PetscFunctionBegin; 5374 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 5375 PetscCall(PetscObjectTypeCompare((PetscObject)mat, matype, &sametype)); 5376 if (sametype) PetscFunctionReturn(0); 5377 5378 PetscCall(PetscFunctionListFind(MatSeqAIJList, matype, &r)); 5379 PetscCheck(r, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown Mat type given: %s", matype); 5380 PetscCall((*r)(mat, matype, MAT_INPLACE_MATRIX, &mat)); 5381 PetscFunctionReturn(0); 5382 } 5383 5384 /*@C 5385 MatSeqAIJRegister - - Adds a new sub-matrix type for sequential `MATSEQAIJ` matrices 5386 5387 Not Collective 5388 5389 Input Parameters: 5390 + name - name of a new user-defined matrix type, for example `MATSEQAIJCRL` 5391 - function - routine to convert to subtype 5392 5393 Notes: 5394 `MatSeqAIJRegister()` may be called multiple times to add several user-defined solvers. 5395 5396 Then, your matrix can be chosen with the procedural interface at runtime via the option 5397 $ -mat_seqaij_type my_mat 5398 5399 Level: advanced 5400 5401 .seealso: `MatSeqAIJRegisterAll()` 5402 @*/ 5403 PetscErrorCode MatSeqAIJRegister(const char sname[], PetscErrorCode (*function)(Mat, MatType, MatReuse, Mat *)) 5404 { 5405 PetscFunctionBegin; 5406 PetscCall(MatInitializePackage()); 5407 PetscCall(PetscFunctionListAdd(&MatSeqAIJList, sname, function)); 5408 PetscFunctionReturn(0); 5409 } 5410 5411 PetscBool MatSeqAIJRegisterAllCalled = PETSC_FALSE; 5412 5413 /*@C 5414 MatSeqAIJRegisterAll - Registers all of the matrix subtypes of `MATSSEQAIJ` 5415 5416 Not Collective 5417 5418 Level: advanced 5419 5420 .seealso: `MatRegisterAll()`, `MatSeqAIJRegister()` 5421 @*/ 5422 PetscErrorCode MatSeqAIJRegisterAll(void) 5423 { 5424 PetscFunctionBegin; 5425 if (MatSeqAIJRegisterAllCalled) PetscFunctionReturn(0); 5426 MatSeqAIJRegisterAllCalled = PETSC_TRUE; 5427 5428 PetscCall(MatSeqAIJRegister(MATSEQAIJCRL, MatConvert_SeqAIJ_SeqAIJCRL)); 5429 PetscCall(MatSeqAIJRegister(MATSEQAIJPERM, MatConvert_SeqAIJ_SeqAIJPERM)); 5430 PetscCall(MatSeqAIJRegister(MATSEQAIJSELL, MatConvert_SeqAIJ_SeqAIJSELL)); 5431 #if defined(PETSC_HAVE_MKL_SPARSE) 5432 PetscCall(MatSeqAIJRegister(MATSEQAIJMKL, MatConvert_SeqAIJ_SeqAIJMKL)); 5433 #endif 5434 #if defined(PETSC_HAVE_CUDA) 5435 PetscCall(MatSeqAIJRegister(MATSEQAIJCUSPARSE, MatConvert_SeqAIJ_SeqAIJCUSPARSE)); 5436 #endif 5437 #if defined(PETSC_HAVE_HIP) 5438 PetscCall(MatSeqAIJRegister(MATSEQAIJHIPSPARSE, MatConvert_SeqAIJ_SeqAIJHIPSPARSE)); 5439 #endif 5440 #if defined(PETSC_HAVE_KOKKOS_KERNELS) 5441 PetscCall(MatSeqAIJRegister(MATSEQAIJKOKKOS, MatConvert_SeqAIJ_SeqAIJKokkos)); 5442 #endif 5443 #if defined(PETSC_HAVE_VIENNACL) && defined(PETSC_HAVE_VIENNACL_NO_CUDA) 5444 PetscCall(MatSeqAIJRegister(MATMPIAIJVIENNACL, MatConvert_SeqAIJ_SeqAIJViennaCL)); 5445 #endif 5446 PetscFunctionReturn(0); 5447 } 5448 5449 /* 5450 Special version for direct calls from Fortran 5451 */ 5452 #include <petsc/private/fortranimpl.h> 5453 #if defined(PETSC_HAVE_FORTRAN_CAPS) 5454 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 5455 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 5456 #define matsetvaluesseqaij_ matsetvaluesseqaij 5457 #endif 5458 5459 /* Change these macros so can be used in void function */ 5460 5461 /* Change these macros so can be used in void function */ 5462 /* Identical to PetscCallVoid, except it assigns to *_ierr */ 5463 #undef PetscCall 5464 #define PetscCall(...) \ 5465 do { \ 5466 PetscErrorCode ierr_msv_mpiaij = __VA_ARGS__; \ 5467 if (PetscUnlikely(ierr_msv_mpiaij)) { \ 5468 *_ierr = PetscError(PETSC_COMM_SELF, __LINE__, PETSC_FUNCTION_NAME, __FILE__, ierr_msv_mpiaij, PETSC_ERROR_REPEAT, " "); \ 5469 return; \ 5470 } \ 5471 } while (0) 5472 5473 #undef SETERRQ 5474 #define SETERRQ(comm, ierr, ...) \ 5475 do { \ 5476 *_ierr = PetscError(comm, __LINE__, PETSC_FUNCTION_NAME, __FILE__, ierr, PETSC_ERROR_INITIAL, __VA_ARGS__); \ 5477 return; \ 5478 } while (0) 5479 5480 PETSC_EXTERN void matsetvaluesseqaij_(Mat *AA, PetscInt *mm, const PetscInt im[], PetscInt *nn, const PetscInt in[], const PetscScalar v[], InsertMode *isis, PetscErrorCode *_ierr) 5481 { 5482 Mat A = *AA; 5483 PetscInt m = *mm, n = *nn; 5484 InsertMode is = *isis; 5485 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 5486 PetscInt *rp, k, low, high, t, ii, row, nrow, i, col, l, rmax, N; 5487 PetscInt *imax, *ai, *ailen; 5488 PetscInt *aj, nonew = a->nonew, lastcol = -1; 5489 MatScalar *ap, value, *aa; 5490 PetscBool ignorezeroentries = a->ignorezeroentries; 5491 PetscBool roworiented = a->roworiented; 5492 5493 PetscFunctionBegin; 5494 MatCheckPreallocated(A, 1); 5495 imax = a->imax; 5496 ai = a->i; 5497 ailen = a->ilen; 5498 aj = a->j; 5499 aa = a->a; 5500 5501 for (k = 0; k < m; k++) { /* loop over added rows */ 5502 row = im[k]; 5503 if (row < 0) continue; 5504 PetscCheck(row < A->rmap->n, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Row too large"); 5505 rp = aj + ai[row]; 5506 ap = aa + ai[row]; 5507 rmax = imax[row]; 5508 nrow = ailen[row]; 5509 low = 0; 5510 high = nrow; 5511 for (l = 0; l < n; l++) { /* loop over added columns */ 5512 if (in[l] < 0) continue; 5513 PetscCheck(in[l] < A->cmap->n, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Column too large"); 5514 col = in[l]; 5515 if (roworiented) value = v[l + k * n]; 5516 else value = v[k + l * m]; 5517 5518 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 5519 5520 if (col <= lastcol) low = 0; 5521 else high = nrow; 5522 lastcol = col; 5523 while (high - low > 5) { 5524 t = (low + high) / 2; 5525 if (rp[t] > col) high = t; 5526 else low = t; 5527 } 5528 for (i = low; i < high; i++) { 5529 if (rp[i] > col) break; 5530 if (rp[i] == col) { 5531 if (is == ADD_VALUES) ap[i] += value; 5532 else ap[i] = value; 5533 goto noinsert; 5534 } 5535 } 5536 if (value == 0.0 && ignorezeroentries) goto noinsert; 5537 if (nonew == 1) goto noinsert; 5538 PetscCheck(nonew != -1, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Inserting a new nonzero in the matrix"); 5539 MatSeqXAIJReallocateAIJ(A, A->rmap->n, 1, nrow, row, col, rmax, aa, ai, aj, rp, ap, imax, nonew, MatScalar); 5540 N = nrow++ - 1; 5541 a->nz++; 5542 high++; 5543 /* shift up all the later entries in this row */ 5544 for (ii = N; ii >= i; ii--) { 5545 rp[ii + 1] = rp[ii]; 5546 ap[ii + 1] = ap[ii]; 5547 } 5548 rp[i] = col; 5549 ap[i] = value; 5550 A->nonzerostate++; 5551 noinsert:; 5552 low = i + 1; 5553 } 5554 ailen[row] = nrow; 5555 } 5556 PetscFunctionReturnVoid(); 5557 } 5558 /* Undefining these here since they were redefined from their original definition above! No 5559 * other PETSc functions should be defined past this point, as it is impossible to recover the 5560 * original definitions */ 5561 #undef PetscCall 5562 #undef SETERRQ 5563