Actual source code: aijbaij.c
1: #include <../src/mat/impls/baij/seq/baij.h>
3: PETSC_INTERN PetscErrorCode MatConvert_SeqBAIJ_SeqAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
4: {
5: Mat B;
6: Mat_SeqAIJ *b;
7: PetscBool roworiented;
8: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
9: PetscInt bs = A->rmap->bs, *ai = a->i, *aj = a->j, n = A->rmap->N / bs, i, j, k;
10: PetscInt *rowlengths, *rows, *cols, maxlen = 0, ncols;
11: MatScalar *aa = a->a;
13: PetscFunctionBegin;
14: if (reuse == MAT_REUSE_MATRIX) {
15: B = *newmat;
16: for (i = 0; i < n; i++) maxlen = PetscMax(maxlen, (ai[i + 1] - ai[i]));
17: } else {
18: PetscCall(PetscMalloc1(n * bs, &rowlengths));
19: for (i = 0; i < n; i++) {
20: maxlen = PetscMax(maxlen, (ai[i + 1] - ai[i]));
21: for (j = 0; j < bs; j++) rowlengths[i * bs + j] = bs * (ai[i + 1] - ai[i]);
22: }
23: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
24: PetscCall(MatSetType(B, MATSEQAIJ));
25: PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
26: PetscCall(MatSetBlockSizes(B, A->rmap->bs, A->cmap->bs));
27: PetscCall(MatSeqAIJSetPreallocation(B, 0, rowlengths));
28: PetscCall(PetscFree(rowlengths));
29: }
30: b = (Mat_SeqAIJ *)B->data;
31: roworiented = b->roworiented;
33: PetscCall(MatSetOption(B, MAT_ROW_ORIENTED, PETSC_FALSE));
34: PetscCall(PetscMalloc1(bs, &rows));
35: PetscCall(PetscMalloc1(bs * maxlen, &cols));
36: for (i = 0; i < n; i++) {
37: for (j = 0; j < bs; j++) rows[j] = i * bs + j;
38: ncols = ai[i + 1] - ai[i];
39: for (k = 0; k < ncols; k++) {
40: for (j = 0; j < bs; j++) cols[k * bs + j] = bs * (*aj) + j;
41: aj++;
42: }
43: PetscCall(MatSetValues(B, bs, rows, bs * ncols, cols, aa, INSERT_VALUES));
44: aa = PetscSafePointerPlusOffset(aa, ncols * bs * bs);
45: }
46: PetscCall(PetscFree(cols));
47: PetscCall(PetscFree(rows));
48: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
49: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
50: PetscCall(MatSetOption(B, MAT_ROW_ORIENTED, roworiented));
52: if (reuse == MAT_INPLACE_MATRIX) {
53: PetscCall(MatHeaderReplace(A, &B));
54: } else *newmat = B;
55: PetscFunctionReturn(PETSC_SUCCESS);
56: }
58: #include <../src/mat/impls/aij/seq/aij.h>
60: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqBAIJ_Preallocate(Mat A, PetscInt **nnz)
61: {
62: Mat_SeqAIJ *Aa = (Mat_SeqAIJ *)A->data;
63: PetscInt m, n, bs = PetscAbs(A->rmap->bs);
64: const PetscInt *ai = Aa->i, *aj = Aa->j;
66: PetscFunctionBegin;
67: PetscCall(MatGetSize(A, &m, &n));
69: if (bs == 1) {
70: PetscCall(PetscMalloc1(m, nnz));
71: for (PetscInt i = 0; i < m; i++) (*nnz)[i] = ai[i + 1] - ai[i];
72: } else {
73: PetscHSetIJ ht;
74: PetscInt k;
75: PetscHashIJKey key;
76: PetscBool missing;
78: PetscCall(PetscHSetIJCreate(&ht));
79: PetscCall(PetscCalloc1(m / bs, nnz));
80: for (PetscInt i = 0; i < m; i++) {
81: key.i = i / bs;
82: for (k = ai[i]; k < ai[i + 1]; k++) {
83: key.j = aj[k] / bs;
84: PetscCall(PetscHSetIJQueryAdd(ht, key, &missing));
85: if (missing) (*nnz)[key.i]++;
86: }
87: }
88: PetscCall(PetscHSetIJDestroy(&ht));
89: }
90: PetscFunctionReturn(PETSC_SUCCESS);
91: }
93: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqBAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
94: {
95: Mat B;
96: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
97: Mat_SeqBAIJ *b;
98: PetscInt m = A->rmap->N, n = A->cmap->N, *rowlengths, bs = PetscAbs(A->rmap->bs);
100: PetscFunctionBegin;
101: if (reuse != MAT_REUSE_MATRIX) {
102: PetscCall(MatConvert_SeqAIJ_SeqBAIJ_Preallocate(A, &rowlengths));
103: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
104: PetscCall(MatSetSizes(B, m, n, m, n));
105: PetscCall(MatSetType(B, MATSEQBAIJ));
106: PetscCall(MatSeqBAIJSetPreallocation(B, bs, 0, rowlengths));
107: PetscCall(PetscFree(rowlengths));
108: } else B = *newmat;
110: if (bs == 1) {
111: b = (Mat_SeqBAIJ *)B->data;
113: PetscCall(PetscArraycpy(b->i, a->i, m + 1));
114: PetscCall(PetscArraycpy(b->ilen, a->ilen, m));
115: PetscCall(PetscArraycpy(b->j, a->j, a->nz));
116: PetscCall(PetscArraycpy(b->a, a->a, a->nz));
118: PetscCall(MatSetOption(B, MAT_ROW_ORIENTED, PETSC_TRUE));
119: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
120: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
121: } else {
122: /* reuse may not be equal to MAT_REUSE_MATRIX, but the basic converter will reallocate or replace newmat if this value is not used */
123: /* if reuse is equal to MAT_INITIAL_MATRIX, it has been appropriately preallocated before */
124: /* MAT_INPLACE_MATRIX, it will be replaced with MatHeaderReplace below */
125: PetscCall(MatConvert_Basic(A, newtype, MAT_REUSE_MATRIX, &B));
126: }
128: if (reuse == MAT_INPLACE_MATRIX) {
129: PetscCall(MatHeaderReplace(A, &B));
130: } else *newmat = B;
131: PetscFunctionReturn(PETSC_SUCCESS);
132: }