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