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