1 /*
2 Defines projective product routines where A is a SeqAIJ matrix
3 C = P^T * A * P
4 */
5
6 #include <../src/mat/impls/aij/seq/aij.h> /*I "petscmat.h" I*/
7 #include <../src/mat/utils/freespace.h>
8 #include <petscbt.h>
9 #include <petsctime.h>
10
11 #if defined(PETSC_HAVE_HYPRE)
12 PETSC_INTERN PetscErrorCode MatPtAPSymbolic_AIJ_AIJ_wHYPRE(Mat, Mat, PetscReal, Mat);
13 #endif
14
MatProductSymbolic_PtAP_SeqAIJ_SeqAIJ(Mat C)15 PetscErrorCode MatProductSymbolic_PtAP_SeqAIJ_SeqAIJ(Mat C)
16 {
17 Mat_Product *product = C->product;
18 Mat A = product->A, P = product->B;
19 MatProductAlgorithm alg = product->alg;
20 PetscReal fill = product->fill;
21 PetscBool flg;
22 Mat Pt;
23
24 PetscFunctionBegin;
25 /* "scalable" */
26 PetscCall(PetscStrcmp(alg, "scalable", &flg));
27 if (flg) {
28 PetscCall(MatPtAPSymbolic_SeqAIJ_SeqAIJ_SparseAxpy(A, P, fill, C));
29 C->ops->productnumeric = MatProductNumeric_PtAP;
30 PetscFunctionReturn(PETSC_SUCCESS);
31 }
32
33 /* "rap" */
34 PetscCall(PetscStrcmp(alg, "rap", &flg));
35 if (flg) {
36 MatProductCtx_MatTransMatMult *atb;
37
38 PetscCall(PetscNew(&atb));
39 PetscCall(MatTranspose(P, MAT_INITIAL_MATRIX, &Pt));
40 PetscCall(MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ(Pt, A, P, fill, C));
41
42 atb->At = Pt;
43 atb->data = C->product->data;
44 atb->destroy = C->product->destroy;
45 C->product->data = atb;
46 C->product->destroy = MatProductCtxDestroy_SeqAIJ_MatTransMatMult;
47 C->ops->ptapnumeric = MatPtAPNumeric_SeqAIJ_SeqAIJ;
48 C->ops->productnumeric = MatProductNumeric_PtAP;
49 PetscFunctionReturn(PETSC_SUCCESS);
50 }
51
52 /* hypre */
53 #if defined(PETSC_HAVE_HYPRE)
54 PetscCall(PetscStrcmp(alg, "hypre", &flg));
55 if (flg) {
56 PetscCall(MatPtAPSymbolic_AIJ_AIJ_wHYPRE(A, P, fill, C));
57 PetscFunctionReturn(PETSC_SUCCESS);
58 }
59 #endif
60
61 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "MatProductType is not supported");
62 }
63
MatPtAPSymbolic_SeqAIJ_SeqAIJ_SparseAxpy(Mat A,Mat P,PetscReal fill,Mat C)64 PetscErrorCode MatPtAPSymbolic_SeqAIJ_SeqAIJ_SparseAxpy(Mat A, Mat P, PetscReal fill, Mat C)
65 {
66 PetscFreeSpaceList free_space = NULL, current_space = NULL;
67 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *p = (Mat_SeqAIJ *)P->data, *c;
68 PetscInt *pti, *ptj, *ptJ, *ai = a->i, *aj = a->j, *ajj, *pi = p->i, *pj = p->j, *pjj;
69 PetscInt *ci, *cj, *ptadenserow, *ptasparserow, *ptaj, nspacedouble = 0;
70 PetscInt an = A->cmap->N, am = A->rmap->N, pn = P->cmap->N, pm = P->rmap->N;
71 PetscInt i, j, k, ptnzi, arow, anzj, ptanzi, prow, pnzj, cnzi, nlnk, *lnk;
72 MatScalar *ca;
73 PetscBT lnkbt;
74 PetscReal afill;
75
76 PetscFunctionBegin;
77 /* Get ij structure of P^T */
78 PetscCall(MatGetSymbolicTranspose_SeqAIJ(P, &pti, &ptj));
79 ptJ = ptj;
80
81 /* Allocate ci array, arrays for fill computation and */
82 /* free space for accumulating nonzero column info */
83 PetscCall(PetscMalloc1(pn + 1, &ci));
84 ci[0] = 0;
85
86 PetscCall(PetscCalloc1(2 * an + 1, &ptadenserow));
87 ptasparserow = ptadenserow + an;
88
89 /* create and initialize a linked list */
90 nlnk = pn + 1;
91 PetscCall(PetscLLCreate(pn, pn, nlnk, lnk, lnkbt));
92
93 /* Set initial free space to be fill*(nnz(A)+ nnz(P)) */
94 PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(fill, PetscIntSumTruncate(ai[am], pi[pm])), &free_space));
95 current_space = free_space;
96
97 /* Determine symbolic info for each row of C: */
98 for (i = 0; i < pn; i++) {
99 ptnzi = pti[i + 1] - pti[i];
100 ptanzi = 0;
101 /* Determine symbolic row of PtA: */
102 for (j = 0; j < ptnzi; j++) {
103 arow = *ptJ++;
104 anzj = ai[arow + 1] - ai[arow];
105 ajj = aj + ai[arow];
106 for (k = 0; k < anzj; k++) {
107 if (!ptadenserow[ajj[k]]) {
108 ptadenserow[ajj[k]] = -1;
109 ptasparserow[ptanzi++] = ajj[k];
110 }
111 }
112 }
113 /* Using symbolic info for row of PtA, determine symbolic info for row of C: */
114 ptaj = ptasparserow;
115 cnzi = 0;
116 for (j = 0; j < ptanzi; j++) {
117 prow = *ptaj++;
118 pnzj = pi[prow + 1] - pi[prow];
119 pjj = pj + pi[prow];
120 /* add non-zero cols of P into the sorted linked list lnk */
121 PetscCall(PetscLLAddSorted(pnzj, pjj, pn, &nlnk, lnk, lnkbt));
122 cnzi += nlnk;
123 }
124
125 /* If free space is not available, make more free space */
126 /* Double the amount of total space in the list */
127 if (current_space->local_remaining < cnzi) {
128 PetscCall(PetscFreeSpaceGet(PetscIntSumTruncate(cnzi, current_space->total_array_size), ¤t_space));
129 nspacedouble++;
130 }
131
132 /* Copy data into free space, and zero out denserows */
133 PetscCall(PetscLLClean(pn, pn, cnzi, lnk, current_space->array, lnkbt));
134
135 current_space->array += cnzi;
136 current_space->local_used += cnzi;
137 current_space->local_remaining -= cnzi;
138
139 for (j = 0; j < ptanzi; j++) ptadenserow[ptasparserow[j]] = 0;
140
141 /* Aside: Perhaps we should save the pta info for the numerical factorization. */
142 /* For now, we will recompute what is needed. */
143 ci[i + 1] = ci[i] + cnzi;
144 }
145 /* nnz is now stored in ci[ptm], column indices are in the list of free space */
146 /* Allocate space for cj, initialize cj, and */
147 /* destroy list of free space and other temporary array(s) */
148 PetscCall(PetscMalloc1(ci[pn], &cj));
149 PetscCall(PetscFreeSpaceContiguous(&free_space, cj));
150 PetscCall(PetscFree(ptadenserow));
151 PetscCall(PetscLLDestroy(lnk, lnkbt));
152
153 PetscCall(PetscCalloc1(ci[pn], &ca));
154
155 /* put together the new matrix */
156 PetscCall(MatSetSeqAIJWithArrays_private(PetscObjectComm((PetscObject)A), pn, pn, ci, cj, ca, ((PetscObject)A)->type_name, C));
157 PetscCall(MatSetBlockSizes(C, P->cmap->bs, P->cmap->bs));
158
159 /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
160 /* Since these are PETSc arrays, change flags to free them as necessary. */
161 c = (Mat_SeqAIJ *)((C)->data);
162 c->free_a = PETSC_TRUE;
163 c->free_ij = PETSC_TRUE;
164 c->nonew = 0;
165
166 C->ops->ptapnumeric = MatPtAPNumeric_SeqAIJ_SeqAIJ_SparseAxpy;
167
168 /* set MatInfo */
169 afill = (PetscReal)ci[pn] / (ai[am] + pi[pm] + 1.e-5);
170 if (afill < 1.0) afill = 1.0;
171 C->info.mallocs = nspacedouble;
172 C->info.fill_ratio_given = fill;
173 C->info.fill_ratio_needed = afill;
174
175 /* Clean up. */
176 PetscCall(MatRestoreSymbolicTranspose_SeqAIJ(P, &pti, &ptj));
177 #if defined(PETSC_USE_INFO)
178 if (ci[pn] != 0) {
179 PetscCall(PetscInfo(C, "Reallocs %" PetscInt_FMT "; Fill ratio: given %g needed %g.\n", nspacedouble, (double)fill, (double)afill));
180 PetscCall(PetscInfo(C, "Use MatPtAP(A,P,MatReuse,%g,&C) for best performance.\n", (double)afill));
181 } else {
182 PetscCall(PetscInfo(C, "Empty matrix product\n"));
183 }
184 #endif
185 PetscFunctionReturn(PETSC_SUCCESS);
186 }
187
MatPtAPNumeric_SeqAIJ_SeqAIJ_SparseAxpy(Mat A,Mat P,Mat C)188 PetscErrorCode MatPtAPNumeric_SeqAIJ_SeqAIJ_SparseAxpy(Mat A, Mat P, Mat C)
189 {
190 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
191 Mat_SeqAIJ *p = (Mat_SeqAIJ *)P->data;
192 Mat_SeqAIJ *c = (Mat_SeqAIJ *)C->data;
193 PetscInt *ai = a->i, *aj = a->j, *apj, *apjdense, *pi = p->i, *pj = p->j, *pJ = p->j, *pjj;
194 PetscInt *ci = c->i, *cj = c->j, *cjj;
195 PetscInt am = A->rmap->N, cn = C->cmap->N, cm = C->rmap->N;
196 PetscInt i, j, k, anzi, pnzi, apnzj, nextap, pnzj, prow, crow;
197 MatScalar *aa, *apa, *pa, *pA, *paj, *ca, *caj;
198
199 PetscFunctionBegin;
200 /* Allocate temporary array for storage of one row of A*P (cn: non-scalable) */
201 PetscCall(PetscCalloc2(cn, &apa, cn, &apjdense));
202 PetscCall(PetscMalloc1(cn, &apj));
203 /* trigger CPU copies if needed and flag CPU mask for C */
204 #if defined(PETSC_HAVE_DEVICE)
205 {
206 const PetscScalar *unused;
207 PetscCall(MatSeqAIJGetArrayRead(A, &unused));
208 PetscCall(MatSeqAIJRestoreArrayRead(A, &unused));
209 PetscCall(MatSeqAIJGetArrayRead(P, &unused));
210 PetscCall(MatSeqAIJRestoreArrayRead(P, &unused));
211 if (C->offloadmask != PETSC_OFFLOAD_UNALLOCATED) C->offloadmask = PETSC_OFFLOAD_CPU;
212 }
213 #endif
214 aa = a->a;
215 pa = p->a;
216 pA = p->a;
217 ca = c->a;
218
219 /* Clear old values in C */
220 PetscCall(PetscArrayzero(ca, ci[cm]));
221
222 for (i = 0; i < am; i++) {
223 /* Form sparse row of A*P */
224 anzi = ai[i + 1] - ai[i];
225 apnzj = 0;
226 for (j = 0; j < anzi; j++) {
227 prow = *aj++;
228 pnzj = pi[prow + 1] - pi[prow];
229 pjj = pj + pi[prow];
230 paj = pa + pi[prow];
231 for (k = 0; k < pnzj; k++) {
232 if (!apjdense[pjj[k]]) {
233 apjdense[pjj[k]] = -1;
234 apj[apnzj++] = pjj[k];
235 }
236 apa[pjj[k]] += (*aa) * paj[k];
237 }
238 PetscCall(PetscLogFlops(2.0 * pnzj));
239 aa++;
240 }
241
242 /* Sort the j index array for quick sparse axpy. */
243 /* Note: a array does not need sorting as it is in dense storage locations. */
244 PetscCall(PetscSortInt(apnzj, apj));
245
246 /* Compute P^T*A*P using outer product (P^T)[:,j]*(A*P)[j,:]. */
247 pnzi = pi[i + 1] - pi[i];
248 for (j = 0; j < pnzi; j++) {
249 nextap = 0;
250 crow = *pJ++;
251 cjj = cj + ci[crow];
252 caj = ca + ci[crow];
253 /* Perform sparse axpy operation. Note cjj includes apj. */
254 for (k = 0; nextap < apnzj; k++) {
255 PetscAssert(k < ci[crow + 1] - ci[crow], PETSC_COMM_SELF, PETSC_ERR_PLIB, "k too large k %" PetscInt_FMT ", crow %" PetscInt_FMT, k, crow);
256 if (cjj[k] == apj[nextap]) caj[k] += (*pA) * apa[apj[nextap++]];
257 }
258 PetscCall(PetscLogFlops(2.0 * apnzj));
259 pA++;
260 }
261
262 /* Zero the current row info for A*P */
263 for (j = 0; j < apnzj; j++) {
264 apa[apj[j]] = 0.;
265 apjdense[apj[j]] = 0;
266 }
267 }
268
269 /* Assemble the final matrix and clean up */
270 PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
271 PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
272
273 PetscCall(PetscFree2(apa, apjdense));
274 PetscCall(PetscFree(apj));
275 PetscFunctionReturn(PETSC_SUCCESS);
276 }
277
MatPtAPNumeric_SeqAIJ_SeqAIJ(Mat A,Mat P,Mat C)278 PetscErrorCode MatPtAPNumeric_SeqAIJ_SeqAIJ(Mat A, Mat P, Mat C)
279 {
280 MatProductCtx_MatTransMatMult *atb;
281
282 PetscFunctionBegin;
283 MatCheckProduct(C, 3);
284 atb = (MatProductCtx_MatTransMatMult *)C->product->data;
285 PetscCheck(atb, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Missing data structure");
286 PetscCall(MatTranspose(P, MAT_REUSE_MATRIX, &atb->At));
287 PetscCheck(C->ops->matmultnumeric, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Missing numeric operation");
288 /* when using rap, MatMatMatMultSymbolic used a different data */
289 if (atb->data) C->product->data = atb->data;
290 PetscCall((*C->ops->matmatmultnumeric)(atb->At, A, P, C));
291 C->product->data = atb;
292 PetscFunctionReturn(PETSC_SUCCESS);
293 }
294