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, m, n, *nidx, isz, val; 2813 const PetscInt *idx; 2814 PetscInt start, end, *ai, *aj; 2815 PetscBT table; 2816 2817 PetscFunctionBegin; 2818 m = A->rmap->n; 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 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 2837 for (j = 0; j < n; ++j) { 2838 if (!PetscBTLookupSet(table, idx[j])) nidx[isz++] = idx[j]; 2839 } 2840 PetscCall(ISRestoreIndices(is[i], &idx)); 2841 PetscCall(ISDestroy(&is[i])); 2842 2843 k = 0; 2844 for (j = 0; j < ov; j++) { /* for each overlap */ 2845 n = isz; 2846 for (; k < n; k++) { /* do only those rows in nidx[k], which are not done yet */ 2847 row = nidx[k]; 2848 start = ai[row]; 2849 end = ai[row + 1]; 2850 for (l = start; l < end; l++) { 2851 val = aj[l]; 2852 if (!PetscBTLookupSet(table, val)) nidx[isz++] = val; 2853 } 2854 } 2855 } 2856 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, isz, nidx, PETSC_COPY_VALUES, (is + i))); 2857 } 2858 PetscCall(PetscBTDestroy(&table)); 2859 PetscCall(PetscFree(nidx)); 2860 PetscFunctionReturn(0); 2861 } 2862 2863 /* -------------------------------------------------------------- */ 2864 PetscErrorCode MatPermute_SeqAIJ(Mat A, IS rowp, IS colp, Mat *B) 2865 { 2866 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 2867 PetscInt i, nz = 0, m = A->rmap->n, n = A->cmap->n; 2868 const PetscInt *row, *col; 2869 PetscInt *cnew, j, *lens; 2870 IS icolp, irowp; 2871 PetscInt *cwork = NULL; 2872 PetscScalar *vwork = NULL; 2873 2874 PetscFunctionBegin; 2875 PetscCall(ISInvertPermutation(rowp, PETSC_DECIDE, &irowp)); 2876 PetscCall(ISGetIndices(irowp, &row)); 2877 PetscCall(ISInvertPermutation(colp, PETSC_DECIDE, &icolp)); 2878 PetscCall(ISGetIndices(icolp, &col)); 2879 2880 /* determine lengths of permuted rows */ 2881 PetscCall(PetscMalloc1(m + 1, &lens)); 2882 for (i = 0; i < m; i++) lens[row[i]] = a->i[i + 1] - a->i[i]; 2883 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), B)); 2884 PetscCall(MatSetSizes(*B, m, n, m, n)); 2885 PetscCall(MatSetBlockSizesFromMats(*B, A, A)); 2886 PetscCall(MatSetType(*B, ((PetscObject)A)->type_name)); 2887 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(*B, 0, lens)); 2888 PetscCall(PetscFree(lens)); 2889 2890 PetscCall(PetscMalloc1(n, &cnew)); 2891 for (i = 0; i < m; i++) { 2892 PetscCall(MatGetRow_SeqAIJ(A, i, &nz, &cwork, &vwork)); 2893 for (j = 0; j < nz; j++) cnew[j] = col[cwork[j]]; 2894 PetscCall(MatSetValues_SeqAIJ(*B, 1, &row[i], nz, cnew, vwork, INSERT_VALUES)); 2895 PetscCall(MatRestoreRow_SeqAIJ(A, i, &nz, &cwork, &vwork)); 2896 } 2897 PetscCall(PetscFree(cnew)); 2898 2899 (*B)->assembled = PETSC_FALSE; 2900 2901 #if defined(PETSC_HAVE_DEVICE) 2902 PetscCall(MatBindToCPU(*B, A->boundtocpu)); 2903 #endif 2904 PetscCall(MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY)); 2905 PetscCall(MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY)); 2906 PetscCall(ISRestoreIndices(irowp, &row)); 2907 PetscCall(ISRestoreIndices(icolp, &col)); 2908 PetscCall(ISDestroy(&irowp)); 2909 PetscCall(ISDestroy(&icolp)); 2910 if (rowp == colp) PetscCall(MatPropagateSymmetryOptions(A, *B)); 2911 PetscFunctionReturn(0); 2912 } 2913 2914 PetscErrorCode MatCopy_SeqAIJ(Mat A, Mat B, MatStructure str) 2915 { 2916 PetscFunctionBegin; 2917 /* If the two matrices have the same copy implementation, use fast copy. */ 2918 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 2919 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 2920 Mat_SeqAIJ *b = (Mat_SeqAIJ *)B->data; 2921 const PetscScalar *aa; 2922 2923 PetscCall(MatSeqAIJGetArrayRead(A, &aa)); 2924 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]); 2925 PetscCall(PetscArraycpy(b->a, aa, a->i[A->rmap->n])); 2926 PetscCall(PetscObjectStateIncrease((PetscObject)B)); 2927 PetscCall(MatSeqAIJRestoreArrayRead(A, &aa)); 2928 } else { 2929 PetscCall(MatCopy_Basic(A, B, str)); 2930 } 2931 PetscFunctionReturn(0); 2932 } 2933 2934 PetscErrorCode MatSetUp_SeqAIJ(Mat A) 2935 { 2936 PetscFunctionBegin; 2937 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(A, PETSC_DEFAULT, NULL)); 2938 PetscFunctionReturn(0); 2939 } 2940 2941 PETSC_INTERN PetscErrorCode MatSeqAIJGetArray_SeqAIJ(Mat A, PetscScalar *array[]) 2942 { 2943 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 2944 2945 PetscFunctionBegin; 2946 *array = a->a; 2947 PetscFunctionReturn(0); 2948 } 2949 2950 PETSC_INTERN PetscErrorCode MatSeqAIJRestoreArray_SeqAIJ(Mat A, PetscScalar *array[]) 2951 { 2952 PetscFunctionBegin; 2953 *array = NULL; 2954 PetscFunctionReturn(0); 2955 } 2956 2957 /* 2958 Computes the number of nonzeros per row needed for preallocation when X and Y 2959 have different nonzero structure. 2960 */ 2961 PetscErrorCode MatAXPYGetPreallocation_SeqX_private(PetscInt m, const PetscInt *xi, const PetscInt *xj, const PetscInt *yi, const PetscInt *yj, PetscInt *nnz) 2962 { 2963 PetscInt i, j, k, nzx, nzy; 2964 2965 PetscFunctionBegin; 2966 /* Set the number of nonzeros in the new matrix */ 2967 for (i = 0; i < m; i++) { 2968 const PetscInt *xjj = xj + xi[i], *yjj = yj + yi[i]; 2969 nzx = xi[i + 1] - xi[i]; 2970 nzy = yi[i + 1] - yi[i]; 2971 nnz[i] = 0; 2972 for (j = 0, k = 0; j < nzx; j++) { /* Point in X */ 2973 for (; k < nzy && yjj[k] < xjj[j]; k++) nnz[i]++; /* Catch up to X */ 2974 if (k < nzy && yjj[k] == xjj[j]) k++; /* Skip duplicate */ 2975 nnz[i]++; 2976 } 2977 for (; k < nzy; k++) nnz[i]++; 2978 } 2979 PetscFunctionReturn(0); 2980 } 2981 2982 PetscErrorCode MatAXPYGetPreallocation_SeqAIJ(Mat Y, Mat X, PetscInt *nnz) 2983 { 2984 PetscInt m = Y->rmap->N; 2985 Mat_SeqAIJ *x = (Mat_SeqAIJ *)X->data; 2986 Mat_SeqAIJ *y = (Mat_SeqAIJ *)Y->data; 2987 2988 PetscFunctionBegin; 2989 /* Set the number of nonzeros in the new matrix */ 2990 PetscCall(MatAXPYGetPreallocation_SeqX_private(m, x->i, x->j, y->i, y->j, nnz)); 2991 PetscFunctionReturn(0); 2992 } 2993 2994 PetscErrorCode MatAXPY_SeqAIJ(Mat Y, PetscScalar a, Mat X, MatStructure str) 2995 { 2996 Mat_SeqAIJ *x = (Mat_SeqAIJ *)X->data, *y = (Mat_SeqAIJ *)Y->data; 2997 2998 PetscFunctionBegin; 2999 if (str == UNKNOWN_NONZERO_PATTERN || (PetscDefined(USE_DEBUG) && str == SAME_NONZERO_PATTERN)) { 3000 PetscBool e = x->nz == y->nz ? PETSC_TRUE : PETSC_FALSE; 3001 if (e) { 3002 PetscCall(PetscArraycmp(x->i, y->i, Y->rmap->n + 1, &e)); 3003 if (e) { 3004 PetscCall(PetscArraycmp(x->j, y->j, y->nz, &e)); 3005 if (e) str = SAME_NONZERO_PATTERN; 3006 } 3007 } 3008 if (!e) PetscCheck(str != SAME_NONZERO_PATTERN, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "MatStructure is not SAME_NONZERO_PATTERN"); 3009 } 3010 if (str == SAME_NONZERO_PATTERN) { 3011 const PetscScalar *xa; 3012 PetscScalar *ya, alpha = a; 3013 PetscBLASInt one = 1, bnz; 3014 3015 PetscCall(PetscBLASIntCast(x->nz, &bnz)); 3016 PetscCall(MatSeqAIJGetArray(Y, &ya)); 3017 PetscCall(MatSeqAIJGetArrayRead(X, &xa)); 3018 PetscCallBLAS("BLASaxpy", BLASaxpy_(&bnz, &alpha, xa, &one, ya, &one)); 3019 PetscCall(MatSeqAIJRestoreArrayRead(X, &xa)); 3020 PetscCall(MatSeqAIJRestoreArray(Y, &ya)); 3021 PetscCall(PetscLogFlops(2.0 * bnz)); 3022 PetscCall(MatSeqAIJInvalidateDiagonal(Y)); 3023 PetscCall(PetscObjectStateIncrease((PetscObject)Y)); 3024 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 3025 PetscCall(MatAXPY_Basic(Y, a, X, str)); 3026 } else { 3027 Mat B; 3028 PetscInt *nnz; 3029 PetscCall(PetscMalloc1(Y->rmap->N, &nnz)); 3030 PetscCall(MatCreate(PetscObjectComm((PetscObject)Y), &B)); 3031 PetscCall(PetscObjectSetName((PetscObject)B, ((PetscObject)Y)->name)); 3032 PetscCall(MatSetLayouts(B, Y->rmap, Y->cmap)); 3033 PetscCall(MatSetType(B, ((PetscObject)Y)->type_name)); 3034 PetscCall(MatAXPYGetPreallocation_SeqAIJ(Y, X, nnz)); 3035 PetscCall(MatSeqAIJSetPreallocation(B, 0, nnz)); 3036 PetscCall(MatAXPY_BasicWithPreallocation(B, Y, a, X, str)); 3037 PetscCall(MatHeaderMerge(Y, &B)); 3038 PetscCall(MatSeqAIJCheckInode(Y)); 3039 PetscCall(PetscFree(nnz)); 3040 } 3041 PetscFunctionReturn(0); 3042 } 3043 3044 PETSC_INTERN PetscErrorCode MatConjugate_SeqAIJ(Mat mat) 3045 { 3046 #if defined(PETSC_USE_COMPLEX) 3047 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 3048 PetscInt i, nz; 3049 PetscScalar *a; 3050 3051 PetscFunctionBegin; 3052 nz = aij->nz; 3053 PetscCall(MatSeqAIJGetArray(mat, &a)); 3054 for (i = 0; i < nz; i++) a[i] = PetscConj(a[i]); 3055 PetscCall(MatSeqAIJRestoreArray(mat, &a)); 3056 #else 3057 PetscFunctionBegin; 3058 #endif 3059 PetscFunctionReturn(0); 3060 } 3061 3062 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A, Vec v, PetscInt idx[]) 3063 { 3064 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 3065 PetscInt i, j, m = A->rmap->n, *ai, *aj, ncols, n; 3066 PetscReal atmp; 3067 PetscScalar *x; 3068 const MatScalar *aa, *av; 3069 3070 PetscFunctionBegin; 3071 PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix"); 3072 PetscCall(MatSeqAIJGetArrayRead(A, &av)); 3073 aa = av; 3074 ai = a->i; 3075 aj = a->j; 3076 3077 PetscCall(VecSet(v, 0.0)); 3078 PetscCall(VecGetArrayWrite(v, &x)); 3079 PetscCall(VecGetLocalSize(v, &n)); 3080 PetscCheck(n == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector"); 3081 for (i = 0; i < m; i++) { 3082 ncols = ai[1] - ai[0]; 3083 ai++; 3084 for (j = 0; j < ncols; j++) { 3085 atmp = PetscAbsScalar(*aa); 3086 if (PetscAbsScalar(x[i]) < atmp) { 3087 x[i] = atmp; 3088 if (idx) idx[i] = *aj; 3089 } 3090 aa++; 3091 aj++; 3092 } 3093 } 3094 PetscCall(VecRestoreArrayWrite(v, &x)); 3095 PetscCall(MatSeqAIJRestoreArrayRead(A, &av)); 3096 PetscFunctionReturn(0); 3097 } 3098 3099 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A, Vec v, PetscInt idx[]) 3100 { 3101 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 3102 PetscInt i, j, m = A->rmap->n, *ai, *aj, ncols, n; 3103 PetscScalar *x; 3104 const MatScalar *aa, *av; 3105 3106 PetscFunctionBegin; 3107 PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix"); 3108 PetscCall(MatSeqAIJGetArrayRead(A, &av)); 3109 aa = av; 3110 ai = a->i; 3111 aj = a->j; 3112 3113 PetscCall(VecSet(v, 0.0)); 3114 PetscCall(VecGetArrayWrite(v, &x)); 3115 PetscCall(VecGetLocalSize(v, &n)); 3116 PetscCheck(n == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector"); 3117 for (i = 0; i < m; i++) { 3118 ncols = ai[1] - ai[0]; 3119 ai++; 3120 if (ncols == A->cmap->n) { /* row is dense */ 3121 x[i] = *aa; 3122 if (idx) idx[i] = 0; 3123 } else { /* row is sparse so already KNOW maximum is 0.0 or higher */ 3124 x[i] = 0.0; 3125 if (idx) { 3126 for (j = 0; j < ncols; j++) { /* find first implicit 0.0 in the row */ 3127 if (aj[j] > j) { 3128 idx[i] = j; 3129 break; 3130 } 3131 } 3132 /* in case first implicit 0.0 in the row occurs at ncols-th column */ 3133 if (j == ncols && j < A->cmap->n) idx[i] = j; 3134 } 3135 } 3136 for (j = 0; j < ncols; j++) { 3137 if (PetscRealPart(x[i]) < PetscRealPart(*aa)) { 3138 x[i] = *aa; 3139 if (idx) idx[i] = *aj; 3140 } 3141 aa++; 3142 aj++; 3143 } 3144 } 3145 PetscCall(VecRestoreArrayWrite(v, &x)); 3146 PetscCall(MatSeqAIJRestoreArrayRead(A, &av)); 3147 PetscFunctionReturn(0); 3148 } 3149 3150 PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A, Vec v, PetscInt idx[]) 3151 { 3152 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 3153 PetscInt i, j, m = A->rmap->n, *ai, *aj, ncols, n; 3154 PetscScalar *x; 3155 const MatScalar *aa, *av; 3156 3157 PetscFunctionBegin; 3158 PetscCall(MatSeqAIJGetArrayRead(A, &av)); 3159 aa = av; 3160 ai = a->i; 3161 aj = a->j; 3162 3163 PetscCall(VecSet(v, 0.0)); 3164 PetscCall(VecGetArrayWrite(v, &x)); 3165 PetscCall(VecGetLocalSize(v, &n)); 3166 PetscCheck(n == m, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector, %" PetscInt_FMT " vs. %" PetscInt_FMT " rows", m, n); 3167 for (i = 0; i < m; i++) { 3168 ncols = ai[1] - ai[0]; 3169 ai++; 3170 if (ncols == A->cmap->n) { /* row is dense */ 3171 x[i] = *aa; 3172 if (idx) idx[i] = 0; 3173 } else { /* row is sparse so already KNOW minimum is 0.0 or higher */ 3174 x[i] = 0.0; 3175 if (idx) { /* find first implicit 0.0 in the row */ 3176 for (j = 0; j < ncols; j++) { 3177 if (aj[j] > j) { 3178 idx[i] = j; 3179 break; 3180 } 3181 } 3182 /* in case first implicit 0.0 in the row occurs at ncols-th column */ 3183 if (j == ncols && j < A->cmap->n) idx[i] = j; 3184 } 3185 } 3186 for (j = 0; j < ncols; j++) { 3187 if (PetscAbsScalar(x[i]) > PetscAbsScalar(*aa)) { 3188 x[i] = *aa; 3189 if (idx) idx[i] = *aj; 3190 } 3191 aa++; 3192 aj++; 3193 } 3194 } 3195 PetscCall(VecRestoreArrayWrite(v, &x)); 3196 PetscCall(MatSeqAIJRestoreArrayRead(A, &av)); 3197 PetscFunctionReturn(0); 3198 } 3199 3200 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A, Vec v, PetscInt idx[]) 3201 { 3202 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 3203 PetscInt i, j, m = A->rmap->n, ncols, n; 3204 const PetscInt *ai, *aj; 3205 PetscScalar *x; 3206 const MatScalar *aa, *av; 3207 3208 PetscFunctionBegin; 3209 PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix"); 3210 PetscCall(MatSeqAIJGetArrayRead(A, &av)); 3211 aa = av; 3212 ai = a->i; 3213 aj = a->j; 3214 3215 PetscCall(VecSet(v, 0.0)); 3216 PetscCall(VecGetArrayWrite(v, &x)); 3217 PetscCall(VecGetLocalSize(v, &n)); 3218 PetscCheck(n == m, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector"); 3219 for (i = 0; i < m; i++) { 3220 ncols = ai[1] - ai[0]; 3221 ai++; 3222 if (ncols == A->cmap->n) { /* row is dense */ 3223 x[i] = *aa; 3224 if (idx) idx[i] = 0; 3225 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 3226 x[i] = 0.0; 3227 if (idx) { /* find first implicit 0.0 in the row */ 3228 for (j = 0; j < ncols; j++) { 3229 if (aj[j] > j) { 3230 idx[i] = j; 3231 break; 3232 } 3233 } 3234 /* in case first implicit 0.0 in the row occurs at ncols-th column */ 3235 if (j == ncols && j < A->cmap->n) idx[i] = j; 3236 } 3237 } 3238 for (j = 0; j < ncols; j++) { 3239 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) { 3240 x[i] = *aa; 3241 if (idx) idx[i] = *aj; 3242 } 3243 aa++; 3244 aj++; 3245 } 3246 } 3247 PetscCall(VecRestoreArrayWrite(v, &x)); 3248 PetscCall(MatSeqAIJRestoreArrayRead(A, &av)); 3249 PetscFunctionReturn(0); 3250 } 3251 3252 PetscErrorCode MatInvertBlockDiagonal_SeqAIJ(Mat A, const PetscScalar **values) 3253 { 3254 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 3255 PetscInt i, bs = PetscAbs(A->rmap->bs), mbs = A->rmap->n / bs, ipvt[5], bs2 = bs * bs, *v_pivots, ij[7], *IJ, j; 3256 MatScalar *diag, work[25], *v_work; 3257 const PetscReal shift = 0.0; 3258 PetscBool allowzeropivot, zeropivotdetected = PETSC_FALSE; 3259 3260 PetscFunctionBegin; 3261 allowzeropivot = PetscNot(A->erroriffailure); 3262 if (a->ibdiagvalid) { 3263 if (values) *values = a->ibdiag; 3264 PetscFunctionReturn(0); 3265 } 3266 PetscCall(MatMarkDiagonal_SeqAIJ(A)); 3267 if (!a->ibdiag) { PetscCall(PetscMalloc1(bs2 * mbs, &a->ibdiag)); } 3268 diag = a->ibdiag; 3269 if (values) *values = a->ibdiag; 3270 /* factor and invert each block */ 3271 switch (bs) { 3272 case 1: 3273 for (i = 0; i < mbs; i++) { 3274 PetscCall(MatGetValues(A, 1, &i, 1, &i, diag + i)); 3275 if (PetscAbsScalar(diag[i] + shift) < PETSC_MACHINE_EPSILON) { 3276 if (allowzeropivot) { 3277 A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3278 A->factorerror_zeropivot_value = PetscAbsScalar(diag[i]); 3279 A->factorerror_zeropivot_row = i; 3280 PetscCall(PetscInfo(A, "Zero pivot, row %" PetscInt_FMT " pivot %g tolerance %g\n", i, (double)PetscAbsScalar(diag[i]), (double)PETSC_MACHINE_EPSILON)); 3281 } 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); 3282 } 3283 diag[i] = (PetscScalar)1.0 / (diag[i] + shift); 3284 } 3285 break; 3286 case 2: 3287 for (i = 0; i < mbs; i++) { 3288 ij[0] = 2 * i; 3289 ij[1] = 2 * i + 1; 3290 PetscCall(MatGetValues(A, 2, ij, 2, ij, diag)); 3291 PetscCall(PetscKernel_A_gets_inverse_A_2(diag, shift, allowzeropivot, &zeropivotdetected)); 3292 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3293 PetscCall(PetscKernel_A_gets_transpose_A_2(diag)); 3294 diag += 4; 3295 } 3296 break; 3297 case 3: 3298 for (i = 0; i < mbs; i++) { 3299 ij[0] = 3 * i; 3300 ij[1] = 3 * i + 1; 3301 ij[2] = 3 * i + 2; 3302 PetscCall(MatGetValues(A, 3, ij, 3, ij, diag)); 3303 PetscCall(PetscKernel_A_gets_inverse_A_3(diag, shift, allowzeropivot, &zeropivotdetected)); 3304 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3305 PetscCall(PetscKernel_A_gets_transpose_A_3(diag)); 3306 diag += 9; 3307 } 3308 break; 3309 case 4: 3310 for (i = 0; i < mbs; i++) { 3311 ij[0] = 4 * i; 3312 ij[1] = 4 * i + 1; 3313 ij[2] = 4 * i + 2; 3314 ij[3] = 4 * i + 3; 3315 PetscCall(MatGetValues(A, 4, ij, 4, ij, diag)); 3316 PetscCall(PetscKernel_A_gets_inverse_A_4(diag, shift, allowzeropivot, &zeropivotdetected)); 3317 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3318 PetscCall(PetscKernel_A_gets_transpose_A_4(diag)); 3319 diag += 16; 3320 } 3321 break; 3322 case 5: 3323 for (i = 0; i < mbs; i++) { 3324 ij[0] = 5 * i; 3325 ij[1] = 5 * i + 1; 3326 ij[2] = 5 * i + 2; 3327 ij[3] = 5 * i + 3; 3328 ij[4] = 5 * i + 4; 3329 PetscCall(MatGetValues(A, 5, ij, 5, ij, diag)); 3330 PetscCall(PetscKernel_A_gets_inverse_A_5(diag, ipvt, work, shift, allowzeropivot, &zeropivotdetected)); 3331 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3332 PetscCall(PetscKernel_A_gets_transpose_A_5(diag)); 3333 diag += 25; 3334 } 3335 break; 3336 case 6: 3337 for (i = 0; i < mbs; i++) { 3338 ij[0] = 6 * i; 3339 ij[1] = 6 * i + 1; 3340 ij[2] = 6 * i + 2; 3341 ij[3] = 6 * i + 3; 3342 ij[4] = 6 * i + 4; 3343 ij[5] = 6 * i + 5; 3344 PetscCall(MatGetValues(A, 6, ij, 6, ij, diag)); 3345 PetscCall(PetscKernel_A_gets_inverse_A_6(diag, shift, allowzeropivot, &zeropivotdetected)); 3346 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3347 PetscCall(PetscKernel_A_gets_transpose_A_6(diag)); 3348 diag += 36; 3349 } 3350 break; 3351 case 7: 3352 for (i = 0; i < mbs; i++) { 3353 ij[0] = 7 * i; 3354 ij[1] = 7 * i + 1; 3355 ij[2] = 7 * i + 2; 3356 ij[3] = 7 * i + 3; 3357 ij[4] = 7 * i + 4; 3358 ij[5] = 7 * i + 5; 3359 ij[5] = 7 * i + 6; 3360 PetscCall(MatGetValues(A, 7, ij, 7, ij, diag)); 3361 PetscCall(PetscKernel_A_gets_inverse_A_7(diag, shift, allowzeropivot, &zeropivotdetected)); 3362 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3363 PetscCall(PetscKernel_A_gets_transpose_A_7(diag)); 3364 diag += 49; 3365 } 3366 break; 3367 default: 3368 PetscCall(PetscMalloc3(bs, &v_work, bs, &v_pivots, bs, &IJ)); 3369 for (i = 0; i < mbs; i++) { 3370 for (j = 0; j < bs; j++) IJ[j] = bs * i + j; 3371 PetscCall(MatGetValues(A, bs, IJ, bs, IJ, diag)); 3372 PetscCall(PetscKernel_A_gets_inverse_A(bs, diag, v_pivots, v_work, allowzeropivot, &zeropivotdetected)); 3373 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3374 PetscCall(PetscKernel_A_gets_transpose_A_N(diag, bs)); 3375 diag += bs2; 3376 } 3377 PetscCall(PetscFree3(v_work, v_pivots, IJ)); 3378 } 3379 a->ibdiagvalid = PETSC_TRUE; 3380 PetscFunctionReturn(0); 3381 } 3382 3383 static PetscErrorCode MatSetRandom_SeqAIJ(Mat x, PetscRandom rctx) 3384 { 3385 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)x->data; 3386 PetscScalar a, *aa; 3387 PetscInt m, n, i, j, col; 3388 3389 PetscFunctionBegin; 3390 if (!x->assembled) { 3391 PetscCall(MatGetSize(x, &m, &n)); 3392 for (i = 0; i < m; i++) { 3393 for (j = 0; j < aij->imax[i]; j++) { 3394 PetscCall(PetscRandomGetValue(rctx, &a)); 3395 col = (PetscInt)(n * PetscRealPart(a)); 3396 PetscCall(MatSetValues(x, 1, &i, 1, &col, &a, ADD_VALUES)); 3397 } 3398 } 3399 } else { 3400 PetscCall(MatSeqAIJGetArrayWrite(x, &aa)); 3401 for (i = 0; i < aij->nz; i++) PetscCall(PetscRandomGetValue(rctx, aa + i)); 3402 PetscCall(MatSeqAIJRestoreArrayWrite(x, &aa)); 3403 } 3404 PetscCall(MatAssemblyBegin(x, MAT_FINAL_ASSEMBLY)); 3405 PetscCall(MatAssemblyEnd(x, MAT_FINAL_ASSEMBLY)); 3406 PetscFunctionReturn(0); 3407 } 3408 3409 /* Like MatSetRandom_SeqAIJ, but do not set values on columns in range of [low, high) */ 3410 PetscErrorCode MatSetRandomSkipColumnRange_SeqAIJ_Private(Mat x, PetscInt low, PetscInt high, PetscRandom rctx) 3411 { 3412 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)x->data; 3413 PetscScalar a; 3414 PetscInt m, n, i, j, col, nskip; 3415 3416 PetscFunctionBegin; 3417 nskip = high - low; 3418 PetscCall(MatGetSize(x, &m, &n)); 3419 n -= nskip; /* shrink number of columns where nonzeros can be set */ 3420 for (i = 0; i < m; i++) { 3421 for (j = 0; j < aij->imax[i]; j++) { 3422 PetscCall(PetscRandomGetValue(rctx, &a)); 3423 col = (PetscInt)(n * PetscRealPart(a)); 3424 if (col >= low) col += nskip; /* shift col rightward to skip the hole */ 3425 PetscCall(MatSetValues(x, 1, &i, 1, &col, &a, ADD_VALUES)); 3426 } 3427 } 3428 PetscCall(MatAssemblyBegin(x, MAT_FINAL_ASSEMBLY)); 3429 PetscCall(MatAssemblyEnd(x, MAT_FINAL_ASSEMBLY)); 3430 PetscFunctionReturn(0); 3431 } 3432 3433 /* -------------------------------------------------------------------*/ 3434 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 3435 MatGetRow_SeqAIJ, 3436 MatRestoreRow_SeqAIJ, 3437 MatMult_SeqAIJ, 3438 /* 4*/ MatMultAdd_SeqAIJ, 3439 MatMultTranspose_SeqAIJ, 3440 MatMultTransposeAdd_SeqAIJ, 3441 NULL, 3442 NULL, 3443 NULL, 3444 /* 10*/ NULL, 3445 MatLUFactor_SeqAIJ, 3446 NULL, 3447 MatSOR_SeqAIJ, 3448 MatTranspose_SeqAIJ, 3449 /*1 5*/ MatGetInfo_SeqAIJ, 3450 MatEqual_SeqAIJ, 3451 MatGetDiagonal_SeqAIJ, 3452 MatDiagonalScale_SeqAIJ, 3453 MatNorm_SeqAIJ, 3454 /* 20*/ NULL, 3455 MatAssemblyEnd_SeqAIJ, 3456 MatSetOption_SeqAIJ, 3457 MatZeroEntries_SeqAIJ, 3458 /* 24*/ MatZeroRows_SeqAIJ, 3459 NULL, 3460 NULL, 3461 NULL, 3462 NULL, 3463 /* 29*/ MatSetUp_SeqAIJ, 3464 NULL, 3465 NULL, 3466 NULL, 3467 NULL, 3468 /* 34*/ MatDuplicate_SeqAIJ, 3469 NULL, 3470 NULL, 3471 MatILUFactor_SeqAIJ, 3472 NULL, 3473 /* 39*/ MatAXPY_SeqAIJ, 3474 MatCreateSubMatrices_SeqAIJ, 3475 MatIncreaseOverlap_SeqAIJ, 3476 MatGetValues_SeqAIJ, 3477 MatCopy_SeqAIJ, 3478 /* 44*/ MatGetRowMax_SeqAIJ, 3479 MatScale_SeqAIJ, 3480 MatShift_SeqAIJ, 3481 MatDiagonalSet_SeqAIJ, 3482 MatZeroRowsColumns_SeqAIJ, 3483 /* 49*/ MatSetRandom_SeqAIJ, 3484 MatGetRowIJ_SeqAIJ, 3485 MatRestoreRowIJ_SeqAIJ, 3486 MatGetColumnIJ_SeqAIJ, 3487 MatRestoreColumnIJ_SeqAIJ, 3488 /* 54*/ MatFDColoringCreate_SeqXAIJ, 3489 NULL, 3490 NULL, 3491 MatPermute_SeqAIJ, 3492 NULL, 3493 /* 59*/ NULL, 3494 MatDestroy_SeqAIJ, 3495 MatView_SeqAIJ, 3496 NULL, 3497 NULL, 3498 /* 64*/ NULL, 3499 MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ, 3500 NULL, 3501 NULL, 3502 NULL, 3503 /* 69*/ MatGetRowMaxAbs_SeqAIJ, 3504 MatGetRowMinAbs_SeqAIJ, 3505 NULL, 3506 NULL, 3507 NULL, 3508 /* 74*/ NULL, 3509 MatFDColoringApply_AIJ, 3510 NULL, 3511 NULL, 3512 NULL, 3513 /* 79*/ MatFindZeroDiagonals_SeqAIJ, 3514 NULL, 3515 NULL, 3516 NULL, 3517 MatLoad_SeqAIJ, 3518 /* 84*/ MatIsSymmetric_SeqAIJ, 3519 MatIsHermitian_SeqAIJ, 3520 NULL, 3521 NULL, 3522 NULL, 3523 /* 89*/ NULL, 3524 NULL, 3525 MatMatMultNumeric_SeqAIJ_SeqAIJ, 3526 NULL, 3527 NULL, 3528 /* 94*/ MatPtAPNumeric_SeqAIJ_SeqAIJ_SparseAxpy, 3529 NULL, 3530 NULL, 3531 MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ, 3532 NULL, 3533 /* 99*/ MatProductSetFromOptions_SeqAIJ, 3534 NULL, 3535 NULL, 3536 MatConjugate_SeqAIJ, 3537 NULL, 3538 /*104*/ MatSetValuesRow_SeqAIJ, 3539 MatRealPart_SeqAIJ, 3540 MatImaginaryPart_SeqAIJ, 3541 NULL, 3542 NULL, 3543 /*109*/ MatMatSolve_SeqAIJ, 3544 NULL, 3545 MatGetRowMin_SeqAIJ, 3546 NULL, 3547 MatMissingDiagonal_SeqAIJ, 3548 /*114*/ NULL, 3549 NULL, 3550 NULL, 3551 NULL, 3552 NULL, 3553 /*119*/ NULL, 3554 NULL, 3555 NULL, 3556 NULL, 3557 MatGetMultiProcBlock_SeqAIJ, 3558 /*124*/ MatFindNonzeroRows_SeqAIJ, 3559 MatGetColumnReductions_SeqAIJ, 3560 MatInvertBlockDiagonal_SeqAIJ, 3561 MatInvertVariableBlockDiagonal_SeqAIJ, 3562 NULL, 3563 /*129*/ NULL, 3564 NULL, 3565 NULL, 3566 MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ, 3567 MatTransposeColoringCreate_SeqAIJ, 3568 /*134*/ MatTransColoringApplySpToDen_SeqAIJ, 3569 MatTransColoringApplyDenToSp_SeqAIJ, 3570 NULL, 3571 NULL, 3572 MatRARtNumeric_SeqAIJ_SeqAIJ, 3573 /*139*/ NULL, 3574 NULL, 3575 NULL, 3576 MatFDColoringSetUp_SeqXAIJ, 3577 MatFindOffBlockDiagonalEntries_SeqAIJ, 3578 MatCreateMPIMatConcatenateSeqMat_SeqAIJ, 3579 /*145*/ MatDestroySubMatrices_SeqAIJ, 3580 NULL, 3581 NULL, 3582 MatCreateGraph_Simple_AIJ, 3583 NULL, 3584 /*150*/ MatTransposeSymbolic_SeqAIJ}; 3585 3586 PetscErrorCode MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat, PetscInt *indices) 3587 { 3588 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 3589 PetscInt i, nz, n; 3590 3591 PetscFunctionBegin; 3592 nz = aij->maxnz; 3593 n = mat->rmap->n; 3594 for (i = 0; i < nz; i++) aij->j[i] = indices[i]; 3595 aij->nz = nz; 3596 for (i = 0; i < n; i++) aij->ilen[i] = aij->imax[i]; 3597 PetscFunctionReturn(0); 3598 } 3599 3600 /* 3601 * Given a sparse matrix with global column indices, compact it by using a local column space. 3602 * The result matrix helps saving memory in other algorithms, such as MatPtAPSymbolic_MPIAIJ_MPIAIJ_scalable() 3603 */ 3604 PetscErrorCode MatSeqAIJCompactOutExtraColumns_SeqAIJ(Mat mat, ISLocalToGlobalMapping *mapping) 3605 { 3606 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 3607 PetscTable gid1_lid1; 3608 PetscTablePosition tpos; 3609 PetscInt gid, lid, i, ec, nz = aij->nz; 3610 PetscInt *garray, *jj = aij->j; 3611 3612 PetscFunctionBegin; 3613 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 3614 PetscValidPointer(mapping, 2); 3615 /* use a table */ 3616 PetscCall(PetscTableCreate(mat->rmap->n, mat->cmap->N + 1, &gid1_lid1)); 3617 ec = 0; 3618 for (i = 0; i < nz; i++) { 3619 PetscInt data, gid1 = jj[i] + 1; 3620 PetscCall(PetscTableFind(gid1_lid1, gid1, &data)); 3621 if (!data) { 3622 /* one based table */ 3623 PetscCall(PetscTableAdd(gid1_lid1, gid1, ++ec, INSERT_VALUES)); 3624 } 3625 } 3626 /* form array of columns we need */ 3627 PetscCall(PetscMalloc1(ec, &garray)); 3628 PetscCall(PetscTableGetHeadPosition(gid1_lid1, &tpos)); 3629 while (tpos) { 3630 PetscCall(PetscTableGetNext(gid1_lid1, &tpos, &gid, &lid)); 3631 gid--; 3632 lid--; 3633 garray[lid] = gid; 3634 } 3635 PetscCall(PetscSortInt(ec, garray)); /* sort, and rebuild */ 3636 PetscCall(PetscTableRemoveAll(gid1_lid1)); 3637 for (i = 0; i < ec; i++) PetscCall(PetscTableAdd(gid1_lid1, garray[i] + 1, i + 1, INSERT_VALUES)); 3638 /* compact out the extra columns in B */ 3639 for (i = 0; i < nz; i++) { 3640 PetscInt gid1 = jj[i] + 1; 3641 PetscCall(PetscTableFind(gid1_lid1, gid1, &lid)); 3642 lid--; 3643 jj[i] = lid; 3644 } 3645 PetscCall(PetscLayoutDestroy(&mat->cmap)); 3646 PetscCall(PetscTableDestroy(&gid1_lid1)); 3647 PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)mat), ec, ec, 1, &mat->cmap)); 3648 PetscCall(ISLocalToGlobalMappingCreate(PETSC_COMM_SELF, mat->cmap->bs, mat->cmap->n, garray, PETSC_OWN_POINTER, mapping)); 3649 PetscCall(ISLocalToGlobalMappingSetType(*mapping, ISLOCALTOGLOBALMAPPINGHASH)); 3650 PetscFunctionReturn(0); 3651 } 3652 3653 /*@ 3654 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 3655 in the matrix. 3656 3657 Input Parameters: 3658 + mat - the `MATSEQAIJ` matrix 3659 - indices - the column indices 3660 3661 Level: advanced 3662 3663 Notes: 3664 This can be called if you have precomputed the nonzero structure of the 3665 matrix and want to provide it to the matrix object to improve the performance 3666 of the `MatSetValues()` operation. 3667 3668 You MUST have set the correct numbers of nonzeros per row in the call to 3669 `MatCreateSeqAIJ()`, and the columns indices MUST be sorted. 3670 3671 MUST be called before any calls to `MatSetValues()` 3672 3673 The indices should start with zero, not one. 3674 3675 @*/ 3676 PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat, PetscInt *indices) 3677 { 3678 PetscFunctionBegin; 3679 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 3680 PetscValidIntPointer(indices, 2); 3681 PetscUseMethod(mat, "MatSeqAIJSetColumnIndices_C", (Mat, PetscInt *), (mat, indices)); 3682 PetscFunctionReturn(0); 3683 } 3684 3685 /* ----------------------------------------------------------------------------------------*/ 3686 3687 PetscErrorCode MatStoreValues_SeqAIJ(Mat mat) 3688 { 3689 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 3690 size_t nz = aij->i[mat->rmap->n]; 3691 3692 PetscFunctionBegin; 3693 PetscCheck(aij->nonew, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3694 3695 /* allocate space for values if not already there */ 3696 if (!aij->saved_values) { PetscCall(PetscMalloc1(nz + 1, &aij->saved_values)); } 3697 3698 /* copy values over */ 3699 PetscCall(PetscArraycpy(aij->saved_values, aij->a, nz)); 3700 PetscFunctionReturn(0); 3701 } 3702 3703 /*@ 3704 MatStoreValues - Stashes a copy of the matrix values; this allows, for 3705 example, reuse of the linear part of a Jacobian, while recomputing the 3706 nonlinear portion. 3707 3708 Collect on mat 3709 3710 Input Parameters: 3711 . mat - the matrix (currently only `MATAIJ` matrices support this option) 3712 3713 Level: advanced 3714 3715 Common Usage, with `SNESSolve()`: 3716 $ Create Jacobian matrix 3717 $ Set linear terms into matrix 3718 $ Apply boundary conditions to matrix, at this time matrix must have 3719 $ final nonzero structure (i.e. setting the nonlinear terms and applying 3720 $ boundary conditions again will not change the nonzero structure 3721 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3722 $ ierr = MatStoreValues(mat); 3723 $ Call SNESSetJacobian() with matrix 3724 $ In your Jacobian routine 3725 $ ierr = MatRetrieveValues(mat); 3726 $ Set nonlinear terms in matrix 3727 3728 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 3729 $ // build linear portion of Jacobian 3730 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3731 $ ierr = MatStoreValues(mat); 3732 $ loop over nonlinear iterations 3733 $ ierr = MatRetrieveValues(mat); 3734 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 3735 $ // call MatAssemblyBegin/End() on matrix 3736 $ Solve linear system with Jacobian 3737 $ endloop 3738 3739 Notes: 3740 Matrix must already be assemblied before calling this routine 3741 Must set the matrix option `MatSetOption`(mat,`MAT_NEW_NONZERO_LOCATIONS`,`PETSC_FALSE`); before 3742 calling this routine. 3743 3744 When this is called multiple times it overwrites the previous set of stored values 3745 and does not allocated additional space. 3746 3747 .seealso: `MatRetrieveValues()` 3748 @*/ 3749 PetscErrorCode MatStoreValues(Mat mat) 3750 { 3751 PetscFunctionBegin; 3752 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 3753 PetscCheck(mat->assembled, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix"); 3754 PetscCheck(!mat->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix"); 3755 PetscUseMethod(mat, "MatStoreValues_C", (Mat), (mat)); 3756 PetscFunctionReturn(0); 3757 } 3758 3759 PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat) 3760 { 3761 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 3762 PetscInt nz = aij->i[mat->rmap->n]; 3763 3764 PetscFunctionBegin; 3765 PetscCheck(aij->nonew, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3766 PetscCheck(aij->saved_values, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call MatStoreValues(A);first"); 3767 /* copy values over */ 3768 PetscCall(PetscArraycpy(aij->a, aij->saved_values, nz)); 3769 PetscFunctionReturn(0); 3770 } 3771 3772 /*@ 3773 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 3774 example, reuse of the linear part of a Jacobian, while recomputing the 3775 nonlinear portion. 3776 3777 Collect on mat 3778 3779 Input Parameters: 3780 . mat - the matrix (currently only `MATAIJ` matrices support this option) 3781 3782 Level: advanced 3783 3784 .seealso: `MatStoreValues()` 3785 @*/ 3786 PetscErrorCode MatRetrieveValues(Mat mat) 3787 { 3788 PetscFunctionBegin; 3789 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 3790 PetscCheck(mat->assembled, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix"); 3791 PetscCheck(!mat->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix"); 3792 PetscUseMethod(mat, "MatRetrieveValues_C", (Mat), (mat)); 3793 PetscFunctionReturn(0); 3794 } 3795 3796 /* --------------------------------------------------------------------------------*/ 3797 /*@C 3798 MatCreateSeqAIJ - Creates a sparse matrix in `MATSEQAIJ` (compressed row) format 3799 (the default parallel PETSc format). For good matrix assembly performance 3800 the user should preallocate the matrix storage by setting the parameter nz 3801 (or the array nnz). By setting these parameters accurately, performance 3802 during matrix assembly can be increased by more than a factor of 50. 3803 3804 Collective 3805 3806 Input Parameters: 3807 + comm - MPI communicator, set to `PETSC_COMM_SELF` 3808 . m - number of rows 3809 . n - number of columns 3810 . nz - number of nonzeros per row (same for all rows) 3811 - nnz - array containing the number of nonzeros in the various rows 3812 (possibly different for each row) or NULL 3813 3814 Output Parameter: 3815 . A - the matrix 3816 3817 It is recommended that one use the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`, 3818 MatXXXXSetPreallocation() paradigm instead of this routine directly. 3819 [MatXXXXSetPreallocation() is, for example, `MatSeqAIJSetPreallocation()`] 3820 3821 Notes: 3822 If nnz is given then nz is ignored 3823 3824 The AIJ format, also called 3825 compressed row storage, is fully compatible with standard Fortran 77 3826 storage. That is, the stored row and column indices can begin at 3827 either one (as in Fortran) or zero. See the users' manual for details. 3828 3829 Specify the preallocated storage with either nz or nnz (not both). 3830 Set nz = `PETSC_DEFAULT` and nnz = NULL for PETSc to control dynamic memory 3831 allocation. For large problems you MUST preallocate memory or you 3832 will get TERRIBLE performance, see the users' manual chapter on matrices. 3833 3834 By default, this format uses inodes (identical nodes) when possible, to 3835 improve numerical efficiency of matrix-vector products and solves. We 3836 search for consecutive rows with the same nonzero structure, thereby 3837 reusing matrix information to achieve increased efficiency. 3838 3839 Options Database Keys: 3840 + -mat_no_inode - Do not use inodes 3841 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3842 3843 Level: intermediate 3844 3845 .seealso: [Sparse Matrix Creation](sec_matsparse), `MatCreate()`, `MatCreateAIJ()`, `MatSetValues()`, `MatSeqAIJSetColumnIndices()`, `MatCreateSeqAIJWithArrays()` 3846 @*/ 3847 PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt nz, const PetscInt nnz[], Mat *A) 3848 { 3849 PetscFunctionBegin; 3850 PetscCall(MatCreate(comm, A)); 3851 PetscCall(MatSetSizes(*A, m, n, m, n)); 3852 PetscCall(MatSetType(*A, MATSEQAIJ)); 3853 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(*A, nz, nnz)); 3854 PetscFunctionReturn(0); 3855 } 3856 3857 /*@C 3858 MatSeqAIJSetPreallocation - For good matrix assembly performance 3859 the user should preallocate the matrix storage by setting the parameter nz 3860 (or the array nnz). By setting these parameters accurately, performance 3861 during matrix assembly can be increased by more than a factor of 50. 3862 3863 Collective 3864 3865 Input Parameters: 3866 + B - The matrix 3867 . nz - number of nonzeros per row (same for all rows) 3868 - nnz - array containing the number of nonzeros in the various rows 3869 (possibly different for each row) or NULL 3870 3871 Notes: 3872 If nnz is given then nz is ignored 3873 3874 The `MATSEQAIJ` format also called 3875 compressed row storage, is fully compatible with standard Fortran 77 3876 storage. That is, the stored row and column indices can begin at 3877 either one (as in Fortran) or zero. See the users' manual for details. 3878 3879 Specify the preallocated storage with either nz or nnz (not both). 3880 Set nz = `PETSC_DEFAULT` and nnz = NULL for PETSc to control dynamic memory 3881 allocation. For large problems you MUST preallocate memory or you 3882 will get TERRIBLE performance, see the users' manual chapter on matrices. 3883 3884 You can call `MatGetInfo()` to get information on how effective the preallocation was; 3885 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 3886 You can also run with the option -info and look for messages with the string 3887 malloc in them to see if additional memory allocation was needed. 3888 3889 Developer Notes: 3890 Use nz of `MAT_SKIP_ALLOCATION` to not allocate any space for the matrix 3891 entries or columns indices 3892 3893 By default, this format uses inodes (identical nodes) when possible, to 3894 improve numerical efficiency of matrix-vector products and solves. We 3895 search for consecutive rows with the same nonzero structure, thereby 3896 reusing matrix information to achieve increased efficiency. 3897 3898 Options Database Keys: 3899 + -mat_no_inode - Do not use inodes 3900 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3901 3902 Level: intermediate 3903 3904 .seealso: `MatCreate()`, `MatCreateAIJ()`, `MatSetValues()`, `MatSeqAIJSetColumnIndices()`, `MatCreateSeqAIJWithArrays()`, `MatGetInfo()`, 3905 `MatSeqAIJSetTotalPreallocation()` 3906 @*/ 3907 PetscErrorCode MatSeqAIJSetPreallocation(Mat B, PetscInt nz, const PetscInt nnz[]) 3908 { 3909 PetscFunctionBegin; 3910 PetscValidHeaderSpecific(B, MAT_CLASSID, 1); 3911 PetscValidType(B, 1); 3912 PetscTryMethod(B, "MatSeqAIJSetPreallocation_C", (Mat, PetscInt, const PetscInt[]), (B, nz, nnz)); 3913 PetscFunctionReturn(0); 3914 } 3915 3916 PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B, PetscInt nz, const PetscInt *nnz) 3917 { 3918 Mat_SeqAIJ *b; 3919 PetscBool skipallocation = PETSC_FALSE, realalloc = PETSC_FALSE; 3920 PetscInt i; 3921 3922 PetscFunctionBegin; 3923 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 3924 if (nz == MAT_SKIP_ALLOCATION) { 3925 skipallocation = PETSC_TRUE; 3926 nz = 0; 3927 } 3928 PetscCall(PetscLayoutSetUp(B->rmap)); 3929 PetscCall(PetscLayoutSetUp(B->cmap)); 3930 3931 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3932 PetscCheck(nz >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "nz cannot be less than 0: value %" PetscInt_FMT, nz); 3933 if (PetscUnlikelyDebug(nnz)) { 3934 for (i = 0; i < B->rmap->n; i++) { 3935 PetscCheck(nnz[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "nnz cannot be less than 0: local row %" PetscInt_FMT " value %" PetscInt_FMT, i, nnz[i]); 3936 PetscCheck(nnz[i] <= B->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "nnz cannot be greater than row length: local row %" PetscInt_FMT " value %" PetscInt_FMT " rowlength %" PetscInt_FMT, i, nnz[i], B->cmap->n); 3937 } 3938 } 3939 3940 B->preallocated = PETSC_TRUE; 3941 3942 b = (Mat_SeqAIJ *)B->data; 3943 3944 if (!skipallocation) { 3945 if (!b->imax) { PetscCall(PetscMalloc1(B->rmap->n, &b->imax)); } 3946 if (!b->ilen) { 3947 /* b->ilen will count nonzeros in each row so far. */ 3948 PetscCall(PetscCalloc1(B->rmap->n, &b->ilen)); 3949 } else { 3950 PetscCall(PetscMemzero(b->ilen, B->rmap->n * sizeof(PetscInt))); 3951 } 3952 if (!b->ipre) { PetscCall(PetscMalloc1(B->rmap->n, &b->ipre)); } 3953 if (!nnz) { 3954 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 3955 else if (nz < 0) nz = 1; 3956 nz = PetscMin(nz, B->cmap->n); 3957 for (i = 0; i < B->rmap->n; i++) b->imax[i] = nz; 3958 nz = nz * B->rmap->n; 3959 } else { 3960 PetscInt64 nz64 = 0; 3961 for (i = 0; i < B->rmap->n; i++) { 3962 b->imax[i] = nnz[i]; 3963 nz64 += nnz[i]; 3964 } 3965 PetscCall(PetscIntCast(nz64, &nz)); 3966 } 3967 3968 /* allocate the matrix space */ 3969 /* FIXME: should B's old memory be unlogged? */ 3970 PetscCall(MatSeqXAIJFreeAIJ(B, &b->a, &b->j, &b->i)); 3971 if (B->structure_only) { 3972 PetscCall(PetscMalloc1(nz, &b->j)); 3973 PetscCall(PetscMalloc1(B->rmap->n + 1, &b->i)); 3974 } else { 3975 PetscCall(PetscMalloc3(nz, &b->a, nz, &b->j, B->rmap->n + 1, &b->i)); 3976 } 3977 b->i[0] = 0; 3978 for (i = 1; i < B->rmap->n + 1; i++) b->i[i] = b->i[i - 1] + b->imax[i - 1]; 3979 if (B->structure_only) { 3980 b->singlemalloc = PETSC_FALSE; 3981 b->free_a = PETSC_FALSE; 3982 } else { 3983 b->singlemalloc = PETSC_TRUE; 3984 b->free_a = PETSC_TRUE; 3985 } 3986 b->free_ij = PETSC_TRUE; 3987 } else { 3988 b->free_a = PETSC_FALSE; 3989 b->free_ij = PETSC_FALSE; 3990 } 3991 3992 if (b->ipre && nnz != b->ipre && b->imax) { 3993 /* reserve user-requested sparsity */ 3994 PetscCall(PetscArraycpy(b->ipre, b->imax, B->rmap->n)); 3995 } 3996 3997 b->nz = 0; 3998 b->maxnz = nz; 3999 B->info.nz_unneeded = (double)b->maxnz; 4000 if (realalloc) PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE)); 4001 B->was_assembled = PETSC_FALSE; 4002 B->assembled = PETSC_FALSE; 4003 /* We simply deem preallocation has changed nonzero state. Updating the state 4004 will give clients (like AIJKokkos) a chance to know something has happened. 4005 */ 4006 B->nonzerostate++; 4007 PetscFunctionReturn(0); 4008 } 4009 4010 PetscErrorCode MatResetPreallocation_SeqAIJ(Mat A) 4011 { 4012 Mat_SeqAIJ *a; 4013 PetscInt i; 4014 4015 PetscFunctionBegin; 4016 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 4017 4018 /* Check local size. If zero, then return */ 4019 if (!A->rmap->n) PetscFunctionReturn(0); 4020 4021 a = (Mat_SeqAIJ *)A->data; 4022 /* if no saved info, we error out */ 4023 PetscCheck(a->ipre, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "No saved preallocation info "); 4024 4025 PetscCheck(a->i && a->j && a->a && a->imax && a->ilen, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Memory info is incomplete, and can not reset preallocation "); 4026 4027 PetscCall(PetscArraycpy(a->imax, a->ipre, A->rmap->n)); 4028 PetscCall(PetscArrayzero(a->ilen, A->rmap->n)); 4029 a->i[0] = 0; 4030 for (i = 1; i < A->rmap->n + 1; i++) a->i[i] = a->i[i - 1] + a->imax[i - 1]; 4031 A->preallocated = PETSC_TRUE; 4032 a->nz = 0; 4033 a->maxnz = a->i[A->rmap->n]; 4034 A->info.nz_unneeded = (double)a->maxnz; 4035 A->was_assembled = PETSC_FALSE; 4036 A->assembled = PETSC_FALSE; 4037 PetscFunctionReturn(0); 4038 } 4039 4040 /*@ 4041 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in `MATSEQAIJ` format. 4042 4043 Input Parameters: 4044 + B - the matrix 4045 . i - the indices into j for the start of each row (starts with zero) 4046 . j - the column indices for each row (starts with zero) these must be sorted for each row 4047 - v - optional values in the matrix 4048 4049 Level: developer 4050 4051 Notes: 4052 The i,j,v values are COPIED with this routine; to avoid the copy use `MatCreateSeqAIJWithArrays()` 4053 4054 This routine may be called multiple times with different nonzero patterns (or the same nonzero pattern). The nonzero 4055 structure will be the union of all the previous nonzero structures. 4056 4057 Developer Notes: 4058 An optimization could be added to the implementation where it checks if the i, and j are identical to the current i and j and 4059 then just copies the v values directly with `PetscMemcpy()`. 4060 4061 This routine could also take a `PetscCopyMode` argument to allow sharing the values instead of always copying them. 4062 4063 .seealso: `MatCreate()`, `MatCreateSeqAIJ()`, `MatSetValues()`, `MatSeqAIJSetPreallocation()`, `MatCreateSeqAIJ()`, `MATSEQAIJ`, `MatResetPreallocation()` 4064 @*/ 4065 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B, const PetscInt i[], const PetscInt j[], const PetscScalar v[]) 4066 { 4067 PetscFunctionBegin; 4068 PetscValidHeaderSpecific(B, MAT_CLASSID, 1); 4069 PetscValidType(B, 1); 4070 PetscTryMethod(B, "MatSeqAIJSetPreallocationCSR_C", (Mat, const PetscInt[], const PetscInt[], const PetscScalar[]), (B, i, j, v)); 4071 PetscFunctionReturn(0); 4072 } 4073 4074 PetscErrorCode MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B, const PetscInt Ii[], const PetscInt J[], const PetscScalar v[]) 4075 { 4076 PetscInt i; 4077 PetscInt m, n; 4078 PetscInt nz; 4079 PetscInt *nnz; 4080 4081 PetscFunctionBegin; 4082 PetscCheck(Ii[0] == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %" PetscInt_FMT, Ii[0]); 4083 4084 PetscCall(PetscLayoutSetUp(B->rmap)); 4085 PetscCall(PetscLayoutSetUp(B->cmap)); 4086 4087 PetscCall(MatGetSize(B, &m, &n)); 4088 PetscCall(PetscMalloc1(m + 1, &nnz)); 4089 for (i = 0; i < m; i++) { 4090 nz = Ii[i + 1] - Ii[i]; 4091 PetscCheck(nz >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local row %" PetscInt_FMT " has a negative number of columns %" PetscInt_FMT, i, nz); 4092 nnz[i] = nz; 4093 } 4094 PetscCall(MatSeqAIJSetPreallocation(B, 0, nnz)); 4095 PetscCall(PetscFree(nnz)); 4096 4097 for (i = 0; i < m; i++) PetscCall(MatSetValues_SeqAIJ(B, 1, &i, Ii[i + 1] - Ii[i], J + Ii[i], v ? v + Ii[i] : NULL, INSERT_VALUES)); 4098 4099 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 4100 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 4101 4102 PetscCall(MatSetOption(B, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE)); 4103 PetscFunctionReturn(0); 4104 } 4105 4106 /*@ 4107 MatSeqAIJKron - Computes C, the Kronecker product of A and B. 4108 4109 Input Parameters: 4110 + A - left-hand side matrix 4111 . B - right-hand side matrix 4112 - reuse - either `MAT_INITIAL_MATRIX` or `MAT_REUSE_MATRIX` 4113 4114 Output Parameter: 4115 . C - Kronecker product of A and B 4116 4117 Level: intermediate 4118 4119 Note: 4120 `MAT_REUSE_MATRIX` can only be used when the nonzero structure of the product matrix has not changed from that last call to `MatSeqAIJKron()`. 4121 4122 .seealso: `MatCreateSeqAIJ()`, `MATSEQAIJ`, `MATKAIJ`, `MatReuse` 4123 @*/ 4124 PetscErrorCode MatSeqAIJKron(Mat A, Mat B, MatReuse reuse, Mat *C) 4125 { 4126 PetscFunctionBegin; 4127 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 4128 PetscValidType(A, 1); 4129 PetscValidHeaderSpecific(B, MAT_CLASSID, 2); 4130 PetscValidType(B, 2); 4131 PetscValidPointer(C, 4); 4132 if (reuse == MAT_REUSE_MATRIX) { 4133 PetscValidHeaderSpecific(*C, MAT_CLASSID, 4); 4134 PetscValidType(*C, 4); 4135 } 4136 PetscTryMethod(A, "MatSeqAIJKron_C", (Mat, Mat, MatReuse, Mat *), (A, B, reuse, C)); 4137 PetscFunctionReturn(0); 4138 } 4139 4140 PetscErrorCode MatSeqAIJKron_SeqAIJ(Mat A, Mat B, MatReuse reuse, Mat *C) 4141 { 4142 Mat newmat; 4143 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 4144 Mat_SeqAIJ *b = (Mat_SeqAIJ *)B->data; 4145 PetscScalar *v; 4146 const PetscScalar *aa, *ba; 4147 PetscInt *i, *j, m, n, p, q, nnz = 0, am = A->rmap->n, bm = B->rmap->n, an = A->cmap->n, bn = B->cmap->n; 4148 PetscBool flg; 4149 4150 PetscFunctionBegin; 4151 PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix"); 4152 PetscCheck(A->assembled, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix"); 4153 PetscCheck(!B->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix"); 4154 PetscCheck(B->assembled, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix"); 4155 PetscCall(PetscObjectTypeCompare((PetscObject)B, MATSEQAIJ, &flg)); 4156 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "MatType %s", ((PetscObject)B)->type_name); 4157 PetscCheck(reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX, PETSC_COMM_SELF, PETSC_ERR_SUP, "MatReuse %d", (int)reuse); 4158 if (reuse == MAT_INITIAL_MATRIX) { 4159 PetscCall(PetscMalloc2(am * bm + 1, &i, a->i[am] * b->i[bm], &j)); 4160 PetscCall(MatCreate(PETSC_COMM_SELF, &newmat)); 4161 PetscCall(MatSetSizes(newmat, am * bm, an * bn, am * bm, an * bn)); 4162 PetscCall(MatSetType(newmat, MATAIJ)); 4163 i[0] = 0; 4164 for (m = 0; m < am; ++m) { 4165 for (p = 0; p < bm; ++p) { 4166 i[m * bm + p + 1] = i[m * bm + p] + (a->i[m + 1] - a->i[m]) * (b->i[p + 1] - b->i[p]); 4167 for (n = a->i[m]; n < a->i[m + 1]; ++n) { 4168 for (q = b->i[p]; q < b->i[p + 1]; ++q) j[nnz++] = a->j[n] * bn + b->j[q]; 4169 } 4170 } 4171 } 4172 PetscCall(MatSeqAIJSetPreallocationCSR(newmat, i, j, NULL)); 4173 *C = newmat; 4174 PetscCall(PetscFree2(i, j)); 4175 nnz = 0; 4176 } 4177 PetscCall(MatSeqAIJGetArray(*C, &v)); 4178 PetscCall(MatSeqAIJGetArrayRead(A, &aa)); 4179 PetscCall(MatSeqAIJGetArrayRead(B, &ba)); 4180 for (m = 0; m < am; ++m) { 4181 for (p = 0; p < bm; ++p) { 4182 for (n = a->i[m]; n < a->i[m + 1]; ++n) { 4183 for (q = b->i[p]; q < b->i[p + 1]; ++q) v[nnz++] = aa[n] * ba[q]; 4184 } 4185 } 4186 } 4187 PetscCall(MatSeqAIJRestoreArray(*C, &v)); 4188 PetscCall(MatSeqAIJRestoreArrayRead(A, &aa)); 4189 PetscCall(MatSeqAIJRestoreArrayRead(B, &ba)); 4190 PetscFunctionReturn(0); 4191 } 4192 4193 #include <../src/mat/impls/dense/seq/dense.h> 4194 #include <petsc/private/kernels/petscaxpy.h> 4195 4196 /* 4197 Computes (B'*A')' since computing B*A directly is untenable 4198 4199 n p p 4200 [ ] [ ] [ ] 4201 m [ A ] * n [ B ] = m [ C ] 4202 [ ] [ ] [ ] 4203 4204 */ 4205 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A, Mat B, Mat C) 4206 { 4207 Mat_SeqDense *sub_a = (Mat_SeqDense *)A->data; 4208 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ *)B->data; 4209 Mat_SeqDense *sub_c = (Mat_SeqDense *)C->data; 4210 PetscInt i, j, n, m, q, p; 4211 const PetscInt *ii, *idx; 4212 const PetscScalar *b, *a, *a_q; 4213 PetscScalar *c, *c_q; 4214 PetscInt clda = sub_c->lda; 4215 PetscInt alda = sub_a->lda; 4216 4217 PetscFunctionBegin; 4218 m = A->rmap->n; 4219 n = A->cmap->n; 4220 p = B->cmap->n; 4221 a = sub_a->v; 4222 b = sub_b->a; 4223 c = sub_c->v; 4224 if (clda == m) { 4225 PetscCall(PetscArrayzero(c, m * p)); 4226 } else { 4227 for (j = 0; j < p; j++) 4228 for (i = 0; i < m; i++) c[j * clda + i] = 0.0; 4229 } 4230 ii = sub_b->i; 4231 idx = sub_b->j; 4232 for (i = 0; i < n; i++) { 4233 q = ii[i + 1] - ii[i]; 4234 while (q-- > 0) { 4235 c_q = c + clda * (*idx); 4236 a_q = a + alda * i; 4237 PetscKernelAXPY(c_q, *b, a_q, m); 4238 idx++; 4239 b++; 4240 } 4241 } 4242 PetscFunctionReturn(0); 4243 } 4244 4245 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A, Mat B, PetscReal fill, Mat C) 4246 { 4247 PetscInt m = A->rmap->n, n = B->cmap->n; 4248 PetscBool cisdense; 4249 4250 PetscFunctionBegin; 4251 PetscCheck(A->cmap->n == B->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "A->cmap->n %" PetscInt_FMT " != B->rmap->n %" PetscInt_FMT, A->cmap->n, B->rmap->n); 4252 PetscCall(MatSetSizes(C, m, n, m, n)); 4253 PetscCall(MatSetBlockSizesFromMats(C, A, B)); 4254 PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSECUDA, "")); 4255 if (!cisdense) PetscCall(MatSetType(C, MATDENSE)); 4256 PetscCall(MatSetUp(C)); 4257 4258 C->ops->matmultnumeric = MatMatMultNumeric_SeqDense_SeqAIJ; 4259 PetscFunctionReturn(0); 4260 } 4261 4262 /* ----------------------------------------------------------------*/ 4263 /*MC 4264 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 4265 based on compressed sparse row format. 4266 4267 Options Database Keys: 4268 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 4269 4270 Level: beginner 4271 4272 Notes: 4273 `MatSetValues()` may be called for this matrix type with a NULL argument for the numerical values, 4274 in this case the values associated with the rows and columns one passes in are set to zero 4275 in the matrix 4276 4277 `MatSetOptions`(,`MAT_STRUCTURE_ONLY`,`PETSC_TRUE`) may be called for this matrix type. In this no 4278 space is allocated for the nonzero entries and any entries passed with `MatSetValues()` are ignored 4279 4280 Developer Note: 4281 It would be nice if all matrix formats supported passing NULL in for the numerical values 4282 4283 .seealso: `MatCreateSeqAIJ()`, `MatSetFromOptions()`, `MatSetType()`, `MatCreate()`, `MatType`, `MATSELL`, `MATSEQSELL`, `MATMPISELL` 4284 M*/ 4285 4286 /*MC 4287 MATAIJ - MATAIJ = "aij" - A matrix type to be used for sparse matrices. 4288 4289 This matrix type is identical to `MATSEQAIJ` when constructed with a single process communicator, 4290 and `MATMPIAIJ` otherwise. As a result, for single process communicators, 4291 `MatSeqAIJSetPreallocation()` is supported, and similarly `MatMPIAIJSetPreallocation()` is supported 4292 for communicators controlling multiple processes. It is recommended that you call both of 4293 the above preallocation routines for simplicity. 4294 4295 Options Database Keys: 4296 . -mat_type aij - sets the matrix type to "aij" during a call to `MatSetFromOptions()` 4297 4298 Note: 4299 Subclasses include `MATAIJCUSPARSE`, `MATAIJPERM`, `MATAIJSELL`, `MATAIJMKL`, `MATAIJCRL`, and also automatically switches over to use inodes when 4300 enough exist. 4301 4302 Level: beginner 4303 4304 .seealso: `MatCreateAIJ()`, `MatCreateSeqAIJ()`, `MATSEQAIJ`, `MATMPIAIJ`, `MATSELL`, `MATSEQSELL`, `MATMPISELL` 4305 M*/ 4306 4307 /*MC 4308 MATAIJCRL - MATAIJCRL = "aijcrl" - A matrix type to be used for sparse matrices. 4309 4310 This matrix type is identical to `MATSEQAIJCRL` when constructed with a single process communicator, 4311 and `MATMPIAIJCRL` otherwise. As a result, for single process communicators, 4312 `MatSeqAIJSetPreallocation()` is supported, and similarly `MatMPIAIJSetPreallocation()` is supported 4313 for communicators controlling multiple processes. It is recommended that you call both of 4314 the above preallocation routines for simplicity. 4315 4316 Options Database Keys: 4317 . -mat_type aijcrl - sets the matrix type to "aijcrl" during a call to `MatSetFromOptions()` 4318 4319 Level: beginner 4320 4321 .seealso: `MatCreateMPIAIJCRL`, `MATSEQAIJCRL`, `MATMPIAIJCRL`, `MATSEQAIJCRL`, `MATMPIAIJCRL` 4322 M*/ 4323 4324 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCRL(Mat, MatType, MatReuse, Mat *); 4325 #if defined(PETSC_HAVE_ELEMENTAL) 4326 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_Elemental(Mat, MatType, MatReuse, Mat *); 4327 #endif 4328 #if defined(PETSC_HAVE_SCALAPACK) 4329 PETSC_INTERN PetscErrorCode MatConvert_AIJ_ScaLAPACK(Mat, MatType, MatReuse, Mat *); 4330 #endif 4331 #if defined(PETSC_HAVE_HYPRE) 4332 PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A, MatType, MatReuse, Mat *); 4333 #endif 4334 4335 PETSC_EXTERN PetscErrorCode MatConvert_SeqAIJ_SeqSELL(Mat, MatType, MatReuse, Mat *); 4336 PETSC_INTERN PetscErrorCode MatConvert_XAIJ_IS(Mat, MatType, MatReuse, Mat *); 4337 PETSC_INTERN PetscErrorCode MatProductSetFromOptions_IS_XAIJ(Mat); 4338 4339 /*@C 4340 MatSeqAIJGetArray - gives read/write access to the array where the data for a `MATSEQAIJ` matrix is stored 4341 4342 Not Collective 4343 4344 Input Parameter: 4345 . mat - a `MATSEQAIJ` matrix 4346 4347 Output Parameter: 4348 . array - pointer to the data 4349 4350 Level: intermediate 4351 4352 .seealso: `MatSeqAIJRestoreArray()`, `MatSeqAIJGetArrayF90()` 4353 @*/ 4354 PetscErrorCode MatSeqAIJGetArray(Mat A, PetscScalar **array) 4355 { 4356 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)A->data; 4357 4358 PetscFunctionBegin; 4359 if (aij->ops->getarray) { 4360 PetscCall((*aij->ops->getarray)(A, array)); 4361 } else { 4362 *array = aij->a; 4363 } 4364 PetscFunctionReturn(0); 4365 } 4366 4367 /*@C 4368 MatSeqAIJRestoreArray - returns access to the array where the data for a `MATSEQAIJ` matrix is stored obtained by `MatSeqAIJGetArray()` 4369 4370 Not Collective 4371 4372 Input Parameters: 4373 + mat - a `MATSEQAIJ` matrix 4374 - array - pointer to the data 4375 4376 Level: intermediate 4377 4378 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArrayF90()` 4379 @*/ 4380 PetscErrorCode MatSeqAIJRestoreArray(Mat A, PetscScalar **array) 4381 { 4382 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)A->data; 4383 4384 PetscFunctionBegin; 4385 if (aij->ops->restorearray) { 4386 PetscCall((*aij->ops->restorearray)(A, array)); 4387 } else { 4388 *array = NULL; 4389 } 4390 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 4391 PetscCall(PetscObjectStateIncrease((PetscObject)A)); 4392 PetscFunctionReturn(0); 4393 } 4394 4395 /*@C 4396 MatSeqAIJGetArrayRead - gives read-only access to the array where the data for a `MATSEQAIJ` matrix is stored 4397 4398 Not Collective 4399 4400 Input Parameter: 4401 . mat - a `MATSEQAIJ` matrix 4402 4403 Output Parameter: 4404 . array - pointer to the data 4405 4406 Level: intermediate 4407 4408 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArrayRead()` 4409 @*/ 4410 PetscErrorCode MatSeqAIJGetArrayRead(Mat A, const PetscScalar **array) 4411 { 4412 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)A->data; 4413 4414 PetscFunctionBegin; 4415 if (aij->ops->getarrayread) { 4416 PetscCall((*aij->ops->getarrayread)(A, array)); 4417 } else { 4418 *array = aij->a; 4419 } 4420 PetscFunctionReturn(0); 4421 } 4422 4423 /*@C 4424 MatSeqAIJRestoreArrayRead - restore the read-only access array obtained from `MatSeqAIJGetArrayRead()` 4425 4426 Not Collective 4427 4428 Input Parameter: 4429 . mat - a `MATSEQAIJ` matrix 4430 4431 Output Parameter: 4432 . array - pointer to the data 4433 4434 Level: intermediate 4435 4436 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJGetArrayRead()` 4437 @*/ 4438 PetscErrorCode MatSeqAIJRestoreArrayRead(Mat A, const PetscScalar **array) 4439 { 4440 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)A->data; 4441 4442 PetscFunctionBegin; 4443 if (aij->ops->restorearrayread) { 4444 PetscCall((*aij->ops->restorearrayread)(A, array)); 4445 } else { 4446 *array = NULL; 4447 } 4448 PetscFunctionReturn(0); 4449 } 4450 4451 /*@C 4452 MatSeqAIJGetArrayWrite - gives write-only access to the array where the data for a `MATSEQAIJ` matrix is stored 4453 4454 Not Collective 4455 4456 Input Parameter: 4457 . mat - a `MATSEQAIJ` matrix 4458 4459 Output Parameter: 4460 . array - pointer to the data 4461 4462 Level: intermediate 4463 4464 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArrayRead()` 4465 @*/ 4466 PetscErrorCode MatSeqAIJGetArrayWrite(Mat A, PetscScalar **array) 4467 { 4468 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)A->data; 4469 4470 PetscFunctionBegin; 4471 if (aij->ops->getarraywrite) { 4472 PetscCall((*aij->ops->getarraywrite)(A, array)); 4473 } else { 4474 *array = aij->a; 4475 } 4476 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 4477 PetscCall(PetscObjectStateIncrease((PetscObject)A)); 4478 PetscFunctionReturn(0); 4479 } 4480 4481 /*@C 4482 MatSeqAIJRestoreArrayWrite - restore the read-only access array obtained from MatSeqAIJGetArrayRead 4483 4484 Not Collective 4485 4486 Input Parameter: 4487 . mat - a MATSEQAIJ matrix 4488 4489 Output Parameter: 4490 . array - pointer to the data 4491 4492 Level: intermediate 4493 4494 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJGetArrayRead()` 4495 @*/ 4496 PetscErrorCode MatSeqAIJRestoreArrayWrite(Mat A, PetscScalar **array) 4497 { 4498 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)A->data; 4499 4500 PetscFunctionBegin; 4501 if (aij->ops->restorearraywrite) { 4502 PetscCall((*aij->ops->restorearraywrite)(A, array)); 4503 } else { 4504 *array = NULL; 4505 } 4506 PetscFunctionReturn(0); 4507 } 4508 4509 /*@C 4510 MatSeqAIJGetCSRAndMemType - Get the CSR arrays and the memory type of the `MATSEQAIJ` matrix 4511 4512 Not Collective 4513 4514 Input Parameter: 4515 . mat - a matrix of type `MATSEQAIJ` or its subclasses 4516 4517 Output Parameters: 4518 + i - row map array of the matrix 4519 . j - column index array of the matrix 4520 . a - data array of the matrix 4521 - memtype - memory type of the arrays 4522 4523 Notes: 4524 Any of the output parameters can be NULL, in which case the corresponding value is not returned. 4525 If mat is a device matrix, the arrays are on the device. Otherwise, they are on the host. 4526 4527 One can call this routine on a preallocated but not assembled matrix to just get the memory of the CSR underneath the matrix. 4528 If the matrix is assembled, the data array 'a' is guaranteed to have the latest values of the matrix. 4529 4530 Level: Developer 4531 4532 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJGetArrayRead()` 4533 @*/ 4534 PetscErrorCode MatSeqAIJGetCSRAndMemType(Mat mat, const PetscInt **i, const PetscInt **j, PetscScalar **a, PetscMemType *mtype) 4535 { 4536 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 4537 4538 PetscFunctionBegin; 4539 PetscCheck(mat->preallocated, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "matrix is not preallocated"); 4540 if (aij->ops->getcsrandmemtype) { 4541 PetscCall((*aij->ops->getcsrandmemtype)(mat, i, j, a, mtype)); 4542 } else { 4543 if (i) *i = aij->i; 4544 if (j) *j = aij->j; 4545 if (a) *a = aij->a; 4546 if (mtype) *mtype = PETSC_MEMTYPE_HOST; 4547 } 4548 PetscFunctionReturn(0); 4549 } 4550 4551 /*@C 4552 MatSeqAIJGetMaxRowNonzeros - returns the maximum number of nonzeros in any row 4553 4554 Not Collective 4555 4556 Input Parameter: 4557 . mat - a `MATSEQAIJ` matrix 4558 4559 Output Parameter: 4560 . nz - the maximum number of nonzeros in any row 4561 4562 Level: intermediate 4563 4564 .seealso: `MatSeqAIJRestoreArray()`, `MatSeqAIJGetArrayF90()` 4565 @*/ 4566 PetscErrorCode MatSeqAIJGetMaxRowNonzeros(Mat A, PetscInt *nz) 4567 { 4568 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)A->data; 4569 4570 PetscFunctionBegin; 4571 *nz = aij->rmax; 4572 PetscFunctionReturn(0); 4573 } 4574 4575 PetscErrorCode MatSetPreallocationCOO_SeqAIJ(Mat mat, PetscCount coo_n, PetscInt coo_i[], PetscInt coo_j[]) 4576 { 4577 MPI_Comm comm; 4578 PetscInt *i, *j; 4579 PetscInt M, N, row; 4580 PetscCount k, p, q, nneg, nnz, start, end; /* Index the coo array, so use PetscCount as their type */ 4581 PetscInt *Ai; /* Change to PetscCount once we use it for row pointers */ 4582 PetscInt *Aj; 4583 PetscScalar *Aa; 4584 Mat_SeqAIJ *seqaij = (Mat_SeqAIJ *)(mat->data); 4585 MatType rtype; 4586 PetscCount *perm, *jmap; 4587 4588 PetscFunctionBegin; 4589 PetscCall(MatResetPreallocationCOO_SeqAIJ(mat)); 4590 PetscCall(PetscObjectGetComm((PetscObject)mat, &comm)); 4591 PetscCall(MatGetSize(mat, &M, &N)); 4592 i = coo_i; 4593 j = coo_j; 4594 PetscCall(PetscMalloc1(coo_n, &perm)); 4595 for (k = 0; k < coo_n; k++) { /* Ignore entries with negative row or col indices */ 4596 if (j[k] < 0) i[k] = -1; 4597 perm[k] = k; 4598 } 4599 4600 /* Sort by row */ 4601 PetscCall(PetscSortIntWithIntCountArrayPair(coo_n, i, j, perm)); 4602 for (k = 0; k < coo_n; k++) { 4603 if (i[k] >= 0) break; 4604 } /* Advance k to the first row with a non-negative index */ 4605 nneg = k; 4606 PetscCall(PetscMalloc1(coo_n - nneg + 1, &jmap)); /* +1 to make a CSR-like data structure. jmap[i] originally is the number of repeats for i-th nonzero */ 4607 nnz = 0; /* Total number of unique nonzeros to be counted */ 4608 jmap++; /* Inc jmap by 1 for convinience */ 4609 4610 PetscCall(PetscCalloc1(M + 1, &Ai)); /* CSR of A */ 4611 PetscCall(PetscMalloc1(coo_n - nneg, &Aj)); /* We have at most coo_n-nneg unique nonzeros */ 4612 4613 /* In each row, sort by column, then unique column indices to get row length */ 4614 Ai++; /* Inc by 1 for convinience */ 4615 q = 0; /* q-th unique nonzero, with q starting from 0 */ 4616 while (k < coo_n) { 4617 row = i[k]; 4618 start = k; /* [start,end) indices for this row */ 4619 while (k < coo_n && i[k] == row) k++; 4620 end = k; 4621 PetscCall(PetscSortIntWithCountArray(end - start, j + start, perm + start)); 4622 /* Find number of unique col entries in this row */ 4623 Aj[q] = j[start]; /* Log the first nonzero in this row */ 4624 jmap[q] = 1; /* Number of repeats of this nozero entry */ 4625 Ai[row] = 1; 4626 nnz++; 4627 4628 for (p = start + 1; p < end; p++) { /* Scan remaining nonzero in this row */ 4629 if (j[p] != j[p - 1]) { /* Meet a new nonzero */ 4630 q++; 4631 jmap[q] = 1; 4632 Aj[q] = j[p]; 4633 Ai[row]++; 4634 nnz++; 4635 } else { 4636 jmap[q]++; 4637 } 4638 } 4639 q++; /* Move to next row and thus next unique nonzero */ 4640 } 4641 4642 Ai--; /* Back to the beginning of Ai[] */ 4643 for (k = 0; k < M; k++) Ai[k + 1] += Ai[k]; 4644 jmap--; /* Back to the beginning of jmap[] */ 4645 jmap[0] = 0; 4646 for (k = 0; k < nnz; k++) jmap[k + 1] += jmap[k]; 4647 if (nnz < coo_n - nneg) { /* Realloc with actual number of unique nonzeros */ 4648 PetscCount *jmap_new; 4649 PetscInt *Aj_new; 4650 4651 PetscCall(PetscMalloc1(nnz + 1, &jmap_new)); 4652 PetscCall(PetscArraycpy(jmap_new, jmap, nnz + 1)); 4653 PetscCall(PetscFree(jmap)); 4654 jmap = jmap_new; 4655 4656 PetscCall(PetscMalloc1(nnz, &Aj_new)); 4657 PetscCall(PetscArraycpy(Aj_new, Aj, nnz)); 4658 PetscCall(PetscFree(Aj)); 4659 Aj = Aj_new; 4660 } 4661 4662 if (nneg) { /* Discard heading entries with negative indices in perm[], as we'll access it from index 0 in MatSetValuesCOO */ 4663 PetscCount *perm_new; 4664 4665 PetscCall(PetscMalloc1(coo_n - nneg, &perm_new)); 4666 PetscCall(PetscArraycpy(perm_new, perm + nneg, coo_n - nneg)); 4667 PetscCall(PetscFree(perm)); 4668 perm = perm_new; 4669 } 4670 4671 PetscCall(MatGetRootType_Private(mat, &rtype)); 4672 PetscCall(PetscCalloc1(nnz, &Aa)); /* Zero the matrix */ 4673 PetscCall(MatSetSeqAIJWithArrays_private(PETSC_COMM_SELF, M, N, Ai, Aj, Aa, rtype, mat)); 4674 4675 seqaij->singlemalloc = PETSC_FALSE; /* Ai, Aj and Aa are not allocated in one big malloc */ 4676 seqaij->free_a = seqaij->free_ij = PETSC_TRUE; /* Let newmat own Ai, Aj and Aa */ 4677 /* Record COO fields */ 4678 seqaij->coo_n = coo_n; 4679 seqaij->Atot = coo_n - nneg; /* Annz is seqaij->nz, so no need to record that again */ 4680 seqaij->jmap = jmap; /* of length nnz+1 */ 4681 seqaij->perm = perm; 4682 PetscFunctionReturn(0); 4683 } 4684 4685 static PetscErrorCode MatSetValuesCOO_SeqAIJ(Mat A, const PetscScalar v[], InsertMode imode) 4686 { 4687 Mat_SeqAIJ *aseq = (Mat_SeqAIJ *)A->data; 4688 PetscCount i, j, Annz = aseq->nz; 4689 PetscCount *perm = aseq->perm, *jmap = aseq->jmap; 4690 PetscScalar *Aa; 4691 4692 PetscFunctionBegin; 4693 PetscCall(MatSeqAIJGetArray(A, &Aa)); 4694 for (i = 0; i < Annz; i++) { 4695 PetscScalar sum = 0.0; 4696 for (j = jmap[i]; j < jmap[i + 1]; j++) sum += v[perm[j]]; 4697 Aa[i] = (imode == INSERT_VALUES ? 0.0 : Aa[i]) + sum; 4698 } 4699 PetscCall(MatSeqAIJRestoreArray(A, &Aa)); 4700 PetscFunctionReturn(0); 4701 } 4702 4703 #if defined(PETSC_HAVE_CUDA) 4704 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCUSPARSE(Mat, MatType, MatReuse, Mat *); 4705 #endif 4706 #if defined(PETSC_HAVE_KOKKOS_KERNELS) 4707 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJKokkos(Mat, MatType, MatReuse, Mat *); 4708 #endif 4709 4710 PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJ(Mat B) 4711 { 4712 Mat_SeqAIJ *b; 4713 PetscMPIInt size; 4714 4715 PetscFunctionBegin; 4716 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B), &size)); 4717 PetscCheck(size <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Comm must be of size 1"); 4718 4719 PetscCall(PetscNew(&b)); 4720 4721 B->data = (void *)b; 4722 4723 PetscCall(PetscMemcpy(B->ops, &MatOps_Values, sizeof(struct _MatOps))); 4724 if (B->sortedfull) B->ops->setvalues = MatSetValues_SeqAIJ_SortedFull; 4725 4726 b->row = NULL; 4727 b->col = NULL; 4728 b->icol = NULL; 4729 b->reallocs = 0; 4730 b->ignorezeroentries = PETSC_FALSE; 4731 b->roworiented = PETSC_TRUE; 4732 b->nonew = 0; 4733 b->diag = NULL; 4734 b->solve_work = NULL; 4735 B->spptr = NULL; 4736 b->saved_values = NULL; 4737 b->idiag = NULL; 4738 b->mdiag = NULL; 4739 b->ssor_work = NULL; 4740 b->omega = 1.0; 4741 b->fshift = 0.0; 4742 b->idiagvalid = PETSC_FALSE; 4743 b->ibdiagvalid = PETSC_FALSE; 4744 b->keepnonzeropattern = PETSC_FALSE; 4745 4746 PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATSEQAIJ)); 4747 #if defined(PETSC_HAVE_MATLAB) 4748 PetscCall(PetscObjectComposeFunction((PetscObject)B, "PetscMatlabEnginePut_C", MatlabEnginePut_SeqAIJ)); 4749 PetscCall(PetscObjectComposeFunction((PetscObject)B, "PetscMatlabEngineGet_C", MatlabEngineGet_SeqAIJ)); 4750 #endif 4751 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqAIJSetColumnIndices_C", MatSeqAIJSetColumnIndices_SeqAIJ)); 4752 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatStoreValues_C", MatStoreValues_SeqAIJ)); 4753 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatRetrieveValues_C", MatRetrieveValues_SeqAIJ)); 4754 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqsbaij_C", MatConvert_SeqAIJ_SeqSBAIJ)); 4755 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqbaij_C", MatConvert_SeqAIJ_SeqBAIJ)); 4756 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqaijperm_C", MatConvert_SeqAIJ_SeqAIJPERM)); 4757 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqaijsell_C", MatConvert_SeqAIJ_SeqAIJSELL)); 4758 #if defined(PETSC_HAVE_MKL_SPARSE) 4759 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqaijmkl_C", MatConvert_SeqAIJ_SeqAIJMKL)); 4760 #endif 4761 #if defined(PETSC_HAVE_CUDA) 4762 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqaijcusparse_C", MatConvert_SeqAIJ_SeqAIJCUSPARSE)); 4763 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaijcusparse_seqaij_C", MatProductSetFromOptions_SeqAIJ)); 4764 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaij_seqaijcusparse_C", MatProductSetFromOptions_SeqAIJ)); 4765 #endif 4766 #if defined(PETSC_HAVE_KOKKOS_KERNELS) 4767 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqaijkokkos_C", MatConvert_SeqAIJ_SeqAIJKokkos)); 4768 #endif 4769 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqaijcrl_C", MatConvert_SeqAIJ_SeqAIJCRL)); 4770 #if defined(PETSC_HAVE_ELEMENTAL) 4771 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_elemental_C", MatConvert_SeqAIJ_Elemental)); 4772 #endif 4773 #if defined(PETSC_HAVE_SCALAPACK) 4774 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_scalapack_C", MatConvert_AIJ_ScaLAPACK)); 4775 #endif 4776 #if defined(PETSC_HAVE_HYPRE) 4777 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_hypre_C", MatConvert_AIJ_HYPRE)); 4778 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_transpose_seqaij_seqaij_C", MatProductSetFromOptions_Transpose_AIJ_AIJ)); 4779 #endif 4780 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqdense_C", MatConvert_SeqAIJ_SeqDense)); 4781 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_seqsell_C", MatConvert_SeqAIJ_SeqSELL)); 4782 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqaij_is_C", MatConvert_XAIJ_IS)); 4783 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatIsTranspose_C", MatIsTranspose_SeqAIJ)); 4784 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatIsHermitianTranspose_C", MatIsTranspose_SeqAIJ)); 4785 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqAIJSetPreallocation_C", MatSeqAIJSetPreallocation_SeqAIJ)); 4786 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatResetPreallocation_C", MatResetPreallocation_SeqAIJ)); 4787 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqAIJSetPreallocationCSR_C", MatSeqAIJSetPreallocationCSR_SeqAIJ)); 4788 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatReorderForNonzeroDiagonal_C", MatReorderForNonzeroDiagonal_SeqAIJ)); 4789 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_is_seqaij_C", MatProductSetFromOptions_IS_XAIJ)); 4790 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdense_seqaij_C", MatProductSetFromOptions_SeqDense_SeqAIJ)); 4791 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaij_seqaij_C", MatProductSetFromOptions_SeqAIJ)); 4792 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqAIJKron_C", MatSeqAIJKron_SeqAIJ)); 4793 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSetPreallocationCOO_C", MatSetPreallocationCOO_SeqAIJ)); 4794 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSetValuesCOO_C", MatSetValuesCOO_SeqAIJ)); 4795 PetscCall(MatCreate_SeqAIJ_Inode(B)); 4796 PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATSEQAIJ)); 4797 PetscCall(MatSeqAIJSetTypeFromOptions(B)); /* this allows changing the matrix subtype to say MATSEQAIJPERM */ 4798 PetscFunctionReturn(0); 4799 } 4800 4801 /* 4802 Given a matrix generated with MatGetFactor() duplicates all the information in A into C 4803 */ 4804 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C, Mat A, MatDuplicateOption cpvalues, PetscBool mallocmatspace) 4805 { 4806 Mat_SeqAIJ *c = (Mat_SeqAIJ *)C->data, *a = (Mat_SeqAIJ *)A->data; 4807 PetscInt m = A->rmap->n, i; 4808 4809 PetscFunctionBegin; 4810 PetscCheck(A->assembled || cpvalues == MAT_DO_NOT_COPY_VALUES, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cannot duplicate unassembled matrix"); 4811 4812 C->factortype = A->factortype; 4813 c->row = NULL; 4814 c->col = NULL; 4815 c->icol = NULL; 4816 c->reallocs = 0; 4817 4818 C->assembled = A->assembled; 4819 C->preallocated = A->preallocated; 4820 4821 if (A->preallocated) { 4822 PetscCall(PetscLayoutReference(A->rmap, &C->rmap)); 4823 PetscCall(PetscLayoutReference(A->cmap, &C->cmap)); 4824 4825 PetscCall(PetscMalloc1(m, &c->imax)); 4826 PetscCall(PetscMemcpy(c->imax, a->imax, m * sizeof(PetscInt))); 4827 PetscCall(PetscMalloc1(m, &c->ilen)); 4828 PetscCall(PetscMemcpy(c->ilen, a->ilen, m * sizeof(PetscInt))); 4829 4830 /* allocate the matrix space */ 4831 if (mallocmatspace) { 4832 PetscCall(PetscMalloc3(a->i[m], &c->a, a->i[m], &c->j, m + 1, &c->i)); 4833 4834 c->singlemalloc = PETSC_TRUE; 4835 4836 PetscCall(PetscArraycpy(c->i, a->i, m + 1)); 4837 if (m > 0) { 4838 PetscCall(PetscArraycpy(c->j, a->j, a->i[m])); 4839 if (cpvalues == MAT_COPY_VALUES) { 4840 const PetscScalar *aa; 4841 4842 PetscCall(MatSeqAIJGetArrayRead(A, &aa)); 4843 PetscCall(PetscArraycpy(c->a, aa, a->i[m])); 4844 PetscCall(MatSeqAIJGetArrayRead(A, &aa)); 4845 } else { 4846 PetscCall(PetscArrayzero(c->a, a->i[m])); 4847 } 4848 } 4849 } 4850 4851 c->ignorezeroentries = a->ignorezeroentries; 4852 c->roworiented = a->roworiented; 4853 c->nonew = a->nonew; 4854 if (a->diag) { 4855 PetscCall(PetscMalloc1(m + 1, &c->diag)); 4856 PetscCall(PetscMemcpy(c->diag, a->diag, m * sizeof(PetscInt))); 4857 } else c->diag = NULL; 4858 4859 c->solve_work = NULL; 4860 c->saved_values = NULL; 4861 c->idiag = NULL; 4862 c->ssor_work = NULL; 4863 c->keepnonzeropattern = a->keepnonzeropattern; 4864 c->free_a = PETSC_TRUE; 4865 c->free_ij = PETSC_TRUE; 4866 4867 c->rmax = a->rmax; 4868 c->nz = a->nz; 4869 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 4870 4871 c->compressedrow.use = a->compressedrow.use; 4872 c->compressedrow.nrows = a->compressedrow.nrows; 4873 if (a->compressedrow.use) { 4874 i = a->compressedrow.nrows; 4875 PetscCall(PetscMalloc2(i + 1, &c->compressedrow.i, i, &c->compressedrow.rindex)); 4876 PetscCall(PetscArraycpy(c->compressedrow.i, a->compressedrow.i, i + 1)); 4877 PetscCall(PetscArraycpy(c->compressedrow.rindex, a->compressedrow.rindex, i)); 4878 } else { 4879 c->compressedrow.use = PETSC_FALSE; 4880 c->compressedrow.i = NULL; 4881 c->compressedrow.rindex = NULL; 4882 } 4883 c->nonzerorowcnt = a->nonzerorowcnt; 4884 C->nonzerostate = A->nonzerostate; 4885 4886 PetscCall(MatDuplicate_SeqAIJ_Inode(A, cpvalues, &C)); 4887 } 4888 PetscCall(PetscFunctionListDuplicate(((PetscObject)A)->qlist, &((PetscObject)C)->qlist)); 4889 PetscFunctionReturn(0); 4890 } 4891 4892 PetscErrorCode MatDuplicate_SeqAIJ(Mat A, MatDuplicateOption cpvalues, Mat *B) 4893 { 4894 PetscFunctionBegin; 4895 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), B)); 4896 PetscCall(MatSetSizes(*B, A->rmap->n, A->cmap->n, A->rmap->n, A->cmap->n)); 4897 if (!(A->rmap->n % A->rmap->bs) && !(A->cmap->n % A->cmap->bs)) PetscCall(MatSetBlockSizesFromMats(*B, A, A)); 4898 PetscCall(MatSetType(*B, ((PetscObject)A)->type_name)); 4899 PetscCall(MatDuplicateNoCreate_SeqAIJ(*B, A, cpvalues, PETSC_TRUE)); 4900 PetscFunctionReturn(0); 4901 } 4902 4903 PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer) 4904 { 4905 PetscBool isbinary, ishdf5; 4906 4907 PetscFunctionBegin; 4908 PetscValidHeaderSpecific(newMat, MAT_CLASSID, 1); 4909 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 4910 /* force binary viewer to load .info file if it has not yet done so */ 4911 PetscCall(PetscViewerSetUp(viewer)); 4912 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary)); 4913 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5)); 4914 if (isbinary) { 4915 PetscCall(MatLoad_SeqAIJ_Binary(newMat, viewer)); 4916 } else if (ishdf5) { 4917 #if defined(PETSC_HAVE_HDF5) 4918 PetscCall(MatLoad_AIJ_HDF5(newMat, viewer)); 4919 #else 4920 SETERRQ(PetscObjectComm((PetscObject)newMat), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5"); 4921 #endif 4922 } else { 4923 SETERRQ(PetscObjectComm((PetscObject)newMat), PETSC_ERR_SUP, "Viewer type %s not yet supported for reading %s matrices", ((PetscObject)viewer)->type_name, ((PetscObject)newMat)->type_name); 4924 } 4925 PetscFunctionReturn(0); 4926 } 4927 4928 PetscErrorCode MatLoad_SeqAIJ_Binary(Mat mat, PetscViewer viewer) 4929 { 4930 Mat_SeqAIJ *a = (Mat_SeqAIJ *)mat->data; 4931 PetscInt header[4], *rowlens, M, N, nz, sum, rows, cols, i; 4932 4933 PetscFunctionBegin; 4934 PetscCall(PetscViewerSetUp(viewer)); 4935 4936 /* read in matrix header */ 4937 PetscCall(PetscViewerBinaryRead(viewer, header, 4, NULL, PETSC_INT)); 4938 PetscCheck(header[0] == MAT_FILE_CLASSID, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Not a matrix object in file"); 4939 M = header[1]; 4940 N = header[2]; 4941 nz = header[3]; 4942 PetscCheck(M >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Matrix row size (%" PetscInt_FMT ") in file is negative", M); 4943 PetscCheck(N >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Matrix column size (%" PetscInt_FMT ") in file is negative", N); 4944 PetscCheck(nz >= 0, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Matrix stored in special format on disk, cannot load as SeqAIJ"); 4945 4946 /* set block sizes from the viewer's .info file */ 4947 PetscCall(MatLoad_Binary_BlockSizes(mat, viewer)); 4948 /* set local and global sizes if not set already */ 4949 if (mat->rmap->n < 0) mat->rmap->n = M; 4950 if (mat->cmap->n < 0) mat->cmap->n = N; 4951 if (mat->rmap->N < 0) mat->rmap->N = M; 4952 if (mat->cmap->N < 0) mat->cmap->N = N; 4953 PetscCall(PetscLayoutSetUp(mat->rmap)); 4954 PetscCall(PetscLayoutSetUp(mat->cmap)); 4955 4956 /* check if the matrix sizes are correct */ 4957 PetscCall(MatGetSize(mat, &rows, &cols)); 4958 PetscCheck(M == rows && N == cols, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Matrix in file of different sizes (%" PetscInt_FMT ", %" PetscInt_FMT ") than the input matrix (%" PetscInt_FMT ", %" PetscInt_FMT ")", M, N, rows, cols); 4959 4960 /* read in row lengths */ 4961 PetscCall(PetscMalloc1(M, &rowlens)); 4962 PetscCall(PetscViewerBinaryRead(viewer, rowlens, M, NULL, PETSC_INT)); 4963 /* check if sum(rowlens) is same as nz */ 4964 sum = 0; 4965 for (i = 0; i < M; i++) sum += rowlens[i]; 4966 PetscCheck(sum == nz, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Inconsistent matrix data in file: nonzeros = %" PetscInt_FMT ", sum-row-lengths = %" PetscInt_FMT, nz, sum); 4967 /* preallocate and check sizes */ 4968 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(mat, 0, rowlens)); 4969 PetscCall(MatGetSize(mat, &rows, &cols)); 4970 PetscCheck(M == rows && N == cols, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Matrix in file of different length (%" PetscInt_FMT ", %" PetscInt_FMT ") than the input matrix (%" PetscInt_FMT ", %" PetscInt_FMT ")", M, N, rows, cols); 4971 /* store row lengths */ 4972 PetscCall(PetscArraycpy(a->ilen, rowlens, M)); 4973 PetscCall(PetscFree(rowlens)); 4974 4975 /* fill in "i" row pointers */ 4976 a->i[0] = 0; 4977 for (i = 0; i < M; i++) a->i[i + 1] = a->i[i] + a->ilen[i]; 4978 /* read in "j" column indices */ 4979 PetscCall(PetscViewerBinaryRead(viewer, a->j, nz, NULL, PETSC_INT)); 4980 /* read in "a" nonzero values */ 4981 PetscCall(PetscViewerBinaryRead(viewer, a->a, nz, NULL, PETSC_SCALAR)); 4982 4983 PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY)); 4984 PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY)); 4985 PetscFunctionReturn(0); 4986 } 4987 4988 PetscErrorCode MatEqual_SeqAIJ(Mat A, Mat B, PetscBool *flg) 4989 { 4990 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data; 4991 const PetscScalar *aa, *ba; 4992 #if defined(PETSC_USE_COMPLEX) 4993 PetscInt k; 4994 #endif 4995 4996 PetscFunctionBegin; 4997 /* If the matrix dimensions are not equal,or no of nonzeros */ 4998 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) || (a->nz != b->nz)) { 4999 *flg = PETSC_FALSE; 5000 PetscFunctionReturn(0); 5001 } 5002 5003 /* if the a->i are the same */ 5004 PetscCall(PetscArraycmp(a->i, b->i, A->rmap->n + 1, flg)); 5005 if (!*flg) PetscFunctionReturn(0); 5006 5007 /* if a->j are the same */ 5008 PetscCall(PetscArraycmp(a->j, b->j, a->nz, flg)); 5009 if (!*flg) PetscFunctionReturn(0); 5010 5011 PetscCall(MatSeqAIJGetArrayRead(A, &aa)); 5012 PetscCall(MatSeqAIJGetArrayRead(B, &ba)); 5013 /* if a->a are the same */ 5014 #if defined(PETSC_USE_COMPLEX) 5015 for (k = 0; k < a->nz; k++) { 5016 if (PetscRealPart(aa[k]) != PetscRealPart(ba[k]) || PetscImaginaryPart(aa[k]) != PetscImaginaryPart(ba[k])) { 5017 *flg = PETSC_FALSE; 5018 PetscFunctionReturn(0); 5019 } 5020 } 5021 #else 5022 PetscCall(PetscArraycmp(aa, ba, a->nz, flg)); 5023 #endif 5024 PetscCall(MatSeqAIJRestoreArrayRead(A, &aa)); 5025 PetscCall(MatSeqAIJRestoreArrayRead(B, &ba)); 5026 PetscFunctionReturn(0); 5027 } 5028 5029 /*@ 5030 MatCreateSeqAIJWithArrays - Creates an sequential `MATSEQAIJ` matrix using matrix elements (in CSR format) 5031 provided by the user. 5032 5033 Collective 5034 5035 Input Parameters: 5036 + comm - must be an MPI communicator of size 1 5037 . m - number of rows 5038 . n - number of columns 5039 . i - row indices; that is i[0] = 0, i[row] = i[row-1] + number of elements in that row of the matrix 5040 . j - column indices 5041 - a - matrix values 5042 5043 Output Parameter: 5044 . mat - the matrix 5045 5046 Level: intermediate 5047 5048 Notes: 5049 The i, j, and a arrays are not copied by this routine, the user must free these arrays 5050 once the matrix is destroyed and not before 5051 5052 You cannot set new nonzero locations into this matrix, that will generate an error. 5053 5054 The i and j indices are 0 based 5055 5056 The format which is used for the sparse matrix input, is equivalent to a 5057 row-major ordering.. i.e for the following matrix, the input data expected is 5058 as shown 5059 5060 $ 1 0 0 5061 $ 2 0 3 5062 $ 4 5 6 5063 $ 5064 $ i = {0,1,3,6} [size = nrow+1 = 3+1] 5065 $ j = {0,0,2,0,1,2} [size = 6]; values must be sorted for each row 5066 $ v = {1,2,3,4,5,6} [size = 6] 5067 5068 .seealso: `MatCreate()`, `MatCreateAIJ()`, `MatCreateSeqAIJ()`, `MatCreateMPIAIJWithArrays()`, `MatMPIAIJSetPreallocationCSR()` 5069 @*/ 5070 PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt i[], PetscInt j[], PetscScalar a[], Mat *mat) 5071 { 5072 PetscInt ii; 5073 Mat_SeqAIJ *aij; 5074 PetscInt jj; 5075 5076 PetscFunctionBegin; 5077 PetscCheck(m <= 0 || i[0] == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "i (row indices) must start with 0"); 5078 PetscCall(MatCreate(comm, mat)); 5079 PetscCall(MatSetSizes(*mat, m, n, m, n)); 5080 /* PetscCall(MatSetBlockSizes(*mat,,)); */ 5081 PetscCall(MatSetType(*mat, MATSEQAIJ)); 5082 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(*mat, MAT_SKIP_ALLOCATION, NULL)); 5083 aij = (Mat_SeqAIJ *)(*mat)->data; 5084 PetscCall(PetscMalloc1(m, &aij->imax)); 5085 PetscCall(PetscMalloc1(m, &aij->ilen)); 5086 5087 aij->i = i; 5088 aij->j = j; 5089 aij->a = a; 5090 aij->singlemalloc = PETSC_FALSE; 5091 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 5092 aij->free_a = PETSC_FALSE; 5093 aij->free_ij = PETSC_FALSE; 5094 5095 for (ii = 0, aij->nonzerorowcnt = 0, aij->rmax = 0; ii < m; ii++) { 5096 aij->ilen[ii] = aij->imax[ii] = i[ii + 1] - i[ii]; 5097 if (PetscDefined(USE_DEBUG)) { 5098 PetscCheck(i[ii + 1] - i[ii] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative row length in i (row indices) row = %" PetscInt_FMT " length = %" PetscInt_FMT, ii, i[ii + 1] - i[ii]); 5099 for (jj = i[ii] + 1; jj < i[ii + 1]; jj++) { 5100 PetscCheck(j[jj] >= j[jj - 1], PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column entry number %" PetscInt_FMT " (actual column %" PetscInt_FMT ") in row %" PetscInt_FMT " is not sorted", jj - i[ii], j[jj], ii); 5101 PetscCheck(j[jj] != j[jj - 1], PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column entry number %" PetscInt_FMT " (actual column %" PetscInt_FMT ") in row %" PetscInt_FMT " is identical to previous entry", jj - i[ii], j[jj], ii); 5102 } 5103 } 5104 } 5105 if (PetscDefined(USE_DEBUG)) { 5106 for (ii = 0; ii < aij->i[m]; ii++) { 5107 PetscCheck(j[ii] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative column index at location = %" PetscInt_FMT " index = %" PetscInt_FMT, ii, j[ii]); 5108 PetscCheck(j[ii] <= n - 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column index to large at location = %" PetscInt_FMT " index = %" PetscInt_FMT, ii, j[ii]); 5109 } 5110 } 5111 5112 PetscCall(MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY)); 5113 PetscCall(MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY)); 5114 PetscFunctionReturn(0); 5115 } 5116 5117 /*@ 5118 MatCreateSeqAIJFromTriple - Creates an sequential `MATSEQAIJ` matrix using matrix elements (in COO format) 5119 provided by the user. 5120 5121 Collective 5122 5123 Input Parameters: 5124 + comm - must be an MPI communicator of size 1 5125 . m - number of rows 5126 . n - number of columns 5127 . i - row indices 5128 . j - column indices 5129 . a - matrix values 5130 . nz - number of nonzeros 5131 - idx - if the i and j indices start with 1 use `PETSC_TRUE` otherwise use `PETSC_FALSE` 5132 5133 Output Parameter: 5134 . mat - the matrix 5135 5136 Level: intermediate 5137 5138 Example: 5139 For the following matrix, the input data expected is as shown (using 0 based indexing) 5140 .vb 5141 1 0 0 5142 2 0 3 5143 4 5 6 5144 5145 i = {0,1,1,2,2,2} 5146 j = {0,0,2,0,1,2} 5147 v = {1,2,3,4,5,6} 5148 .ve 5149 5150 .seealso: `MatCreate()`, `MatCreateAIJ()`, `MatCreateSeqAIJ()`, `MatCreateSeqAIJWithArrays()`, `MatMPIAIJSetPreallocationCSR()`, `MatSetValuesCOO()` 5151 @*/ 5152 PetscErrorCode MatCreateSeqAIJFromTriple(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt i[], PetscInt j[], PetscScalar a[], Mat *mat, PetscInt nz, PetscBool idx) 5153 { 5154 PetscInt ii, *nnz, one = 1, row, col; 5155 5156 PetscFunctionBegin; 5157 PetscCall(PetscCalloc1(m, &nnz)); 5158 for (ii = 0; ii < nz; ii++) nnz[i[ii] - !!idx] += 1; 5159 PetscCall(MatCreate(comm, mat)); 5160 PetscCall(MatSetSizes(*mat, m, n, m, n)); 5161 PetscCall(MatSetType(*mat, MATSEQAIJ)); 5162 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(*mat, 0, nnz)); 5163 for (ii = 0; ii < nz; ii++) { 5164 if (idx) { 5165 row = i[ii] - 1; 5166 col = j[ii] - 1; 5167 } else { 5168 row = i[ii]; 5169 col = j[ii]; 5170 } 5171 PetscCall(MatSetValues(*mat, one, &row, one, &col, &a[ii], ADD_VALUES)); 5172 } 5173 PetscCall(MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY)); 5174 PetscCall(MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY)); 5175 PetscCall(PetscFree(nnz)); 5176 PetscFunctionReturn(0); 5177 } 5178 5179 PetscErrorCode MatSeqAIJInvalidateDiagonal(Mat A) 5180 { 5181 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 5182 5183 PetscFunctionBegin; 5184 a->idiagvalid = PETSC_FALSE; 5185 a->ibdiagvalid = PETSC_FALSE; 5186 5187 PetscCall(MatSeqAIJInvalidateDiagonal_Inode(A)); 5188 PetscFunctionReturn(0); 5189 } 5190 5191 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqAIJ(MPI_Comm comm, Mat inmat, PetscInt n, MatReuse scall, Mat *outmat) 5192 { 5193 PetscFunctionBegin; 5194 PetscCall(MatCreateMPIMatConcatenateSeqMat_MPIAIJ(comm, inmat, n, scall, outmat)); 5195 PetscFunctionReturn(0); 5196 } 5197 5198 /* 5199 Permute A into C's *local* index space using rowemb,colemb. 5200 The embedding are supposed to be injections and the above implies that the range of rowemb is a subset 5201 of [0,m), colemb is in [0,n). 5202 If pattern == DIFFERENT_NONZERO_PATTERN, C is preallocated according to A. 5203 */ 5204 PetscErrorCode MatSetSeqMat_SeqAIJ(Mat C, IS rowemb, IS colemb, MatStructure pattern, Mat B) 5205 { 5206 /* If making this function public, change the error returned in this function away from _PLIB. */ 5207 Mat_SeqAIJ *Baij; 5208 PetscBool seqaij; 5209 PetscInt m, n, *nz, i, j, count; 5210 PetscScalar v; 5211 const PetscInt *rowindices, *colindices; 5212 5213 PetscFunctionBegin; 5214 if (!B) PetscFunctionReturn(0); 5215 /* Check to make sure the target matrix (and embeddings) are compatible with C and each other. */ 5216 PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATSEQAIJ, &seqaij)); 5217 PetscCheck(seqaij, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Input matrix is of wrong type"); 5218 if (rowemb) { 5219 PetscCall(ISGetLocalSize(rowemb, &m)); 5220 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); 5221 } else { 5222 PetscCheck(C->rmap->n == B->rmap->n, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Input matrix is row-incompatible with the target matrix"); 5223 } 5224 if (colemb) { 5225 PetscCall(ISGetLocalSize(colemb, &n)); 5226 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); 5227 } else { 5228 PetscCheck(C->cmap->n == B->cmap->n, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Input matrix is col-incompatible with the target matrix"); 5229 } 5230 5231 Baij = (Mat_SeqAIJ *)(B->data); 5232 if (pattern == DIFFERENT_NONZERO_PATTERN) { 5233 PetscCall(PetscMalloc1(B->rmap->n, &nz)); 5234 for (i = 0; i < B->rmap->n; i++) nz[i] = Baij->i[i + 1] - Baij->i[i]; 5235 PetscCall(MatSeqAIJSetPreallocation(C, 0, nz)); 5236 PetscCall(PetscFree(nz)); 5237 } 5238 if (pattern == SUBSET_NONZERO_PATTERN) PetscCall(MatZeroEntries(C)); 5239 count = 0; 5240 rowindices = NULL; 5241 colindices = NULL; 5242 if (rowemb) PetscCall(ISGetIndices(rowemb, &rowindices)); 5243 if (colemb) PetscCall(ISGetIndices(colemb, &colindices)); 5244 for (i = 0; i < B->rmap->n; i++) { 5245 PetscInt row; 5246 row = i; 5247 if (rowindices) row = rowindices[i]; 5248 for (j = Baij->i[i]; j < Baij->i[i + 1]; j++) { 5249 PetscInt col; 5250 col = Baij->j[count]; 5251 if (colindices) col = colindices[col]; 5252 v = Baij->a[count]; 5253 PetscCall(MatSetValues(C, 1, &row, 1, &col, &v, INSERT_VALUES)); 5254 ++count; 5255 } 5256 } 5257 /* FIXME: set C's nonzerostate correctly. */ 5258 /* Assembly for C is necessary. */ 5259 C->preallocated = PETSC_TRUE; 5260 C->assembled = PETSC_TRUE; 5261 C->was_assembled = PETSC_FALSE; 5262 PetscFunctionReturn(0); 5263 } 5264 5265 PetscFunctionList MatSeqAIJList = NULL; 5266 5267 /*@C 5268 MatSeqAIJSetType - Converts a `MATSEQAIJ` matrix to a subtype 5269 5270 Collective on mat 5271 5272 Input Parameters: 5273 + mat - the matrix object 5274 - matype - matrix type 5275 5276 Options Database Key: 5277 . -mat_seqai_type <method> - for example seqaijcrl 5278 5279 Level: intermediate 5280 5281 .seealso: `PCSetType()`, `VecSetType()`, `MatCreate()`, `MatType`, `Mat` 5282 @*/ 5283 PetscErrorCode MatSeqAIJSetType(Mat mat, MatType matype) 5284 { 5285 PetscBool sametype; 5286 PetscErrorCode (*r)(Mat, MatType, MatReuse, Mat *); 5287 5288 PetscFunctionBegin; 5289 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 5290 PetscCall(PetscObjectTypeCompare((PetscObject)mat, matype, &sametype)); 5291 if (sametype) PetscFunctionReturn(0); 5292 5293 PetscCall(PetscFunctionListFind(MatSeqAIJList, matype, &r)); 5294 PetscCheck(r, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown Mat type given: %s", matype); 5295 PetscCall((*r)(mat, matype, MAT_INPLACE_MATRIX, &mat)); 5296 PetscFunctionReturn(0); 5297 } 5298 5299 /*@C 5300 MatSeqAIJRegister - - Adds a new sub-matrix type for sequential `MATSEQAIJ` matrices 5301 5302 Not Collective 5303 5304 Input Parameters: 5305 + name - name of a new user-defined matrix type, for example `MATSEQAIJCRL` 5306 - function - routine to convert to subtype 5307 5308 Notes: 5309 `MatSeqAIJRegister()` may be called multiple times to add several user-defined solvers. 5310 5311 Then, your matrix can be chosen with the procedural interface at runtime via the option 5312 $ -mat_seqaij_type my_mat 5313 5314 Level: advanced 5315 5316 .seealso: `MatSeqAIJRegisterAll()` 5317 @*/ 5318 PetscErrorCode MatSeqAIJRegister(const char sname[], PetscErrorCode (*function)(Mat, MatType, MatReuse, Mat *)) 5319 { 5320 PetscFunctionBegin; 5321 PetscCall(MatInitializePackage()); 5322 PetscCall(PetscFunctionListAdd(&MatSeqAIJList, sname, function)); 5323 PetscFunctionReturn(0); 5324 } 5325 5326 PetscBool MatSeqAIJRegisterAllCalled = PETSC_FALSE; 5327 5328 /*@C 5329 MatSeqAIJRegisterAll - Registers all of the matrix subtypes of `MATSSEQAIJ` 5330 5331 Not Collective 5332 5333 Level: advanced 5334 5335 .seealso: `MatRegisterAll()`, `MatSeqAIJRegister()` 5336 @*/ 5337 PetscErrorCode MatSeqAIJRegisterAll(void) 5338 { 5339 PetscFunctionBegin; 5340 if (MatSeqAIJRegisterAllCalled) PetscFunctionReturn(0); 5341 MatSeqAIJRegisterAllCalled = PETSC_TRUE; 5342 5343 PetscCall(MatSeqAIJRegister(MATSEQAIJCRL, MatConvert_SeqAIJ_SeqAIJCRL)); 5344 PetscCall(MatSeqAIJRegister(MATSEQAIJPERM, MatConvert_SeqAIJ_SeqAIJPERM)); 5345 PetscCall(MatSeqAIJRegister(MATSEQAIJSELL, MatConvert_SeqAIJ_SeqAIJSELL)); 5346 #if defined(PETSC_HAVE_MKL_SPARSE) 5347 PetscCall(MatSeqAIJRegister(MATSEQAIJMKL, MatConvert_SeqAIJ_SeqAIJMKL)); 5348 #endif 5349 #if defined(PETSC_HAVE_CUDA) 5350 PetscCall(MatSeqAIJRegister(MATSEQAIJCUSPARSE, MatConvert_SeqAIJ_SeqAIJCUSPARSE)); 5351 #endif 5352 #if defined(PETSC_HAVE_KOKKOS_KERNELS) 5353 PetscCall(MatSeqAIJRegister(MATSEQAIJKOKKOS, MatConvert_SeqAIJ_SeqAIJKokkos)); 5354 #endif 5355 #if defined(PETSC_HAVE_VIENNACL) && defined(PETSC_HAVE_VIENNACL_NO_CUDA) 5356 PetscCall(MatSeqAIJRegister(MATMPIAIJVIENNACL, MatConvert_SeqAIJ_SeqAIJViennaCL)); 5357 #endif 5358 PetscFunctionReturn(0); 5359 } 5360 5361 /* 5362 Special version for direct calls from Fortran 5363 */ 5364 #include <petsc/private/fortranimpl.h> 5365 #if defined(PETSC_HAVE_FORTRAN_CAPS) 5366 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 5367 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 5368 #define matsetvaluesseqaij_ matsetvaluesseqaij 5369 #endif 5370 5371 /* Change these macros so can be used in void function */ 5372 5373 /* Change these macros so can be used in void function */ 5374 /* Identical to PetscCallVoid, except it assigns to *_ierr */ 5375 #undef PetscCall 5376 #define PetscCall(...) \ 5377 do { \ 5378 PetscErrorCode ierr_msv_mpiaij = __VA_ARGS__; \ 5379 if (PetscUnlikely(ierr_msv_mpiaij)) { \ 5380 *_ierr = PetscError(PETSC_COMM_SELF, __LINE__, PETSC_FUNCTION_NAME, __FILE__, ierr_msv_mpiaij, PETSC_ERROR_REPEAT, " "); \ 5381 return; \ 5382 } \ 5383 } while (0) 5384 5385 #undef SETERRQ 5386 #define SETERRQ(comm, ierr, ...) \ 5387 do { \ 5388 *_ierr = PetscError(comm, __LINE__, PETSC_FUNCTION_NAME, __FILE__, ierr, PETSC_ERROR_INITIAL, __VA_ARGS__); \ 5389 return; \ 5390 } while (0) 5391 5392 PETSC_EXTERN void matsetvaluesseqaij_(Mat *AA, PetscInt *mm, const PetscInt im[], PetscInt *nn, const PetscInt in[], const PetscScalar v[], InsertMode *isis, PetscErrorCode *_ierr) 5393 { 5394 Mat A = *AA; 5395 PetscInt m = *mm, n = *nn; 5396 InsertMode is = *isis; 5397 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 5398 PetscInt *rp, k, low, high, t, ii, row, nrow, i, col, l, rmax, N; 5399 PetscInt *imax, *ai, *ailen; 5400 PetscInt *aj, nonew = a->nonew, lastcol = -1; 5401 MatScalar *ap, value, *aa; 5402 PetscBool ignorezeroentries = a->ignorezeroentries; 5403 PetscBool roworiented = a->roworiented; 5404 5405 PetscFunctionBegin; 5406 MatCheckPreallocated(A, 1); 5407 imax = a->imax; 5408 ai = a->i; 5409 ailen = a->ilen; 5410 aj = a->j; 5411 aa = a->a; 5412 5413 for (k = 0; k < m; k++) { /* loop over added rows */ 5414 row = im[k]; 5415 if (row < 0) continue; 5416 PetscCheck(row < A->rmap->n, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Row too large"); 5417 rp = aj + ai[row]; 5418 ap = aa + ai[row]; 5419 rmax = imax[row]; 5420 nrow = ailen[row]; 5421 low = 0; 5422 high = nrow; 5423 for (l = 0; l < n; l++) { /* loop over added columns */ 5424 if (in[l] < 0) continue; 5425 PetscCheck(in[l] < A->cmap->n, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Column too large"); 5426 col = in[l]; 5427 if (roworiented) value = v[l + k * n]; 5428 else value = v[k + l * m]; 5429 5430 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 5431 5432 if (col <= lastcol) low = 0; 5433 else high = nrow; 5434 lastcol = col; 5435 while (high - low > 5) { 5436 t = (low + high) / 2; 5437 if (rp[t] > col) high = t; 5438 else low = t; 5439 } 5440 for (i = low; i < high; i++) { 5441 if (rp[i] > col) break; 5442 if (rp[i] == col) { 5443 if (is == ADD_VALUES) ap[i] += value; 5444 else ap[i] = value; 5445 goto noinsert; 5446 } 5447 } 5448 if (value == 0.0 && ignorezeroentries) goto noinsert; 5449 if (nonew == 1) goto noinsert; 5450 PetscCheck(nonew != -1, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Inserting a new nonzero in the matrix"); 5451 MatSeqXAIJReallocateAIJ(A, A->rmap->n, 1, nrow, row, col, rmax, aa, ai, aj, rp, ap, imax, nonew, MatScalar); 5452 N = nrow++ - 1; 5453 a->nz++; 5454 high++; 5455 /* shift up all the later entries in this row */ 5456 for (ii = N; ii >= i; ii--) { 5457 rp[ii + 1] = rp[ii]; 5458 ap[ii + 1] = ap[ii]; 5459 } 5460 rp[i] = col; 5461 ap[i] = value; 5462 A->nonzerostate++; 5463 noinsert:; 5464 low = i + 1; 5465 } 5466 ailen[row] = nrow; 5467 } 5468 PetscFunctionReturnVoid(); 5469 } 5470 /* Undefining these here since they were redefined from their original definition above! No 5471 * other PETSc functions should be defined past this point, as it is impossible to recover the 5472 * original definitions */ 5473 #undef PetscCall 5474 #undef SETERRQ 5475