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:   /* Allocate space for symbolic transpose info and work array */
 20:   PetscCalloc1(an + 1, &ati);
 21:   PetscMalloc1(ai[am], &atj);

 23:   /* Walk through aj and count ## of non-zeros in each row of A^T. */
 24:   /* Note: offset by 1 for fast conversion into csr format. */
 25:   for (i = 0; i < ai[am]; i++) ati[aj[i] + 1] += 1;
 26:   /* Form ati for csr format of A^T. */
 27:   for (i = 0; i < an; i++) ati[i + 1] += ati[i];

 29:   /* Copy ati into atfill so we have locations of the next free space in atj */
 30:   PetscMalloc1(an, &atfill);
 31:   PetscArraycpy(atfill, ati, an);

 33:   /* Walk through A row-wise and mark nonzero entries of A^T. */
 34:   for (i = 0; i < am; i++) {
 35:     anzj = ai[i + 1] - ai[i];
 36:     for (j = 0; j < anzj; j++) {
 37:       atj[atfill[*aj]] = i;
 38:       atfill[*aj++] += 1;
 39:     }
 40:   }
 41:   PetscFree(atfill);

 43:   MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A), an, am, ati, atj, NULL, &At);
 44:   MatSetBlockSizes(At, PetscAbs(A->cmap->bs), PetscAbs(A->rmap->bs));
 45:   MatSetType(At, ((PetscObject)A)->type_name);
 46:   at          = (Mat_SeqAIJ *)At->data;
 47:   at->free_a  = PETSC_FALSE;
 48:   at->free_ij = PETSC_TRUE;
 49:   at->nonew   = 0;
 50:   at->maxnz   = ati[an];
 51:   *B          = At;
 52:   return 0;
 53: }

 55: PetscErrorCode MatTranspose_SeqAIJ(Mat A, MatReuse reuse, Mat *B)
 56: {
 57:   PetscInt         i, j, anzj;
 58:   Mat              At;
 59:   Mat_SeqAIJ      *a  = (Mat_SeqAIJ *)A->data, *at;
 60:   PetscInt         an = A->cmap->N, am = A->rmap->N;
 61:   PetscInt        *ati, *atj, *atfill, *ai = a->i, *aj = a->j;
 62:   MatScalar       *ata;
 63:   const MatScalar *aa, *av;
 64:   PetscContainer   rB;
 65:   MatParentState  *rb;
 66:   PetscBool        nonzerochange = PETSC_FALSE;

 68:   if (reuse == MAT_REUSE_MATRIX) {
 69:     PetscObjectQuery((PetscObject)*B, "MatTransposeParent", (PetscObject *)&rB);
 71:     PetscContainerGetPointer(rB, (void **)&rb);
 72:     if (rb->nonzerostate != A->nonzerostate) nonzerochange = PETSC_TRUE;
 73:   }

 75:   MatSeqAIJGetArrayRead(A, &av);
 76:   aa = av;
 77:   if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX || nonzerochange) {
 78:     /* Allocate space for symbolic transpose info and work array */
 79:     PetscCalloc1(an + 1, &ati);
 80:     PetscMalloc1(ai[am], &atj);
 81:     /* Walk through aj and count ## of non-zeros in each row of A^T. */
 82:     /* Note: offset by 1 for fast conversion into csr format. */
 83:     for (i = 0; i < ai[am]; i++) ati[aj[i] + 1] += 1;
 84:     /* Form ati for csr format of A^T. */
 85:     for (i = 0; i < an; i++) ati[i + 1] += ati[i];
 86:     PetscMalloc1(ai[am], &ata);
 87:   } else {
 88:     Mat_SeqAIJ *sub_B = (Mat_SeqAIJ *)(*B)->data;
 89:     ati               = sub_B->i;
 90:     atj               = sub_B->j;
 91:     ata               = sub_B->a;
 92:     At                = *B;
 93:   }

 95:   /* Copy ati into atfill so we have locations of the next free space in atj */
 96:   PetscMalloc1(an, &atfill);
 97:   PetscArraycpy(atfill, ati, an);

 99:   /* Walk through A row-wise and mark nonzero entries of A^T. */
