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