Actual source code: symtranspose.c
petsc-3.13.6 2020-09-29
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 is the new
11: standard since it is much faster than MatTranspose_AIJ.
13: */
15: #include <../src/mat/impls/aij/seq/aij.h>
18: PetscErrorCode MatGetSymbolicTranspose_SeqAIJ(Mat A,PetscInt *Ati[],PetscInt *Atj[])
19: {
21: PetscInt i,j,anzj;
22: Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data;
23: PetscInt an=A->cmap->N,am=A->rmap->N;
24: PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j;
27: PetscInfo(A,"Getting Symbolic Transpose.\n");
29: /* Set up timers */
30: PetscLogEventBegin(MAT_Getsymtranspose,A,0,0,0);
32: /* Allocate space for symbolic transpose info and work array */
33: PetscCalloc1(an+1,&ati);
34: PetscMalloc1(ai[am],&atj);
35: PetscMalloc1(an,&atfill);
37: /* Walk through aj and count ## of non-zeros in each row of A^T. */
38: /* Note: offset by 1 for fast conversion into csr format. */
39: for (i=0;i<ai[am];i++) {
40: ati[aj[i]+1] += 1;
41: }
42: /* Form ati for csr format of A^T. */
43: for (i=0;i<an;i++) {
44: ati[i+1] += ati[i];
45: }
47: /* Copy ati into atfill so we have locations of the next free space in atj */
48: PetscArraycpy(atfill,ati,an);
50: /* Walk through A row-wise and mark nonzero entries of A^T. */
51: for (i=0; i<am; i++) {
52: anzj = ai[i+1] - ai[i];
53: for (j=0; j<anzj; j++) {
54: atj[atfill[*aj]] = i;
55: atfill[*aj++] += 1;
56: }
57: }
59: /* Clean up temporary space and complete requests. */
60: PetscFree(atfill);
61: *Ati = ati;
62: *Atj = atj;
64: PetscLogEventEnd(MAT_Getsymtranspose,A,0,0,0);
65: return(0);
66: }
67: /*
68: MatGetSymbolicTransposeReduced_SeqAIJ() - Get symbolic matrix structure of submatrix A[rstart:rend,:],
69: modified from MatGetSymbolicTranspose_SeqAIJ()
70: */
71: PetscErrorCode MatGetSymbolicTransposeReduced_SeqAIJ(Mat A,PetscInt rstart,PetscInt rend,PetscInt *Ati[],PetscInt *Atj[])
72: {
74: PetscInt i,j,anzj;
75: Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data;
76: PetscInt an=A->cmap->N;
77: PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j;
80: PetscInfo(A,"Getting Symbolic Transpose\n");
81: PetscLogEventBegin(MAT_Getsymtransreduced,A,0,0,0);
83: /* Allocate space for symbolic transpose info and work array */
84: PetscCalloc1(an+1,&ati);
85: anzj = ai[rend] - ai[rstart];
86: PetscMalloc1(anzj+1,&atj);
87: PetscMalloc1(an+1,&atfill);
89: /* Walk through aj and count ## of non-zeros in each row of A^T. */
90: /* Note: offset by 1 for fast conversion into csr format. */
91: for (i=ai[rstart]; i<ai[rend]; i++) {
92: ati[aj[i]+1] += 1;
93: }
94: /* Form ati for csr format of A^T. */
95: for (i=0;i<an;i++) {
96: ati[i+1] += ati[i];
97: }
99: /* Copy ati into atfill so we have locations of the next free space in atj */
100: PetscArraycpy(atfill,ati,an);
102: /* Walk through A row-wise and mark nonzero entries of A^T. */
103: aj = aj + ai[rstart];
104: for (i=rstart; i<rend; i++) {
105: anzj = ai[i+1] - ai[i];
106: for (j=0; j<anzj; j++) {
107: atj[atfill[*aj]] = i-rstart;
108: atfill[*aj++] += 1;
109: }
110: }
112: /* Clean up temporary space and complete requests. */
113: PetscFree(atfill);
114: *Ati = ati;
115: *Atj = atj;
117: PetscLogEventEnd(MAT_Getsymtransreduced,A,0,0,0);
118: return(0);
119: }
121: PetscErrorCode MatTranspose_SeqAIJ(Mat A,MatReuse reuse,Mat *B)
122: {
124: PetscInt i,j,anzj;
125: Mat At;
126: Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*at;
127: PetscInt an=A->cmap->N,am=A->rmap->N;
128: PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j;
129: MatScalar *ata,*aa=a->a;
132: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) {
133: /* Allocate space for symbolic transpose info and work array */
134: PetscCalloc1(an+1,&ati);
135: PetscMalloc1(ai[am],&atj);
136: PetscMalloc1(ai[am],&ata);
137: /* Walk through aj and count ## of non-zeros in each row of A^T. */
138: /* Note: offset by 1 for fast conversion into csr format. */
139: for (i=0;i<ai[am];i++) {
140: ati[aj[i]+1] += 1; /* count ## of non-zeros for row aj[i] of A^T */
141: }
142: /* Form ati for csr format of A^T. */
143: for (i=0;i<an;i++) {
144: ati[i+1] += ati[i];
145: }
146: } else { /* This segment is called by MatTranspose_MPIAIJ(...,MAT_INITIAL_MATRIX,..) directly! */
147: Mat_SeqAIJ *sub_B = (Mat_SeqAIJ*) (*B)->data;
148: ati = sub_B->i;
149: atj = sub_B->j;
150: ata = sub_B->a;
151: At = *B;
152: }
154: /* Copy ati into atfill so we have locations of the next free space in atj */
155: PetscMalloc1(an,&atfill);
156: PetscArraycpy(atfill,ati,an);
158: /* Walk through A row-wise and mark nonzero entries of A^T. */
159: for (i=0;i<am;i++) {
160: anzj = ai[i+1] - ai[i];
161: for (j=0;j<anzj;j++) {
162: atj[atfill[*aj]] = i;
163: ata[atfill[*aj]] = *aa++;
164: atfill[*aj++] += 1;
165: }
166: }
168: /* Clean up temporary space and complete requests. */
169: PetscFree(atfill);
170: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) {
171: MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),an,am,ati,atj,ata,&At);
172: MatSetBlockSizes(At,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));
174: at = (Mat_SeqAIJ*)(At->data);
175: at->free_a = PETSC_TRUE;
176: at->free_ij = PETSC_TRUE;
177: at->nonew = 0;
178: at->maxnz = ati[an];
180: MatSetType(At,((PetscObject)A)->type_name);
181: }
183: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) {
184: *B = At;
185: } else {
186: MatHeaderMerge(A,&At);
187: }
188: return(0);
189: }
191: PetscErrorCode MatRestoreSymbolicTranspose_SeqAIJ(Mat A,PetscInt *ati[],PetscInt *atj[])
192: {
196: PetscInfo(A,"Restoring Symbolic Transpose.\n");
197: PetscFree(*ati);
198: PetscFree(*atj);
199: return(0);
200: }