1 #define PETSCMAT_DLL 2 3 /* 4 Defines symbolic transpose routines for SeqAIJ matrices. 5 6 Currently Get/Restore only allocates/frees memory for holding the 7 (i,j) info for the transpose. Someday, this info could be 8 maintained so successive calls to Get will not recompute the info. 9 10 Also defined is a "faster" implementation of MatTranspose for SeqAIJ 11 matrices which avoids calls to MatSetValues. This routine has not 12 been adopted as the standard yet as it is somewhat untested. 13 14 */ 15 16 #include "src/mat/impls/aij/seq/aij.h" 17 18 static PetscEvent logkey_matgetsymtranspose = 0; 19 static PetscEvent logkey_mattranspose = 0; 20 21 #undef __FUNCT__ 22 #define __FUNCT__ "MatGetSymbolicTranspose_SeqIJ" 23 PetscErrorCode MatGetSymbolicTranspose_SeqAIJ(Mat A,PetscInt *Ati[],PetscInt *Atj[]) 24 { 25 PetscErrorCode ierr; 26 PetscInt i,j,anzj; 27 Mat_SeqAIJ *a=(Mat_SeqAIJ *)A->data; 28 PetscInt an=A->N,am=A->M; 29 PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j; 30 31 PetscFunctionBegin; 32 33 ierr = PetscVerboseInfo((A,"MatGetSymbolicTranspose_SeqAIJ:Getting Symbolic Transpose.\n"));CHKERRQ(ierr); 34 35 /* Set up timers */ 36 if (!logkey_matgetsymtranspose) { 37 ierr = PetscLogEventRegister(&logkey_matgetsymtranspose,"MatGetSymbolicTranspose",MAT_COOKIE);CHKERRQ(ierr); 38 } 39 ierr = PetscLogEventBegin(logkey_matgetsymtranspose,A,0,0,0);CHKERRQ(ierr); 40 41 /* Allocate space for symbolic transpose info and work array */ 42 ierr = PetscMalloc((an+1)*sizeof(PetscInt),&ati);CHKERRQ(ierr); 43 ierr = PetscMalloc(ai[am]*sizeof(PetscInt),&atj);CHKERRQ(ierr); 44 ierr = PetscMalloc(an*sizeof(PetscInt),&atfill);CHKERRQ(ierr); 45 ierr = PetscMemzero(ati,(an+1)*sizeof(PetscInt));CHKERRQ(ierr); 46 47 /* Walk through aj and count ## of non-zeros in each row of A^T. */ 48 /* Note: offset by 1 for fast conversion into csr format. */ 49 for (i=0;i<ai[am];i++) { 50 ati[aj[i]+1] += 1; 51 } 52 /* Form ati for csr format of A^T. */ 53 for (i=0;i<an;i++) { 54 ati[i+1] += ati[i]; 55 } 56 57 /* Copy ati into atfill so we have locations of the next free space in atj */ 58 ierr = PetscMemcpy(atfill,ati,an*sizeof(PetscInt));CHKERRQ(ierr); 59 60 /* Walk through A row-wise and mark nonzero entries of A^T. */ 61 for (i=0;i<am;i++) { 62 anzj = ai[i+1] - ai[i]; 63 for (j=0;j<anzj;j++) { 64 atj[atfill[*aj]] = i; 65 atfill[*aj++] += 1; 66 } 67 } 68 69 /* Clean up temporary space and complete requests. */ 70 ierr = PetscFree(atfill);CHKERRQ(ierr); 71 *Ati = ati; 72 *Atj = atj; 73 74 ierr = PetscLogEventEnd(logkey_matgetsymtranspose,A,0,0,0);CHKERRQ(ierr); 75 PetscFunctionReturn(0); 76 } 77 /* 78 MatGetSymbolicTransposeReduced_SeqAIJ() - Get symbolic matrix structure of submatrix A[rstart:rend,:], 79 modified from MatGetSymbolicTranspose_SeqAIJ() 80 */ 81 static PetscEvent logkey_matgetsymtransreduced = 0; 82 #undef __FUNCT__ 83 #define __FUNCT__ "MatGetSymbolicTransposeReduced_SeqIJ" 84 PetscErrorCode MatGetSymbolicTransposeReduced_SeqAIJ(Mat A,PetscInt rstart,PetscInt rend,PetscInt *Ati[],PetscInt *Atj[]) 85 { 86 PetscErrorCode ierr; 87 PetscInt i,j,anzj; 88 Mat_SeqAIJ *a=(Mat_SeqAIJ *)A->data; 89 PetscInt an=A->N; 90 PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j; 91 92 PetscFunctionBegin; 93 ierr = PetscVerboseInfo((A,"MatGetSymbolicTransposeReduced_SeqAIJ:Getting Symbolic Transpose.\n"));CHKERRQ(ierr); 94 95 /* Set up timers */ 96 if (!logkey_matgetsymtransreduced) { 97 ierr = PetscLogEventRegister(&logkey_matgetsymtransreduced,"MatGetSymbolicTransposeReduced",MAT_COOKIE);CHKERRQ(ierr); 98 } 99 ierr = PetscLogEventBegin(logkey_matgetsymtransreduced,A,0,0,0);CHKERRQ(ierr); 100 101 /* Allocate space for symbolic transpose info and work array */ 102 ierr = PetscMalloc((an+1)*sizeof(PetscInt),&ati);CHKERRQ(ierr); 103 anzj = ai[rend] - ai[rstart]; 104 ierr = PetscMalloc((anzj+1)*sizeof(PetscInt),&atj);CHKERRQ(ierr); 105 ierr = PetscMalloc((an+1)*sizeof(PetscInt),&atfill);CHKERRQ(ierr); 106 ierr = PetscMemzero(ati,(an+1)*sizeof(PetscInt));CHKERRQ(ierr); 107 108 /* Walk through aj and count ## of non-zeros in each row of A^T. */ 109 /* Note: offset by 1 for fast conversion into csr format. */ 110 for (i=ai[rstart]; i<ai[rend]; i++) { 111 ati[aj[i]+1] += 1; 112 } 113 /* Form ati for csr format of A^T. */ 114 for (i=0;i<an;i++) { 115 ati[i+1] += ati[i]; 116 } 117 118 /* Copy ati into atfill so we have locations of the next free space in atj */ 119 ierr = PetscMemcpy(atfill,ati,an*sizeof(PetscInt));CHKERRQ(ierr); 120 121 /* Walk through A row-wise and mark nonzero entries of A^T. */ 122 aj = aj + ai[rstart]; 123 for (i=rstart; i<rend; i++) { 124 anzj = ai[i+1] - ai[i]; 125 for (j=0;j<anzj;j++) { 126 atj[atfill[*aj]] = i-rstart; 127 atfill[*aj++] += 1; 128 } 129 } 130 131 /* Clean up temporary space and complete requests. */ 132 ierr = PetscFree(atfill);CHKERRQ(ierr); 133 *Ati = ati; 134 *Atj = atj; 135 136 ierr = PetscLogEventEnd(logkey_matgetsymtransreduced,A,0,0,0);CHKERRQ(ierr); 137 PetscFunctionReturn(0); 138 } 139 140 #undef __FUNCT__ 141 #define __FUNCT__ "MatTranspose_SeqIJ_FAST" 142 PetscErrorCode MatTranspose_SeqAIJ_FAST(Mat A,Mat *B) 143 { 144 PetscErrorCode ierr; 145 PetscInt i,j,anzj; 146 Mat At; 147 Mat_SeqAIJ *a=(Mat_SeqAIJ *)A->data,*at; 148 PetscInt an=A->N,am=A->M; 149 PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j; 150 MatScalar *ata,*aa=a->a; 151 152 PetscFunctionBegin; 153 154 /* Set up timers */ 155 if (!logkey_mattranspose) { 156 ierr = PetscLogEventRegister(&logkey_mattranspose,"MatTranspose_SeqAIJ_FAST",MAT_COOKIE);CHKERRQ(ierr); 157 } 158 ierr = PetscLogEventBegin(logkey_mattranspose,A,0,0,0);CHKERRQ(ierr); 159 160 /* Allocate space for symbolic transpose info and work array */ 161 ierr = PetscMalloc((an+1)*sizeof(PetscInt),&ati);CHKERRQ(ierr); 162 ierr = PetscMalloc(ai[am]*sizeof(PetscInt),&atj);CHKERRQ(ierr); 163 ierr = PetscMalloc(ai[am]*sizeof(MatScalar),&ata);CHKERRQ(ierr); 164 ierr = PetscMalloc(an*sizeof(PetscInt),&atfill);CHKERRQ(ierr); 165 ierr = PetscMemzero(ati,(an+1)*sizeof(PetscInt));CHKERRQ(ierr); 166 /* Walk through aj and count ## of non-zeros in each row of A^T. */ 167 /* Note: offset by 1 for fast conversion into csr format. */ 168 for (i=0;i<ai[am];i++) { 169 ati[aj[i]+1] += 1; 170 } 171 /* Form ati for csr format of A^T. */ 172 for (i=0;i<an;i++) { 173 ati[i+1] += ati[i]; 174 } 175 176 /* Copy ati into atfill so we have locations of the next free space in atj */ 177 ierr = PetscMemcpy(atfill,ati,an*sizeof(PetscInt));CHKERRQ(ierr); 178 179 /* Walk through A row-wise and mark nonzero entries of A^T. */ 180 for (i=0;i<am;i++) { 181 anzj = ai[i+1] - ai[i]; 182 for (j=0;j<anzj;j++) { 183 atj[atfill[*aj]] = i; 184 ata[atfill[*aj]] = *aa++; 185 atfill[*aj++] += 1; 186 } 187 } 188 189 /* Clean up temporary space and complete requests. */ 190 ierr = PetscFree(atfill);CHKERRQ(ierr); 191 ierr = MatCreateSeqAIJWithArrays(A->comm,an,am,ati,atj,ata,&At);CHKERRQ(ierr); 192 at = (Mat_SeqAIJ *)(At->data); 193 at->freedata = PETSC_TRUE; 194 at->nonew = 0; 195 if (B) { 196 *B = At; 197 } else { 198 ierr = MatHeaderCopy(A,At); 199 } 200 ierr = PetscLogEventEnd(logkey_mattranspose,A,0,0,0);CHKERRQ(ierr); 201 PetscFunctionReturn(0); 202 } 203 204 #undef __FUNCT__ 205 #define __FUNCT__ "MatRestoreSymbolicTranspose_SeqAIJ" 206 PetscErrorCode MatRestoreSymbolicTranspose_SeqAIJ(Mat A,PetscInt *ati[],PetscInt *atj[]) 207 { 208 PetscErrorCode ierr; 209 210 PetscFunctionBegin; 211 ierr = PetscVerboseInfo((A,"MatRestoreSymbolicTranspose_SeqAIJ:Restoring Symbolic Transpose.\n"));CHKERRQ(ierr); 212 ierr = PetscFree(*ati);CHKERRQ(ierr); 213 ati = PETSC_NULL; 214 ierr = PetscFree(*atj);CHKERRQ(ierr); 215 atj = PETSC_NULL; 216 PetscFunctionReturn(0); 217 } 218 219