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