Actual source code: symtranspose.c

petsc-master 2019-11-16
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 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: }