Actual source code: symtranspose.c


  2: /*
  3:   Defines transpose routines for SeqAIJ matrices.
  4: */

  6: #include <../src/mat/impls/aij/seq/aij.h>

  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;

 19:   PetscFunctionBegin;
 20:   /* Allocate space for symbolic transpose info and work array */
 21:   PetscCall(PetscCalloc1(an + 1, &ati));
 22:   PetscCall(PetscMalloc1(ai[am], &atj));

 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];

 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));

 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));

 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: }

 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;

 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:   }

 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:   }

 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));

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)));

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:   }

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: }

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];

146:   PetscFunctionBegin;
147:   PetscCall(PetscLogEventBegin(MAT_Getsymtransreduced, A, 0, 0, 0));

149:   /* Allocate space for symbolic transpose info and work array */
150:   PetscCall(PetscCalloc1(an + 1, &ati));
151:   PetscCall(PetscMalloc1(am + 1, &atj));

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];

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));

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;

176:   PetscCall(PetscLogEventEnd(MAT_Getsymtransreduced, A, 0, 0, 0));
177:   PetscFunctionReturn(PETSC_SUCCESS);
178: }

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: }

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: }