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 is the new 11 standard since it is much faster than MatTranspose_AIJ. 12 13 */ 14 15 #include <../src/mat/impls/aij/seq/aij.h> 16 17 18 PetscErrorCode MatGetSymbolicTranspose_SeqAIJ(Mat A,PetscInt *Ati[],PetscInt *Atj[]) 19 { 20 PetscErrorCode ierr; 21 PetscInt i,j,anzj; 22 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data; 23 PetscInt an=A->cmap->N,am=A->rmap->N; 24 PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j; 25 26 PetscFunctionBegin; 27 ierr = PetscInfo(A,"Getting Symbolic Transpose.\n");CHKERRQ(ierr); 28 29 /* Set up timers */ 30 ierr = PetscLogEventBegin(MAT_Getsymtranspose,A,0,0,0);CHKERRQ(ierr); 31 32 /* Allocate space for symbolic transpose info and work array */ 33 ierr = PetscCalloc1(an+1,&ati);CHKERRQ(ierr); 34 ierr = PetscMalloc1(ai[am],&atj);CHKERRQ(ierr); 35 ierr = PetscMalloc1(an,&atfill);CHKERRQ(ierr); 36 37 /* Walk through aj and count ## of non-zeros in each row of A^T. */ 38 /* Note: offset by 1 for fast conversion into csr format. */ 39 for (i=0;i<ai[am];i++) { 40 ati[aj[i]+1] += 1; 41 } 42 /* Form ati for csr format of A^T. */ 43 for (i=0;i<an;i++) { 44 ati[i+1] += ati[i]; 45 } 46 47 /* Copy ati into atfill so we have locations of the next free space in atj */ 48 ierr = PetscArraycpy(atfill,ati,an);CHKERRQ(ierr); 49 50 /* Walk through A row-wise and mark nonzero entries of A^T. */ 51 for (i=0; i<am; i++) { 52 anzj = ai[i+1] - ai[i]; 53 for (j=0; j<anzj; j++) { 54 atj[atfill[*aj]] = i; 55 atfill[*aj++] += 1; 56 } 57 } 58 59 /* Clean up temporary space and complete requests. */ 60 ierr = PetscFree(atfill);CHKERRQ(ierr); 61 *Ati = ati; 62 *Atj = atj; 63 64 ierr = PetscLogEventEnd(MAT_Getsymtranspose,A,0,0,0);CHKERRQ(ierr); 65 PetscFunctionReturn(0); 66 } 67 /* 68 MatGetSymbolicTransposeReduced_SeqAIJ() - Get symbolic matrix structure of submatrix A[rstart:rend,:], 69 modified from MatGetSymbolicTranspose_SeqAIJ() 70 */ 71 PetscErrorCode MatGetSymbolicTransposeReduced_SeqAIJ(Mat A,PetscInt rstart,PetscInt rend,PetscInt *Ati[],PetscInt *Atj[]) 72 { 73 PetscErrorCode ierr; 74 PetscInt i,j,anzj; 75 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data; 76 PetscInt an=A->cmap->N; 77 PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j; 78 79 PetscFunctionBegin; 80 ierr = PetscInfo(A,"Getting Symbolic Transpose\n");CHKERRQ(ierr); 81 ierr = PetscLogEventBegin(MAT_Getsymtransreduced,A,0,0,0);CHKERRQ(ierr); 82 83 /* Allocate space for symbolic transpose info and work array */ 84 ierr = PetscCalloc1(an+1,&ati);CHKERRQ(ierr); 85 anzj = ai[rend] - ai[rstart]; 86 ierr = PetscMalloc1(anzj+1,&atj);CHKERRQ(ierr); 87 ierr = PetscMalloc1(an+1,&atfill);CHKERRQ(ierr); 88 89 /* Walk through aj and count ## of non-zeros in each row of A^T. */ 90 /* Note: offset by 1 for fast conversion into csr format. */ 91 for (i=ai[rstart]; i<ai[rend]; i++) { 92 ati[aj[i]+1] += 1; 93 } 94 /* Form ati for csr format of A^T. */ 95 for (i=0;i<an;i++) { 96 ati[i+1] += ati[i]; 97 } 98 99 /* Copy ati into atfill so we have locations of the next free space in atj */ 100 ierr = PetscArraycpy(atfill,ati,an);CHKERRQ(ierr); 101 102 /* Walk through A row-wise and mark nonzero entries of A^T. */ 103 aj = aj + ai[rstart]; 104 for (i=rstart; i<rend; i++) { 105 anzj = ai[i+1] - ai[i]; 106 for (j=0; j<anzj; j++) { 107 atj[atfill[*aj]] = i-rstart; 108 atfill[*aj++] += 1; 109 } 110 } 111 112 /* Clean up temporary space and complete requests. */ 113 ierr = PetscFree(atfill);CHKERRQ(ierr); 114 *Ati = ati; 115 *Atj = atj; 116 117 ierr = PetscLogEventEnd(MAT_Getsymtransreduced,A,0,0,0);CHKERRQ(ierr); 118 PetscFunctionReturn(0); 119 } 120 121 PetscErrorCode MatTranspose_SeqAIJ(Mat A,MatReuse reuse,Mat *B) 122 { 123 PetscErrorCode ierr; 124 PetscInt i,j,anzj; 125 Mat At; 126 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*at; 127 PetscInt an=A->cmap->N,am=A->rmap->N; 128 PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j; 129 MatScalar *ata; 130 const MatScalar *aa,*av; 131 132 PetscFunctionBegin; 133 ierr = MatSeqAIJGetArrayRead(A,&av);CHKERRQ(ierr); 134 aa = av; 135 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) { 136 /* Allocate space for symbolic transpose info and work array */ 137 ierr = PetscCalloc1(an+1,&ati);CHKERRQ(ierr); 138 ierr = PetscMalloc1(ai[am],&atj);CHKERRQ(ierr); 139 ierr = PetscMalloc1(ai[am],&ata);CHKERRQ(ierr); 140 /* Walk through aj and count ## of non-zeros in each row of A^T. */ 141 /* Note: offset by 1 for fast conversion into csr format. */ 142 for (i=0;i<ai[am];i++) { 143 ati[aj[i]+1] += 1; /* count ## of non-zeros for row aj[i] of A^T */ 144 } 145 /* Form ati for csr format of A^T. */ 146 for (i=0;i<an;i++) { 147 ati[i+1] += ati[i]; 148 } 149 } else { /* This segment is called by MatTranspose_MPIAIJ(...,MAT_INITIAL_MATRIX,..) directly! */ 150 Mat_SeqAIJ *sub_B = (Mat_SeqAIJ*) (*B)->data; 151 ati = sub_B->i; 152 atj = sub_B->j; 153 ata = sub_B->a; 154 At = *B; 155 } 156 157 /* Copy ati into atfill so we have locations of the next free space in atj */ 158 ierr = PetscMalloc1(an,&atfill);CHKERRQ(ierr); 159 ierr = PetscArraycpy(atfill,ati,an);CHKERRQ(ierr); 160 161 /* Walk through A row-wise and mark nonzero entries of A^T. */ 162 for (i=0;i<am;i++) { 163 anzj = ai[i+1] - ai[i]; 164 for (j=0;j<anzj;j++) { 165 atj[atfill[*aj]] = i; 166 ata[atfill[*aj]] = *aa++; 167 atfill[*aj++] += 1; 168 } 169 } 170 ierr = MatSeqAIJRestoreArrayRead(A,&av);CHKERRQ(ierr); 171 172 /* Clean up temporary space and complete requests. */ 173 ierr = PetscFree(atfill);CHKERRQ(ierr); 174 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) { 175 ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),an,am,ati,atj,ata,&At);CHKERRQ(ierr); 176 ierr = MatSetBlockSizes(At,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));CHKERRQ(ierr); 177 178 at = (Mat_SeqAIJ*)(At->data); 179 at->free_a = PETSC_TRUE; 180 at->free_ij = PETSC_TRUE; 181 at->nonew = 0; 182 at->maxnz = ati[an]; 183 184 ierr = MatSetType(At,((PetscObject)A)->type_name);CHKERRQ(ierr); 185 } 186 187 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) { 188 *B = At; 189 } else { 190 ierr = MatHeaderMerge(A,&At);CHKERRQ(ierr); 191 } 192 #if defined(PETSC_HAVE_DEVICE) 193 (*B)->offloadmask = PETSC_OFFLOAD_CPU; 194 #endif 195 PetscFunctionReturn(0); 196 } 197 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