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,int *Ati[],int *Atj[]) 23 { 24 PetscErrorCode ierr; 25 int i,j,anzj; 26 Mat_SeqAIJ *a=(Mat_SeqAIJ *)A->data; 27 int an=A->N,am=A->M; 28 int *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(int),&ati);CHKERRQ(ierr); 42 ierr = PetscMalloc(ai[am]*sizeof(int),&atj);CHKERRQ(ierr); 43 ierr = PetscMalloc(an*sizeof(int),&atfill);CHKERRQ(ierr); 44 ierr = PetscMemzero(ati,(an+1)*sizeof(int));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(int));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 int i,j,anzj; 83 Mat At; 84 Mat_SeqAIJ *a=(Mat_SeqAIJ *)A->data,*at; 85 int an=A->N,am=A->M; 86 int *ati,*atj,*atfill,*ai=a->i,*aj=a->j; 87 MatScalar *ata,*aa=a->a; 88 PetscFunctionBegin; 89 90 /* Set up timers */ 91 if (!logkey_mattranspose) { 92 ierr = PetscLogEventRegister(&logkey_mattranspose,"MatTranspose_SeqAIJ_FAST",MAT_COOKIE);CHKERRQ(ierr); 93 } 94 ierr = PetscLogEventBegin(logkey_mattranspose,A,0,0,0);CHKERRQ(ierr); 95 96 /* Allocate space for symbolic transpose info and work array */ 97 ierr = PetscMalloc((an+1)*sizeof(int),&ati);CHKERRQ(ierr); 98 ierr = PetscMalloc(ai[am]*sizeof(int),&atj);CHKERRQ(ierr); 99 ierr = PetscMalloc(ai[am]*sizeof(MatScalar),&ata);CHKERRQ(ierr); 100 ierr = PetscMalloc(an*sizeof(int),&atfill);CHKERRQ(ierr); 101 ierr = PetscMemzero(ati,(an+1)*sizeof(int));CHKERRQ(ierr); 102 /* Walk through aj and count ## of non-zeros in each row of A^T. */ 103 /* Note: offset by 1 for fast conversion into csr format. */ 104 for (i=0;i<ai[am];i++) { 105 ati[aj[i]+1] += 1; 106 } 107 /* Form ati for csr format of A^T. */ 108 for (i=0;i<an;i++) { 109 ati[i+1] += ati[i]; 110 } 111 112 /* Copy ati into atfill so we have locations of the next free space in atj */ 113 ierr = PetscMemcpy(atfill,ati,an*sizeof(int));CHKERRQ(ierr); 114 115 /* Walk through A row-wise and mark nonzero entries of A^T. */ 116 for (i=0;i<am;i++) { 117 anzj = ai[i+1] - ai[i]; 118 for (j=0;j<anzj;j++) { 119 atj[atfill[*aj]] = i; 120 ata[atfill[*aj]] = *aa++; 121 atfill[*aj++] += 1; 122 } 123 } 124 125 /* Clean up temporary space and complete requests. */ 126 ierr = PetscFree(atfill);CHKERRQ(ierr); 127 ierr = MatCreateSeqAIJWithArrays(A->comm,an,am,ati,atj,ata,&At);CHKERRQ(ierr); 128 at = (Mat_SeqAIJ *)(At->data); 129 at->freedata = PETSC_TRUE; 130 at->nonew = 0; 131 if (B) { 132 *B = At; 133 } else { 134 ierr = MatHeaderCopy(A,At); 135 } 136 ierr = PetscLogEventEnd(logkey_mattranspose,A,0,0,0);CHKERRQ(ierr); 137 PetscFunctionReturn(0); 138 } 139 140 #undef __FUNCT__ 141 #define __FUNCT__ "MatRestoreSymbolicTranspose_SeqAIJ" 142 PetscErrorCode MatRestoreSymbolicTranspose_SeqAIJ(Mat A,int *ati[],int *atj[]) { 143 PetscErrorCode ierr; 144 145 PetscFunctionBegin; 146 ierr = PetscLogInfo(A,"Restoring Symbolic Transpose.\n");CHKERRQ(ierr); 147 ierr = PetscFree(*ati);CHKERRQ(ierr); 148 ati = PETSC_NULL; 149 ierr = PetscFree(*atj);CHKERRQ(ierr); 150 atj = PETSC_NULL; 151 PetscFunctionReturn(0); 152 } 153 154