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