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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 197 } 198