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