Actual source code: aijbaij.c
2: #include <../src/mat/impls/baij/seq/baij.h>
4: PETSC_INTERN PetscErrorCode MatConvert_SeqBAIJ_SeqAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
5: {
6: Mat B;
7: Mat_SeqAIJ *b;
8: PetscBool roworiented;
9: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
10: PetscInt bs = A->rmap->bs, *ai = a->i, *aj = a->j, n = A->rmap->N / bs, i, j, k;
11: PetscInt *rowlengths, *rows, *cols, maxlen = 0, ncols;
12: MatScalar *aa = a->a;
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: 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: MatCreate(PetscObjectComm((PetscObject)A), &B);
24: MatSetType(B, MATSEQAIJ);
25: MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N);
26: MatSetBlockSizes(B, A->rmap->bs, A->cmap->bs);
27: MatSeqAIJSetPreallocation(B, 0, rowlengths);
28: PetscFree(rowlengths);
29: }
30: b = (Mat_SeqAIJ *)B->data;
31: roworiented = b->roworiented;
33: MatSetOption(B, MAT_ROW_ORIENTED, PETSC_FALSE);
34: PetscMalloc1(bs, &rows);
35: 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: MatSetValues(B, bs, rows, bs * ncols, cols, aa, INSERT_VALUES);
44: aa += ncols * bs * bs;
45: }
46: PetscFree(cols);
47: PetscFree(rows);
48: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
49: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
50: MatSetOption(B, MAT_ROW_ORIENTED, roworiented);
52: if (reuse == MAT_INPLACE_MATRIX) {
53: MatHeaderReplace(A, &B);
54: } else *newmat = B;
55: return 0;
56: }
58: #include <../src/mat/impls/aij/seq/aij.h>
60: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqBAIJ_Preallocate(Mat A, PetscInt **nnz)
61: {
62: PetscInt m, n, bs = PetscAbs(A->rmap->bs);
63: PetscInt *ii;
64: Mat_SeqAIJ *Aa = (Mat_SeqAIJ *)A->data;
66: MatGetSize(A, &m, &n);
67: PetscCalloc1(m / bs, nnz);
68: PetscMalloc1(bs, &ii);
70: /* loop over all potential blocks in the matrix to see if the potential block has a nonzero */
71: for (PetscInt i = 0; i < m / bs; i++) {
72: for (PetscInt k = 0; k < bs; k++) ii[k] = Aa->i[i * bs + k];
73: for (PetscInt j = 0; j < n / bs; j++) {
74: for (PetscInt k = 0; k < bs; k++) {
75: for (; ii[k] < Aa->i[i * bs + k + 1] && Aa->j[ii[k]] < (j + 1) * bs; ii[k]++) {
76: if (Aa->j[ii[k]] >= j * bs) {
77: /* found a nonzero in the potential block so allocate for that block and jump to check the next potential block */
78: (*nnz)[i]++;
79: goto theend;
80: }
81: }
82: }
83: theend:;
84: }
85: }
86: PetscFree(ii);
87: return 0;
88: }
90: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqBAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
91: {
92: Mat B;
93: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
94: Mat_SeqBAIJ *b;
95: PetscInt m = A->rmap->N, n = A->cmap->N, *rowlengths, bs = PetscAbs(A->rmap->bs);
97: if (reuse != MAT_REUSE_MATRIX) {
98: MatConvert_SeqAIJ_SeqBAIJ_Preallocate(A, &rowlengths);
99: MatCreate(PetscObjectComm((PetscObject)A), &B);
100: MatSetSizes(B, m, n, m, n);
101: MatSetType(B, MATSEQBAIJ);
102: MatSeqBAIJSetPreallocation(B, bs, 0, rowlengths);
103: PetscFree(rowlengths);
104: } else B = *newmat;
106: if (bs == 1) {
107: b = (Mat_SeqBAIJ *)(B->data);
109: PetscArraycpy(b->i, a->i, m + 1);
110: PetscArraycpy(b->ilen, a->ilen, m);
111: PetscArraycpy(b->j, a->j, a->nz);
112: PetscArraycpy(b->a, a->a, a->nz);
114: MatSetOption(B, MAT_ROW_ORIENTED, PETSC_TRUE);
115: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
116: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
117: } else {
118: /* reuse may not be equal to MAT_REUSE_MATRIX, but the basic converter will reallocate or replace newmat if this value is not used */
119: /* if reuse is equal to MAT_INITIAL_MATRIX, it has been appropriately preallocated before */
120: /* MAT_INPLACE_MATRIX, it will be replaced with MatHeaderReplace below */
121: MatConvert_Basic(A, newtype, MAT_REUSE_MATRIX, &B);
122: }
124: if (reuse == MAT_INPLACE_MATRIX) {
125: MatHeaderReplace(A, &B);
126: } else *newmat = B;
127: return 0;
128: }