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