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(PetscTableDestroy((PetscTable *)&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(PetscTableDestroy((PetscTable *)&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 PetscTable gid1_lid1; 3608 PetscTablePosition 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(PetscTableCreate(mat->rmap->n, mat->cmap->N + 1, &gid1_lid1)); 3617 ec = 0; 3618 for (i = 0; i < nz; i++) { 3619 PetscInt data, gid1 = jj[i] + 1; 3620 PetscCall(PetscTableFind(gid1_lid1, gid1, &data)); 3621 if (!data) { 3622 /* one based table */ 3623 PetscCall(PetscTableAdd(gid1_lid1, gid1, ++ec, INSERT_VALUES)); 3624 } 3625 } 3626 /* form array of columns we need */ 3627 PetscCall(PetscMalloc1(ec, &garray)); 3628 PetscCall(PetscTableGetHeadPosition(gid1_lid1, &tpos)); 3629 while (tpos) { 3630 PetscCall(PetscTableGetNext(gid1_lid1, &tpos, &gid, &lid)); 3631 gid--; 3632 lid--; 3633 garray[lid] = gid; 3634 } 3635 PetscCall(PetscSortInt(ec, garray)); /* sort, and rebuild */ 3636 PetscCall(PetscTableRemoveAll(gid1_lid1)); 3637 for (i = 0; i < ec; i++) PetscCall(PetscTableAdd(gid1_lid1, garray[i] + 1, i + 1, INSERT_VALUES)); 3638 /* compact out the extra columns in B */ 3639 for (i = 0; i < nz; i++) { 3640 PetscInt gid1 = jj[i] + 1; 3641 PetscCall(PetscTableFind(gid1_lid1, gid1, &lid)); 3642 lid--; 3643 jj[i] = lid; 3644 } 3645 PetscCall(PetscLayoutDestroy(&mat->cmap)); 3646 PetscCall(PetscTableDestroy(&gid1_lid1)); 3647 PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)mat), ec, ec, 1, &mat->cmap)); 3648 PetscCall(ISLocalToGlobalMappingCreate(PETSC_COMM_SELF, mat->cmap->bs, mat->cmap->n, garray, PETSC_OWN_POINTER, mapping)); 3649 PetscCall(ISLocalToGlobalMappingSetType(*mapping, ISLOCALTOGLOBALMAPPINGHASH)); 3650 PetscFunctionReturn(0); 3651 } 3652 3653 /*@ 3654 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 3655 in the matrix. 3656 3657 Input Parameters: 3658 + mat - the `MATSEQAIJ` matrix 3659 - indices - the column indices 3660 3661 Level: advanced 3662 3663 Notes: 3664 This can be called if you have precomputed the nonzero structure of the 3665 matrix and want to provide it to the matrix object to improve the performance 3666 of the `MatSetValues()` operation. 3667 3668 You MUST have set the correct numbers of nonzeros per row in the call to 3669 `MatCreateSeqAIJ()`, and the columns indices MUST be sorted. 3670 3671 MUST be called before any calls to `MatSetValues()` 3672 3673 The indices should start with zero, not one. 3674 3675 @*/ 3676 PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat, PetscInt *indices) 3677 { 3678 PetscFunctionBegin; 3679 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 3680 PetscValidIntPointer(indices, 2); 3681 PetscUseMethod(mat, "MatSeqAIJSetColumnIndices_C", (Mat, PetscInt *), (mat, indices)); 3682 PetscFunctionReturn(0); 3683 } 3684 3685 /* ----------------------------------------------------------------------------------------*/ 3686 3687 PetscErrorCode MatStoreValues_SeqAIJ(Mat mat) 3688 { 3689 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 3690 size_t nz = aij->i[mat->rmap->n]; 3691 3692 PetscFunctionBegin; 3693 PetscCheck(aij->nonew, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3694 3695 /* allocate space for values if not already there */ 3696 if (!aij->saved_values) { PetscCall(PetscMalloc1(nz + 1, &aij->saved_values)); } 3697 3698 /* copy values over */ 3699 PetscCall(PetscArraycpy(aij->saved_values, aij->a, nz)); 3700 PetscFunctionReturn(0); 3701 } 3702 3703 /*@ 3704 MatStoreValues - Stashes a copy of the matrix values; this allows, for 3705 example, reuse of the linear part of a Jacobian, while recomputing the 3706 nonlinear portion. 3707 3708 Collect on mat 3709 3710 Input Parameters: 3711 . mat - the matrix (currently only `MATAIJ` matrices support this option) 3712 3713 Level: advanced 3714 3715 Common Usage, with `SNESSolve()`: 3716 $ Create Jacobian matrix 3717 $ Set linear terms into matrix 3718 $ Apply boundary conditions to matrix, at this time matrix must have 3719 $ final nonzero structure (i.e. setting the nonlinear terms and applying 3720 $ boundary conditions again will not change the nonzero structure 3721 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3722 $ ierr = MatStoreValues(mat); 3723 $ Call SNESSetJacobian() with matrix 3724 $ In your Jacobian routine 3725 $ ierr = MatRetrieveValues(mat); 3726 $ Set nonlinear terms in matrix 3727 3728 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 3729 $ // build linear portion of Jacobian 3730 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3731 $ ierr = MatStoreValues(mat); 3732 $ loop over nonlinear iterations 3733 $ ierr = MatRetrieveValues(mat); 3734 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 3735 $ // call MatAssemblyBegin/End() on matrix 3736 $ Solve linear system with Jacobian 3737 $ endloop 3738 3739 Notes: 3740 Matrix must already be assemblied before calling this routine 3741 Must set the matrix option `MatSetOption`(mat,`MAT_NEW_NONZERO_LOCATIONS`,`PETSC_FALSE`); before 3742 calling this routine. 3743 3744 When this is called multiple times it overwrites the previous set of stored values 3745 and does not allocated additional space. 3746 3747 .seealso: `MatRetrieveValues()` 3748 @*/ 3749 PetscErrorCode MatStoreValues(Mat mat) 3750 { 3751 PetscFunctionBegin; 3752 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 3753 PetscCheck(mat->assembled, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix"); 3754 PetscCheck(!mat->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix"); 3755 PetscUseMethod(mat, "MatStoreValues_C", (Mat), (mat)); 3756 PetscFunctionReturn(0); 3757 } 3758 3759 PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat) 3760 { 3761 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 3762 PetscInt nz = aij->i[mat->rmap->n]; 3763 3764 PetscFunctionBegin; 3765 PetscCheck(aij->nonew, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3766 PetscCheck(aij->saved_values, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call MatStoreValues(A);first"); 3767 /* copy values over */ 3768 PetscCall(PetscArraycpy(aij->a, aij->saved_values, nz)); 3769 PetscFunctionReturn(0); 3770 } 3771 3772 /*@ 3773 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 3774 example, reuse of the linear part of a Jacobian, while recomputing the 3775 nonlinear portion. 3776 3777 Collect on mat 3778 3779 Input Parameters: 3780 . mat - the matrix (currently only `MATAIJ` matrices support this option) 3781 3782 Level: advanced 3783 3784 .seealso: `MatStoreValues()` 3785 @*/ 3786 PetscErrorCode MatRetrieveValues(Mat mat) 3787 { 3788 PetscFunctionBegin; 3789 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 3790 PetscCheck(mat->assembled, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix"); 3791 PetscCheck(!mat->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix"); 3792 PetscUseMethod(mat, "MatRetrieveValues_C", (Mat), (mat)); 3793 PetscFunctionReturn(0); 3794 } 3795 3796 /* --------------------------------------------------------------------------------*/ 3797 /*@C 3798 MatCreateSeqAIJ - Creates a sparse matrix in `MATSEQAIJ` (compressed row) format 3799 (the default parallel PETSc format). For good matrix assembly performance 3800 the user should preallocate the matrix storage by setting the parameter nz 3801 (or the array nnz). By setting these parameters accurately, performance 3802 during matrix assembly can be increased by more than a factor of 50. 3803 3804 Collective 3805 3806 Input Parameters: 3807 + comm - MPI communicator, set to `PETSC_COMM_SELF` 3808 . m - number of rows 3809 . n - number of columns 3810 . nz - number of nonzeros per row (same for all rows) 3811 - nnz - array containing the number of nonzeros in the various rows 3812 (possibly different for each row) or NULL 3813 3814 Output Parameter: 3815 . A - the matrix 3816 3817 It is recommended that one use the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`, 3818 MatXXXXSetPreallocation() paradigm instead of this routine directly. 3819 [MatXXXXSetPreallocation() is, for example, `MatSeqAIJSetPreallocation()`] 3820 3821 Notes: 3822 If nnz is given then nz is ignored 3823 3824 The AIJ format, also called 3825 compressed row storage, is fully compatible with standard Fortran 77 3826 storage. That is, the stored row and column indices can begin at 3827 either one (as in Fortran) or zero. See the users' manual for details. 3828 3829 Specify the preallocated storage with either nz or nnz (not both). 3830 Set nz = `PETSC_DEFAULT` and nnz = NULL for PETSc to control dynamic memory 3831 allocation. For large problems you MUST preallocate memory or you 3832 will get TERRIBLE performance, see the users' manual chapter on matrices. 3833 3834 By default, this format uses inodes (identical nodes) when possible, to 3835 improve numerical efficiency of matrix-vector products and solves. We 3836 search for consecutive rows with the same nonzero structure, thereby 3837 reusing matrix information to achieve increased efficiency. 3838 3839 Options Database Keys: 3840 + -mat_no_inode - Do not use inodes 3841 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3842 3843 Level: intermediate 3844 3845 .seealso: [Sparse Matrix Creation](sec_matsparse), `MatCreate()`, `MatCreateAIJ()`, `MatSetValues()`, `MatSeqAIJSetColumnIndices()`, `MatCreateSeqAIJWithArrays()` 3846 @*/ 3847 PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt nz, const PetscInt nnz[], Mat *A) 3848 { 3849 PetscFunctionBegin; 3850 PetscCall(MatCreate(comm, A)); 3851 PetscCall(MatSetSizes(*A, m, n, m, n)); 3852 PetscCall(MatSetType(*A, MATSEQAIJ)); 3853 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(*A, nz, nnz)); 3854 PetscFunctionReturn(0); 3855 } 3856 3857 /*@C 3858 MatSeqAIJSetPreallocation - For good matrix assembly performance 3859 the user should preallocate the matrix storage by setting the parameter nz 3860 (or the array nnz). By setting these parameters accurately, performance 3861 during matrix assembly can be increased by more than a factor of 50. 3862 3863 Collective 3864 3865 Input Parameters: 3866 + B - The matrix 3867 . nz - number of nonzeros per row (same for all rows) 3868 - nnz - array containing the number of nonzeros in the various rows 3869 (possibly different for each row) or NULL 3870 3871 Notes: 3872 If nnz is given then nz is ignored 3873 3874 The `MATSEQAIJ` format also called 3875 compressed row storage, is fully compatible with standard Fortran 77 3876 storage. That is, the stored row and column indices can begin at 3877 either one (as in Fortran) or zero. See the users' manual for details. 3878 3879 Specify the preallocated storage with either nz or nnz (not both). 3880 Set nz = `PETSC_DEFAULT` and nnz = NULL for PETSc to control dynamic memory 3881 allocation. For large problems you MUST preallocate memory or you 3882 will get TERRIBLE performance, see the users' manual chapter on matrices. 3883 3884 You can call `MatGetInfo()` to get information on how effective the preallocation was; 3885 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 3886 You can also run with the option -info and look for messages with the string 3887 malloc in them to see if additional memory allocation was needed. 3888 3889 Developer Notes: 3890 Use nz of `MAT_SKIP_ALLOCATION` to not allocate any space for the matrix 3891 entries or columns indices 3892 3893 By default, this format uses inodes (identical nodes) when possible, to 3894 improve numerical efficiency of matrix-vector products and solves. We 3895 search for consecutive rows with the same nonzero structure, thereby 3896 reusing matrix information to achieve increased efficiency. 3897 3898 Options Database Keys: 3899 + -mat_no_inode - Do not use inodes 3900 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3901 3902 Level: intermediate 3903 3904 .seealso: `MatCreate()`, `MatCreateAIJ()`, `MatSetValues()`, `MatSeqAIJSetColumnIndices()`, `MatCreateSeqAIJWithArrays()`, `MatGetInfo()`, 3905 `MatSeqAIJSetTotalPreallocation()` 3906 @*/ 3907 PetscErrorCode MatSeqAIJSetPreallocation(Mat B, PetscInt nz, const PetscInt nnz[]) 3908 { 3909 PetscFunctionBegin; 3910 PetscValidHeaderSpecific(B, MAT_CLASSID, 1); 3911 PetscValidType(B, 1); 3912 PetscTryMethod(B, "MatSeqAIJSetPreallocation_C", (Mat, PetscInt, const PetscInt[]), (B, nz, nnz)); 3913 PetscFunctionReturn(0); 3914 } 3915 3916 PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B, PetscInt nz, const PetscInt *nnz) 3917 { 3918 Mat_SeqAIJ *b; 3919 PetscBool skipallocation = PETSC_FALSE, realalloc = PETSC_FALSE; 3920 PetscInt i; 3921 3922 PetscFunctionBegin; 3923 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 3924 if (nz == MAT_SKIP_ALLOCATION) { 3925 skipallocation = PETSC_TRUE; 3926 nz = 0; 3927 } 3928 PetscCall(PetscLayoutSetUp(B->rmap)); 3929 PetscCall(PetscLayoutSetUp(B->cmap)); 3930 3931 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3932 PetscCheck(nz >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "nz cannot be less than 0: value %" PetscInt_FMT, nz); 3933 if (PetscUnlikelyDebug(nnz)) { 3934 for (i = 0; i < B->rmap->n; i++) { 3935 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]); 3936 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); 3937 } 3938 } 3939 3940 B->preallocated = PETSC_TRUE; 3941 3942 b = (Mat_SeqAIJ *)B->data; 3943 3944 if (!skipallocation) { 3945 if (!b->imax) { PetscCall(PetscMalloc1(B->rmap->n, &b->imax)); } 3946 if (!b->ilen) { 3947 /* b->ilen will count nonzeros in each row so far. */ 3948 PetscCall(PetscCalloc1(B->rmap->n, &b->ilen)); 3949 } else { 3950 PetscCall(PetscMemzero(b->ilen, B->rmap->n * sizeof(PetscInt))); 3951 } 3952 if (!b->ipre) { PetscCall(PetscMalloc1(B->rmap->n, &b->ipre)); } 3953 if (!nnz) { 3954 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 3955 else if (nz < 0) nz = 1; 3956 nz = PetscMin(nz, B->cmap->n); 3957 for (i = 0; i < B->rmap->n; i++) b->imax[i] = nz; 3958 nz = nz * B->rmap->n; 3959 } else { 3960 PetscInt64 nz64 = 0; 3961 for (i = 0; i < B->rmap->n; i++) { 3962 b->imax[i] = nnz[i]; 3963 nz64 += nnz[i]; 3964 } 3965 PetscCall(PetscIntCast(nz64, &nz)); 3966 } 3967 3968 /* allocate the matrix space */ 3969 /* FIXME: should B's old memory be unlogged? */ 3970 PetscCall(MatSeqXAIJFreeAIJ(B, &b->a, &b->j, &b->i)); 3971 if (B->structure_only) { 3972 PetscCall(PetscMalloc1(nz, &b->j)); 3973 PetscCall(PetscMalloc1(B->rmap->n + 1, &b->i)); 3974 } else { 3975 PetscCall(PetscMalloc3(nz, &b->a, nz, &b->j, B->rmap->n + 1, &b->i)); 3976 } 3977 b->i[0] = 0; 3978 for (i = 1; i < B->rmap->n + 1; i++) b->i[i] = b->i[i - 1] + b->imax[i - 1]; 3979 if (B->structure_only) { 3980 b->singlemalloc = PETSC_FALSE; 3981 b->free_a = PETSC_FALSE; 3982 } else { 3983 b->singlemalloc = PETSC_TRUE; 3984 b->free_a = PETSC_TRUE; 3985 } 3986 b->free_ij = PETSC_TRUE; 3987 } else { 3988 b->free_a = PETSC_FALSE; 3989 b->free_ij = PETSC_FALSE; 3990 } 3991 3992 if (b->ipre && nnz != b->ipre && b->imax) { 3993 /* reserve user-requested sparsity */ 3994 PetscCall(PetscArraycpy(b->ipre, b->imax, B->rmap->n)); 3995 } 3996 3997 b->nz = 0; 3998 b->maxnz = nz; 3999 B->info.nz_unneeded = (double)b->maxnz; 4000 if (realalloc) PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE)); 4001 B->was_assembled = PETSC_FALSE; 4002 B->assembled = PETSC_FALSE; 4003 /* We simply deem preallocation has changed nonzero state. Updating the state 4004 will give clients (like AIJKokkos) a chance to know something has happened. 4005 */ 4006 B->nonzerostate++; 4007 PetscFunctionReturn(0); 4008 } 4009 4010 PetscErrorCode MatResetPreallocation_SeqAIJ(Mat A) 4011 { 4012 Mat_SeqAIJ *a; 4013 PetscInt i; 4014 4015 PetscFunctionBegin; 4016 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 4017 4018 /* Check local size. If zero, then return */ 4019 if (!A->rmap->n) PetscFunctionReturn(0); 4020 4021 a = (Mat_SeqAIJ *)A->data; 4022 /* if no saved info, we error out */ 4023 PetscCheck(a->ipre, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "No saved preallocation info "); 4024 4025 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 "); 4026 4027 PetscCall(PetscArraycpy(a->imax, a->ipre, A->rmap->n)); 4028 PetscCall(PetscArrayzero(a->ilen, A->rmap->n)); 4029 a->i[0] = 0; 4030 for (i = 1; i < A->rmap->n + 1; i++) a->i[i] = a->i[i - 1] + a->imax[i - 1]; 4031 A->preallocated = PETSC_TRUE; 4032 a->nz = 0; 4033 a->maxnz = a->i[A->rmap->n]; 4034 A->info.nz_unneeded = (double)a->maxnz; 4035 A->was_assembled = PETSC_FALSE; 4036 A->assembled = PETSC_FALSE; 4037 PetscFunctionReturn(0); 4038 } 4039 4040 /*@ 4041 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in `MATSEQAIJ` format. 4042 4043 Input Parameters: 4044 + B - the matrix 4045 . i - the indices into j for the start of each row (starts with zero) 4046 . j - the column indices for each row (starts with zero) these must be sorted for each row 4047 - v - optional values in the matrix 4048 4049 Level: developer 4050 4051 Notes: 4052 The i,j,v values are COPIED with this routine; to avoid the copy use `MatCreateSeqAIJWithArrays()` 4053 4054 This routine may be called multiple times with different nonzero patterns (or the same nonzero pattern). The nonzero 4055 structure will be the union of all the previous nonzero structures. 4056 4057 Developer Notes: 4058 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 4059 then just copies the v values directly with `PetscMemcpy()`. 4060 4061 This routine could also take a `PetscCopyMode` argument to allow sharing the values instead of always copying them. 4062 4063 .seealso: `MatCreate()`, `MatCreateSeqAIJ()`, `MatSetValues()`, `MatSeqAIJSetPreallocation()`, `MatCreateSeqAIJ()`, `MATSEQAIJ`, `MatResetPreallocation()` 4064 @*/ 4065 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B, const PetscInt i[], const PetscInt j[], const PetscScalar v[]) 4066 { 4067 PetscFunctionBegin; 4068 PetscValidHeaderSpecific(B, MAT_CLASSID, 1); 4069 PetscValidType(B, 1); 4070 PetscTryMethod(B, "MatSeqAIJSetPreallocationCSR_C", (Mat, const PetscInt[], const PetscInt[], const PetscScalar[]), (B, i, j, v)); 4071 PetscFunctionReturn(0); 4072 } 4073 4074 PetscErrorCode MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B, const PetscInt Ii[], const PetscInt J[], const PetscScalar v[]) 4075 { 4076 PetscInt i; 4077 PetscInt m, n; 4078 PetscInt nz; 4079 PetscInt *nnz; 4080 4081 PetscFunctionBegin; 4082 PetscCheck(Ii[0] == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %" PetscInt_FMT, Ii[0]); 4083 4084 PetscCall(PetscLayoutSetUp(B->rmap)); 4085 PetscCall(PetscLayoutSetUp(B->cmap)); 4086 4087 PetscCall(MatGetSize(B, &m, &n)); 4088 PetscCall(PetscMalloc1(m + 1, &nnz)); 4089 for (i = 0; i < m; i++) { 4090 nz = Ii[i + 1] - Ii[i]; 4091 PetscCheck(nz >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local row %" PetscInt_FMT " has a negative number of columns %" PetscInt_FMT, i, nz); 4092 nnz[i] = nz; 4093 } 4094 PetscCall(MatSeqAIJSetPreallocation(B, 0, nnz)); 4095 PetscCall(PetscFree(nnz)); 4096 4097 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)); 4098 4099 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 4100 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 4101 4102 PetscCall(MatSetOption(B, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE)); 4103 PetscFunctionReturn(0); 4104 } 4105 4106 /*@ 4107 MatSeqAIJKron - Computes C, the Kronecker product of A and B. 4108 4109 Input Parameters: 4110 + A - left-hand side matrix 4111 . B - right-hand side matrix 4112 - reuse - either `MAT_INITIAL_MATRIX` or `MAT_REUSE_MATRIX` 4113 4114 Output Parameter: 4115 . C - Kronecker product of A and B 4116 4117 Level: intermediate 4118 4119 Note: 4120 `MAT_REUSE_MATRIX` can only be used when the nonzero structure of the product matrix has not changed from that last call to `MatSeqAIJKron()`. 4121 4122 .seealso: `MatCreateSeqAIJ()`, `MATSEQAIJ`, `MATKAIJ`, `MatReuse` 4123 @*/ 4124 PetscErrorCode MatSeqAIJKron(Mat A, Mat B, MatReuse reuse, Mat *C) 4125 { 4126 PetscFunctionBegin; 4127 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 4128 PetscValidType(A, 1); 4129 PetscValidHeaderSpecific(B, MAT_CLASSID, 2); 4130 PetscValidType(B, 2); 4131 PetscValidPointer(C, 4); 4132 if (reuse == MAT_REUSE_MATRIX) { 4133 PetscValidHeaderSpecific(*C, MAT_CLASSID, 4); 4134 PetscValidType(*C, 4); 4135 } 4136 PetscTryMethod(A, "MatSeqAIJKron_C", (Mat, Mat, MatReuse, Mat *), (A, B, reuse, C)); 4137 PetscFunctionReturn(0); 4138 } 4139 4140 PetscErrorCode MatSeqAIJKron_SeqAIJ(Mat A, Mat B, MatReuse reuse, Mat *C) 4141 { 4142 Mat newmat; 4143 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 4144 Mat_SeqAIJ *b = (Mat_SeqAIJ *)B->data; 4145 PetscScalar *v; 4146 const PetscScalar *aa, *ba; 4147 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; 4148 PetscBool flg; 4149 4150 PetscFunctionBegin; 4151 PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix"); 4152 PetscCheck(A->assembled, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix"); 4153 PetscCheck(!B->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix"); 4154 PetscCheck(B->assembled, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix"); 4155 PetscCall(PetscObjectTypeCompare((PetscObject)B, MATSEQAIJ, &flg)); 4156 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "MatType %s", ((PetscObject)B)->type_name); 4157 PetscCheck(reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX, PETSC_COMM_SELF, PETSC_ERR_SUP, "MatReuse %d", (int)reuse); 4158 if (reuse == MAT_INITIAL_MATRIX) { 4159 PetscCall(PetscMalloc2(am * bm + 1, &i, a->i[am] * b->i[bm], &j)); 4160 PetscCall(MatCreate(PETSC_COMM_SELF, &newmat)); 4161 PetscCall(MatSetSizes(newmat, am * bm, an * bn, am * bm, an * bn)); 4162 PetscCall(MatSetType(newmat, MATAIJ)); 4163 i[0] = 0; 4164 for (m = 0; m < am; ++m) { 4165 for (p = 0; p < bm; ++p) { 4166 i[m * bm + p + 1] = i[m * bm + p] + (a->i[m + 1] - a->i[m]) * (b->i[p + 1] - b->i[p]); 4167 for (n = a->i[m]; n < a->i[m + 1]; ++n) { 4168 for (q = b->i[p]; q < b->i[p + 1]; ++q) j[nnz++] = a->j[n] * bn + b->j[q]; 4169 } 4170 } 4171 } 4172 PetscCall(MatSeqAIJSetPreallocationCSR(newmat, i, j, NULL)); 4173 *C = newmat; 4174 PetscCall(PetscFree2(i, j)); 4175 nnz = 0; 4176 } 4177 PetscCall(MatSeqAIJGetArray(*C, &v)); 4178 PetscCall(MatSeqAIJGetArrayRead(A, &aa)); 4179 PetscCall(MatSeqAIJGetArrayRead(B, &ba)); 4180 for (m = 0; m < am; ++m) { 4181 for (p = 0; p < bm; ++p) { 4182 for (n = a->i[m]; n < a->i[m + 1]; ++n) { 4183 for (q = b->i[p]; q < b->i[p + 1]; ++q) v[nnz++] = aa[n] * ba[q]; 4184 } 4185 } 4186 } 4187 PetscCall(MatSeqAIJRestoreArray(*C, &v)); 4188 PetscCall(MatSeqAIJRestoreArrayRead(A, &aa)); 4189 PetscCall(MatSeqAIJRestoreArrayRead(B, &ba)); 4190 PetscFunctionReturn(0); 4191 } 4192 4193 #include <../src/mat/impls/dense/seq/dense.h> 4194 #include <petsc/private/kernels/petscaxpy.h> 4195 4196 /* 4197 Computes (B'*A')' since computing B*A directly is untenable 4198 4199 n p p 4200 [ ] [ ] [ ] 4201 m [ A ] * n [ B ] = m [ C ] 4202 [ ] [ ] [ ] 4203 4204 */ 4205 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A, Mat B, Mat C) 4206 { 4207 Mat_SeqDense *sub_a = (Mat_SeqDense *)A->data; 4208 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ *)B->data; 4209 Mat_SeqDense *sub_c = (Mat_SeqDense *)C->data; 4210 PetscInt i, j, n, m, q, p; 4211 const PetscInt *ii, *idx; 4212 const PetscScalar *b, *a, *a_q; 4213 PetscScalar *c, *c_q; 4214 PetscInt clda = sub_c->lda; 4215 PetscInt alda = sub_a->lda; 4216 4217 PetscFunctionBegin; 4218 m = A->rmap->n; 4219 n = A->cmap->n; 4220 p = B->cmap->n; 4221 a = sub_a->v; 4222 b = sub_b->a; 4223 c = sub_c->v; 4224 if (clda == m) { 4225 PetscCall(PetscArrayzero(c, m * p)); 4226 } else { 4227 for (j = 0; j < p; j++) 4228 for (i = 0; i < m; i++) c[j * clda + i] = 0.0; 4229 } 4230 ii = sub_b->i; 4231 idx = sub_b->j; 4232 for (i = 0; i < n; i++) { 4233 q = ii[i + 1] - ii[i]; 4234 while (q-- > 0) { 4235 c_q = c + clda * (*idx); 4236 a_q = a + alda * i; 4237 PetscKernelAXPY(c_q, *b, a_q, m); 4238 idx++; 4239 b++; 4240 } 4241 } 4242 PetscFunctionReturn(0); 4243 } 4244 4245 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A, Mat B, PetscReal fill, Mat C) 4246 { 4247 PetscInt m = A->rmap->n, n = B->cmap->n; 4248 PetscBool cisdense; 4249 4250 PetscFunctionBegin; 4251 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); 4252 PetscCall(MatSetSizes(C, m, n, m, n)); 4253 PetscCall(MatSetBlockSizesFromMats(C, A, B)); 4254 PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSECUDA, "")); 4255 if (!cisdense) PetscCall(MatSetType(C, MATDENSE)); 4256 PetscCall(MatSetUp(C)); 4257 4258 C->ops->matmultnumeric = MatMatMultNumeric_SeqDense_SeqAIJ; 4259 PetscFunctionReturn(0); 4260 } 4261 4262 /* ----------------------------------------------------------------*/ 4263 /*MC 4264 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 4265 based on compressed sparse row format. 4266 4267 Options Database Keys: 4268 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 4269 4270 Level: beginner 4271 4272 Notes: 4273 `MatSetValues()` may be called for this matrix type with a NULL argument for the numerical values, 4274 in this case the values associated with the rows and columns one passes in are set to zero 4275 in the matrix 4276 4277 `MatSetOptions`(,`MAT_STRUCTURE_ONLY`,`PETSC_TRUE`) may be called for this matrix type. In this no 4278 space is allocated for the nonzero entries and any entries passed with `MatSetValues()` are ignored 4279 4280 Developer Note: 4281 It would be nice if all matrix formats supported passing NULL in for the numerical values 4282 4283 .seealso: `MatCreateSeqAIJ()`, `MatSetFromOptions()`, `MatSetType()`, `MatCreate()`, `MatType`, `MATSELL`, `MATSEQSELL`, `MATMPISELL` 4284 M*/ 4285 4286 /*MC 4287 MATAIJ - MATAIJ = "aij" - A matrix type to be used for sparse matrices. 4288 4289 This matrix type is identical to `MATSEQAIJ` when constructed with a single process communicator, 4290 and `MATMPIAIJ` otherwise. As a result, for single process communicators, 4291 `MatSeqAIJSetPreallocation()` is supported, and similarly `MatMPIAIJSetPreallocation()` is supported 4292 for communicators controlling multiple processes. It is recommended that you call both of 4293 the above preallocation routines for simplicity. 4294 4295 Options Database Keys: 4296 . -mat_type aij - sets the matrix type to "aij" during a call to `MatSetFromOptions()` 4297 4298 Note: 4299 Subclasses include `MATAIJCUSPARSE`, `MATAIJPERM`, `MATAIJSELL`, `MATAIJMKL`, `MATAIJCRL`, and also automatically switches over to use inodes when 4300 enough exist. 4301 4302 Level: beginner 4303 4304 .seealso: `MatCreateAIJ()`, `MatCreateSeqAIJ()`, `MATSEQAIJ`, `MATMPIAIJ`, `MATSELL`, `MATSEQSELL`, `MATMPISELL` 4305 M*/ 4306 4307 /*MC 4308 MATAIJCRL - MATAIJCRL = "aijcrl" - A matrix type to be used for sparse matrices. 4309 4310 This matrix type is identical to `MATSEQAIJCRL` when constructed with a single process communicator, 4311 and `MATMPIAIJCRL` otherwise. As a result, for single process communicators, 4312 `MatSeqAIJSetPreallocation()` is supported, and similarly `MatMPIAIJSetPreallocation()` is supported 4313 for communicators controlling multiple processes. It is recommended that you call both of 4314 the above preallocation routines for simplicity. 4315 4316 Options Database Keys: 4317 . -mat_type aijcrl - sets the matrix type to "aijcrl" during a call to `MatSetFromOptions()` 4318 4319 Level: beginner 4320 4321 .seealso: `MatCreateMPIAIJCRL`, `MATSEQAIJCRL`, `MATMPIAIJCRL`, `MATSEQAIJCRL`, `MATMPIAIJCRL` 4322 M*/ 4323 4324 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCRL(Mat, MatType, MatReuse, Mat *); 4325 #if defined(PETSC_HAVE_ELEMENTAL) 4326 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_Elemental(Mat, MatType, MatReuse, Mat *); 4327 #endif 4328 #if defined(PETSC_HAVE_SCALAPACK) 4329 PETSC_INTERN PetscErrorCode MatConvert_AIJ_ScaLAPACK(Mat, MatType, MatReuse, Mat *); 4330 #endif 4331 #if defined(PETSC_HAVE_HYPRE) 4332 PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A, MatType, MatReuse, Mat *); 4333 #endif 4334 4335 PETSC_EXTERN PetscErrorCode MatConvert_SeqAIJ_SeqSELL(Mat, MatType, MatReuse, Mat *); 4336 PETSC_INTERN PetscErrorCode MatConvert_XAIJ_IS(Mat, MatType, MatReuse, Mat *); 4337 PETSC_INTERN PetscErrorCode MatProductSetFromOptions_IS_XAIJ(Mat); 4338 4339 /*@C 4340 MatSeqAIJGetArray - gives read/write access to the array where the data for a `MATSEQAIJ` matrix is stored 4341 4342 Not Collective 4343 4344 Input Parameter: 4345 . mat - a `MATSEQAIJ` matrix 4346 4347 Output Parameter: 4348 . array - pointer to the data 4349 4350 Level: intermediate 4351 4352 .seealso: `MatSeqAIJRestoreArray()`, `MatSeqAIJGetArrayF90()` 4353 @*/ 4354 PetscErrorCode MatSeqAIJGetArray(Mat A, PetscScalar **array) 4355 { 4356 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)A->data; 4357 4358 PetscFunctionBegin; 4359 if (aij->ops->getarray) { 4360 PetscCall((*aij->ops->getarray)(A, array)); 4361 } else { 4362 *array = aij->a; 4363 } 4364 PetscFunctionReturn(0); 4365 } 4366 4367 /*@C 4368 MatSeqAIJRestoreArray - returns access to the array where the data for a `MATSEQAIJ` matrix is stored obtained by `MatSeqAIJGetArray()` 4369 4370 Not Collective 4371 4372 Input Parameters: 4373 + mat - a `MATSEQAIJ` matrix 4374 - array - pointer to the data 4375 4376 Level: intermediate 4377 4378 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArrayF90()` 4379 @*/ 4380 PetscErrorCode MatSeqAIJRestoreArray(Mat A, PetscScalar **array) 4381 { 4382 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)A->data; 4383 4384 PetscFunctionBegin; 4385 if (aij->ops->restorearray) { 4386 PetscCall((*aij->ops->restorearray)(A, array)); 4387 } else { 4388 *array = NULL; 4389 } 4390 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 4391 PetscCall(PetscObjectStateIncrease((PetscObject)A)); 4392 PetscFunctionReturn(0); 4393 } 4394 4395 /*@C 4396 MatSeqAIJGetArrayRead - gives read-only access to the array where the data for a `MATSEQAIJ` matrix is stored 4397 4398 Not Collective 4399 4400 Input Parameter: 4401 . mat - a `MATSEQAIJ` matrix 4402 4403 Output Parameter: 4404 . array - pointer to the data 4405 4406 Level: intermediate 4407 4408 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArrayRead()` 4409 @*/ 4410 PetscErrorCode MatSeqAIJGetArrayRead(Mat A, const PetscScalar **array) 4411 { 4412 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)A->data; 4413 4414 PetscFunctionBegin; 4415 if (aij->ops->getarrayread) { 4416 PetscCall((*aij->ops->getarrayread)(A, array)); 4417 } else { 4418 *array = aij->a; 4419 } 4420 PetscFunctionReturn(0); 4421 } 4422 4423 /*@C 4424 MatSeqAIJRestoreArrayRead - restore the read-only access array obtained from `MatSeqAIJGetArrayRead()` 4425 4426 Not Collective 4427 4428 Input Parameter: 4429 . mat - a `MATSEQAIJ` matrix 4430 4431 Output Parameter: 4432 . array - pointer to the data 4433 4434 Level: intermediate 4435 4436 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJGetArrayRead()` 4437 @*/ 4438 PetscErrorCode MatSeqAIJRestoreArrayRead(Mat A, const PetscScalar **array) 4439 { 4440 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)A->data; 4441 4442 PetscFunctionBegin; 4443 if (aij->ops->restorearrayread) { 4444 PetscCall((*aij->ops->restorearrayread)(A, array)); 4445 } else { 4446 *array = NULL; 4447 } 4448 PetscFunctionReturn(0); 4449 } 4450 4451 /*@C 4452 MatSeqAIJGetArrayWrite - gives write-only access to the array where the data for a `MATSEQAIJ` matrix is stored 4453 4454 Not Collective 4455 4456 Input Parameter: 4457 . mat - a `MATSEQAIJ` matrix 4458 4459 Output Parameter: 4460 . array - pointer to the data 4461 4462 Level: intermediate 4463 4464 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArrayRead()` 4465 @*/ 4466 PetscErrorCode MatSeqAIJGetArrayWrite(Mat A, PetscScalar **array) 4467 { 4468 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)A->data; 4469 4470 PetscFunctionBegin; 4471 if (aij->ops->getarraywrite) { 4472 PetscCall((*aij->ops->getarraywrite)(A, array)); 4473 } else { 4474 *array = aij->a; 4475 } 4476 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 4477 PetscCall(PetscObjectStateIncrease((PetscObject)A)); 4478 PetscFunctionReturn(0); 4479 } 4480 4481 /*@C 4482 MatSeqAIJRestoreArrayWrite - restore the read-only access array obtained from MatSeqAIJGetArrayRead 4483 4484 Not Collective 4485 4486 Input Parameter: 4487 . mat - a MATSEQAIJ matrix 4488 4489 Output Parameter: 4490 . array - pointer to the data 4491 4492 Level: intermediate 4493 4494 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJGetArrayRead()` 4495 @*/ 4496 PetscErrorCode MatSeqAIJRestoreArrayWrite(Mat A, PetscScalar **array) 4497 { 4498 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)A->data; 4499 4500 PetscFunctionBegin; 4501 if (aij->ops->restorearraywrite) { 4502 PetscCall((*aij->ops->restorearraywrite)(A, array)); 4503 } else { 4504 *array = NULL; 4505 } 4506 PetscFunctionReturn(0); 4507 } 4508 4509 /*@C 4510 MatSeqAIJGetCSRAndMemType - Get the CSR arrays and the memory type of the `MATSEQAIJ` matrix 4511 4512 Not Collective 4513 4514 Input Parameter: 4515 . mat - a matrix of type `MATSEQAIJ` or its subclasses 4516 4517 Output Parameters: 4518 + i - row map array of the matrix 4519 . j - column index array of the matrix 4520 . a - data array of the matrix 4521 - memtype - memory type of the arrays 4522 4523 Notes: 4524 Any of the output parameters can be NULL, in which case the corresponding value is not returned. 4525 If mat is a device matrix, the arrays are on the device. Otherwise, they are on the host. 4526 4527 One can call this routine on a preallocated but not assembled matrix to just get the memory of the CSR underneath the matrix. 4528 If the matrix is assembled, the data array 'a' is guaranteed to have the latest values of the matrix. 4529 4530 Level: Developer 4531 4532 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJGetArrayRead()` 4533 @*/ 4534 PetscErrorCode MatSeqAIJGetCSRAndMemType(Mat mat, const PetscInt **i, const PetscInt **j, PetscScalar **a, PetscMemType *mtype) 4535 { 4536 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 4537 4538 PetscFunctionBegin; 4539 PetscCheck(mat->preallocated, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "matrix is not preallocated"); 4540 if (aij->ops->getcsrandmemtype) { 4541 PetscCall((*aij->ops->getcsrandmemtype)(mat, i, j, a, mtype)); 4542 } else { 4543 if (i) *i = aij->i; 4544 if (j) *j = aij->j; 4545 if (a) *a = aij->a; 4546 if (mtype) *mtype = PETSC_MEMTYPE_HOST; 4547 } 4548 PetscFunctionReturn(0); 4549 } 4550 4551 /*@C 4552 MatSeqAIJGetMaxRowNonzeros - returns the maximum number of nonzeros in any row 4553 4554 Not Collective 4555 4556 Input Parameter: 4557 . mat - a `MATSEQAIJ` matrix 4558 4559 Output Parameter: 4560 . nz - the maximum number of nonzeros in any row 4561 4562 Level: intermediate 4563 4564 .seealso: `MatSeqAIJRestoreArray()`, `MatSeqAIJGetArrayF90()` 4565 @*/ 4566 PetscErrorCode MatSeqAIJGetMaxRowNonzeros(Mat A, PetscInt *nz) 4567 { 4568 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)A->data; 4569 4570 PetscFunctionBegin; 4571 *nz = aij->rmax; 4572 PetscFunctionReturn(0); 4573 } 4574 4575 PetscErrorCode MatSetPreallocationCOO_SeqAIJ(Mat mat, PetscCount coo_n, PetscInt coo_i[], PetscInt coo_j[]) 4576 { 4577 MPI_Comm comm; 4578 PetscInt *i, *j; 4579 PetscInt M, N, row; 4580 PetscCount k, p, q, nneg, nnz, start, end; /* Index the coo array, so use PetscCount as their type */ 4581 PetscInt *Ai; /* Change to PetscCount once we use it for row pointers */ 4582 PetscInt *Aj; 4583 PetscScalar *Aa; 4584 Mat_SeqAIJ *seqaij = (Mat_SeqAIJ *)(mat->data); 4585 MatType rtype; 4586 PetscCount *perm, *jmap; 4587 4588 PetscFunctionBegin; 4589 PetscCall(MatResetPreallocationCOO_SeqAIJ(mat)); 4590 PetscCall(PetscObjectGetComm((PetscObject)mat, &comm)); 4591 PetscCall(MatGetSize(mat, &M, &N)); 4592 i = coo_i; 4593 j = coo_j; 4594 PetscCall(PetscMalloc1(coo_n, &perm)); 4595 for (k = 0; k < coo_n; k++) { /* Ignore entries with negative row or col indices */ 4596 if (j[k] < 0) i[k] = -1; 4597 perm[k] = k; 4598 } 4599 4600 /* Sort by row */ 4601 PetscCall(PetscSortIntWithIntCountArrayPair(coo_n, i, j, perm)); 4602 for (k = 0; k < coo_n; k++) { 4603 if (i[k] >= 0) break; 4604 } /* Advance k to the first row with a non-negative index */ 4605 nneg = k; 4606 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 */ 4607 nnz = 0; /* Total number of unique nonzeros to be counted */ 4608 jmap++; /* Inc jmap by 1 for convinience */ 4609 4610 PetscCall(PetscCalloc1(M + 1, &Ai)); /* CSR of A */ 4611 PetscCall(PetscMalloc1(coo_n - nneg, &Aj)); /* We have at most coo_n-nneg unique nonzeros */ 4612 4613 /* In each row, sort by column, then unique column indices to get row length */ 4614 Ai++; /* Inc by 1 for convinience */ 4615 q = 0; /* q-th unique nonzero, with q starting from 0 */ 4616 while (k < coo_n) { 4617 row = i[k]; 4618 start = k; /* [start,end) indices for this row */ 4619 while (k < coo_n && i[k] == row) k++; 4620 end = k; 4621 PetscCall(PetscSortIntWithCountArray(end - start, j + start, perm + start)); 4622 /* Find number of unique col entries in this row */ 4623 Aj[q] = j[start]; /* Log the first nonzero in this row */ 4624 jmap[q] = 1; /* Number of repeats of this nozero entry */ 4625 Ai[row] = 1; 4626 nnz++; 4627 4628 for (p = start + 1; p < end; p++) { /* Scan remaining nonzero in this row */ 4629 if (j[p] != j[p - 1]) { /* Meet a new nonzero */ 4630 q++; 4631 jmap[q] = 1; 4632 Aj[q] = j[p]; 4633 Ai[row]++; 4634 nnz++; 4635 } else { 4636 jmap[q]++; 4637 } 4638 } 4639 q++; /* Move to next row and thus next unique nonzero */ 4640 } 4641 4642 Ai--; /* Back to the beginning of Ai[] */ 4643 for (k = 0; k < M; k++) Ai[k + 1] += Ai[k]; 4644 jmap--; /* Back to the beginning of jmap[] */ 4645 jmap[0] = 0; 4646 for (k = 0; k < nnz; k++) jmap[k + 1] += jmap[k]; 4647 if (nnz < coo_n - nneg) { /* Realloc with actual number of unique nonzeros */ 4648 PetscCount *jmap_new; 4649 PetscInt *Aj_new; 4650 4651 PetscCall(PetscMalloc1(nnz + 1, &jmap_new)); 4652 PetscCall(PetscArraycpy(jmap_new, jmap, nnz + 1)); 4653 PetscCall(PetscFree(jmap)); 4654 jmap = jmap_new; 4655 4656 PetscCall(PetscMalloc1(nnz, &Aj_new)); 4657 PetscCall(PetscArraycpy(Aj_new, Aj, nnz)); 4658 PetscCall(PetscFree(Aj)); 4659 Aj = Aj_new; 4660 } 4661 4662 if (nneg) { /* Discard heading entries with negative indices in perm[], as we'll access it from index 0 in MatSetValuesCOO */ 4663 PetscCount *perm_new; 4664 4665 PetscCall(PetscMalloc1(coo_n - nneg, &perm_new)); 4666 PetscCall(PetscArraycpy(perm_new, perm + nneg, coo_n - nneg)); 4667 PetscCall(PetscFree(perm)); 4668 perm = perm_new; 4669 } 4670 4671 PetscCall(MatGetRootType_Private(mat, &rtype)); 4672 PetscCall(PetscCalloc1(nnz, &Aa)); /* Zero the matrix */ 4673 PetscCall(MatSetSeqAIJWithArrays_private(PETSC_COMM_SELF, M, N, Ai, Aj, Aa, rtype, mat)); 4674 4675 seqaij->singlemalloc = PETSC_FALSE; /* Ai, Aj and Aa are not allocated in one big malloc */ 4676 seqaij->free_a = seqaij->free_ij = PETSC_TRUE; /* Let newmat own Ai, Aj and Aa */ 4677 /* Record COO fields */ 4678 seqaij->coo_n = coo_n; 4679 seqaij->Atot = coo_n - nneg; /* Annz is seqaij->nz, so no need to record that again */ 4680 seqaij->jmap = jmap; /* of length nnz+1 */ 4681 seqaij->perm = perm; 4682 PetscFunctionReturn(0); 4683 } 4684 4685 static PetscErrorCode MatSetValuesCOO_SeqAIJ(Mat A, const PetscScalar v[], InsertMode imode) 4686 { 4687 Mat_SeqAIJ *aseq = (Mat_SeqAIJ *)A->data; 4688 PetscCount i, j, Annz = aseq->nz; 4689 PetscCount *perm = aseq->perm, *jmap = aseq->jmap; 4690 PetscScalar *Aa; 4691 4692 PetscFunctionBegin; 4693 PetscCall(MatSeqAIJGetArray(A, &Aa)); 4694 for (i = 0; i < Annz; i++) { 4695 PetscScalar sum = 0.0; 4696 for (j = jmap[i]; j < jmap[i + 1]; j++) sum += v[perm[j]]; 4697 Aa[i] = (imode == INSERT_VALUES ? 0.0 : Aa[i]) + sum; 4698 } 4699 PetscCall(MatSeqAIJRestoreArray(A, &Aa)); 4700 PetscFunctionReturn(0); 4701 } 4702 4703 #if defined(PETSC_HAVE_CUDA) 4704 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCUSPARSE(Mat, MatType, MatReuse, Mat *); 4705 #endif 4706 #if defined(PETSC_HAVE_KOKKOS_KERNELS) 4707 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJKokkos(Mat, MatType, MatReuse, Mat *); 4708 #endif 4709 4710 PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJ(Mat B) 4711 { 4712 Mat_SeqAIJ *b; 4713 PetscMPIInt size; 4714 4715 PetscFunctionBegin; 4716 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B), &size)); 4717 PetscCheck(size <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Comm must be of size 1"); 4718 4719 PetscCall(PetscNew(&b)); 4720 4721 B->data = (void *)b; 4722 4723 PetscCall(PetscMemcpy(B->ops, &MatOps_Values, sizeof(struct _MatOps))); 4724 if (B->sortedfull) B->ops->setvalues = MatSetValues_SeqAIJ_SortedFull; 4725 4726 b->row = NULL; 4727 b->col = NULL; 4728 b->icol = NULL; 4729 b->reallocs = 0; 4730 b->ignorezeroentries = PETSC_FALSE; 4731 b->roworiented = PETSC_TRUE; 4732 b->nonew = 0; 4733 b->diag = NULL; 4734 b->solve_work = NULL; 4735 B->spptr = NULL; 4736 b->saved_values = NULL; 4737 b->idiag = NULL; 4738 b->mdiag = NULL; 4739 b->ssor_work = NULL; 4740 b->omega = 1.0; 4741 b->fshift = 0.0; 4742 b->idiagvalid = PETSC_FALSE; 4743 b->ibdiagvalid = PETSC_FALSE; 4744 b->keepnonzeropattern = PETSC_FALSE; 4745 4746 PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATSEQAIJ)); 4747 #if defined(PETSC_HAVE_MATLAB) 4748 PetscCall(PetscObjectComposeFunction((PetscObject)B, "PetscMatlabEnginePut_C", MatlabEnginePut_SeqAIJ)); 4749 PetscCall(PetscObjectComposeFunction((PetscObject)B, "PetscMatlabEngineGet_C", MatlabEngineGet_SeqAIJ)); 4750 #endif 4751 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqAIJSetColumnIndices_C", MatSeqAIJSetColumnIndices_SeqAIJ)); 4752 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatStoreValues_C", MatStoreValues_SeqAIJ)); 4753 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatRetrieveValues_C", MatRetrieveValues_SeqAIJ)); 4754 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqsbaij_C", MatConvert_SeqAIJ_SeqSBAIJ)); 4755 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqbaij_C", MatConvert_SeqAIJ_SeqBAIJ)); 4756 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqaijperm_C", MatConvert_SeqAIJ_SeqAIJPERM)); 4757 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqaijsell_C", MatConvert_SeqAIJ_SeqAIJSELL)); 4758 #if defined(PETSC_HAVE_MKL_SPARSE) 4759 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqaijmkl_C", MatConvert_SeqAIJ_SeqAIJMKL)); 4760 #endif 4761 #if defined(PETSC_HAVE_CUDA) 4762 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqaijcusparse_C", MatConvert_SeqAIJ_SeqAIJCUSPARSE)); 4763 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaijcusparse_seqaij_C", MatProductSetFromOptions_SeqAIJ)); 4764 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaij_seqaijcusparse_C", MatProductSetFromOptions_SeqAIJ)); 4765 #endif 4766 #if defined(PETSC_HAVE_KOKKOS_KERNELS) 4767 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqaijkokkos_C", MatConvert_SeqAIJ_SeqAIJKokkos)); 4768 #endif 4769 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqaijcrl_C", MatConvert_SeqAIJ_SeqAIJCRL)); 4770 #if defined(PETSC_HAVE_ELEMENTAL) 4771 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_elemental_C", MatConvert_SeqAIJ_Elemental)); 4772 #endif 4773 #if defined(PETSC_HAVE_SCALAPACK) 4774 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_scalapack_C", MatConvert_AIJ_ScaLAPACK)); 4775 #endif 4776 #if defined(PETSC_HAVE_HYPRE) 4777 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_hypre_C", MatConvert_AIJ_HYPRE)); 4778 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_transpose_seqaij_seqaij_C", MatProductSetFromOptions_Transpose_AIJ_AIJ)); 4779 #endif 4780 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqdense_C", MatConvert_SeqAIJ_SeqDense)); 4781 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqsell_C", MatConvert_SeqAIJ_SeqSELL)); 4782 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_is_C", MatConvert_XAIJ_IS)); 4783 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatIsTranspose_C", MatIsTranspose_SeqAIJ)); 4784 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatIsHermitianTranspose_C", MatIsTranspose_SeqAIJ)); 4785 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqAIJSetPreallocation_C", MatSeqAIJSetPreallocation_SeqAIJ)); 4786 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatResetPreallocation_C", MatResetPreallocation_SeqAIJ)); 4787 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqAIJSetPreallocationCSR_C", MatSeqAIJSetPreallocationCSR_SeqAIJ)); 4788 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatReorderForNonzeroDiagonal_C", MatReorderForNonzeroDiagonal_SeqAIJ)); 4789 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_is_seqaij_C", MatProductSetFromOptions_IS_XAIJ)); 4790 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdense_seqaij_C", MatProductSetFromOptions_SeqDense_SeqAIJ)); 4791 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaij_seqaij_C", MatProductSetFromOptions_SeqAIJ)); 4792 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqAIJKron_C", MatSeqAIJKron_SeqAIJ)); 4793 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSetPreallocationCOO_C", MatSetPreallocationCOO_SeqAIJ)); 4794 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSetValuesCOO_C", MatSetValuesCOO_SeqAIJ)); 4795 PetscCall(MatCreate_SeqAIJ_Inode(B)); 4796 PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATSEQAIJ)); 4797 PetscCall(MatSeqAIJSetTypeFromOptions(B)); /* this allows changing the matrix subtype to say MATSEQAIJPERM */ 4798 PetscFunctionReturn(0); 4799 } 4800 4801 /* 4802 Given a matrix generated with MatGetFactor() duplicates all the information in A into C 4803 */ 4804 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C, Mat A, MatDuplicateOption cpvalues, PetscBool mallocmatspace) 4805 { 4806 Mat_SeqAIJ *c = (Mat_SeqAIJ *)C->data, *a = (Mat_SeqAIJ *)A->data; 4807 PetscInt m = A->rmap->n, i; 4808 4809 PetscFunctionBegin; 4810 PetscCheck(A->assembled || cpvalues == MAT_DO_NOT_COPY_VALUES, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cannot duplicate unassembled matrix"); 4811 4812 C->factortype = A->factortype; 4813 c->row = NULL; 4814 c->col = NULL; 4815 c->icol = NULL; 4816 c->reallocs = 0; 4817 4818 C->assembled = A->assembled; 4819 C->preallocated = A->preallocated; 4820 4821 if (A->preallocated) { 4822 PetscCall(PetscLayoutReference(A->rmap, &C->rmap)); 4823 PetscCall(PetscLayoutReference(A->cmap, &C->cmap)); 4824 4825 PetscCall(PetscMalloc1(m, &c->imax)); 4826 PetscCall(PetscMemcpy(c->imax, a->imax, m * sizeof(PetscInt))); 4827 PetscCall(PetscMalloc1(m, &c->ilen)); 4828 PetscCall(PetscMemcpy(c->ilen, a->ilen, m * sizeof(PetscInt))); 4829 4830 /* allocate the matrix space */ 4831 if (mallocmatspace) { 4832 PetscCall(PetscMalloc3(a->i[m], &c->a, a->i[m], &c->j, m + 1, &c->i)); 4833 4834 c->singlemalloc = PETSC_TRUE; 4835 4836 PetscCall(PetscArraycpy(c->i, a->i, m + 1)); 4837 if (m > 0) { 4838 PetscCall(PetscArraycpy(c->j, a->j, a->i[m])); 4839 if (cpvalues == MAT_COPY_VALUES) { 4840 const PetscScalar *aa; 4841 4842 PetscCall(MatSeqAIJGetArrayRead(A, &aa)); 4843 PetscCall(PetscArraycpy(c->a, aa, a->i[m])); 4844 PetscCall(MatSeqAIJGetArrayRead(A, &aa)); 4845 } else { 4846 PetscCall(PetscArrayzero(c->a, a->i[m])); 4847 } 4848 } 4849 } 4850 4851 c->ignorezeroentries = a->ignorezeroentries; 4852 c->roworiented = a->roworiented; 4853 c->nonew = a->nonew; 4854 if (a->diag) { 4855 PetscCall(PetscMalloc1(m + 1, &c->diag)); 4856 PetscCall(PetscMemcpy(c->diag, a->diag, m * sizeof(PetscInt))); 4857 } else c->diag = NULL; 4858 4859 c->solve_work = NULL; 4860 c->saved_values = NULL; 4861 c->idiag = NULL; 4862 c->ssor_work = NULL; 4863 c->keepnonzeropattern = a->keepnonzeropattern; 4864 c->free_a = PETSC_TRUE; 4865 c->free_ij = PETSC_TRUE; 4866 4867 c->rmax = a->rmax; 4868 c->nz = a->nz; 4869 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 4870 4871 c->compressedrow.use = a->compressedrow.use; 4872 c->compressedrow.nrows = a->compressedrow.nrows; 4873 if (a->compressedrow.use) { 4874 i = a->compressedrow.nrows; 4875 PetscCall(PetscMalloc2(i + 1, &c->compressedrow.i, i, &c->compressedrow.rindex)); 4876 PetscCall(PetscArraycpy(c->compressedrow.i, a->compressedrow.i, i + 1)); 4877 PetscCall(PetscArraycpy(c->compressedrow.rindex, a->compressedrow.rindex, i)); 4878 } else { 4879 c->compressedrow.use = PETSC_FALSE; 4880 c->compressedrow.i = NULL; 4881 c->compressedrow.rindex = NULL; 4882 } 4883 c->nonzerorowcnt = a->nonzerorowcnt; 4884 C->nonzerostate = A->nonzerostate; 4885 4886 PetscCall(MatDuplicate_SeqAIJ_Inode(A, cpvalues, &C)); 4887 } 4888 PetscCall(PetscFunctionListDuplicate(((PetscObject)A)->qlist, &((PetscObject)C)->qlist)); 4889 PetscFunctionReturn(0); 4890 } 4891 4892 PetscErrorCode MatDuplicate_SeqAIJ(Mat A, MatDuplicateOption cpvalues, Mat *B) 4893 { 4894 PetscFunctionBegin; 4895 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), B)); 4896 PetscCall(MatSetSizes(*B, A->rmap->n, A->cmap->n, A->rmap->n, A->cmap->n)); 4897 if (!(A->rmap->n % A->rmap->bs) && !(A->cmap->n % A->cmap->bs)) PetscCall(MatSetBlockSizesFromMats(*B, A, A)); 4898 PetscCall(MatSetType(*B, ((PetscObject)A)->type_name)); 4899 PetscCall(MatDuplicateNoCreate_SeqAIJ(*B, A, cpvalues, PETSC_TRUE)); 4900 PetscFunctionReturn(0); 4901 } 4902 4903 PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer) 4904 { 4905 PetscBool isbinary, ishdf5; 4906 4907 PetscFunctionBegin; 4908 PetscValidHeaderSpecific(newMat, MAT_CLASSID, 1); 4909 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 4910 /* force binary viewer to load .info file if it has not yet done so */ 4911 PetscCall(PetscViewerSetUp(viewer)); 4912 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary)); 4913 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5)); 4914 if (isbinary) { 4915 PetscCall(MatLoad_SeqAIJ_Binary(newMat, viewer)); 4916 } else if (ishdf5) { 4917 #if defined(PETSC_HAVE_HDF5) 4918 PetscCall(MatLoad_AIJ_HDF5(newMat, viewer)); 4919 #else 4920 SETERRQ(PetscObjectComm((PetscObject)newMat), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5"); 4921 #endif 4922 } else { 4923 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); 4924 } 4925 PetscFunctionReturn(0); 4926 } 4927 4928 PetscErrorCode MatLoad_SeqAIJ_Binary(Mat mat, PetscViewer viewer) 4929 { 4930 Mat_SeqAIJ *a = (Mat_SeqAIJ *)mat->data; 4931 PetscInt header[4], *rowlens, M, N, nz, sum, rows, cols, i; 4932 4933 PetscFunctionBegin; 4934 PetscCall(PetscViewerSetUp(viewer)); 4935 4936 /* read in matrix header */ 4937 PetscCall(PetscViewerBinaryRead(viewer, header, 4, NULL, PETSC_INT)); 4938 PetscCheck(header[0] == MAT_FILE_CLASSID, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Not a matrix object in file"); 4939 M = header[1]; 4940 N = header[2]; 4941 nz = header[3]; 4942 PetscCheck(M >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Matrix row size (%" PetscInt_FMT ") in file is negative", M); 4943 PetscCheck(N >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Matrix column size (%" PetscInt_FMT ") in file is negative", N); 4944 PetscCheck(nz >= 0, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Matrix stored in special format on disk, cannot load as SeqAIJ"); 4945 4946 /* set block sizes from the viewer's .info file */ 4947 PetscCall(MatLoad_Binary_BlockSizes(mat, viewer)); 4948 /* set local and global sizes if not set already */ 4949 if (mat->rmap->n < 0) mat->rmap->n = M; 4950 if (mat->cmap->n < 0) mat->cmap->n = N; 4951 if (mat->rmap->N < 0) mat->rmap->N = M; 4952 if (mat->cmap->N < 0) mat->cmap->N = N; 4953 PetscCall(PetscLayoutSetUp(mat->rmap)); 4954 PetscCall(PetscLayoutSetUp(mat->cmap)); 4955 4956 /* check if the matrix sizes are correct */ 4957 PetscCall(MatGetSize(mat, &rows, &cols)); 4958 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); 4959 4960 /* read in row lengths */ 4961 PetscCall(PetscMalloc1(M, &rowlens)); 4962 PetscCall(PetscViewerBinaryRead(viewer, rowlens, M, NULL, PETSC_INT)); 4963 /* check if sum(rowlens) is same as nz */ 4964 sum = 0; 4965 for (i = 0; i < M; i++) sum += rowlens[i]; 4966 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); 4967 /* preallocate and check sizes */ 4968 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(mat, 0, rowlens)); 4969 PetscCall(MatGetSize(mat, &rows, &cols)); 4970 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); 4971 /* store row lengths */ 4972 PetscCall(PetscArraycpy(a->ilen, rowlens, M)); 4973 PetscCall(PetscFree(rowlens)); 4974 4975 /* fill in "i" row pointers */ 4976 a->i[0] = 0; 4977 for (i = 0; i < M; i++) a->i[i + 1] = a->i[i] + a->ilen[i]; 4978 /* read in "j" column indices */ 4979 PetscCall(PetscViewerBinaryRead(viewer, a->j, nz, NULL, PETSC_INT)); 4980 /* read in "a" nonzero values */ 4981 PetscCall(PetscViewerBinaryRead(viewer, a->a, nz, NULL, PETSC_SCALAR)); 4982 4983 PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY)); 4984 PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY)); 4985 PetscFunctionReturn(0); 4986 } 4987 4988 PetscErrorCode MatEqual_SeqAIJ(Mat A, Mat B, PetscBool *flg) 4989 { 4990 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data; 4991 const PetscScalar *aa, *ba; 4992 #if defined(PETSC_USE_COMPLEX) 4993 PetscInt k; 4994 #endif 4995 4996 PetscFunctionBegin; 4997 /* If the matrix dimensions are not equal,or no of nonzeros */ 4998 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) || (a->nz != b->nz)) { 4999 *flg = PETSC_FALSE; 5000 PetscFunctionReturn(0); 5001 } 5002 5003 /* if the a->i are the same */ 5004 PetscCall(PetscArraycmp(a->i, b->i, A->rmap->n + 1, flg)); 5005 if (!*flg) PetscFunctionReturn(0); 5006 5007 /* if a->j are the same */ 5008 PetscCall(PetscArraycmp(a->j, b->j, a->nz, flg)); 5009 if (!*flg) PetscFunctionReturn(0); 5010 5011 PetscCall(MatSeqAIJGetArrayRead(A, &aa)); 5012 PetscCall(MatSeqAIJGetArrayRead(B, &ba)); 5013 /* if a->a are the same */ 5014 #if defined(PETSC_USE_COMPLEX) 5015 for (k = 0; k < a->nz; k++) { 5016 if (PetscRealPart(aa[k]) != PetscRealPart(ba[k]) || PetscImaginaryPart(aa[k]) != PetscImaginaryPart(ba[k])) { 5017 *flg = PETSC_FALSE; 5018 PetscFunctionReturn(0); 5019 } 5020 } 5021 #else 5022 PetscCall(PetscArraycmp(aa, ba, a->nz, flg)); 5023 #endif 5024 PetscCall(MatSeqAIJRestoreArrayRead(A, &aa)); 5025 PetscCall(MatSeqAIJRestoreArrayRead(B, &ba)); 5026 PetscFunctionReturn(0); 5027 } 5028 5029 /*@ 5030 MatCreateSeqAIJWithArrays - Creates an sequential `MATSEQAIJ` matrix using matrix elements (in CSR format) 5031 provided by the user. 5032 5033 Collective 5034 5035 Input Parameters: 5036 + comm - must be an MPI communicator of size 1 5037 . m - number of rows 5038 . n - number of columns 5039 . i - row indices; that is i[0] = 0, i[row] = i[row-1] + number of elements in that row of the matrix 5040 . j - column indices 5041 - a - matrix values 5042 5043 Output Parameter: 5044 . mat - the matrix 5045 5046 Level: intermediate 5047 5048 Notes: 5049 The i, j, and a arrays are not copied by this routine, the user must free these arrays 5050 once the matrix is destroyed and not before 5051 5052 You cannot set new nonzero locations into this matrix, that will generate an error. 5053 5054 The i and j indices are 0 based 5055 5056 The format which is used for the sparse matrix input, is equivalent to a 5057 row-major ordering.. i.e for the following matrix, the input data expected is 5058 as shown 5059 5060 $ 1 0 0 5061 $ 2 0 3 5062 $ 4 5 6 5063 $ 5064 $ i = {0,1,3,6} [size = nrow+1 = 3+1] 5065 $ j = {0,0,2,0,1,2} [size = 6]; values must be sorted for each row 5066 $ v = {1,2,3,4,5,6} [size = 6] 5067 5068 .seealso: `MatCreate()`, `MatCreateAIJ()`, `MatCreateSeqAIJ()`, `MatCreateMPIAIJWithArrays()`, `MatMPIAIJSetPreallocationCSR()` 5069 @*/ 5070 PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt i[], PetscInt j[], PetscScalar a[], Mat *mat) 5071 { 5072 PetscInt ii; 5073 Mat_SeqAIJ *aij; 5074 PetscInt jj; 5075 5076 PetscFunctionBegin; 5077 PetscCheck(m <= 0 || i[0] == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "i (row indices) must start with 0"); 5078 PetscCall(MatCreate(comm, mat)); 5079 PetscCall(MatSetSizes(*mat, m, n, m, n)); 5080 /* PetscCall(MatSetBlockSizes(*mat,,)); */ 5081 PetscCall(MatSetType(*mat, MATSEQAIJ)); 5082 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(*mat, MAT_SKIP_ALLOCATION, NULL)); 5083 aij = (Mat_SeqAIJ *)(*mat)->data; 5084 PetscCall(PetscMalloc1(m, &aij->imax)); 5085 PetscCall(PetscMalloc1(m, &aij->ilen)); 5086 5087 aij->i = i; 5088 aij->j = j; 5089 aij->a = a; 5090 aij->singlemalloc = PETSC_FALSE; 5091 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 5092 aij->free_a = PETSC_FALSE; 5093 aij->free_ij = PETSC_FALSE; 5094 5095 for (ii = 0, aij->nonzerorowcnt = 0, aij->rmax = 0; ii < m; ii++) { 5096 aij->ilen[ii] = aij->imax[ii] = i[ii + 1] - i[ii]; 5097 if (PetscDefined(USE_DEBUG)) { 5098 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]); 5099 for (jj = i[ii] + 1; jj < i[ii + 1]; jj++) { 5100 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); 5101 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); 5102 } 5103 } 5104 } 5105 if (PetscDefined(USE_DEBUG)) { 5106 for (ii = 0; ii < aij->i[m]; ii++) { 5107 PetscCheck(j[ii] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative column index at location = %" PetscInt_FMT " index = %" PetscInt_FMT, ii, j[ii]); 5108 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]); 5109 } 5110 } 5111 5112 PetscCall(MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY)); 5113 PetscCall(MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY)); 5114 PetscFunctionReturn(0); 5115 } 5116 5117 /*@ 5118 MatCreateSeqAIJFromTriple - Creates an sequential `MATSEQAIJ` matrix using matrix elements (in COO format) 5119 provided by the user. 5120 5121 Collective 5122 5123 Input Parameters: 5124 + comm - must be an MPI communicator of size 1 5125 . m - number of rows 5126 . n - number of columns 5127 . i - row indices 5128 . j - column indices 5129 . a - matrix values 5130 . nz - number of nonzeros 5131 - idx - if the i and j indices start with 1 use `PETSC_TRUE` otherwise use `PETSC_FALSE` 5132 5133 Output Parameter: 5134 . mat - the matrix 5135 5136 Level: intermediate 5137 5138 Example: 5139 For the following matrix, the input data expected is as shown (using 0 based indexing) 5140 .vb 5141 1 0 0 5142 2 0 3 5143 4 5 6 5144 5145 i = {0,1,1,2,2,2} 5146 j = {0,0,2,0,1,2} 5147 v = {1,2,3,4,5,6} 5148 .ve 5149 Notes: 5150 Instead of using this function, users should also consider `MatSetPreallocationCOO()` and `MatSetValuesCOO()`, which allow repeated or remote entries, 5151 and are particularly useful in iterative applications. 5152 5153 .seealso: `MatCreate()`, `MatCreateAIJ()`, `MatCreateSeqAIJ()`, `MatCreateSeqAIJWithArrays()`, `MatMPIAIJSetPreallocationCSR()`, `MatSetValuesCOO()`, `MatSetPreallocationCOO()` 5154 @*/ 5155 PetscErrorCode MatCreateSeqAIJFromTriple(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt i[], PetscInt j[], PetscScalar a[], Mat *mat, PetscInt nz, PetscBool idx) 5156 { 5157 PetscInt ii, *nnz, one = 1, row, col; 5158 5159 PetscFunctionBegin; 5160 PetscCall(PetscCalloc1(m, &nnz)); 5161 for (ii = 0; ii < nz; ii++) nnz[i[ii] - !!idx] += 1; 5162 PetscCall(MatCreate(comm, mat)); 5163 PetscCall(MatSetSizes(*mat, m, n, m, n)); 5164 PetscCall(MatSetType(*mat, MATSEQAIJ)); 5165 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(*mat, 0, nnz)); 5166 for (ii = 0; ii < nz; ii++) { 5167 if (idx) { 5168 row = i[ii] - 1; 5169 col = j[ii] - 1; 5170 } else { 5171 row = i[ii]; 5172 col = j[ii]; 5173 } 5174 PetscCall(MatSetValues(*mat, one, &row, one, &col, &a[ii], ADD_VALUES)); 5175 } 5176 PetscCall(MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY)); 5177 PetscCall(MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY)); 5178 PetscCall(PetscFree(nnz)); 5179 PetscFunctionReturn(0); 5180 } 5181 5182 PetscErrorCode MatSeqAIJInvalidateDiagonal(Mat A) 5183 { 5184 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 5185 5186 PetscFunctionBegin; 5187 a->idiagvalid = PETSC_FALSE; 5188 a->ibdiagvalid = PETSC_FALSE; 5189 5190 PetscCall(MatSeqAIJInvalidateDiagonal_Inode(A)); 5191 PetscFunctionReturn(0); 5192 } 5193 5194 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqAIJ(MPI_Comm comm, Mat inmat, PetscInt n, MatReuse scall, Mat *outmat) 5195 { 5196 PetscFunctionBegin; 5197 PetscCall(MatCreateMPIMatConcatenateSeqMat_MPIAIJ(comm, inmat, n, scall, outmat)); 5198 PetscFunctionReturn(0); 5199 } 5200 5201 /* 5202 Permute A into C's *local* index space using rowemb,colemb. 5203 The embedding are supposed to be injections and the above implies that the range of rowemb is a subset 5204 of [0,m), colemb is in [0,n). 5205 If pattern == DIFFERENT_NONZERO_PATTERN, C is preallocated according to A. 5206 */ 5207 PetscErrorCode MatSetSeqMat_SeqAIJ(Mat C, IS rowemb, IS colemb, MatStructure pattern, Mat B) 5208 { 5209 /* If making this function public, change the error returned in this function away from _PLIB. */ 5210 Mat_SeqAIJ *Baij; 5211 PetscBool seqaij; 5212 PetscInt m, n, *nz, i, j, count; 5213 PetscScalar v; 5214 const PetscInt *rowindices, *colindices; 5215 5216 PetscFunctionBegin; 5217 if (!B) PetscFunctionReturn(0); 5218 /* Check to make sure the target matrix (and embeddings) are compatible with C and each other. */ 5219 PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATSEQAIJ, &seqaij)); 5220 PetscCheck(seqaij, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Input matrix is of wrong type"); 5221 if (rowemb) { 5222 PetscCall(ISGetLocalSize(rowemb, &m)); 5223 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); 5224 } else { 5225 PetscCheck(C->rmap->n == B->rmap->n, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Input matrix is row-incompatible with the target matrix"); 5226 } 5227 if (colemb) { 5228 PetscCall(ISGetLocalSize(colemb, &n)); 5229 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); 5230 } else { 5231 PetscCheck(C->cmap->n == B->cmap->n, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Input matrix is col-incompatible with the target matrix"); 5232 } 5233 5234 Baij = (Mat_SeqAIJ *)(B->data); 5235 if (pattern == DIFFERENT_NONZERO_PATTERN) { 5236 PetscCall(PetscMalloc1(B->rmap->n, &nz)); 5237 for (i = 0; i < B->rmap->n; i++) nz[i] = Baij->i[i + 1] - Baij->i[i]; 5238 PetscCall(MatSeqAIJSetPreallocation(C, 0, nz)); 5239 PetscCall(PetscFree(nz)); 5240 } 5241 if (pattern == SUBSET_NONZERO_PATTERN) PetscCall(MatZeroEntries(C)); 5242 count = 0; 5243 rowindices = NULL; 5244 colindices = NULL; 5245 if (rowemb) PetscCall(ISGetIndices(rowemb, &rowindices)); 5246 if (colemb) PetscCall(ISGetIndices(colemb, &colindices)); 5247 for (i = 0; i < B->rmap->n; i++) { 5248 PetscInt row; 5249 row = i; 5250 if (rowindices) row = rowindices[i]; 5251 for (j = Baij->i[i]; j < Baij->i[i + 1]; j++) { 5252 PetscInt col; 5253 col = Baij->j[count]; 5254 if (colindices) col = colindices[col]; 5255 v = Baij->a[count]; 5256 PetscCall(MatSetValues(C, 1, &row, 1, &col, &v, INSERT_VALUES)); 5257 ++count; 5258 } 5259 } 5260 /* FIXME: set C's nonzerostate correctly. */ 5261 /* Assembly for C is necessary. */ 5262 C->preallocated = PETSC_TRUE; 5263 C->assembled = PETSC_TRUE; 5264 C->was_assembled = PETSC_FALSE; 5265 PetscFunctionReturn(0); 5266 } 5267 5268 PetscErrorCode MatEliminateZeros_SeqAIJ(Mat A) 5269 { 5270 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 5271 MatScalar *aa = a->a; 5272 PetscInt m = A->rmap->n, fshift = 0, fshift_prev = 0, i, k; 5273 PetscInt *ailen = a->ilen, *imax = a->imax, *ai = a->i, *aj = a->j, rmax = 0; 5274 5275 PetscFunctionBegin; 5276 PetscCheck(A->assembled, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cannot eliminate zeros for unassembled matrix"); 5277 if (m) rmax = ailen[0]; /* determine row with most nonzeros */ 5278 for (i = 1; i <= m; i++) { 5279 /* move each nonzero entry back by the amount of zero slots (fshift) before it*/ 5280 for (k = ai[i - 1]; k < ai[i]; k++) { 5281 if (aa[k] == 0 && aj[k] != i - 1) fshift++; 5282 else { 5283 if (aa[k] == 0 && aj[k] == i - 1) PetscCall(PetscInfo(A, "Keep the diagonal zero at row %" PetscInt_FMT "\n", i - 1)); 5284 aa[k - fshift] = aa[k]; 5285 aj[k - fshift] = aj[k]; 5286 } 5287 } 5288 ai[i - 1] -= fshift_prev; // safe to update ai[i-1] now since it will not be used in the next iteration 5289 fshift_prev = fshift; 5290 /* reset ilen and imax for each row */ 5291 ailen[i - 1] = imax[i - 1] = ai[i] - fshift - ai[i - 1]; 5292 a->nonzerorowcnt += ((ai[i] - fshift - ai[i - 1]) > 0); 5293 rmax = PetscMax(rmax, ailen[i - 1]); 5294 } 5295 if (m) { 5296 ai[m] -= fshift; 5297 a->nz = ai[m]; 5298 } 5299 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)); 5300 A->nonzerostate -= fshift; 5301 A->info.nz_unneeded += (PetscReal)fshift; 5302 a->rmax = rmax; 5303 if (a->inode.use && a->inode.checked) PetscCall(MatSeqAIJCheckInode(A)); 5304 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 5305 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 5306 PetscFunctionReturn(0); 5307 } 5308 5309 PetscFunctionList MatSeqAIJList = NULL; 5310 5311 /*@C 5312 MatSeqAIJSetType - Converts a `MATSEQAIJ` matrix to a subtype 5313 5314 Collective on mat 5315 5316 Input Parameters: 5317 + mat - the matrix object 5318 - matype - matrix type 5319 5320 Options Database Key: 5321 . -mat_seqai_type <method> - for example seqaijcrl 5322 5323 Level: intermediate 5324 5325 .seealso: `PCSetType()`, `VecSetType()`, `MatCreate()`, `MatType`, `Mat` 5326 @*/ 5327 PetscErrorCode MatSeqAIJSetType(Mat mat, MatType matype) 5328 { 5329 PetscBool sametype; 5330 PetscErrorCode (*r)(Mat, MatType, MatReuse, Mat *); 5331 5332 PetscFunctionBegin; 5333 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 5334 PetscCall(PetscObjectTypeCompare((PetscObject)mat, matype, &sametype)); 5335 if (sametype) PetscFunctionReturn(0); 5336 5337 PetscCall(PetscFunctionListFind(MatSeqAIJList, matype, &r)); 5338 PetscCheck(r, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown Mat type given: %s", matype); 5339 PetscCall((*r)(mat, matype, MAT_INPLACE_MATRIX, &mat)); 5340 PetscFunctionReturn(0); 5341 } 5342 5343 /*@C 5344 MatSeqAIJRegister - - Adds a new sub-matrix type for sequential `MATSEQAIJ` matrices 5345 5346 Not Collective 5347 5348 Input Parameters: 5349 + name - name of a new user-defined matrix type, for example `MATSEQAIJCRL` 5350 - function - routine to convert to subtype 5351 5352 Notes: 5353 `MatSeqAIJRegister()` may be called multiple times to add several user-defined solvers. 5354 5355 Then, your matrix can be chosen with the procedural interface at runtime via the option 5356 $ -mat_seqaij_type my_mat 5357 5358 Level: advanced 5359 5360 .seealso: `MatSeqAIJRegisterAll()` 5361 @*/ 5362 PetscErrorCode MatSeqAIJRegister(const char sname[], PetscErrorCode (*function)(Mat, MatType, MatReuse, Mat *)) 5363 { 5364 PetscFunctionBegin; 5365 PetscCall(MatInitializePackage()); 5366 PetscCall(PetscFunctionListAdd(&MatSeqAIJList, sname, function)); 5367 PetscFunctionReturn(0); 5368 } 5369 5370 PetscBool MatSeqAIJRegisterAllCalled = PETSC_FALSE; 5371 5372 /*@C 5373 MatSeqAIJRegisterAll - Registers all of the matrix subtypes of `MATSSEQAIJ` 5374 5375 Not Collective 5376 5377 Level: advanced 5378 5379 .seealso: `MatRegisterAll()`, `MatSeqAIJRegister()` 5380 @*/ 5381 PetscErrorCode MatSeqAIJRegisterAll(void) 5382 { 5383 PetscFunctionBegin; 5384 if (MatSeqAIJRegisterAllCalled) PetscFunctionReturn(0); 5385 MatSeqAIJRegisterAllCalled = PETSC_TRUE; 5386 5387 PetscCall(MatSeqAIJRegister(MATSEQAIJCRL, MatConvert_SeqAIJ_SeqAIJCRL)); 5388 PetscCall(MatSeqAIJRegister(MATSEQAIJPERM, MatConvert_SeqAIJ_SeqAIJPERM)); 5389 PetscCall(MatSeqAIJRegister(MATSEQAIJSELL, MatConvert_SeqAIJ_SeqAIJSELL)); 5390 #if defined(PETSC_HAVE_MKL_SPARSE) 5391 PetscCall(MatSeqAIJRegister(MATSEQAIJMKL, MatConvert_SeqAIJ_SeqAIJMKL)); 5392 #endif 5393 #if defined(PETSC_HAVE_CUDA) 5394 PetscCall(MatSeqAIJRegister(MATSEQAIJCUSPARSE, MatConvert_SeqAIJ_SeqAIJCUSPARSE)); 5395 #endif 5396 #if defined(PETSC_HAVE_KOKKOS_KERNELS) 5397 PetscCall(MatSeqAIJRegister(MATSEQAIJKOKKOS, MatConvert_SeqAIJ_SeqAIJKokkos)); 5398 #endif 5399 #if defined(PETSC_HAVE_VIENNACL) && defined(PETSC_HAVE_VIENNACL_NO_CUDA) 5400 PetscCall(MatSeqAIJRegister(MATMPIAIJVIENNACL, MatConvert_SeqAIJ_SeqAIJViennaCL)); 5401 #endif 5402 PetscFunctionReturn(0); 5403 } 5404 5405 /* 5406 Special version for direct calls from Fortran 5407 */ 5408 #include <petsc/private/fortranimpl.h> 5409 #if defined(PETSC_HAVE_FORTRAN_CAPS) 5410 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 5411 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 5412 #define matsetvaluesseqaij_ matsetvaluesseqaij 5413 #endif 5414 5415 /* Change these macros so can be used in void function */ 5416 5417 /* Change these macros so can be used in void function */ 5418 /* Identical to PetscCallVoid, except it assigns to *_ierr */ 5419 #undef PetscCall 5420 #define PetscCall(...) \ 5421 do { \ 5422 PetscErrorCode ierr_msv_mpiaij = __VA_ARGS__; \ 5423 if (PetscUnlikely(ierr_msv_mpiaij)) { \ 5424 *_ierr = PetscError(PETSC_COMM_SELF, __LINE__, PETSC_FUNCTION_NAME, __FILE__, ierr_msv_mpiaij, PETSC_ERROR_REPEAT, " "); \ 5425 return; \ 5426 } \ 5427 } while (0) 5428 5429 #undef SETERRQ 5430 #define SETERRQ(comm, ierr, ...) \ 5431 do { \ 5432 *_ierr = PetscError(comm, __LINE__, PETSC_FUNCTION_NAME, __FILE__, ierr, PETSC_ERROR_INITIAL, __VA_ARGS__); \ 5433 return; \ 5434 } while (0) 5435 5436 PETSC_EXTERN void matsetvaluesseqaij_(Mat *AA, PetscInt *mm, const PetscInt im[], PetscInt *nn, const PetscInt in[], const PetscScalar v[], InsertMode *isis, PetscErrorCode *_ierr) 5437 { 5438 Mat A = *AA; 5439 PetscInt m = *mm, n = *nn; 5440 InsertMode is = *isis; 5441 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 5442 PetscInt *rp, k, low, high, t, ii, row, nrow, i, col, l, rmax, N; 5443 PetscInt *imax, *ai, *ailen; 5444 PetscInt *aj, nonew = a->nonew, lastcol = -1; 5445 MatScalar *ap, value, *aa; 5446 PetscBool ignorezeroentries = a->ignorezeroentries; 5447 PetscBool roworiented = a->roworiented; 5448 5449 PetscFunctionBegin; 5450 MatCheckPreallocated(A, 1); 5451 imax = a->imax; 5452 ai = a->i; 5453 ailen = a->ilen; 5454 aj = a->j; 5455 aa = a->a; 5456 5457 for (k = 0; k < m; k++) { /* loop over added rows */ 5458 row = im[k]; 5459 if (row < 0) continue; 5460 PetscCheck(row < A->rmap->n, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Row too large"); 5461 rp = aj + ai[row]; 5462 ap = aa + ai[row]; 5463 rmax = imax[row]; 5464 nrow = ailen[row]; 5465 low = 0; 5466 high = nrow; 5467 for (l = 0; l < n; l++) { /* loop over added columns */ 5468 if (in[l] < 0) continue; 5469 PetscCheck(in[l] < A->cmap->n, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Column too large"); 5470 col = in[l]; 5471 if (roworiented) value = v[l + k * n]; 5472 else value = v[k + l * m]; 5473 5474 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 5475 5476 if (col <= lastcol) low = 0; 5477 else high = nrow; 5478 lastcol = col; 5479 while (high - low > 5) { 5480 t = (low + high) / 2; 5481 if (rp[t] > col) high = t; 5482 else low = t; 5483 } 5484 for (i = low; i < high; i++) { 5485 if (rp[i] > col) break; 5486 if (rp[i] == col) { 5487 if (is == ADD_VALUES) ap[i] += value; 5488 else ap[i] = value; 5489 goto noinsert; 5490 } 5491 } 5492 if (value == 0.0 && ignorezeroentries) goto noinsert; 5493 if (nonew == 1) goto noinsert; 5494 PetscCheck(nonew != -1, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Inserting a new nonzero in the matrix"); 5495 MatSeqXAIJReallocateAIJ(A, A->rmap->n, 1, nrow, row, col, rmax, aa, ai, aj, rp, ap, imax, nonew, MatScalar); 5496 N = nrow++ - 1; 5497 a->nz++; 5498 high++; 5499 /* shift up all the later entries in this row */ 5500 for (ii = N; ii >= i; ii--) { 5501 rp[ii + 1] = rp[ii]; 5502 ap[ii + 1] = ap[ii]; 5503 } 5504 rp[i] = col; 5505 ap[i] = value; 5506 A->nonzerostate++; 5507 noinsert:; 5508 low = i + 1; 5509 } 5510 ailen[row] = nrow; 5511 } 5512 PetscFunctionReturnVoid(); 5513 } 5514 /* Undefining these here since they were redefined from their original definition above! No 5515 * other PETSc functions should be defined past this point, as it is impossible to recover the 5516 * original definitions */ 5517 #undef PetscCall 5518 #undef SETERRQ 5519