Actual source code: symtranspose.c
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>
17: PetscErrorCode MatGetSymbolicTranspose_SeqAIJ(Mat A,PetscInt *Ati[],PetscInt *Atj[])
18: {
19: PetscInt i,j,anzj;
20: Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data;
21: PetscInt an=A->cmap->N,am=A->rmap->N;
22: PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j;
24: PetscInfo(A,"Getting Symbolic Transpose.\n");
26: /* Set up timers */
27: PetscLogEventBegin(MAT_Getsymtranspose,A,0,0,0);
29: /* Allocate space for symbolic transpose info and work array */
30: PetscCalloc1(an+1,&ati);
31: PetscMalloc1(ai[am],&atj);
32: PetscMalloc1(an,&atfill);
34: /* Walk through aj and count ## of non-zeros in each row of A^T. */
35: /* Note: offset by 1 for fast conversion into csr format. */
36: for (i=0;i<ai[am];i++) {
37: ati[aj[i]+1] += 1;
38: }
39: /* Form ati for csr format of A^T. */
40: for (i=0;i<an;i++) {
41: ati[i+1] += ati[i];
42: }
44: /* Copy ati into atfill so we have locations of the next free space in atj */
45: PetscArraycpy(atfill,ati,an);
47: /* Walk through A row-wise and mark nonzero entries of A^T. */
48: for (i=0; i<am; i++) {
49: anzj = ai[i+1] - ai[i];
50: for (j=0; j<anzj; j++) {
51: atj[atfill[*aj]] = i;
52: atfill[*aj++] += 1;
53: }
54: }
56: /* Clean up temporary space and complete requests. */
57: PetscFree(atfill);
58: *Ati = ati;
59: *Atj = atj;
61: PetscLogEventEnd(MAT_Getsymtranspose,A,0,0,0);
62: return 0;
63: }
64: /*
65: MatGetSymbolicTransposeReduced_SeqAIJ() - Get symbolic matrix structure of submatrix A[rstart:rend,:],
66: modified from MatGetSymbolicTranspose_SeqAIJ()
67: */
68: PetscErrorCode MatGetSymbolicTransposeReduced_SeqAIJ(Mat A,PetscInt rstart,PetscInt rend,PetscInt *Ati[],PetscInt *Atj[])
69: {
70: PetscInt i,j,anzj;
71: Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data;
72: PetscInt an=A->cmap->N;
73: PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j;
75: PetscInfo(A,"Getting Symbolic Transpose\n");
76: PetscLogEventBegin(MAT_Getsymtransreduced,A,0,0,0);
78: /* Allocate space for symbolic transpose info and work array */
79: PetscCalloc1(an+1,&ati);
80: anzj = ai[rend] - ai[rstart];
81: PetscMalloc1(anzj+1,&atj);
82: PetscMalloc1(an+1,&atfill);
84: /* Walk through aj and count ## of non-zeros in each row of A^T. */
85: /* Note: offset by 1 for fast conversion into csr format. */
86: for (i=ai[rstart]; i<ai[rend]; i++) {
87: ati[aj[i]+1] += 1;
88: }
89: /* Form ati for csr format of A^T. */
90: for (i=0;i<an;i++) {
91: ati[i+1] += ati[i];
92: }
94: /* Copy ati into atfill so we have locations of the next free space in atj */
95: PetscArraycpy(atfill,ati,an);
97: /* Walk through A row-wise and mark nonzero entries of A^T. */
98: aj = aj + ai[rstart];
99: for (i=rstart; i<rend; i++) {
100: anzj = ai[i+1] - ai[i];
101: for (j=0; j<anzj; j++) {
102: atj[atfill[*aj]] = i-rstart;
103: atfill[*aj++] += 1;
104: }
105: }
107: /* Clean up temporary space and complete requests. */
108: PetscFree(atfill);
109: *Ati = ati;
110: *Atj = atj;
112: PetscLogEventEnd(MAT_Getsymtransreduced,A,0,0,0);
113: return 0;
114: }
116: PetscErrorCode MatTranspose_SeqAIJ(Mat A,MatReuse reuse,Mat *B)
117: {
118: PetscInt i,j,anzj;
119: Mat At;
120: Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*at;
121: PetscInt an=A->cmap->N,am=A->rmap->N;
122: PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j;
123: MatScalar *ata;
124: const MatScalar *aa,*av;
126: MatSeqAIJGetArrayRead(A,&av);
127: aa = av;
128: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) {
129: /* Allocate space for symbolic transpose info and work array */
130: PetscCalloc1(an+1,&ati);
131: PetscMalloc1(ai[am],&atj);
132: PetscMalloc1(ai[am],&ata);
133: /* Walk through aj and count ## of non-zeros in each row of A^T. */
134: /* Note: offset by 1 for fast conversion into csr format. */
135: for (i=0;i<ai[am];i++) {
136: ati[aj[i]+1] += 1; /* count ## of non-zeros for row aj[i] of A^T */
137: }
138: /* Form ati for csr format of A^T. */
139: for (i=0;i<an;i++) {
140: ati[i+1] += ati[i];
141: }
142: } else { /* This segment is called by MatTranspose_MPIAIJ(...,MAT_INITIAL_MATRIX,..) directly! */
143: Mat_SeqAIJ *sub_B = (Mat_SeqAIJ*) (*B)->data;
144: ati = sub_B->i;
145: atj = sub_B->j;
146: ata = sub_B->a;
147: At = *B;
148: }
150: /* Copy ati into atfill so we have locations of the next free space in atj */
151: PetscMalloc1(an,&atfill);
152: PetscArraycpy(atfill,ati,an);
154: /* Walk through A row-wise and mark nonzero entries of A^T. */
155: for (i=0;i<am;i++) {
156: anzj = ai[i+1] - ai[i];
157: for (j=0;j<anzj;j++) {
158: atj[atfill[*aj]] = i;
159: ata[atfill[*aj]] = *aa++;
160: atfill[*aj++] += 1;
161: }
162: }
163: MatSeqAIJRestoreArrayRead(A,&av);
165: /* Clean up temporary space and complete requests. */
166: PetscFree(atfill);
167: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) {
168: MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),an,am,ati,atj,ata,&At);
169: MatSetBlockSizes(At,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));
171: at = (Mat_SeqAIJ*)(At->data);
172: at->free_a = PETSC_TRUE;
173: at->free_ij = PETSC_TRUE;
174: at->nonew = 0;
175: at->maxnz = ati[an];
177: MatSetType(At,((PetscObject)A)->type_name);
178: }
180: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) {
181: *B = At;
182: } else {
183: MatHeaderMerge(A,&At);
184: }
185: return 0;
186: }
188: PetscErrorCode MatRestoreSymbolicTranspose_SeqAIJ(Mat A,PetscInt *ati[],PetscInt *atj[])
189: {
190: PetscInfo(A,"Restoring Symbolic Transpose.\n");
191: PetscFree(*ati);
192: PetscFree(*atj);
193: return 0;
194: }