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