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