1 2 /* 3 Defines transpose routines for SeqAIJ matrices. 4 */ 5 6 #include <../src/mat/impls/aij/seq/aij.h> 7 8 /* 9 The symbolic and full transpose versions share several similar code blocks but the macros to reuse the code would be confusing and ugly 10 */ 11 PetscErrorCode MatTransposeSymbolic_SeqAIJ(Mat A,Mat *B) 12 { 13 PetscInt i,j,anzj; 14 Mat At; 15 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*at; 16 PetscInt an=A->cmap->N,am=A->rmap->N; 17 PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j; 18 19 PetscFunctionBegin; 20 /* Allocate space for symbolic transpose info and work array */ 21 PetscCall(PetscCalloc1(an+1,&ati)); 22 PetscCall(PetscMalloc1(ai[am],&atj)); 23 24 /* Walk through aj and count ## of non-zeros in each row of A^T. */ 25 /* Note: offset by 1 for fast conversion into csr format. */ 26 for (i=0;i<ai[am];i++) ati[aj[i]+1] += 1; 27 /* Form ati for csr format of A^T. */ 28 for (i=0;i<an;i++) ati[i+1] += ati[i]; 29 30 /* Copy ati into atfill so we have locations of the next free space in atj */ 31 PetscCall(PetscMalloc1(an,&atfill)); 32 PetscCall(PetscArraycpy(atfill,ati,an)); 33 34 /* Walk through A row-wise and mark nonzero entries of A^T. */ 35 for (i=0;i<am;i++) { 36 anzj = ai[i+1] - ai[i]; 37 for (j=0;j<anzj;j++) { 38 atj[atfill[*aj]] = i; 39 atfill[*aj++] += 1; 40 } 41 } 42 PetscCall(PetscFree(atfill)); 43 44 PetscCall(MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),an,am,ati,atj,NULL,&At)); 45 PetscCall(MatSetBlockSizes(At,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs))); 46 PetscCall(MatSetType(At,((PetscObject)A)->type_name)); 47 at = (Mat_SeqAIJ*)At->data; 48 at->free_a = PETSC_FALSE; 49 at->free_ij = PETSC_TRUE; 50 at->nonew = 0; 51 at->maxnz = ati[an]; 52 *B = At; 53 PetscFunctionReturn(0); 54 } 55 56 PetscErrorCode MatTranspose_SeqAIJ(Mat A,MatReuse reuse,Mat *B) 57 { 58 PetscInt i,j,anzj; 59 Mat At; 60 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*at; 61 PetscInt an=A->cmap->N,am=A->rmap->N; 62 PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j; 63 MatScalar *ata; 64 const MatScalar *aa,*av; 65 PetscContainer rB; 66 MatParentState *rb; 67 PetscBool nonzerochange = PETSC_FALSE; 68 69 PetscFunctionBegin; 70 if (reuse == MAT_REUSE_MATRIX) { 71 PetscCall(PetscObjectQuery((PetscObject)*B,"MatTransposeParent",(PetscObject*)&rB)); 72 PetscCheck(rB,PetscObjectComm((PetscObject)*B),PETSC_ERR_ARG_WRONG,"Reuse matrix used was not generated from call to MatTranspose()"); 73 PetscCall(PetscContainerGetPointer(rB,(void**)&rb)); 74 if (rb->nonzerostate != A->nonzerostate) nonzerochange = PETSC_TRUE; 75 } 76 77 PetscCall(MatSeqAIJGetArrayRead(A,&av)); 78 aa = av; 79 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX || nonzerochange) { 80 /* Allocate space for symbolic transpose info and work array */ 81 PetscCall(PetscCalloc1(an+1,&ati)); 82 PetscCall(PetscMalloc1(ai[am],&atj)); 83 /* Walk through aj and count ## of non-zeros in each row of A^T. */ 84 /* Note: offset by 1 for fast conversion into csr format. */ 85 for (i=0;i<ai[am];i++) ati[aj[i]+1] += 1; 86 /* Form ati for csr format of A^T. */ 87 for (i=0;i<an;i++) ati[i+1] += ati[i]; 88 PetscCall(PetscMalloc1(ai[am],&ata)); 89 } else { 90 Mat_SeqAIJ *sub_B = (Mat_SeqAIJ*) (*B)->data; 91 ati = sub_B->i; 92 atj = sub_B->j; 93 ata = sub_B->a; 94 At = *B; 95 } 96 97 /* Copy ati into atfill so we have locations of the next free space in atj */ 98 PetscCall(PetscMalloc1(an,&atfill)); 99 PetscCall(PetscArraycpy(atfill,ati,an)); 100 101 /* Walk through A row-wise and mark nonzero entries of A^T. */ 102 for (i=0;i<am;i++) { 103 anzj = ai[i+1] - ai[i]; 104 for (j=0;j<anzj;j++) { 105 atj[atfill[*aj]] = i; 106 ata[atfill[*aj]] = *aa++; 107 atfill[*aj++] += 1; 108 } 109 } 110 PetscCall(PetscFree(atfill)); 111 PetscCall(MatSeqAIJRestoreArrayRead(A,&av)); 112 if (reuse == MAT_REUSE_MATRIX) PetscCall(PetscObjectStateIncrease((PetscObject)(*B))); 113 114 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX || nonzerochange) { 115 PetscCall(MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),an,am,ati,atj,ata,&At)); 116 PetscCall(MatSetBlockSizes(At,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs))); 117 PetscCall(MatSetType(At,((PetscObject)A)->type_name)); 118 at = (Mat_SeqAIJ*)(At->data); 119 at->free_a = PETSC_TRUE; 120 at->free_ij = PETSC_TRUE; 121 at->nonew = 0; 122 at->maxnz = ati[an]; 123 } 124 125 if (reuse == MAT_INITIAL_MATRIX || (reuse == MAT_REUSE_MATRIX && !nonzerochange)) { 126 *B = At; 127 } else if (nonzerochange) { 128 PetscCall(MatHeaderMerge(*B,&At)); 129 PetscCall(MatTransposeSetPrecursor(A,*B)); 130 } else if (reuse == MAT_INPLACE_MATRIX) { 131 PetscCall(MatHeaderMerge(A,&At)); 132 } 133 PetscFunctionReturn(0); 134 } 135 136 /* 137 Get symbolic matrix structure of a submatrix of A, A[rstart:rend,:], 138 */ 139 PetscErrorCode MatGetSymbolicTransposeReduced_SeqAIJ(Mat A,PetscInt rstart,PetscInt rend,PetscInt *Ati[],PetscInt *Atj[]) 140 { 141 PetscInt i,j,anzj; 142 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data; 143 PetscInt an=A->cmap->N; 144 PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j, am=ai[rend] - ai[rstart]; 145 146 PetscFunctionBegin; 147 PetscCall(PetscLogEventBegin(MAT_Getsymtransreduced,A,0,0,0)); 148 149 /* Allocate space for symbolic transpose info and work array */ 150 PetscCall(PetscCalloc1(an+1,&ati)); 151 PetscCall(PetscMalloc1(am+1,&atj)); 152 153 /* Walk through aj and count ## of non-zeros in each row of A^T. */ 154 /* Note: offset by 1 for fast conversion into csr format. */ 155 for (i=ai[rstart]; i<ai[rend]; i++) ati[aj[i]+1] += 1; 156 /* Form ati for csr format of A^T. */ 157 for (i=0;i<an;i++) ati[i+1] += ati[i]; 158 159 /* Copy ati into atfill so we have locations of the next free space in atj */ 160 PetscCall(PetscMalloc1(an+1,&atfill)); 161 PetscCall(PetscArraycpy(atfill,ati,an)); 162 163 /* Walk through A row-wise and mark nonzero entries of A^T. */ 164 aj = aj + ai[rstart]; 165 for (i=rstart; i<rend; i++) { 166 anzj = ai[i+1] - ai[i]; 167 for (j=0; j<anzj; j++) { 168 atj[atfill[*aj]] = i-rstart; 169 atfill[*aj++] += 1; 170 } 171 } 172 PetscCall(PetscFree(atfill)); 173 *Ati = ati; 174 *Atj = atj; 175 176 PetscCall(PetscLogEventEnd(MAT_Getsymtransreduced,A,0,0,0)); 177 PetscFunctionReturn(0); 178 } 179 180 /* 181 Returns the i and j arrays for a symbolic transpose, this is used internally within SeqAIJ code when the full 182 symbolic matrix (which can be obtained with MatTransposeSymbolic() is not needed. MatRestoreSymbolicTranspose_SeqAIJ() should be used to free the arrays. 183 */ 184 PetscErrorCode MatGetSymbolicTranspose_SeqAIJ(Mat A,PetscInt *Ati[],PetscInt *Atj[]) 185 { 186 PetscFunctionBegin; 187 PetscCall(MatGetSymbolicTransposeReduced_SeqAIJ(A,0,A->rmap->N,Ati,Atj)); 188 PetscFunctionReturn(0); 189 } 190 191 PetscErrorCode MatRestoreSymbolicTranspose_SeqAIJ(Mat A,PetscInt *ati[],PetscInt *atj[]) 192 { 193 PetscFunctionBegin; 194 PetscCall(PetscFree(*ati)); 195 PetscCall(PetscFree(*atj)); 196 PetscFunctionReturn(0); 197 } 198 199