100:   for (i = 0; i < am; i++) {
101:     anzj = ai[i + 1] - ai[i];
102:     for (j = 0; j < anzj; j++) {
103:       atj[atfill[*aj]] = i;
104:       ata[atfill[*aj]] = *aa++;
105:       atfill[*aj++] += 1;
106:     }
107:   }
108:   PetscFree(atfill);
109:   MatSeqAIJRestoreArrayRead(A, &av);
110:   if (reuse == MAT_REUSE_MATRIX) PetscObjectStateIncrease((PetscObject)(*B));

112:   if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX || nonzerochange) {
113:     MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A), an, am, ati, atj, ata, &At);
114:     MatSetBlockSizes(At, PetscAbs(A->cmap->bs), PetscAbs(A->rmap->bs));
115:     MatSetType(At, ((PetscObject)A)->type_name);
116:     at          = (Mat_SeqAIJ *)(At->data);
117:     at->free_a  = PETSC_TRUE;
118:     at->free_ij = PETSC_TRUE;
119:     at->nonew   = 0;
120:     at->maxnz   = ati[an];
121:   }

123:   if (reuse == MAT_INITIAL_MATRIX || (reuse == MAT_REUSE_MATRIX && !nonzerochange)) {
124:     *B = At;
125:   } else if (nonzerochange) {
126:     MatHeaderMerge(*B, &At);
127:     MatTransposeSetPrecursor(A, *B);
128:   } else if (reuse == MAT_INPLACE_MATRIX) {
129:     MatHeaderMerge(A, &At);
130:   }
131:   return 0;
132: }

134: /*
135:    Get symbolic matrix structure of a submatrix of A, A[rstart:rend,:],
136: */
137: PetscErrorCode MatGetSymbolicTransposeReduced_SeqAIJ(Mat A, PetscInt rstart, PetscInt rend, PetscInt *Ati[], PetscInt *Atj[])
138: {
139:   PetscInt    i, j, anzj;
140:   Mat_SeqAIJ *a  = (Mat_SeqAIJ *)A->data;
141:   PetscInt    an = A->cmap->N;
142:   PetscInt   *ati, *atj, *atfill, *ai = a->i, *aj = a->j, am = ai[rend] - ai[rstart];

144:   PetscLogEventBegin(MAT_Getsymtransreduced, A, 0, 0, 0);

146:   /* Allocate space for symbolic transpose info and work array */
147:   PetscCalloc1(an + 1, &ati);
148:   PetscMalloc1(am + 1, &atj);

150:   /* Walk through aj and count ## of non-zeros in each row of A^T. */
151:   /* Note: offset by 1 for fast conversion into csr format. */
152:   for (i = ai[rstart]; i < ai[rend]; i++) ati[aj[i] + 1] += 1;
153:   /* Form ati for csr format of A^T. */
154:   for (i = 0; i < an; i++) ati[i + 1] += ati[i];

156:   /* Copy ati into atfill so we have locations of the next free space in atj */
157:   PetscMalloc1(an + 1, &atfill);
158:   PetscArraycpy(atfill, ati, an);

160:   /* Walk through A row-wise and mark nonzero entries of A^T. */
161:   aj = aj + ai[rstart];
162:   for (i = rstart; i < rend; i++) {
163:     anzj = ai[i + 1] - ai[i];
164:     for (j = 0; j < anzj; j++) {
165:       atj[atfill[*aj]] = i - rstart;
166:       atfill[*aj++] += 1;
167:     }
168:   }
169:   PetscFree(atfill);
170:   *Ati = ati;
171:   *Atj = atj;

173:   PetscLogEventEnd(MAT_Getsymtransreduced, A, 0, 0, 0);
174:   return 0;
175: }

177: /*
178:     Returns the i and j arrays for a symbolic transpose, this is used internally within SeqAIJ code when the full
179:     symbolic matrix (which can be obtained with MatTransposeSymbolic() is not needed. MatRestoreSymbolicTranspose_SeqAIJ() should be used to free the arrays.
180: */
181: PetscErrorCode MatGetSymbolicTranspose_SeqAIJ(Mat A, PetscInt *Ati[], PetscInt *Atj[])
182: {
183:   MatGetSymbolicTransposeReduced_SeqAIJ(A, 0, A->rmap->N, Ati, Atj);
184:   return 0;
185: }

187: PetscErrorCode MatRestoreSymbolicTranspose_SeqAIJ(Mat A, PetscInt *ati[], PetscInt *atj[])
188: {
189:   PetscFree(*ati);
190:   PetscFree(*atj);
191:   return 0;
192: }