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