Actual source code: symtranspose.c

      1: #define PETSCMAT_DLL

      3: /*
      4:   Defines symbolic transpose routines for SeqAIJ matrices.

      6:   Currently Get/Restore only allocates/frees memory for holding the
      7:   (i,j) info for the transpose.  Someday, this info could be
      8:   maintained so successive calls to Get will not recompute the info.

     10:   Also defined is a "faster" implementation of MatTranspose for SeqAIJ
     11:   matrices which avoids calls to MatSetValues.  This routine has not
     12:   been adopted as the standard yet as it is somewhat untested.

     14: */

     16:  #include src/mat/impls/aij/seq/aij.h

     18: static PetscEvent logkey_matgetsymtranspose    = 0;
     19: static PetscEvent logkey_mattranspose          = 0;

     23: PetscErrorCode MatGetSymbolicTranspose_SeqAIJ(Mat A,PetscInt *Ati[],PetscInt *Atj[])
     24: {
     26:   PetscInt       i,j,anzj;
     27:   Mat_SeqAIJ     *a=(Mat_SeqAIJ *)A->data;
     28:   PetscInt       an=A->N,am=A->M;
     29:   PetscInt       *ati,*atj,*atfill,*ai=a->i,*aj=a->j;

     33:   PetscLogInfo((A,"MatGetSymbolicTranspose_SeqAIJ:Getting Symbolic Transpose.\n"));

     35:   /* Set up timers */
     36:   if (!logkey_matgetsymtranspose) {
     37:     PetscLogEventRegister(&logkey_matgetsymtranspose,"MatGetSymbolicTranspose",MAT_COOKIE);
     38:   }
     39:   PetscLogEventBegin(logkey_matgetsymtranspose,A,0,0,0);

     41:   /* Allocate space for symbolic transpose info and work array */
     42:   PetscMalloc((an+1)*sizeof(PetscInt),&ati);
     43:   PetscMalloc(ai[am]*sizeof(PetscInt),&atj);
     44:   PetscMalloc(an*sizeof(PetscInt),&atfill);
     45:   PetscMemzero(ati,(an+1)*sizeof(PetscInt));

     47:   /* Walk through aj and count ## of non-zeros in each row of A^T. */
     48:   /* Note: offset by 1 for fast conversion into csr format. */
     49:   for (i=0;i<ai[am];i++) {
     50:     ati[aj[i]+1] += 1;
     51:   }
     52:   /* Form ati for csr format of A^T. */
     53:   for (i=0;i<an;i++) {
     54:     ati[i+1] += ati[i];
     55:   }

     57:   /* Copy ati into atfill so we have locations of the next free space in atj */
     58:   PetscMemcpy(atfill,ati,an*sizeof(PetscInt));

     60:   /* Walk through A row-wise and mark nonzero entries of A^T. */
     61:   for (i=0;i<am;i++) {
     62:     anzj = ai[i+1] - ai[i];
     63:     for (j=0;j<anzj;j++) {
     64:       atj[atfill[*aj]] = i;
     65:       atfill[*aj++]   += 1;
     66:     }
     67:   }

     69:   /* Clean up temporary space and complete requests. */
     70:   PetscFree(atfill);
     71:   *Ati = ati;
     72:   *Atj = atj;

     74:   PetscLogEventEnd(logkey_matgetsymtranspose,A,0,0,0);
     75:   return(0);
     76: }
     77: /*
     78:   MatGetSymbolicTransposeReduced_SeqAIJ() - Get symbolic matrix structure of submatrix A[rstart:rend,:],
     79:      modified from MatGetSymbolicTranspose_SeqAIJ()
     80: */
     81: static PetscEvent logkey_matgetsymtransreduced = 0;
     84: PetscErrorCode MatGetSymbolicTransposeReduced_SeqAIJ(Mat A,PetscInt rstart,PetscInt rend,PetscInt *Ati[],PetscInt *Atj[])
     85: {
     87:   PetscInt       i,j,anzj;
     88:   Mat_SeqAIJ     *a=(Mat_SeqAIJ *)A->data;
     89:   PetscInt       an=A->N;
     90:   PetscInt       *ati,*atj,*atfill,*ai=a->i,*aj=a->j;

     93:   PetscLogInfo((A,"MatGetSymbolicTransposeReduced_SeqAIJ:Getting Symbolic Transpose.\n"));

     95:   /* Set up timers */
     96:   if (!logkey_matgetsymtransreduced) {
     97:     PetscLogEventRegister(&logkey_matgetsymtransreduced,"MatGetSymbolicTransposeReduced",MAT_COOKIE);
     98:   }
     99:   PetscLogEventBegin(logkey_matgetsymtransreduced,A,0,0,0);

    101:   /* Allocate space for symbolic transpose info and work array */
    102:   PetscMalloc((an+1)*sizeof(PetscInt),&ati);
    103:   anzj = ai[rend] - ai[rstart];
    104:   PetscMalloc((anzj+1)*sizeof(PetscInt),&atj);
    105:   PetscMalloc((an+1)*sizeof(PetscInt),&atfill);
    106:   PetscMemzero(ati,(an+1)*sizeof(PetscInt));

    108:   /* Walk through aj and count ## of non-zeros in each row of A^T. */
    109:   /* Note: offset by 1 for fast conversion into csr format. */
    110:   for (i=ai[rstart]; i<ai[rend]; i++) {
    111:     ati[aj[i]+1] += 1;
    112:   }
    113:   /* Form ati for csr format of A^T. */
    114:   for (i=0;i<an;i++) {
    115:     ati[i+1] += ati[i];
    116:   }

    118:   /* Copy ati into atfill so we have locations of the next free space in atj */
    119:   PetscMemcpy(atfill,ati,an*sizeof(PetscInt));

    121:   /* Walk through A row-wise and mark nonzero entries of A^T. */
    122:   aj = aj + ai[rstart];
    123:   for (i=rstart; i<rend; i++) {
    124:     anzj = ai[i+1] - ai[i];
    125:     for (j=0;j<anzj;j++) {
    126:       atj[atfill[*aj]] = i-rstart;
    127:       atfill[*aj++]   += 1;
    128:     }
    129:   }

    131:   /* Clean up temporary space and complete requests. */
    132:   PetscFree(atfill);
    133:   *Ati = ati;
    134:   *Atj = atj;

    136:   PetscLogEventEnd(logkey_matgetsymtransreduced,A,0,0,0);
    137:   return(0);
    138: }

    142: PetscErrorCode MatTranspose_SeqAIJ_FAST(Mat A,Mat *B)
    143: {
    145:   PetscInt       i,j,anzj;
    146:   Mat            At;
    147:   Mat_SeqAIJ     *a=(Mat_SeqAIJ *)A->data,*at;
    148:   PetscInt       an=A->N,am=A->M;
    149:   PetscInt       *ati,*atj,*atfill,*ai=a->i,*aj=a->j;
    150:   MatScalar      *ata,*aa=a->a;

    154:   /* Set up timers */
    155:   if (!logkey_mattranspose) {
    156:     PetscLogEventRegister(&logkey_mattranspose,"MatTranspose_SeqAIJ_FAST",MAT_COOKIE);
    157:   }
    158:   PetscLogEventBegin(logkey_mattranspose,A,0,0,0);

    160:   /* Allocate space for symbolic transpose info and work array */
    161:   PetscMalloc((an+1)*sizeof(PetscInt),&ati);
    162:   PetscMalloc(ai[am]*sizeof(PetscInt),&atj);
    163:   PetscMalloc(ai[am]*sizeof(MatScalar),&ata);
    164:   PetscMalloc(an*sizeof(PetscInt),&atfill);
    165:   PetscMemzero(ati,(an+1)*sizeof(PetscInt));
    166:   /* Walk through aj and count ## of non-zeros in each row of A^T. */
    167:   /* Note: offset by 1 for fast conversion into csr format. */
    168:   for (i=0;i<ai[am];i++) {
    169:     ati[aj[i]+1] += 1;
    170:   }
    171:   /* Form ati for csr format of A^T. */
    172:   for (i=0;i<an;i++) {
    173:     ati[i+1] += ati[i];
    174:   }

    176:   /* Copy ati into atfill so we have locations of the next free space in atj */
    177:   PetscMemcpy(atfill,ati,an*sizeof(PetscInt));

    179:   /* Walk through A row-wise and mark nonzero entries of A^T. */
    180:   for (i=0;i<am;i++) {
    181:     anzj = ai[i+1] - ai[i];
    182:     for (j=0;j<anzj;j++) {
    183:       atj[atfill[*aj]] = i;
    184:       ata[atfill[*aj]] = *aa++;
    185:       atfill[*aj++]   += 1;
    186:     }
    187:   }

    189:   /* Clean up temporary space and complete requests. */
    190:   PetscFree(atfill);
    191:   MatCreateSeqAIJWithArrays(A->comm,an,am,ati,atj,ata,&At);
    192:   at   = (Mat_SeqAIJ *)(At->data);
    193:   at->freedata = PETSC_TRUE;
    194:   at->nonew    = 0;
    195:   if (B) {
    196:     *B = At;
    197:   } else {
    198:     MatHeaderCopy(A,At);
    199:   }
    200:   PetscLogEventEnd(logkey_mattranspose,A,0,0,0);
    201:   return(0);
    202: }

    206: PetscErrorCode MatRestoreSymbolicTranspose_SeqAIJ(Mat A,PetscInt *ati[],PetscInt *atj[])
    207: {

    211:   PetscLogInfo((A,"MatRestoreSymbolicTranspose_SeqAIJ:Restoring Symbolic Transpose.\n"));
    212:   PetscFree(*ati);
    213:   ati  = PETSC_NULL;
    214:   PetscFree(*atj);
    215:   atj  = PETSC_NULL;
    216:   return(0);
    217: }