Actual source code: mpiaijsbaij.c
2: #include <../src/mat/impls/sbaij/mpi/mpisbaij.h>
3: #include <../src/mat/impls/aij/mpi/mpiaij.h>
4: #include <petsc/private/matimpl.h>
5: #include <petscmat.h>
7: PETSC_INTERN PetscErrorCode MatConvert_MPIAIJ_MPISBAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
8: {
9: Mat M;
10: Mat_MPIAIJ *mpimat = (Mat_MPIAIJ *)A->data;
11: Mat_SeqAIJ *Aa = (Mat_SeqAIJ *)mpimat->A->data, *Ba = (Mat_SeqAIJ *)mpimat->B->data;
12: PetscInt *d_nnz, *o_nnz;
13: PetscInt i, j, nz;
14: PetscInt m, n, lm, ln;
15: PetscInt rstart, rend, bs = PetscAbs(A->rmap->bs);
16: const PetscScalar *vwork;
17: const PetscInt *cwork;
20: if (reuse != MAT_REUSE_MATRIX) {
21: MatGetSize(A, &m, &n);
22: MatGetLocalSize(A, &lm, &ln);
23: PetscMalloc2(lm / bs, &d_nnz, lm / bs, &o_nnz);
25: MatMarkDiagonal_SeqAIJ(mpimat->A);
26: for (i = 0; i < lm / bs; i++) {
27: if (Aa->i[i * bs + 1] == Aa->diag[i * bs]) { /* misses diagonal entry */
28: d_nnz[i] = (Aa->i[i * bs + 1] - Aa->i[i * bs]) / bs;
29: } else {
30: d_nnz[i] = (Aa->i[i * bs + 1] - Aa->diag[i * bs]) / bs;
31: }
32: o_nnz[i] = (Ba->i[i * bs + 1] - Ba->i[i * bs]) / bs;
33: }
35: MatCreate(PetscObjectComm((PetscObject)A), &M);
36: MatSetSizes(M, lm, ln, m, n);
37: MatSetType(M, MATMPISBAIJ);
38: MatSeqSBAIJSetPreallocation(M, bs, 0, d_nnz);
39: MatMPISBAIJSetPreallocation(M, bs, 0, d_nnz, 0, o_nnz);
40: PetscFree2(d_nnz, o_nnz);
41: } else M = *newmat;
43: if (bs == 1) {
44: MatGetOwnershipRange(A, &rstart, &rend);
45: for (i = rstart; i < rend; i++) {
46: MatGetRow(A, i, &nz, &cwork, &vwork);
47: if (nz) {
48: j = 0;
49: while (cwork[j] < i) {
50: j++;
51: nz--;
52: }
53: MatSetValues(M, 1, &i, nz, cwork + j, vwork + j, INSERT_VALUES);
54: }
55: MatRestoreRow(A, i, &nz, &cwork, &vwork);
56: }
57: MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
58: MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);
59: } else {
60: MatSetOption(M, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE);
61: /* reuse may not be equal to MAT_REUSE_MATRIX, but the basic converter will reallocate or replace newmat if this value is not used */
62: /* if reuse is equal to MAT_INITIAL_MATRIX, it has been appropriately preallocated before */
63: /* MAT_INPLACE_MATRIX, it will be replaced with MatHeaderReplace below */
64: MatConvert_Basic(A, newtype, MAT_REUSE_MATRIX, &M);
65: }
67: if (reuse == MAT_INPLACE_MATRIX) {
68: MatHeaderReplace(A, &M);
69: } else *newmat = M;
70: return 0;
71: }
72: /* contributed by Dahai Guo <dhguo@ncsa.uiuc.edu> April 2011 */
73: PETSC_INTERN PetscErrorCode MatConvert_MPIBAIJ_MPISBAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
74: {
75: Mat M;
76: Mat_MPIBAIJ *mpimat = (Mat_MPIBAIJ *)A->data;
77: Mat_SeqBAIJ *Aa = (Mat_SeqBAIJ *)mpimat->A->data, *Ba = (Mat_SeqBAIJ *)mpimat->B->data;
78: PetscInt *d_nnz, *o_nnz;
79: PetscInt i, nz;
80: PetscInt m, n, lm, ln;
81: PetscInt rstart, rend;
82: const PetscScalar *vwork;
83: const PetscInt *cwork;
84: PetscInt bs = A->rmap->bs;
86: if (reuse != MAT_REUSE_MATRIX) {
87: MatGetSize(A, &m, &n);
88: MatGetLocalSize(A, &lm, &ln);
89: PetscMalloc2(lm / bs, &d_nnz, lm / bs, &o_nnz);
91: MatMarkDiagonal_SeqBAIJ(mpimat->A);
92: for (i = 0; i < lm / bs; i++) {
93: d_nnz[i] = Aa->i[i + 1] - Aa->diag[i];
94: o_nnz[i] = Ba->i[i + 1] - Ba->i[i];
95: }
97: MatCreate(PetscObjectComm((PetscObject)A), &M);
98: MatSetSizes(M, lm, ln, m, n);
99: MatSetType(M, MATMPISBAIJ);
100: MatSeqSBAIJSetPreallocation(M, bs, 0, d_nnz);
101: MatMPISBAIJSetPreallocation(M, bs, 0, d_nnz, 0, o_nnz);
103: PetscFree2(d_nnz, o_nnz);
104: } else M = *newmat;
106: MatGetOwnershipRange(A, &rstart, &rend);
107: MatSetOption(M, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE);
108: for (i = rstart; i < rend; i++) {
109: MatGetRow(A, i, &nz, &cwork, &vwork);
110: MatSetValues(M, 1, &i, nz, cwork, vwork, INSERT_VALUES);
111: MatRestoreRow(A, i, &nz, &cwork, &vwork);
112: }
113: MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
114: MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);
116: if (reuse == MAT_INPLACE_MATRIX) {
117: MatHeaderReplace(A, &M);
118: } else *newmat = M;
119: return 0;
120: }