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 int logkey_matgetsymtranspose = 0; 17 static int logkey_mattranspose = 0; 18 19 20 #undef __FUNCT__ 21 #define __FUNCT__ "MatGetSymbolicTranspose_SeqIJ" 22 int MatGetSymbolicTranspose_SeqAIJ(Mat A,int *Ati[],int *Atj[]) { 23 int ierr,i,j,anzj; 24 Mat_SeqAIJ *a=(Mat_SeqAIJ *)A->data; 25 int an=A->N,am=A->M; 26 int *ati,*atj,*atfill,*ai=a->i,*aj=a->j; 27 28 PetscFunctionBegin; 29 30 ierr = PetscLogInfo(A,"Getting Symbolic Transpose.\n");CHKERRQ(ierr); 31 32 /* Set up timers */ 33 if (!logkey_matgetsymtranspose) { 34 ierr = PetscLogEventRegister(&logkey_matgetsymtranspose,"MatGetSymbolicTranspose",MAT_COOKIE);CHKERRQ(ierr); 35 } 36 ierr = PetscLogEventBegin(logkey_matgetsymtranspose,A,0,0,0);CHKERRQ(ierr); 37 38 /* Allocate space for symbolic transpose info and work array */ 39 ierr = PetscMalloc((an+1)*sizeof(int),&ati);CHKERRQ(ierr); 40 ierr = PetscMalloc(ai[am]*sizeof(int),&atj);CHKERRQ(ierr); 41 ierr = PetscMalloc(an*sizeof(int),&atfill);CHKERRQ(ierr); 42 ierr = PetscMemzero(ati,(an+1)*sizeof(int));CHKERRQ(ierr); 43 44 /* Walk through aj and count ## of non-zeros in each row of A^T. */ 45 /* Note: offset by 1 for fast conversion into csr format. */ 46 for (i=0;i<ai[am];i++) { 47 ati[aj[i]+1] += 1; 48 } 49 /* Form ati for csr format of A^T. */ 50 for (i=0;i<an;i++) { 51 ati[i+1] += ati[i]; 52 } 53 54 /* Copy ati into atfill so we have locations of the next free space in atj */ 55 ierr = PetscMemcpy(atfill,ati,an*sizeof(int));CHKERRQ(ierr); 56 57 /* Walk through A row-wise and mark nonzero entries of A^T. */ 58 for (i=0;i<am;i++) { 59 anzj = ai[i+1] - ai[i]; 60 for (j=0;j<anzj;j++) { 61 atj[atfill[*aj]] = i; 62 atfill[*aj++] += 1; 63 } 64 } 65 66 /* Clean up temporary space and complete requests. */ 67 ierr = PetscFree(atfill);CHKERRQ(ierr); 68 *Ati = ati; 69 *Atj = atj; 70 71 ierr = PetscLogEventEnd(logkey_matgetsymtranspose,A,0,0,0);CHKERRQ(ierr); 72 PetscFunctionReturn(0); 73 } 74 75 #undef __FUNCT__ 76 #define __FUNCT__ "MatTranspose_SeqIJ_FAST" 77 int MatTranspose_SeqAIJ_FAST(Mat A,Mat *B) { 78 int ierr,i,j,anzj; 79 Mat At; 80 Mat_SeqAIJ *a=(Mat_SeqAIJ *)A->data,*at; 81 int an=A->N,am=A->M; 82 int *ati,*atj,*atfill,*ai=a->i,*aj=a->j; 83 MatScalar *ata,*aa=a->a; 84 PetscFunctionBegin; 85 86 /* Set up timers */ 87 if (!logkey_mattranspose) { 88 ierr = PetscLogEventRegister(&logkey_mattranspose,"MatTranspose_SeqAIJ_FAST",MAT_COOKIE);CHKERRQ(ierr); 89 } 90 ierr = PetscLogEventBegin(logkey_mattranspose,A,0,0,0);CHKERRQ(ierr); 91 92 /* Allocate space for symbolic transpose info and work array */ 93 ierr = PetscMalloc((an+1)*sizeof(int),&ati);CHKERRQ(ierr); 94 ierr = PetscMalloc(ai[am]*sizeof(int),&atj);CHKERRQ(ierr); 95 ierr = PetscMalloc(ai[am]*sizeof(MatScalar),&ata);CHKERRQ(ierr); 96 ierr = PetscMalloc(an*sizeof(int),&atfill);CHKERRQ(ierr); 97 ierr = PetscMemzero(ati,(an+1)*sizeof(int));CHKERRQ(ierr); 98 /* Walk through aj and count ## of non-zeros in each row of A^T. */ 99 /* Note: offset by 1 for fast conversion into csr format. */ 100 for (i=0;i<ai[am];i++) { 101 ati[aj[i]+1] += 1; 102 } 103 /* Form ati for csr format of A^T. */ 104 for (i=0;i<an;i++) { 105 ati[i+1] += ati[i]; 106 } 107 108 /* Copy ati into atfill so we have locations of the next free space in atj */ 109 ierr = PetscMemcpy(atfill,ati,an*sizeof(int));CHKERRQ(ierr); 110 111 /* Walk through A row-wise and mark nonzero entries of A^T. */ 112 for (i=0;i<am;i++) { 113 anzj = ai[i+1] - ai[i]; 114 for (j=0;j<anzj;j++) { 115 atj[atfill[*aj]] = i; 116 ata[atfill[*aj]] = *aa++; 117 atfill[*aj++] += 1; 118 } 119 } 120 121 /* Clean up temporary space and complete requests. */ 122 ierr = PetscFree(atfill);CHKERRQ(ierr); 123 ierr = MatCreateSeqAIJWithArrays(A->comm,an,am,ati,atj,ata,&At);CHKERRQ(ierr); 124 at = (Mat_SeqAIJ *)(At->data); 125 at->freedata = PETSC_TRUE; 126 at->nonew = 0; 127 if (B) { 128 *B = At; 129 } else { 130 ierr = MatHeaderCopy(A,At); 131 } 132 ierr = PetscLogEventEnd(logkey_mattranspose,A,0,0,0);CHKERRQ(ierr); 133 PetscFunctionReturn(0); 134 } 135 136 #undef __FUNCT__ 137 #define __FUNCT__ "MatRestoreSymbolicTranspose_SeqAIJ" 138 int MatRestoreSymbolicTranspose_SeqAIJ(Mat A,int *ati[],int *atj[]) { 139 int ierr; 140 141 PetscFunctionBegin; 142 ierr = PetscLogInfo(A,"Restoring Symbolic Transpose.\n");CHKERRQ(ierr); 143 ierr = PetscFree(*ati);CHKERRQ(ierr); 144 ati = PETSC_NULL; 145 ierr = PetscFree(*atj);CHKERRQ(ierr); 146 atj = PETSC_NULL; 147 PetscFunctionReturn(0); 148 } 149 150