Actual source code: symtranspose.c

petsc-3.9.4 2018-09-11
Report Typos and Errors

  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>


 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:   PetscMemcpy(atfill,ati,an*sizeof(PetscInt));

 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:   PetscMemcpy(atfill,ati,an*sizeof(PetscInt));

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_FAST(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:   PetscLogEventBegin(MAT_Transpose_SeqAIJ,A,0,0,0);

134:   if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) {
135:     /* Allocate space for symbolic transpose info and work array */
136:     PetscCalloc1(an+1,&ati);
137:     PetscMalloc1(ai[am],&atj);
138:     PetscMalloc1(ai[am],&ata);
139:     /* Walk through aj and count ## of non-zeros in each row of A^T. */
140:     /* Note: offset by 1 for fast conversion into csr format. */
141:     for (i=0;i<ai[am];i++) {
142:       ati[aj[i]+1] += 1;
143:     }
144:     /* Form ati for csr format of A^T. */
145:     for (i=0;i<an;i++) {
146:       ati[i+1] += ati[i];
147:     }
148:   } else {
149:     Mat_SeqAIJ *sub_B = (Mat_SeqAIJ*) (*B)->data;
150:     ati = sub_B->i;
151:     atj = sub_B->j;
152:     ata = sub_B->a;
153:     At  = *B;
154:   }

156:   /* Copy ati into atfill so we have locations of the next free space in atj */
157:   PetscMalloc1(an,&atfill);
158:   PetscMemcpy(atfill,ati,an*sizeof(PetscInt));

160:   /* Walk through A row-wise and mark nonzero entries of A^T. */
161:   for (i=0;i<am;i++) {
162:     anzj = ai[i+1] - ai[i];
163:     for (j=0;j<anzj;j++) {
164:       atj[atfill[*aj]] = i;
165:       ata[atfill[*aj]] = *aa++;
166:       atfill[*aj++]   += 1;
167:     }
168:   }

170:   /* Clean up temporary space and complete requests. */
171:   PetscFree(atfill);
172:   if (reuse == MAT_INITIAL_MATRIX) {
173:     MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),an,am,ati,atj,ata,&At);

175:     at          = (Mat_SeqAIJ*)(At->data);
176:     at->free_a  = PETSC_TRUE;
177:     at->free_ij = PETSC_TRUE;
178:     at->nonew   = 0;
179:   }

181:   if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) {
182:     *B = At;
183:   } else {
184:     MatHeaderMerge(A,&At);
185:   }
186:   PetscLogEventEnd(MAT_Transpose_SeqAIJ,A,0,0,0);
187:   return(0);
188: }

190: PetscErrorCode MatRestoreSymbolicTranspose_SeqAIJ(Mat A,PetscInt *ati[],PetscInt *atj[])
191: {

195:   PetscInfo(A,"Restoring Symbolic Transpose.\n");
196:   PetscFree(*ati);
197:   PetscFree(*atj);
198:   return(0);
199: }