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