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: /* Allocate space for symbolic transpose info and work array */
20: PetscCalloc1(an + 1, &ati);
21: PetscMalloc1(ai[am], &atj);
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];
29: /* Copy ati into atfill so we have locations of the next free space in atj */
30: PetscMalloc1(an, &atfill);
31: PetscArraycpy(atfill, ati, an);
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: PetscFree(atfill);
43: MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A), an, am, ati, atj, NULL, &At);
44: MatSetBlockSizes(At, PetscAbs(A->cmap->bs), PetscAbs(A->rmap->bs));
45: 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: return 0;
53: }
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;
68: if (reuse == MAT_REUSE_MATRIX) {
69: PetscObjectQuery((PetscObject)*B, "MatTransposeParent", (PetscObject *)&rB);
71: PetscContainerGetPointer(rB, (void **)&rb);
72: if (rb->nonzerostate != A->nonzerostate) nonzerochange = PETSC_TRUE;
73: }
75: MatSeqAIJGetArrayRead(A, &av);
76: aa = av;
77: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX || nonzerochange) {
78: /* Allocate space for symbolic transpose info and work array */
79: PetscCalloc1(an + 1, &ati);
80: PetscMalloc1(ai[am], &atj);
81: /* Walk through aj and count ## of non-zeros in each row of A^T. */
82: /* Note: offset by 1 for fast conversion into csr format. */
83: for (i = 0; i < ai[am]; i++) ati[aj[i] + 1] += 1;
84: /* Form ati for csr format of A^T. */
85: for (i = 0; i < an; i++) ati[i + 1] += ati[i];
86: PetscMalloc1(ai[am], &ata);
87: } else {
88: Mat_SeqAIJ *sub_B = (Mat_SeqAIJ *)(*B)->data;
89: ati = sub_B->i;
90: atj = sub_B->j;
91: ata = sub_B->a;
92: At = *B;
93: }
95: /* Copy ati into atfill so we have locations of the next free space in atj */
96: PetscMalloc1(an, &atfill);
97: PetscArraycpy(atfill, ati, an);
99: /* Walk through A row-wise and mark nonzero entries of A^T. */
100: for (i = 0; i < am; i++) {
101: anzj = ai[i + 1] - ai[i];
102: for (j = 0; j < anzj; j++) {
103: atj[atfill[*aj]] = i;
104: ata[atfill[*aj]] = *aa++;
105: atfill[*aj++] += 1;
106: }
107: }
108: PetscFree(atfill);
109: MatSeqAIJRestoreArrayRead(A, &av);
110: if (reuse == MAT_REUSE_MATRIX) PetscObjectStateIncrease((PetscObject)(*B));
112: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX || nonzerochange) {
113: MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A), an, am, ati, atj, ata, &At);
114: MatSetBlockSizes(At, PetscAbs(A->cmap->bs), PetscAbs(A->rmap->bs));
115: MatSetType(At, ((PetscObject)A)->type_name);
116: at = (Mat_SeqAIJ *)(At->data);
117: at->free_a = PETSC_TRUE;
118: at->free_ij = PETSC_TRUE;
119: at->nonew = 0;
120: at->maxnz = ati[an];
121: }
123: if (reuse == MAT_INITIAL_MATRIX || (reuse == MAT_REUSE_MATRIX && !nonzerochange)) {
124: *B = At;
125: } else if (nonzerochange) {
126: MatHeaderMerge(*B, &At);
127: MatTransposeSetPrecursor(A, *B);
128: } else if (reuse == MAT_INPLACE_MATRIX) {
129: MatHeaderMerge(A, &At);
130: }
131: return 0;
132: }
134: /*
135: Get symbolic matrix structure of a submatrix of A, A[rstart:rend,:],
136: */
137: PetscErrorCode MatGetSymbolicTransposeReduced_SeqAIJ(Mat A, PetscInt rstart, PetscInt rend, PetscInt *Ati[], PetscInt *Atj[])
138: {
139: PetscInt i, j, anzj;
140: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
141: PetscInt an = A->cmap->N;
142: PetscInt *ati, *atj, *atfill, *ai = a->i, *aj = a->j, am = ai[rend] - ai[rstart];
144: PetscLogEventBegin(MAT_Getsymtransreduced, A, 0, 0, 0);
146: /* Allocate space for symbolic transpose info and work array */
147: PetscCalloc1(an + 1, &ati);
148: PetscMalloc1(am + 1, &atj);
150: /* Walk through aj and count ## of non-zeros in each row of A^T. */
151: /* Note: offset by 1 for fast conversion into csr format. */
152: for (i = ai[rstart]; i < ai[rend]; i++) ati[aj[i] + 1] += 1;
153: /* Form ati for csr format of A^T. */
154: for (i = 0; i < an; i++) ati[i + 1] += ati[i];
156: /* Copy ati into atfill so we have locations of the next free space in atj */
157: PetscMalloc1(an + 1, &atfill);
158: PetscArraycpy(atfill, ati, an);
160: /* Walk through A row-wise and mark nonzero entries of A^T. */
161: aj = aj + ai[rstart];
162: for (i = rstart; i < rend; i++) {
163: anzj = ai[i + 1] - ai[i];
164: for (j = 0; j < anzj; j++) {
165: atj[atfill[*aj]] = i - rstart;
166: atfill[*aj++] += 1;
167: }
168: }
169: PetscFree(atfill);
170: *Ati = ati;
171: *Atj = atj;
173: PetscLogEventEnd(MAT_Getsymtransreduced, A, 0, 0, 0);
174: return 0;
175: }
177: /*
178: Returns the i and j arrays for a symbolic transpose, this is used internally within SeqAIJ code when the full
179: symbolic matrix (which can be obtained with MatTransposeSymbolic() is not needed. MatRestoreSymbolicTranspose_SeqAIJ() should be used to free the arrays.
180: */
181: PetscErrorCode MatGetSymbolicTranspose_SeqAIJ(Mat A, PetscInt *Ati[], PetscInt *Atj[])
182: {
183: MatGetSymbolicTransposeReduced_SeqAIJ(A, 0, A->rmap->N, Ati, Atj);
184: return 0;
185: }
187: PetscErrorCode MatRestoreSymbolicTranspose_SeqAIJ(Mat A, PetscInt *ati[], PetscInt *atj[])
188: {
189: PetscFree(*ati);
190: PetscFree(*atj);
191: return 0;
192: }