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: }