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,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 /* 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 = PetscMalloc(an*sizeof(PetscInt),&atfill);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 161 /* Copy ati into atfill so we have locations of the next free space in atj */ 162 ierr = PetscMemcpy(atfill,ati,an*sizeof(PetscInt));CHKERRQ(ierr); 163 164 /* Walk through A row-wise and mark nonzero entries of A^T. */ 165 for (i=0;i<am;i++) { 166 anzj = ai[i+1] - ai[i]; 167 for (j=0;j<anzj;j++) { 168 atj[atfill[*aj]] = i; 169 ata[atfill[*aj]] = *aa++; 170 atfill[*aj++] += 1; 171 } 172 } 173 174 /* Clean up temporary space and complete requests. */ 175 ierr = PetscFree(atfill);CHKERRQ(ierr); 176 ierr = MatCreateSeqAIJWithArrays(((PetscObject)A)->comm,an,am,ati,atj,ata,&At);CHKERRQ(ierr); 177 at = (Mat_SeqAIJ *)(At->data); 178 at->free_a = PETSC_TRUE; 179 at->free_ij = PETSC_TRUE; 180 at->nonew = 0; 181 if (B) { 182 *B = At; 183 } else { 184 ierr = MatHeaderCopy(A,At); 185 } 186 ierr = PetscLogEventEnd(MAT_Transpose_SeqAIJ,A,0,0,0);CHKERRQ(ierr); 187 PetscFunctionReturn(0); 188 } 189 190 #undef __FUNCT__ 191 #define __FUNCT__ "MatRestoreSymbolicTranspose_SeqAIJ" 192 PetscErrorCode MatRestoreSymbolicTranspose_SeqAIJ(Mat A,PetscInt *ati[],PetscInt *atj[]) 193 { 194 PetscErrorCode ierr; 195 196 PetscFunctionBegin; 197 ierr = PetscInfo(A,"Restoring Symbolic Transpose.\n");CHKERRQ(ierr); 198 ierr = PetscFree(*ati);CHKERRQ(ierr); 199 ati = PETSC_NULL; 200 ierr = PetscFree(*atj);CHKERRQ(ierr); 201 atj = PETSC_NULL; 202 PetscFunctionReturn(0); 203 } 204 205