1 /* 2 Defines symbolic transpose routines for SeqAIJ matrices. 3 4 Currently Get/Restore only allocates/frees memory for holding the 5 (i,j) info for the transpose. Someday, this info could be 6 maintained so successive calls to Get will not recompute the info. 7 8 Also defined is a "faster" implementation of MatTranspose for SeqAIJ 9 matrices which avoids calls to MatSetValues. This routine has not 10 been adopted as the standard yet as it is somewhat untested. 11 12 */ 13 14 #include "src/mat/impls/aij/seq/aij.h" 15 16 static PetscEvent logkey_matgetsymtranspose = 0; 17 static PetscEvent logkey_mattranspose = 0; 18 19 20 #undef __FUNCT__ 21 #define __FUNCT__ "MatGetSymbolicTranspose_SeqIJ" 22 PetscErrorCode MatGetSymbolicTranspose_SeqAIJ(Mat A,PetscInt *Ati[],PetscInt *Atj[]) 23 { 24 PetscErrorCode ierr; 25 PetscInt i,j,anzj; 26 Mat_SeqAIJ *a=(Mat_SeqAIJ *)A->data; 27 PetscInt an=A->N,am=A->M; 28 PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j; 29 30 PetscFunctionBegin; 31 32 ierr = PetscLogInfo(A,"Getting Symbolic Transpose.\n");CHKERRQ(ierr); 33 34 /* Set up timers */ 35 if (!logkey_matgetsymtranspose) { 36 ierr = PetscLogEventRegister(&logkey_matgetsymtranspose,"MatGetSymbolicTranspose",MAT_COOKIE);CHKERRQ(ierr); 37 } 38 ierr = PetscLogEventBegin(logkey_matgetsymtranspose,A,0,0,0);CHKERRQ(ierr); 39 40 /* Allocate space for symbolic transpose info and work array */ 41 ierr = PetscMalloc((an+1)*sizeof(PetscInt),&ati);CHKERRQ(ierr); 42 ierr = PetscMalloc(ai[am]*sizeof(PetscInt),&atj);CHKERRQ(ierr); 43 ierr = PetscMalloc(an*sizeof(PetscInt),&atfill);CHKERRQ(ierr); 44 ierr = PetscMemzero(ati,(an+1)*sizeof(PetscInt));CHKERRQ(ierr); 45 46 /* Walk through aj and count ## of non-zeros in each row of A^T. */ 47 /* Note: offset by 1 for fast conversion into csr format. */ 48 for (i=0;i<ai[am];i++) { 49 ati[aj[i]+1] += 1; 50 } 51 /* Form ati for csr format of A^T. */ 52 for (i=0;i<an;i++) { 53 ati[i+1] += ati[i]; 54 } 55 56 /* Copy ati into atfill so we have locations of the next free space in atj */ 57 ierr = PetscMemcpy(atfill,ati,an*sizeof(PetscInt));CHKERRQ(ierr); 58 59 /* Walk through A row-wise and mark nonzero entries of A^T. */ 60 for (i=0;i<am;i++) { 61 anzj = ai[i+1] - ai[i]; 62 for (j=0;j<anzj;j++) { 63 atj[atfill[*aj]] = i; 64 atfill[*aj++] += 1; 65 } 66 } 67 68 /* Clean up temporary space and complete requests. */ 69 ierr = PetscFree(atfill);CHKERRQ(ierr); 70 *Ati = ati; 71 *Atj = atj; 72 73 ierr = PetscLogEventEnd(logkey_matgetsymtranspose,A,0,0,0);CHKERRQ(ierr); 74 PetscFunctionReturn(0); 75 } 76 77 #undef __FUNCT__ 78 #define __FUNCT__ "MatTranspose_SeqIJ_FAST" 79 PetscErrorCode MatTranspose_SeqAIJ_FAST(Mat A,Mat *B) 80 { 81 PetscErrorCode ierr; 82 PetscInt i,j,anzj; 83 Mat At; 84 Mat_SeqAIJ *a=(Mat_SeqAIJ *)A->data,*at; 85 PetscInt an=A->N,am=A->M; 86 PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j; 87 MatScalar *ata,*aa=a->a; 88 89 PetscFunctionBegin; 90 91 /* Set up timers */ 92 if (!logkey_mattranspose) { 93 ierr = PetscLogEventRegister(&logkey_mattranspose,"MatTranspose_SeqAIJ_FAST",MAT_COOKIE);CHKERRQ(ierr); 94 } 95 ierr = PetscLogEventBegin(logkey_mattranspose,A,0,0,0);CHKERRQ(ierr); 96 97 /* Allocate space for symbolic transpose info and work array */ 98 ierr = PetscMalloc((an+1)*sizeof(PetscInt),&ati);CHKERRQ(ierr); 99 ierr = PetscMalloc(ai[am]*sizeof(PetscInt),&atj);CHKERRQ(ierr); 100 ierr = PetscMalloc(ai[am]*sizeof(MatScalar),&ata);CHKERRQ(ierr); 101 ierr = PetscMalloc(an*sizeof(PetscInt),&atfill);CHKERRQ(ierr); 102 ierr = PetscMemzero(ati,(an+1)*sizeof(PetscInt));CHKERRQ(ierr); 103 /* Walk through aj and count ## of non-zeros in each row of A^T. */ 104 /* Note: offset by 1 for fast conversion into csr format. */ 105 for (i=0;i<ai[am];i++) { 106 ati[aj[i]+1] += 1; 107 } 108 /* Form ati for csr format of A^T. */ 109 for (i=0;i<an;i++) { 110 ati[i+1] += ati[i]; 111 } 112 113 /* Copy ati into atfill so we have locations of the next free space in atj */ 114 ierr = PetscMemcpy(atfill,ati,an*sizeof(PetscInt));CHKERRQ(ierr); 115 116 /* Walk through A row-wise and mark nonzero entries of A^T. */ 117 for (i=0;i<am;i++) { 118 anzj = ai[i+1] - ai[i]; 119 for (j=0;j<anzj;j++) { 120 atj[atfill[*aj]] = i; 121 ata[atfill[*aj]] = *aa++; 122 atfill[*aj++] += 1; 123 } 124 } 125 126 /* Clean up temporary space and complete requests. */ 127 ierr = PetscFree(atfill);CHKERRQ(ierr); 128 ierr = MatCreateSeqAIJWithArrays(A->comm,an,am,ati,atj,ata,&At);CHKERRQ(ierr); 129 at = (Mat_SeqAIJ *)(At->data); 130 at->freedata = PETSC_TRUE; 131 at->nonew = 0; 132 if (B) { 133 *B = At; 134 } else { 135 ierr = MatHeaderCopy(A,At); 136 } 137 ierr = PetscLogEventEnd(logkey_mattranspose,A,0,0,0);CHKERRQ(ierr); 138 PetscFunctionReturn(0); 139 } 140 141 #undef __FUNCT__ 142 #define __FUNCT__ "MatRestoreSymbolicTranspose_SeqAIJ" 143 PetscErrorCode MatRestoreSymbolicTranspose_SeqAIJ(Mat A,PetscInt *ati[],PetscInt *atj[]) 144 { 145 PetscErrorCode ierr; 146 147 PetscFunctionBegin; 148 ierr = PetscLogInfo(A,"Restoring Symbolic Transpose.\n");CHKERRQ(ierr); 149 ierr = PetscFree(*ati);CHKERRQ(ierr); 150 ati = PETSC_NULL; 151 ierr = PetscFree(*atj);CHKERRQ(ierr); 152 atj = PETSC_NULL; 153 PetscFunctionReturn(0); 154 } 155 156