Actual source code: symtranspose.c
petsc-3.6.1 2015-08-06
2: /*
3: Defines symbolic transpose routines for SeqAIJ matrices.
5: Currently Get/Restore only allocates/frees memory for holding the
6: (i,j) info for the transpose. Someday, this info could be
7: maintained so successive calls to Get will not recompute the info.
9: Also defined is a "faster" implementation of MatTranspose for SeqAIJ
10: matrices which avoids calls to MatSetValues. This routine has not
11: been adopted as the standard yet as it is somewhat untested.
13: */
15: #include <../src/mat/impls/aij/seq/aij.h>
20: PetscErrorCode MatGetSymbolicTranspose_SeqAIJ(Mat A,PetscInt *Ati[],PetscInt *Atj[])
21: {
23: PetscInt i,j,anzj;
24: Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data;
25: PetscInt an=A->cmap->N,am=A->rmap->N;
26: PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j;
29: PetscInfo(A,"Getting Symbolic Transpose.\n");
31: /* Set up timers */
32: PetscLogEventBegin(MAT_Getsymtranspose,A,0,0,0);
34: /* Allocate space for symbolic transpose info and work array */
35: PetscCalloc1(an+1,&ati);
36: PetscMalloc1(ai[am],&atj);
37: PetscMalloc1(an,&atfill);
39: /* Walk through aj and count ## of non-zeros in each row of A^T. */
40: /* Note: offset by 1 for fast conversion into csr format. */
41: for (i=0;i<ai[am];i++) {
42: ati[aj[i]+1] += 1;
43: }
44: /* Form ati for csr format of A^T. */
45: for (i=0;i<an;i++) {
46: ati[i+1] += ati[i];
47: }
49: /* Copy ati into atfill so we have locations of the next free space in atj */
50: PetscMemcpy(atfill,ati,an*sizeof(PetscInt));
52: /* Walk through A row-wise and mark nonzero entries of A^T. */
53: for (i=0; i<am; i++) {
54: anzj = ai[i+1] - ai[i];
55: for (j=0; j<anzj; j++) {
56: atj[atfill[*aj]] = i;
57: atfill[*aj++] += 1;
58: }
59: }
61: /* Clean up temporary space and complete requests. */
62: PetscFree(atfill);
63: *Ati = ati;
64: *Atj = atj;
66: PetscLogEventEnd(MAT_Getsymtranspose,A,0,0,0);
67: return(0);
68: }
69: /*
70: MatGetSymbolicTransposeReduced_SeqAIJ() - Get symbolic matrix structure of submatrix A[rstart:rend,:],
71: modified from MatGetSymbolicTranspose_SeqAIJ()
72: */
75: PetscErrorCode MatGetSymbolicTransposeReduced_SeqAIJ(Mat A,PetscInt rstart,PetscInt rend,PetscInt *Ati[],PetscInt *Atj[])
76: {
78: PetscInt i,j,anzj;
79: Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data;
80: PetscInt an=A->cmap->N;
81: PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j;
84: PetscInfo(A,"Getting Symbolic Transpose\n");
85: PetscLogEventBegin(MAT_Getsymtransreduced,A,0,0,0);
87: /* Allocate space for symbolic transpose info and work array */
88: PetscCalloc1(an+1,&ati);
89: anzj = ai[rend] - ai[rstart];
90: PetscMalloc1(anzj+1,&atj);
91: PetscMalloc1(an+1,&atfill);
93: /* Walk through aj and count ## of non-zeros in each row of A^T. */
94: /* Note: offset by 1 for fast conversion into csr format. */
95: for (i=ai[rstart]; i<ai[rend]; i++) {
96: ati[aj[i]+1] += 1;
97: }
98: /* Form ati for csr format of A^T. */
99: for (i=0;i<an;i++) {
100: ati[i+1] += ati[i];
101: }
103: /* Copy ati into atfill so we have locations of the next free space in atj */
104: PetscMemcpy(atfill,ati,an*sizeof(PetscInt));
106: /* Walk through A row-wise and mark nonzero entries of A^T. */
107: aj = aj + ai[rstart];
108: for (i=rstart; i<rend; i++) {
109: anzj = ai[i+1] - ai[i];
110: for (j=0; j<anzj; j++) {
111: atj[atfill[*aj]] = i-rstart;
112: atfill[*aj++] += 1;
113: }
114: }
116: /* Clean up temporary space and complete requests. */
117: PetscFree(atfill);
118: *Ati = ati;
119: *Atj = atj;
121: PetscLogEventEnd(MAT_Getsymtransreduced,A,0,0,0);
122: return(0);
123: }
127: PetscErrorCode MatTranspose_SeqAIJ_FAST(Mat A,MatReuse reuse,Mat *B)
128: {
130: PetscInt i,j,anzj;
131: Mat At;
132: Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*at;
133: PetscInt an=A->cmap->N,am=A->rmap->N;
134: PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j;
135: MatScalar *ata,*aa=a->a;
138: PetscLogEventBegin(MAT_Transpose_SeqAIJ,A,0,0,0);
140: if (reuse == MAT_INITIAL_MATRIX || *B == A) {
141: /* Allocate space for symbolic transpose info and work array */
142: PetscCalloc1(an+1,&ati);
143: PetscMalloc1(ai[am],&atj);
144: PetscMalloc1(ai[am],&ata);
145: /* Walk through aj and count ## of non-zeros in each row of A^T. */
146: /* Note: offset by 1 for fast conversion into csr format. */
147: for (i=0;i<ai[am];i++) {
148: ati[aj[i]+1] += 1;
149: }
150: /* Form ati for csr format of A^T. */
151: for (i=0;i<an;i++) {
152: ati[i+1] += ati[i];
153: }
154: } else {
155: Mat_SeqAIJ *sub_B = (Mat_SeqAIJ*) (*B)->data;
156: ati = sub_B->i;
157: atj = sub_B->j;
158: ata = sub_B->a;
159: At = *B;
160: }
162: /* Copy ati into atfill so we have locations of the next free space in atj */
163: PetscMalloc1(an,&atfill);
164: PetscMemcpy(atfill,ati,an*sizeof(PetscInt));
166: /* Walk through A row-wise and mark nonzero entries of A^T. */
167: for (i=0;i<am;i++) {
168: anzj = ai[i+1] - ai[i];
169: for (j=0;j<anzj;j++) {
170: atj[atfill[*aj]] = i;
171: ata[atfill[*aj]] = *aa++;
172: atfill[*aj++] += 1;
173: }
174: }
176: /* Clean up temporary space and complete requests. */
177: PetscFree(atfill);
178: if (reuse == MAT_INITIAL_MATRIX) {
179: MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),an,am,ati,atj,ata,&At);
181: at = (Mat_SeqAIJ*)(At->data);
182: at->free_a = PETSC_TRUE;
183: at->free_ij = PETSC_TRUE;
184: at->nonew = 0;
185: }
187: if (reuse == MAT_INITIAL_MATRIX || *B != A) {
188: *B = At;
189: } else {
190: MatHeaderMerge(A,At);
191: }
192: PetscLogEventEnd(MAT_Transpose_SeqAIJ,A,0,0,0);
193: return(0);
194: }
198: PetscErrorCode MatRestoreSymbolicTranspose_SeqAIJ(Mat A,PetscInt *ati[],PetscInt *atj[])
199: {
203: PetscInfo(A,"Restoring Symbolic Transpose.\n");
204: PetscFree(*ati);
205: PetscFree(*atj);
206: return(0);
207: }