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