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 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->cmap->N,am=A->rmap->N; 27 PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j; 28 29 PetscFunctionBegin; 30 31 ierr = PetscInfo(A,"Getting Symbolic Transpose.\n");CHKERRQ(ierr); 32 33 /* Set up timers */ 34 ierr = PetscLogEventBegin(MAT_Getsymtranspose,A,0,0,0);CHKERRQ(ierr); 35 36 /* Allocate space for symbolic transpose info and work array */ 37 ierr = PetscMalloc((an+1)*sizeof(PetscInt),&ati);CHKERRQ(ierr); 38 ierr = PetscMalloc(ai[am]*sizeof(PetscInt),&atj);CHKERRQ(ierr); 39 ierr = PetscMalloc(an*sizeof(PetscInt),&atfill);CHKERRQ(ierr); 40 ierr = PetscMemzero(ati,(an+1)*sizeof(PetscInt));CHKERRQ(ierr); 41 42 /* Walk through aj and count ## of non-zeros in each row of A^T. */ 43 /* Note: offset by 1 for fast conversion into csr format. */ 44 for (i=0;i<ai[am];i++) { 45 ati[aj[i]+1] += 1; 46 } 47 /* Form ati for csr format of A^T. */ 48 for (i=0;i<an;i++) { 49 ati[i+1] += ati[i]; 50 } 51 52 /* Copy ati into atfill so we have locations of the next free space in atj */ 53 ierr = PetscMemcpy(atfill,ati,an*sizeof(PetscInt));CHKERRQ(ierr); 54 55 /* Walk through A row-wise and mark nonzero entries of A^T. */ 56 for (i=0;i<am;i++) { 57 anzj = ai[i+1] - ai[i]; 58 for (j=0;j<anzj;j++) { 59 atj[atfill[*aj]] = i; 60 atfill[*aj++] += 1; 61 } 62 } 63 64 /* Clean up temporary space and complete requests. */ 65 ierr = PetscFree(atfill);CHKERRQ(ierr); 66 *Ati = ati; 67 *Atj = atj; 68 69 ierr = PetscLogEventEnd(MAT_Getsymtranspose,A,0,0,0);CHKERRQ(ierr); 70 PetscFunctionReturn(0); 71 } 72 /* 73 MatGetSymbolicTransposeReduced_SeqAIJ() - Get symbolic matrix structure of submatrix A[rstart:rend,:], 74 modified from MatGetSymbolicTranspose_SeqAIJ() 75 */ 76 #undef __FUNCT__ 77 #define __FUNCT__ "MatGetSymbolicTransposeReduced_SeqIJ" 78 PetscErrorCode MatGetSymbolicTransposeReduced_SeqAIJ(Mat A,PetscInt rstart,PetscInt rend,PetscInt *Ati[],PetscInt *Atj[]) 79 { 80 PetscErrorCode ierr; 81 PetscInt i,j,anzj; 82 Mat_SeqAIJ *a=(Mat_SeqAIJ *)A->data; 83 PetscInt an=A->cmap->N; 84 PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j; 85 86 PetscFunctionBegin; 87 ierr = PetscInfo(A,"Getting Symbolic Transpose\n");CHKERRQ(ierr); 88 ierr = PetscLogEventBegin(MAT_Getsymtransreduced,A,0,0,0);CHKERRQ(ierr); 89 90 /* Allocate space for symbolic transpose info and work array */ 91 ierr = PetscMalloc((an+1)*sizeof(PetscInt),&ati);CHKERRQ(ierr); 92 anzj = ai[rend] - ai[rstart]; 93 ierr = PetscMalloc((anzj+1)*sizeof(PetscInt),&atj);CHKERRQ(ierr); 94 ierr = PetscMalloc((an+1)*sizeof(PetscInt),&atfill);CHKERRQ(ierr); 95 ierr = PetscMemzero(ati,(an+1)*sizeof(PetscInt));CHKERRQ(ierr); 96 97 /* Walk through aj and count ## of non-zeros in each row of A^T. */ 98 /* Note: offset by 1 for fast conversion into csr format. */ 99 for (i=ai[rstart]; i<ai[rend]; i++) { 100 ati[aj[i]+1] += 1; 101 } 102 /* Form ati for csr format of A^T. */ 103 for (i=0;i<an;i++) { 104 ati[i+1] += ati[i]; 105 } 106 107 /* Copy ati into atfill so we have locations of the next free space in atj */ 108 ierr = PetscMemcpy(atfill,ati,an*sizeof(PetscInt));CHKERRQ(ierr); 109 110 /* Walk through A row-wise and mark nonzero entries of A^T. */ 111 aj = aj + ai[rstart]; 112 for (i=rstart; i<rend; i++) { 113 anzj = ai[i+1] - ai[i]; 114 for (j=0;j<anzj;j++) { 115 atj[atfill[*aj]] = i-rstart; 116 atfill[*aj++] += 1; 117 } 118 } 119 120 /* Clean up temporary space and complete requests. */ 121 ierr = PetscFree(atfill);CHKERRQ(ierr); 122 *Ati = ati; 123 *Atj = atj; 124 125 ierr = PetscLogEventEnd(MAT_Getsymtransreduced,A,0,0,0);CHKERRQ(ierr); 126 PetscFunctionReturn(0); 127 } 128 129 #undef __FUNCT__ 130 #define __FUNCT__ "MatTranspose_SeqIJ_FAST" 131 PetscErrorCode MatTranspose_SeqAIJ_FAST(Mat A,MatReuse reuse,Mat *B) 132 { 133 PetscErrorCode ierr; 134 PetscInt i,j,anzj; 135 Mat At; 136 Mat_SeqAIJ *a=(Mat_SeqAIJ *)A->data,*at; 137 PetscInt an=A->cmap->N,am=A->rmap->N; 138 PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j; 139 MatScalar *ata,*aa=a->a; 140 141 PetscFunctionBegin; 142 143 ierr = PetscLogEventBegin(MAT_Transpose_SeqAIJ,A,0,0,0);CHKERRQ(ierr); 144 145 if (reuse == MAT_INITIAL_MATRIX || *B == A) { 146 /* Allocate space for symbolic transpose info and work array */ 147 ierr = PetscMalloc((an+1)*sizeof(PetscInt),&ati);CHKERRQ(ierr); 148 ierr = PetscMalloc(ai[am]*sizeof(PetscInt),&atj);CHKERRQ(ierr); 149 ierr = PetscMalloc(ai[am]*sizeof(MatScalar),&ata);CHKERRQ(ierr); 150 ierr = PetscMemzero(ati,(an+1)*sizeof(PetscInt));CHKERRQ(ierr); 151 /* Walk through aj and count ## of non-zeros in each row of A^T. */ 152 /* Note: offset by 1 for fast conversion into csr format. */ 153 for (i=0;i<ai[am];i++) { 154 ati[aj[i]+1] += 1; 155 } 156 /* Form ati for csr format of A^T. */ 157 for (i=0;i<an;i++) { 158 ati[i+1] += ati[i]; 159 } 160 } else { 161 Mat_SeqAIJ *sub_B = (Mat_SeqAIJ*) (*B)->data; 162 ati = sub_B->i; 163 atj = sub_B->j; 164 ata = sub_B->a; 165 At = *B; 166 } 167 168 /* Copy ati into atfill so we have locations of the next free space in atj */ 169 ierr = PetscMalloc(an*sizeof(PetscInt),&atfill);CHKERRQ(ierr); 170 ierr = PetscMemcpy(atfill,ati,an*sizeof(PetscInt));CHKERRQ(ierr); 171 172 /* Walk through A row-wise and mark nonzero entries of A^T. */ 173 for (i=0;i<am;i++) { 174 anzj = ai[i+1] - ai[i]; 175 for (j=0;j<anzj;j++) { 176 atj[atfill[*aj]] = i; 177 ata[atfill[*aj]] = *aa++; 178 atfill[*aj++] += 1; 179 } 180 } 181 182 /* Clean up temporary space and complete requests. */ 183 ierr = PetscFree(atfill);CHKERRQ(ierr); 184 if (reuse == MAT_INITIAL_MATRIX) { 185 ierr = MatCreateSeqAIJWithArrays(((PetscObject)A)->comm,an,am,ati,atj,ata,&At);CHKERRQ(ierr); 186 at = (Mat_SeqAIJ *)(At->data); 187 at->free_a = PETSC_TRUE; 188 at->free_ij = PETSC_TRUE; 189 at->nonew = 0; 190 } 191 192 if (reuse == MAT_INITIAL_MATRIX || *B != A) { 193 *B = At; 194 } else { 195 ierr = MatHeaderCopy(A,At); 196 } 197 ierr = PetscLogEventEnd(MAT_Transpose_SeqAIJ,A,0,0,0);CHKERRQ(ierr); 198 PetscFunctionReturn(0); 199 } 200 201 #undef __FUNCT__ 202 #define __FUNCT__ "MatRestoreSymbolicTranspose_SeqAIJ" 203 PetscErrorCode MatRestoreSymbolicTranspose_SeqAIJ(Mat A,PetscInt *ati[],PetscInt *atj[]) 204 { 205 PetscErrorCode ierr; 206 207 PetscFunctionBegin; 208 ierr = PetscInfo(A,"Restoring Symbolic Transpose.\n");CHKERRQ(ierr); 209 ierr = PetscFree(*ati);CHKERRQ(ierr); 210 ati = PETSC_NULL; 211 ierr = PetscFree(*atj);CHKERRQ(ierr); 212 atj = PETSC_NULL; 213 PetscFunctionReturn(0); 214 } 215 216