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 #undef __FUNCT__ 20 #define __FUNCT__ "MatGetSymbolicTranspose_SeqIJ" 21 PetscErrorCode MatGetSymbolicTranspose_SeqAIJ(Mat A,PetscInt *Ati[],PetscInt *Atj[]) 22 { 23 PetscErrorCode ierr; 24 PetscInt i,j,anzj; 25 Mat_SeqAIJ *a=(Mat_SeqAIJ *)A->data; 26 PetscInt an=A->N,am=A->M; 27 PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j; 28 29 PetscFunctionBegin; 30 31 ierr = PetscLogInfo(A,"Getting Symbolic Transpose.\n");CHKERRQ(ierr); 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(PetscInt),&ati);CHKERRQ(ierr); 41 ierr = PetscMalloc(ai[am]*sizeof(PetscInt),&atj);CHKERRQ(ierr); 42 ierr = PetscMalloc(an*sizeof(PetscInt),&atfill);CHKERRQ(ierr); 43 ierr = PetscMemzero(ati,(an+1)*sizeof(PetscInt));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(PetscInt));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 MatGetSymbolicTransposeReduced_SeqAIJ() - Get symbolic matrix structure of submatrix A[rstart:rend,:], 77 modified from MatGetSymbolicTranspose_SeqAIJ() 78 */ 79 static PetscEvent logkey_matgetsymtransreduced = 0; 80 #undef __FUNCT__ 81 #define __FUNCT__ "MatGetSymbolicTransposeReduced_SeqIJ" 82 PetscErrorCode MatGetSymbolicTransposeReduced_SeqAIJ(Mat A,PetscInt rstart,PetscInt rend,PetscInt *Ati[],PetscInt *Atj[]) 83 { 84 PetscErrorCode ierr; 85 PetscInt i,j,anzj; 86 Mat_SeqAIJ *a=(Mat_SeqAIJ *)A->data; 87 PetscInt an=A->N; 88 PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j; 89 90 PetscFunctionBegin; 91 92 ierr = PetscLogInfo(A,"Getting Symbolic Transpose.\n");CHKERRQ(ierr); 93 94 /* Set up timers */ 95 if (!logkey_matgetsymtransreduced) { 96 ierr = PetscLogEventRegister(&logkey_matgetsymtransreduced,"MatGetSymbolicTransposeReduced",MAT_COOKIE);CHKERRQ(ierr); 97 } 98 ierr = PetscLogEventBegin(logkey_matgetsymtransreduced,A,0,0,0);CHKERRQ(ierr); 99 100 /* Allocate space for symbolic transpose info and work array */ 101 ierr = PetscMalloc((an+1)*sizeof(PetscInt),&ati);CHKERRQ(ierr); 102 anzj = ai[rend] - ai[rstart]; 103 ierr = PetscMalloc((anzj+1)*sizeof(PetscInt),&atj);CHKERRQ(ierr); 104 ierr = PetscMalloc((an+1)*sizeof(PetscInt),&atfill);CHKERRQ(ierr); 105 ierr = PetscMemzero(ati,(an+1)*sizeof(PetscInt));CHKERRQ(ierr); 106 107 /* Walk through aj and count ## of non-zeros in each row of A^T. */ 108 /* Note: offset by 1 for fast conversion into csr format. */ 109 for (i=ai[rstart]; i<ai[rend]; i++) { 110 ati[aj[i]+1] += 1; 111 } 112 /* Form ati for csr format of A^T. */ 113 for (i=0;i<an;i++) { 114 ati[i+1] += ati[i]; 115 } 116 117 /* Copy ati into atfill so we have locations of the next free space in atj */ 118 ierr = PetscMemcpy(atfill,ati,an*sizeof(PetscInt));CHKERRQ(ierr); 119 120 /* Walk through A row-wise and mark nonzero entries of A^T. */ 121 aj = aj + ai[rstart]; 122 for (i=rstart; i<rend; i++) { 123 anzj = ai[i+1] - ai[i]; 124 for (j=0;j<anzj;j++) { 125 atj[atfill[*aj]] = i-rstart; 126 atfill[*aj++] += 1; 127 } 128 } 129 130 /* Clean up temporary space and complete requests. */ 131 ierr = PetscFree(atfill);CHKERRQ(ierr); 132 *Ati = ati; 133 *Atj = atj; 134 135 ierr = PetscLogEventEnd(logkey_matgetsymtransreduced,A,0,0,0);CHKERRQ(ierr); 136 PetscFunctionReturn(0); 137 } 138 139 #undef __FUNCT__ 140 #define __FUNCT__ "MatTranspose_SeqIJ_FAST" 141 PetscErrorCode MatTranspose_SeqAIJ_FAST(Mat A,Mat *B) 142 { 143 PetscErrorCode ierr; 144 PetscInt i,j,anzj; 145 Mat At; 146 Mat_SeqAIJ *a=(Mat_SeqAIJ *)A->data,*at; 147 PetscInt an=A->N,am=A->M; 148 PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j; 149 MatScalar *ata,*aa=a->a; 150 151 PetscFunctionBegin; 152 153 /* Set up timers */ 154 if (!logkey_mattranspose) { 155 ierr = PetscLogEventRegister(&logkey_mattranspose,"MatTranspose_SeqAIJ_FAST",MAT_COOKIE);CHKERRQ(ierr); 156 } 157 ierr = PetscLogEventBegin(logkey_mattranspose,A,0,0,0);CHKERRQ(ierr); 158 159 /* Allocate space for symbolic transpose info and work array */ 160 ierr = PetscMalloc((an+1)*sizeof(PetscInt),&ati);CHKERRQ(ierr); 161 ierr = PetscMalloc(ai[am]*sizeof(PetscInt),&atj);CHKERRQ(ierr); 162 ierr = PetscMalloc(ai[am]*sizeof(MatScalar),&ata);CHKERRQ(ierr); 163 ierr = PetscMalloc(an*sizeof(PetscInt),&atfill);CHKERRQ(ierr); 164 ierr = PetscMemzero(ati,(an+1)*sizeof(PetscInt));CHKERRQ(ierr); 165 /* Walk through aj and count ## of non-zeros in each row of A^T. */ 166 /* Note: offset by 1 for fast conversion into csr format. */ 167 for (i=0;i<ai[am];i++) { 168 ati[aj[i]+1] += 1; 169 } 170 /* Form ati for csr format of A^T. */ 171 for (i=0;i<an;i++) { 172 ati[i+1] += ati[i]; 173 } 174 175 /* Copy ati into atfill so we have locations of the next free space in atj */ 176 ierr = PetscMemcpy(atfill,ati,an*sizeof(PetscInt));CHKERRQ(ierr); 177 178 /* Walk through A row-wise and mark nonzero entries of A^T. */ 179 for (i=0;i<am;i++) { 180 anzj = ai[i+1] - ai[i]; 181 for (j=0;j<anzj;j++) { 182 atj[atfill[*aj]] = i; 183 ata[atfill[*aj]] = *aa++; 184 atfill[*aj++] += 1; 185 } 186 } 187 188 /* Clean up temporary space and complete requests. */ 189 ierr = PetscFree(atfill);CHKERRQ(ierr); 190 ierr = MatCreateSeqAIJWithArrays(A->comm,an,am,ati,atj,ata,&At);CHKERRQ(ierr); 191 at = (Mat_SeqAIJ *)(At->data); 192 at->freedata = PETSC_TRUE; 193 at->nonew = 0; 194 if (B) { 195 *B = At; 196 } else { 197 ierr = MatHeaderCopy(A,At); 198 } 199 ierr = PetscLogEventEnd(logkey_mattranspose,A,0,0,0);CHKERRQ(ierr); 200 PetscFunctionReturn(0); 201 } 202 203 #undef __FUNCT__ 204 #define __FUNCT__ "MatRestoreSymbolicTranspose_SeqAIJ" 205 PetscErrorCode MatRestoreSymbolicTranspose_SeqAIJ(Mat A,PetscInt *ati[],PetscInt *atj[]) 206 { 207 PetscErrorCode ierr; 208 209 PetscFunctionBegin; 210 ierr = PetscLogInfo(A,"Restoring Symbolic Transpose.\n");CHKERRQ(ierr); 211 ierr = PetscFree(*ati);CHKERRQ(ierr); 212 ati = PETSC_NULL; 213 ierr = PetscFree(*atj);CHKERRQ(ierr); 214 atj = PETSC_NULL; 215 PetscFunctionReturn(0); 216 } 217 218