xref: /petsc/src/mat/impls/aij/seq/matptap.c (revision 9ff76d2c1165698255fb7c007af3c9986d2cdf39) 
1 
2 /*
3   Defines projective product routines where A is a SeqAIJ matrix
4           C = P^T * A * P
5 */
6 
7 #include <../src/mat/impls/aij/seq/aij.h>   /*I "petscmat.h" I*/
8 #include <../src/mat/utils/freespace.h>
9 #include <petscbt.h>
10 
11 #undef __FUNCT__
12 #define __FUNCT__ "MatPtAPSymbolic_SeqAIJ"
13 PetscErrorCode MatPtAPSymbolic_SeqAIJ(Mat A,Mat P,PetscReal fill,Mat *C)
14 {
15   PetscErrorCode ierr;
16 
17   PetscFunctionBegin;
18   if (!P->ops->ptapsymbolic_seqaij) SETERRQ2(((PetscObject)A)->comm,PETSC_ERR_SUP,"Not implemented for A=%s and P=%s",((PetscObject)A)->type_name,((PetscObject)P)->type_name);
19   ierr = (*P->ops->ptapsymbolic_seqaij)(A,P,fill,C);CHKERRQ(ierr);
20   PetscFunctionReturn(0);
21 }
22 
23 #undef __FUNCT__
24 #define __FUNCT__ "MatPtAPNumeric_SeqAIJ"
25 PetscErrorCode MatPtAPNumeric_SeqAIJ(Mat A,Mat P,Mat C)
26 {
27   PetscErrorCode ierr;
28 
29   PetscFunctionBegin;
30   if (!P->ops->ptapnumeric_seqaij) SETERRQ2(((PetscObject)A)->comm,PETSC_ERR_SUP,"Not implemented for A=%s and P=%s",((PetscObject)A)->type_name,((PetscObject)P)->type_name);
31   ierr = (*P->ops->ptapnumeric_seqaij)(A,P,C);CHKERRQ(ierr);
32   PetscFunctionReturn(0);
33 }
34 
35 #undef __FUNCT__
36 #define __FUNCT__ "PetscContainerDestroy_Mat_PtAP"
37 PetscErrorCode PetscContainerDestroy_Mat_PtAP(void *ptr)
38 {
39   PetscErrorCode ierr;
40   Mat_PtAP       *ptap=(Mat_PtAP*)ptr;
41 
42   PetscFunctionBegin;
43   ierr = PetscFree(ptap->apa);CHKERRQ(ierr);
44   ierr = MatDestroy(&ptap->AP);CHKERRQ(ierr);
45   ierr = PetscFree(ptap);CHKERRQ(ierr);
46   PetscFunctionReturn(0);
47 }
48 
49 #undef __FUNCT__
50 #define __FUNCT__ "MatDestroy_SeqAIJ_PtAP"
51 PetscErrorCode MatDestroy_SeqAIJ_PtAP(Mat A)
52 {
53   PetscErrorCode ierr;
54   PetscContainer container;
55   Mat_PtAP       *ptap=PETSC_NULL;
56 
57   PetscFunctionBegin;
58   ierr = PetscObjectQuery((PetscObject)A,"Mat_PtAP",(PetscObject *)&container);CHKERRQ(ierr);
59   if (!container) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Container does not exit");
60   ierr = PetscContainerGetPointer(container,(void **)&ptap);CHKERRQ(ierr);
61   A->ops->destroy = ptap->destroy;
62   if (A->ops->destroy) {
63     ierr = (*A->ops->destroy)(A);CHKERRQ(ierr);
64   }
65   ierr = PetscObjectCompose((PetscObject)A,"Mat_PtAP",0);CHKERRQ(ierr);
66   PetscFunctionReturn(0);
67 }
68 
69 #undef __FUNCT__
70 #define __FUNCT__ "MatPtAPSymbolic_SeqAIJ_SeqAIJ"
71 PetscErrorCode MatPtAPSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat P,PetscReal fill,Mat *C)
72 {
73   PetscErrorCode     ierr;
74   PetscFreeSpaceList free_space=PETSC_NULL,current_space=PETSC_NULL;
75   Mat_SeqAIJ         *a = (Mat_SeqAIJ*)A->data,*p = (Mat_SeqAIJ*)P->data,*c;
76   PetscInt           *pti,*ptj,*ptJ,*ai=a->i,*aj=a->j,*ajj,*pi=p->i,*pj=p->j,*pjj;
77   PetscInt           *ci,*cj,*ptadenserow,*ptasparserow,*ptaj,nspacedouble=0;
78   PetscInt           an=A->cmap->N,am=A->rmap->N,pn=P->cmap->N;
79   PetscInt           i,j,k,ptnzi,arow,anzj,ptanzi,prow,pnzj,cnzi,nlnk,*lnk;
80   MatScalar          *ca;
81   PetscBT            lnkbt;
82   Mat_PtAP           *ptap;
83   PetscContainer     container;
84 
85   PetscFunctionBegin;
86   /* Get ij structure of P^T */
87   ierr = MatGetSymbolicTranspose_SeqAIJ(P,&pti,&ptj);CHKERRQ(ierr);
88   ptJ=ptj;
89 
90   /* Allocate ci array, arrays for fill computation and */
91   /* free space for accumulating nonzero column info */
92   ierr = PetscMalloc((pn+1)*sizeof(PetscInt),&ci);CHKERRQ(ierr);
93   ci[0] = 0;
94 
95   ierr = PetscMalloc((2*an+1)*sizeof(PetscInt),&ptadenserow);CHKERRQ(ierr);
96   ierr = PetscMemzero(ptadenserow,(2*an+1)*sizeof(PetscInt));CHKERRQ(ierr);
97   ptasparserow = ptadenserow  + an;
98 
99   /* create and initialize a linked list */
100   nlnk = pn+1;
101   ierr = PetscLLCreate(pn,pn,nlnk,lnk,lnkbt);CHKERRQ(ierr);
102 
103   /* Set initial free space to be fill*nnz(A). */
104   /* This should be reasonable if sparsity of PtAP is similar to that of A. */
105   ierr          = PetscFreeSpaceGet((PetscInt)(fill*ai[am]),&free_space);
106   current_space = free_space;
107 
108   /* Determine symbolic info for each row of C: */
109   for (i=0;i<pn;i++) {
110     ptnzi  = pti[i+1] - pti[i];
111     ptanzi = 0;
112     /* Determine symbolic row of PtA: */
113     for (j=0;j<ptnzi;j++) {
114       arow = *ptJ++;
115       anzj = ai[arow+1] - ai[arow];
116       ajj  = aj + ai[arow];
117       for (k=0;k<anzj;k++) {
118         if (!ptadenserow[ajj[k]]) {
119           ptadenserow[ajj[k]]    = -1;
120           ptasparserow[ptanzi++] = ajj[k];
121         }
122       }
123     }
124     /* Using symbolic info for row of PtA, determine symbolic info for row of C: */
125     ptaj = ptasparserow;
126     cnzi   = 0;
127     for (j=0;j<ptanzi;j++) {
128       prow = *ptaj++;
129       pnzj = pi[prow+1] - pi[prow];
130       pjj  = pj + pi[prow];
131       /* add non-zero cols of P into the sorted linked list lnk */
132       ierr = PetscLLAdd(pnzj,pjj,pn,nlnk,lnk,lnkbt);CHKERRQ(ierr);
133       cnzi += nlnk;
134     }
135 
136     /* If free space is not available, make more free space */
137     /* Double the amount of total space in the list */
138     if (current_space->local_remaining<cnzi) {
139       ierr = PetscFreeSpaceGet(cnzi+current_space->total_array_size,&current_space);CHKERRQ(ierr);
140       nspacedouble++;
141     }
142 
143     /* Copy data into free space, and zero out denserows */
144     ierr = PetscLLClean(pn,pn,cnzi,lnk,current_space->array,lnkbt);CHKERRQ(ierr);
145     current_space->array           += cnzi;
146     current_space->local_used      += cnzi;
147     current_space->local_remaining -= cnzi;
148 
149     for (j=0;j<ptanzi;j++) {
150       ptadenserow[ptasparserow[j]] = 0;
151     }
152     /* Aside: Perhaps we should save the pta info for the numerical factorization. */
153     /*        For now, we will recompute what is needed. */
154     ci[i+1] = ci[i] + cnzi;
155   }
156   /* nnz is now stored in ci[ptm], column indices are in the list of free space */
157   /* Allocate space for cj, initialize cj, and */
158   /* destroy list of free space and other temporary array(s) */
159   ierr = PetscMalloc((ci[pn]+1)*sizeof(PetscInt),&cj);CHKERRQ(ierr);
160   ierr = PetscFreeSpaceContiguous(&free_space,cj);CHKERRQ(ierr);
161   ierr = PetscFree(ptadenserow);CHKERRQ(ierr);
162   ierr = PetscLLDestroy(lnk,lnkbt);CHKERRQ(ierr);
163 
164   /* Allocate space for ca */
165   ierr = PetscMalloc((ci[pn]+1)*sizeof(MatScalar),&ca);CHKERRQ(ierr);
166   ierr = PetscMemzero(ca,(ci[pn]+1)*sizeof(MatScalar));CHKERRQ(ierr);
167 
168   /* put together the new matrix */
169   ierr = MatCreateSeqAIJWithArrays(((PetscObject)A)->comm,pn,pn,ci,cj,ca,C);CHKERRQ(ierr);
170 
171   /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
172   /* Since these are PETSc arrays, change flags to free them as necessary. */
173   c          = (Mat_SeqAIJ *)(*C)->data;
174   c->free_a  = PETSC_TRUE;
175   c->free_ij = PETSC_TRUE;
176   c->nonew   = 0;
177 
178   /* create a supporting struct for reuse by MatPtAPNumeric() */
179   ierr = PetscNew(Mat_PtAP,&ptap);CHKERRQ(ierr);
180 
181   /* attach the supporting struct to C */
182   ierr = PetscContainerCreate(PETSC_COMM_SELF,&container);CHKERRQ(ierr);
183   ierr = PetscContainerSetPointer(container,ptap);CHKERRQ(ierr);
184   ierr = PetscContainerSetUserDestroy(container,PetscContainerDestroy_Mat_PtAP);CHKERRQ(ierr);
185   ierr = PetscObjectCompose((PetscObject)(*C),"Mat_PtAP",(PetscObject)container);CHKERRQ(ierr);
186   ierr = PetscContainerDestroy(&container);CHKERRQ(ierr);
187 
188   ptap->destroy      = (*C)->ops->destroy;
189   (*C)->ops->destroy = MatDestroy_SeqAIJ_PtAP;
190 
191   /* Allocate temporary array for storage of one row of A*P */
192   ierr = PetscMalloc((pn+1)*sizeof(PetscScalar),&ptap->apa);CHKERRQ(ierr);
193 
194   /* Get structure of A*P */
195   ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(A,P,fill,&ptap->AP);CHKERRQ(ierr);
196 
197   /* Clean up. */
198   ierr = MatRestoreSymbolicTranspose_SeqAIJ(P,&pti,&ptj);CHKERRQ(ierr);
199 #if defined(PETSC_USE_INFO)
200   if (ci[pn] != 0) {
201     PetscReal afill = ((PetscReal)ci[pn])/ai[am];
202     if (afill < 1.0) afill = 1.0;
203     ierr = PetscInfo3((*C),"Reallocs %D; Fill ratio: given %G needed %G.\n",nspacedouble,fill,afill);CHKERRQ(ierr);
204     ierr = PetscInfo1((*C),"Use MatPtAP(A,P,MatReuse,%G,&C) for best performance.\n",afill);CHKERRQ(ierr);
205   } else {
206     ierr = PetscInfo((*C),"Empty matrix product\n");CHKERRQ(ierr);
207   }
208 #endif
209   PetscFunctionReturn(0);
210 }
211 
212 #undef __FUNCT__
213 #define __FUNCT__ "MatPtAPNumeric_SeqAIJ_SeqAIJ"
214 PetscErrorCode MatPtAPNumeric_SeqAIJ_SeqAIJ(Mat A,Mat P,Mat C)
215 {
216   PetscErrorCode ierr;
217   Mat_SeqAIJ     *a  = (Mat_SeqAIJ *) A->data;
218   Mat_SeqAIJ     *p  = (Mat_SeqAIJ *) P->data,*ap;
219   Mat_SeqAIJ     *c  = (Mat_SeqAIJ *) C->data;
220   PetscInt       *ai=a->i,*aj=a->j,*apj,*pi=p->i,*pj=p->j,*pcol;
221   PetscInt       *ci=c->i,*cj=c->j,*cjj,cnz;
222   PetscInt       am=A->rmap->N,cn=C->cmap->N,cm=C->rmap->N;
223   PetscInt       i,j,k,anz,apnz,pnz,prow,crow,apcol,nextap;
224   MatScalar      *aa=a->a,*apa,*pa=p->a,*pval,*ca=c->a,*caj;
225   PetscBool      sparse_axpy=PETSC_FALSE;
226   Mat_PtAP       *ptap;
227   PetscContainer container;
228 
229   PetscFunctionBegin;
230   ierr = PetscOptionsGetBool(PETSC_NULL,"-matptap_spaxpy",&sparse_axpy,PETSC_NULL);CHKERRQ(ierr);
231   ierr = PetscObjectQuery((PetscObject)C,"Mat_PtAP",(PetscObject *)&container);CHKERRQ(ierr);
232   if (!container) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Container does not exit");
233   ierr  = PetscContainerGetPointer(container,(void **)&ptap);CHKERRQ(ierr);
234 
235   /* Get temporary array for storage of one row of A*P */
236   apa = ptap->apa;
237   ierr = PetscMemzero(apa,cn*sizeof(MatScalar));CHKERRQ(ierr);
238   ap   = (Mat_SeqAIJ *)(ptap->AP)->data;
239 
240   /* Clear old values in C */
241   ierr = PetscMemzero(ca,ci[cm]*sizeof(MatScalar));CHKERRQ(ierr);
242 
243   for (i=0;i<am;i++) {
244     /* Form sparse row of AP[i,:] = A[i,:]*P */
245     anz  = ai[i+1] - ai[i];
246     apnz = 0;
247     for (j=0; j<anz; j++) {
248       prow = aj[j];
249       pnz  = pi[prow+1] - pi[prow];
250       pcol = pj + pi[prow];
251       pval = pa + pi[prow];
252       for (k=0; k<pnz; k++) {
253         apa[pcol[k]] += aa[j]*pval[k];
254       }
255       ierr = PetscLogFlops(2.0*pnz);CHKERRQ(ierr);
256     }
257     aj += anz; aa += anz;
258 
259     /* Compute P^T*A*P using outer product P[i,:]^T*AP[i,:]. */
260     apj  = ap->j + ap->i[i];
261     apnz = ap->i[i+1] - ap->i[i];
262     pnz  = pi[i+1] - pi[i];
263     pcol = pj + pi[i];
264     pval = pa + pi[i];
265 
266     for (j=0; j<pnz; j++) {
267       crow = pcol[j];
268       cjj  = cj + ci[crow];
269       caj  = ca + ci[crow];
270 
271       if (sparse_axpy){  /* Perform sparse axpy */
272         nextap = 0;
273         apcol = apj[nextap];
274         for (k=0; nextap<apnz; k++) {
275           if (cjj[k] == apcol) {
276             caj[k] += pval[j]*apa[apcol];
277             apcol   = apj[++nextap];
278           }
279         }
280         ierr = PetscLogFlops(2.0*apnz);CHKERRQ(ierr);
281       } else { /* Perform dense axpy */
282         cnz  = ci[crow+1] - ci[crow];
283         for (k=0; k<cnz; k++){
284           caj[k] += pval[j]*apa[cjj[k]];
285         }
286         ierr = PetscLogFlops(2.0*cnz);CHKERRQ(ierr);
287       }
288     }
289 
290     /* Zero the current row info for A*P */
291     for (j=0; j<apnz; j++) {
292       apcol      = apj[j];
293       apa[apcol] = 0.;
294     }
295   }
296 
297   /* Assemble the final matrix and clean up */
298   ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
299   ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
300   PetscFunctionReturn(0);
301 }
302