1 /*
2 Defines transpose routines for SeqAIJ matrices.
3 */
4
5 #include <../src/mat/impls/aij/seq/aij.h>
6
7 /*
8 The symbolic and full transpose versions share several similar code blocks but the macros to reuse the code would be confusing and ugly
9 */
MatTransposeSymbolic_SeqAIJ(Mat A,Mat * B)10 PetscErrorCode MatTransposeSymbolic_SeqAIJ(Mat A, Mat *B)
11 {
12 PetscInt i, j, anzj;
13 Mat At;
14 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *at;
15 PetscInt an = A->cmap->N, am = A->rmap->N;
16 PetscInt *ati, *atj, *atfill, *ai = a->i, *aj = a->j;
17
18 PetscFunctionBegin;
19 /* Allocate space for symbolic transpose info and work array */
20 PetscCall(PetscCalloc1(an + 1, &ati));
21 PetscCall(PetscMalloc1(ai[am], &atj));
22
23 /* Walk through aj and count ## of non-zeros in each row of A^T. */
24 /* Note: offset by 1 for fast conversion into csr format. */
25 for (i = 0; i < ai[am]; i++) ati[aj[i] + 1] += 1;
26 /* Form ati for csr format of A^T. */
27 for (i = 0; i < an; i++) ati[i + 1] += ati[i];
28
29 /* Copy ati into atfill so we have locations of the next free space in atj */
30 PetscCall(PetscMalloc1(an, &atfill));
31 PetscCall(PetscArraycpy(atfill, ati, an));
32
33 /* Walk through A row-wise and mark nonzero entries of A^T. */
34 for (i = 0; i < am; i++) {
35 anzj = ai[i + 1] - ai[i];
36 for (j = 0; j < anzj; j++) {
37 atj[atfill[*aj]] = i;
38 atfill[*aj++] += 1;
39 }
40 }
41 PetscCall(PetscFree(atfill));
42
43 PetscCall(MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A), an, am, ati, atj, NULL, &At));
44 PetscCall(MatSetBlockSizes(At, A->cmap->bs, A->rmap->bs));
45 PetscCall(MatSetType(At, ((PetscObject)A)->type_name));
46 at = (Mat_SeqAIJ *)At->data;
47 at->free_a = PETSC_FALSE;
48 at->free_ij = PETSC_TRUE;
49 at->nonew = 0;
50 at->maxnz = ati[an];
51 *B = At;
52 PetscFunctionReturn(PETSC_SUCCESS);
53 }
54
MatTranspose_SeqAIJ(Mat A,MatReuse reuse,Mat * B)55 PetscErrorCode MatTranspose_SeqAIJ(Mat A, MatReuse reuse, Mat *B)
56 {
57 PetscInt i, j, anzj;
58 Mat At;
59 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *at;
60 PetscInt an = A->cmap->N, am = A->rmap->N;
61 PetscInt *ati, *atj, *atfill, *ai = a->i, *aj = a->j;
62 MatScalar *ata;
63 const MatScalar *aa, *av;
64 PetscContainer rB;
65 MatParentState *rb;
66 PetscBool nonzerochange = PETSC_FALSE;
67
68 PetscFunctionBegin;
69 if (reuse == MAT_REUSE_MATRIX) {
70 PetscCall(PetscObjectQuery((PetscObject)*B, "MatTransposeParent", (PetscObject *)&rB));
71 PetscCheck(rB, PetscObjectComm((PetscObject)*B), PETSC_ERR_ARG_WRONG, "Reuse matrix used was not generated from call to MatTranspose()");
72 PetscCall(PetscContainerGetPointer(rB, &rb));
73 if (rb->nonzerostate != A->nonzerostate) nonzerochange = PETSC_TRUE;
74 }
75
76 PetscCall(MatSeqAIJGetArrayRead(A, &av));
77 aa = av;
78 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX || nonzerochange) {
79 /* Allocate space for symbolic transpose info and work array */
80 PetscCall(PetscCalloc1(an + 1, &ati));
81 PetscCall(PetscMalloc1(ai[am], &atj));
82 /* Walk through aj and count ## of non-zeros in each row of A^T. */
83 /* Note: offset by 1 for fast conversion into csr format. */
84 for (i = 0; i < ai[am]; i++) ati[aj[i] + 1] += 1;
85 /* Form ati for csr format of A^T. */
86 for (i = 0; i < an; i++) ati[i + 1] += ati[i];
87 PetscCall(PetscMalloc1(ai[am], &ata));
88 } else {
89 Mat_SeqAIJ *sub_B = (Mat_SeqAIJ *)(*B)->data;
90 ati = sub_B->i;
91 atj = sub_B->j;
92 ata = sub_B->a;
93 At = *B;
94 }
95
96 /* Copy ati into atfill so we have locations of the next free space in atj */
97 PetscCall(PetscMalloc1(an, &atfill));
98 PetscCall(PetscArraycpy(atfill, ati, an));
99
100 /* Walk through A row-wise and mark nonzero entries of A^T. */
101 if (aa) {
102 for (i = 0; i < am; i++) {
103 anzj = ai[i + 1] - ai[i];
104 for (j = 0; j < anzj; j++) {
105 atj[atfill[*aj]] = i;
106 ata[atfill[*aj]] = *aa++;
107 atfill[*aj++] += 1;
108 }
109 }
110 } else {
111 for (i = 0; i < am; i++) {
112 anzj = ai[i + 1] - ai[i];
113 for (j = 0; j < anzj; j++) {
114 atj[atfill[*aj]] = i;
115 atfill[*aj++] += 1;
116 }
117 }
118 }
119 PetscCall(PetscFree(atfill));
120 PetscCall(MatSeqAIJRestoreArrayRead(A, &av));
121 if (reuse == MAT_REUSE_MATRIX) PetscCall(PetscObjectStateIncrease((PetscObject)*B));
122
123 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX || nonzerochange) {
124 PetscCall(MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A), an, am, ati, atj, ata, &At));
125 PetscCall(MatSetBlockSizes(At, A->cmap->bs, A->rmap->bs));
126 PetscCall(MatSetType(At, ((PetscObject)A)->type_name));
127 at = (Mat_SeqAIJ *)At->data;
128 at->free_a = PETSC_TRUE;
129 at->free_ij = PETSC_TRUE;
130 at->nonew = 0;
131 at->maxnz = ati[an];
132 }
133
134 if (reuse == MAT_INITIAL_MATRIX || (reuse == MAT_REUSE_MATRIX && !nonzerochange)) {
135 *B = At;
136 } else if (nonzerochange) {
137 PetscCall(MatHeaderMerge(*B, &At));
138 PetscCall(MatTransposeSetPrecursor(A, *B));
139 } else if (reuse == MAT_INPLACE_MATRIX) {
140 PetscCall(MatHeaderMerge(A, &At));
141 }
142 PetscFunctionReturn(PETSC_SUCCESS);
143 }
144
145 /*
146 Get symbolic matrix nonzero structure of a submatrix of A, A[rstart:rend,:],
147 */
MatGetSymbolicTransposeReduced_SeqAIJ(Mat A,PetscInt rstart,PetscInt rend,PetscInt * Ati[],PetscInt * Atj[])148 PetscErrorCode MatGetSymbolicTransposeReduced_SeqAIJ(Mat A, PetscInt rstart, PetscInt rend, PetscInt *Ati[], PetscInt *Atj[])
149 {
150 PetscInt i, j, anzj;
151 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
152 PetscInt an = A->cmap->N;
153 PetscInt *ati, *atj, *atfill, *ai = a->i, *aj = a->j, am = ai[rend] - ai[rstart];
154
155 PetscFunctionBegin;
156 PetscCall(PetscLogEventBegin(MAT_Getsymtransreduced, A, 0, 0, 0));
157
158 /* Allocate space for symbolic transpose info and work array */
159 PetscCall(PetscCalloc1(an + 1, &ati));
160 PetscCall(PetscMalloc1(am + 1, &atj));
161
162 /* Walk through aj and count ## of non-zeros in each row of A^T. */
163 /* Note: offset by 1 for fast conversion into csr format. */
164 for (i = ai[rstart]; i < ai[rend]; i++) ati[aj[i] + 1] += 1;
165 /* Form ati for csr format of A^T. */
166 for (i = 0; i < an; i++) ati[i + 1] += ati[i];
167
168 /* Copy ati into atfill so we have locations of the next free space in atj */
169 PetscCall(PetscMalloc1(an + 1, &atfill));
170 PetscCall(PetscArraycpy(atfill, ati, an));
171
172 /* Walk through A row-wise and mark nonzero entries of A^T. */
173 aj = PetscSafePointerPlusOffset(aj, ai[rstart]);
174 for (i = rstart; i < rend; i++) {
175 anzj = ai[i + 1] - ai[i];
176 for (j = 0; j < anzj; j++) {
177 atj[atfill[*aj]] = i - rstart;
178 atfill[*aj++] += 1;
179 }
180 }
181 PetscCall(PetscFree(atfill));
182 *Ati = ati;
183 *Atj = atj;
184
185 PetscCall(PetscLogEventEnd(MAT_Getsymtransreduced, A, 0, 0, 0));
186 PetscFunctionReturn(PETSC_SUCCESS);
187 }
188
189 /*
190 Returns the i and j arrays for a symbolic transpose, this is used internally within SeqAIJ code when the full
191 symbolic matrix (which can be obtained with MatTransposeSymbolic() is not needed. MatRestoreSymbolicTranspose_SeqAIJ() should be used to free the arrays.
192 */
MatGetSymbolicTranspose_SeqAIJ(Mat A,PetscInt * Ati[],PetscInt * Atj[])193 PetscErrorCode MatGetSymbolicTranspose_SeqAIJ(Mat A, PetscInt *Ati[], PetscInt *Atj[])
194 {
195 PetscFunctionBegin;
196 PetscCall(MatGetSymbolicTransposeReduced_SeqAIJ(A, 0, A->rmap->N, Ati, Atj));
197 PetscFunctionReturn(PETSC_SUCCESS);
198 }
199
MatRestoreSymbolicTranspose_SeqAIJ(Mat A,PetscInt * ati[],PetscInt * atj[])200 PetscErrorCode MatRestoreSymbolicTranspose_SeqAIJ(Mat A, PetscInt *ati[], PetscInt *atj[])
201 {
202 PetscFunctionBegin;
203 PetscCall(PetscFree(*ati));
204 PetscCall(PetscFree(*atj));
205 PetscFunctionReturn(PETSC_SUCCESS);
206 }
207