1 #include <petscmat.h>              /*I  "mat.h"  I*/
2 
3 PETSC_INTERN PetscErrorCode MatCreateADA(Mat,Vec, Vec, Mat*);
4 
5 typedef struct {
6   Mat      A;
7   Vec      D1;
8   Vec      D2;
9   Vec      W;
10   Vec      W2;
11   Vec      ADADiag;
12   PetscInt GotDiag;
13 } _p_TaoMatADACtx;
14 typedef  _p_TaoMatADACtx* TaoMatADACtx;
15 
MatMult_ADA(Mat mat,Vec a,Vec y)16 static PetscErrorCode MatMult_ADA(Mat mat,Vec a,Vec y)
17 {
18   TaoMatADACtx   ctx;
19   PetscReal      one = 1.0;
20   PetscErrorCode ierr;
21 
22   PetscFunctionBegin;
23   ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr);
24   ierr = MatMult(ctx->A,a,ctx->W);CHKERRQ(ierr);
25   if (ctx->D1){
26     ierr = VecPointwiseMult(ctx->W,ctx->D1,ctx->W);CHKERRQ(ierr);
27   }
28   ierr = MatMultTranspose(ctx->A,ctx->W,y);CHKERRQ(ierr);
29   if (ctx->D2){
30     ierr = VecPointwiseMult(ctx->W2, ctx->D2, a);CHKERRQ(ierr);
31     ierr = VecAXPY(y, one, ctx->W2);CHKERRQ(ierr);
32   }
33   PetscFunctionReturn(0);
34 }
35 
MatMultTranspose_ADA(Mat mat,Vec a,Vec y)36 static PetscErrorCode MatMultTranspose_ADA(Mat mat,Vec a,Vec y)
37 {
38   PetscErrorCode ierr;
39 
40   PetscFunctionBegin;
41   ierr = MatMult_ADA(mat,a,y);CHKERRQ(ierr);
42   PetscFunctionReturn(0);
43 }
44 
MatDiagonalSet_ADA(Mat M,Vec D,InsertMode mode)45 static PetscErrorCode MatDiagonalSet_ADA(Mat M,Vec D, InsertMode mode)
46 {
47   TaoMatADACtx   ctx;
48   PetscReal      zero=0.0,one = 1.0;
49   PetscErrorCode ierr;
50 
51   PetscFunctionBegin;
52   if (mode == INSERT_VALUES) SETERRQ(PetscObjectComm((PetscObject)M),PETSC_ERR_SUP,"Cannot insert diagonal entries of this matrix type, can only add");
53   ierr = MatShellGetContext(M,(void **)&ctx);CHKERRQ(ierr);
54   if (!ctx->D2){
55     ierr = VecDuplicate(D,&ctx->D2);CHKERRQ(ierr);
56     ierr = VecSet(ctx->D2, zero);CHKERRQ(ierr);
57   }
58   ierr = VecAXPY(ctx->D2, one, D);CHKERRQ(ierr);
59   PetscFunctionReturn(0);
60 }
61 
MatDestroy_ADA(Mat mat)62 static PetscErrorCode MatDestroy_ADA(Mat mat)
63 {
64   PetscErrorCode ierr;
65   TaoMatADACtx   ctx;
66 
67   PetscFunctionBegin;
68   ierr=MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr);
69   ierr=VecDestroy(&ctx->W);CHKERRQ(ierr);
70   ierr=VecDestroy(&ctx->W2);CHKERRQ(ierr);
71   ierr=VecDestroy(&ctx->ADADiag);CHKERRQ(ierr);
72   ierr=MatDestroy(&ctx->A);CHKERRQ(ierr);
73   ierr=VecDestroy(&ctx->D1);CHKERRQ(ierr);
74   ierr=VecDestroy(&ctx->D2);CHKERRQ(ierr);
75   ierr = PetscFree(ctx);CHKERRQ(ierr);
76   PetscFunctionReturn(0);
77 }
78 
MatView_ADA(Mat mat,PetscViewer viewer)79 static PetscErrorCode MatView_ADA(Mat mat,PetscViewer viewer)
80 {
81   PetscFunctionBegin;
82   PetscFunctionReturn(0);
83 }
84 
MatShift_ADA(Mat Y,PetscReal a)85 static PetscErrorCode MatShift_ADA(Mat Y, PetscReal a)
86 {
87   PetscErrorCode ierr;
88   TaoMatADACtx   ctx;
89 
90   PetscFunctionBegin;
91   ierr = MatShellGetContext(Y,(void **)&ctx);CHKERRQ(ierr);
92   ierr = VecShift(ctx->D2,a);CHKERRQ(ierr);
93   PetscFunctionReturn(0);
94 }
95 
MatDuplicate_ADA(Mat mat,MatDuplicateOption op,Mat * M)96 static PetscErrorCode MatDuplicate_ADA(Mat mat,MatDuplicateOption op,Mat *M)
97 {
98   PetscErrorCode    ierr;
99   TaoMatADACtx      ctx;
100   Mat               A2;
101   Vec               D1b=NULL,D2b;
102 
103   PetscFunctionBegin;
104   ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr);
105   ierr = MatDuplicate(ctx->A,op,&A2);CHKERRQ(ierr);
106   if (ctx->D1){
107     ierr = VecDuplicate(ctx->D1,&D1b);CHKERRQ(ierr);
108     ierr = VecCopy(ctx->D1,D1b);CHKERRQ(ierr);
109   }
110   ierr = VecDuplicate(ctx->D2,&D2b);CHKERRQ(ierr);
111   ierr = VecCopy(ctx->D2,D2b);CHKERRQ(ierr);
112   ierr = MatCreateADA(A2,D1b,D2b,M);CHKERRQ(ierr);
113   if (ctx->D1){
114     ierr = PetscObjectDereference((PetscObject)D1b);CHKERRQ(ierr);
115   }
116   ierr = PetscObjectDereference((PetscObject)D2b);CHKERRQ(ierr);
117   ierr = PetscObjectDereference((PetscObject)A2);CHKERRQ(ierr);
118   PetscFunctionReturn(0);
119 }
120 
MatEqual_ADA(Mat A,Mat B,PetscBool * flg)121 static PetscErrorCode MatEqual_ADA(Mat A,Mat B,PetscBool *flg)
122 {
123   PetscErrorCode ierr;
124   TaoMatADACtx   ctx1,ctx2;
125 
126   PetscFunctionBegin;
127   ierr = MatShellGetContext(A,(void **)&ctx1);CHKERRQ(ierr);
128   ierr = MatShellGetContext(B,(void **)&ctx2);CHKERRQ(ierr);
129   ierr = VecEqual(ctx1->D2,ctx2->D2,flg);CHKERRQ(ierr);
130   if (*flg==PETSC_TRUE){
131     ierr = VecEqual(ctx1->D1,ctx2->D1,flg);CHKERRQ(ierr);
132   }
133   if (*flg==PETSC_TRUE){
134     ierr = MatEqual(ctx1->A,ctx2->A,flg);CHKERRQ(ierr);
135   }
136   PetscFunctionReturn(0);
137 }
138 
MatScale_ADA(Mat mat,PetscReal a)139 static PetscErrorCode MatScale_ADA(Mat mat, PetscReal a)
140 {
141   PetscErrorCode ierr;
142   TaoMatADACtx   ctx;
143 
144   PetscFunctionBegin;
145   ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr);
146   ierr = VecScale(ctx->D1,a);CHKERRQ(ierr);
147   if (ctx->D2){
148     ierr = VecScale(ctx->D2,a);CHKERRQ(ierr);
149   }
150   PetscFunctionReturn(0);
151 }
152 
MatTranspose_ADA(Mat mat,MatReuse reuse,Mat * B)153 static PetscErrorCode MatTranspose_ADA(Mat mat,MatReuse reuse,Mat *B)
154 {
155   PetscErrorCode ierr;
156   TaoMatADACtx   ctx;
157 
158   PetscFunctionBegin;
159   ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr);
160   if (reuse == MAT_INITIAL_MATRIX){
161     ierr = MatDuplicate(mat,MAT_COPY_VALUES,B);CHKERRQ(ierr);
162   } else if (reuse == MAT_REUSE_MATRIX){
163     ierr = MatCopy(mat,*B,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
164   } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Does not support inplace transpose");
165   PetscFunctionReturn(0);
166 }
167 
MatADAComputeDiagonal(Mat mat)168 static PetscErrorCode MatADAComputeDiagonal(Mat mat)
169 {
170   PetscErrorCode ierr;
171   PetscInt       i,m,n,low,high;
172   PetscScalar    *dtemp,*dptr;
173   TaoMatADACtx   ctx;
174 
175   PetscFunctionBegin;
176   ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr);
177   ierr = MatGetOwnershipRange(mat, &low, &high);CHKERRQ(ierr);
178   ierr = MatGetSize(mat,&m,&n);CHKERRQ(ierr);
179 
180   ierr = PetscMalloc1(n,&dtemp);CHKERRQ(ierr);
181   for (i=0; i<n; i++){
182     ierr = MatGetColumnVector(ctx->A, ctx->W, i);CHKERRQ(ierr);
183     ierr = VecPointwiseMult(ctx->W,ctx->W,ctx->W);CHKERRQ(ierr);
184     ierr = VecDotBegin(ctx->D1, ctx->W,dtemp+i);CHKERRQ(ierr);
185   }
186   for (i=0; i<n; i++){
187     ierr = VecDotEnd(ctx->D1, ctx->W,dtemp+i);CHKERRQ(ierr);
188   }
189 
190   ierr = VecGetArray(ctx->ADADiag,&dptr);CHKERRQ(ierr);
191   for (i=low; i<high; i++){
192     dptr[i-low]= dtemp[i];
193   }
194   ierr = VecRestoreArray(ctx->ADADiag,&dptr);CHKERRQ(ierr);
195   ierr = PetscFree(dtemp);CHKERRQ(ierr);
196   PetscFunctionReturn(0);
197 }
198 
MatGetDiagonal_ADA(Mat mat,Vec v)199 static PetscErrorCode MatGetDiagonal_ADA(Mat mat,Vec v)
200 {
201   PetscErrorCode  ierr;
202   PetscReal       one=1.0;
203   TaoMatADACtx    ctx;
204 
205   PetscFunctionBegin;
206   ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr);
207   ierr = MatADAComputeDiagonal(mat);CHKERRQ(ierr);
208   ierr = VecCopy(ctx->ADADiag,v);CHKERRQ(ierr);
209   if (ctx->D2){
210     ierr = VecAXPY(v, one, ctx->D2);CHKERRQ(ierr);
211   }
212   PetscFunctionReturn(0);
213 }
214 
MatCreateSubMatrix_ADA(Mat mat,IS isrow,IS iscol,MatReuse cll,Mat * newmat)215 static PetscErrorCode MatCreateSubMatrix_ADA(Mat mat,IS isrow,IS iscol,MatReuse cll, Mat *newmat)
216 {
217   PetscErrorCode    ierr;
218   PetscInt          low,high;
219   IS                ISrow;
220   Vec               D1,D2;
221   Mat               Atemp;
222   TaoMatADACtx      ctx;
223   PetscBool         isequal;
224 
225   PetscFunctionBegin;
226   ierr = ISEqual(isrow,iscol,&isequal);CHKERRQ(ierr);
227   if (!isequal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only for identical column and row indices");
228   ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr);
229 
230   ierr = MatGetOwnershipRange(ctx->A,&low,&high);CHKERRQ(ierr);
231   ierr = ISCreateStride(PetscObjectComm((PetscObject)mat),high-low,low,1,&ISrow);CHKERRQ(ierr);
232   ierr = MatCreateSubMatrix(ctx->A,ISrow,iscol,cll,&Atemp);CHKERRQ(ierr);
233   ierr = ISDestroy(&ISrow);CHKERRQ(ierr);
234 
235   if (ctx->D1){
236     ierr=VecDuplicate(ctx->D1,&D1);CHKERRQ(ierr);
237     ierr=VecCopy(ctx->D1,D1);CHKERRQ(ierr);
238   } else {
239     D1 = NULL;
240   }
241 
242   if (ctx->D2){
243     Vec D2sub;
244 
245     ierr=VecGetSubVector(ctx->D2,isrow,&D2sub);CHKERRQ(ierr);
246     ierr=VecDuplicate(D2sub,&D2);CHKERRQ(ierr);
247     ierr=VecCopy(D2sub,D2);CHKERRQ(ierr);
248     ierr=VecRestoreSubVector(ctx->D2,isrow,&D2sub);CHKERRQ(ierr);
249   } else {
250     D2 = NULL;
251   }
252 
253   ierr = MatCreateADA(Atemp,D1,D2,newmat);CHKERRQ(ierr);
254   ierr = MatShellGetContext(*newmat,(void **)&ctx);CHKERRQ(ierr);
255   ierr = PetscObjectDereference((PetscObject)Atemp);CHKERRQ(ierr);
256   if (ctx->D1){
257     ierr = PetscObjectDereference((PetscObject)D1);CHKERRQ(ierr);
258   }
259   if (ctx->D2){
260     ierr = PetscObjectDereference((PetscObject)D2);CHKERRQ(ierr);
261   }
262   PetscFunctionReturn(0);
263 }
264 
MatCreateSubMatrices_ADA(Mat A,PetscInt n,IS * irow,IS * icol,MatReuse scall,Mat ** B)265 static PetscErrorCode MatCreateSubMatrices_ADA(Mat A,PetscInt n, IS *irow,IS *icol,MatReuse scall,Mat **B)
266 {
267   PetscErrorCode ierr;
268   PetscInt       i;
269 
270   PetscFunctionBegin;
271   if (scall == MAT_INITIAL_MATRIX) {
272     ierr = PetscCalloc1(n+1,B);CHKERRQ(ierr);
273   }
274   for (i=0; i<n; i++) {
275     ierr = MatCreateSubMatrix_ADA(A,irow[i],icol[i],scall,&(*B)[i]);CHKERRQ(ierr);
276   }
277   PetscFunctionReturn(0);
278 }
279 
MatGetColumnVector_ADA(Mat mat,Vec Y,PetscInt col)280 static PetscErrorCode MatGetColumnVector_ADA(Mat mat,Vec Y, PetscInt col)
281 {
282   PetscErrorCode ierr;
283   PetscInt       low,high;
284   PetscScalar    zero=0.0,one=1.0;
285 
286   PetscFunctionBegin;
287   ierr=VecSet(Y, zero);CHKERRQ(ierr);
288   ierr=VecGetOwnershipRange(Y,&low,&high);CHKERRQ(ierr);
289   if (col>=low && col<high){
290     ierr=VecSetValue(Y,col,one,INSERT_VALUES);CHKERRQ(ierr);
291   }
292   ierr=VecAssemblyBegin(Y);CHKERRQ(ierr);
293   ierr=VecAssemblyEnd(Y);CHKERRQ(ierr);
294   ierr=MatMult_ADA(mat,Y,Y);CHKERRQ(ierr);
295   PetscFunctionReturn(0);
296 }
297 
MatConvert_ADA(Mat mat,MatType newtype,Mat * NewMat)298 PETSC_INTERN PetscErrorCode MatConvert_ADA(Mat mat,MatType newtype,Mat *NewMat)
299 {
300   PetscErrorCode ierr;
301   PetscMPIInt    size;
302   PetscBool      sametype, issame, isdense, isseqdense;
303   TaoMatADACtx   ctx;
304 
305   PetscFunctionBegin;
306   ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr);
307   ierr = MPI_Comm_size(PetscObjectComm((PetscObject)mat),&size);CHKERRQ(ierr);
308 
309   ierr = PetscObjectTypeCompare((PetscObject)mat,newtype,&sametype);CHKERRQ(ierr);
310   ierr = PetscObjectTypeCompare((PetscObject)mat,MATSAME,&issame);CHKERRQ(ierr);
311   ierr = PetscObjectTypeCompare((PetscObject)mat,MATMPIDENSE,&isdense);CHKERRQ(ierr);
312   ierr = PetscObjectTypeCompare((PetscObject)mat,MATSEQDENSE,&isseqdense);CHKERRQ(ierr);
313 
314   if (sametype || issame) {
315     ierr=MatDuplicate(mat,MAT_COPY_VALUES,NewMat);CHKERRQ(ierr);
316   } else if (isdense) {
317     PetscInt          i,j,low,high,m,n,M,N;
318     const PetscScalar *dptr;
319     Vec               X;
320 
321     ierr = VecDuplicate(ctx->D2,&X);CHKERRQ(ierr);
322     ierr=MatGetSize(mat,&M,&N);CHKERRQ(ierr);
323     ierr=MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr);
324     ierr = MatCreateDense(PetscObjectComm((PetscObject)mat),m,m,N,N,NULL,NewMat);CHKERRQ(ierr);
325     ierr = MatGetOwnershipRange(*NewMat,&low,&high);CHKERRQ(ierr);
326     for (i=0;i<M;i++){
327       ierr = MatGetColumnVector_ADA(mat,X,i);CHKERRQ(ierr);
328       ierr = VecGetArrayRead(X,&dptr);CHKERRQ(ierr);
329       for (j=0; j<high-low; j++){
330         ierr = MatSetValue(*NewMat,low+j,i,dptr[j],INSERT_VALUES);CHKERRQ(ierr);
331       }
332       ierr = VecRestoreArrayRead(X,&dptr);CHKERRQ(ierr);
333     }
334     ierr = MatAssemblyBegin(*NewMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
335     ierr = MatAssemblyEnd(*NewMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
336     ierr = VecDestroy(&X);CHKERRQ(ierr);
337   } else if (isseqdense && size==1){
338     PetscInt          i,j,low,high,m,n,M,N;
339     const PetscScalar *dptr;
340     Vec               X;
341 
342     ierr = VecDuplicate(ctx->D2,&X);CHKERRQ(ierr);
343     ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr);
344     ierr = MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr);
345     ierr = MatCreateSeqDense(PetscObjectComm((PetscObject)mat),N,N,NULL,NewMat);CHKERRQ(ierr);
346     ierr = MatGetOwnershipRange(*NewMat,&low,&high);CHKERRQ(ierr);
347     for (i=0;i<M;i++){
348       ierr = MatGetColumnVector_ADA(mat,X,i);CHKERRQ(ierr);
349       ierr = VecGetArrayRead(X,&dptr);CHKERRQ(ierr);
350       for (j=0; j<high-low; j++){
351         ierr = MatSetValue(*NewMat,low+j,i,dptr[j],INSERT_VALUES);CHKERRQ(ierr);
352       }
353       ierr = VecRestoreArrayRead(X,&dptr);CHKERRQ(ierr);
354     }
355     ierr = MatAssemblyBegin(*NewMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
356     ierr = MatAssemblyEnd(*NewMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
357     ierr = VecDestroy(&X);CHKERRQ(ierr);
358   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"No support to convert objects to that type");
359   PetscFunctionReturn(0);
360 }
361 
MatNorm_ADA(Mat mat,NormType type,PetscReal * norm)362 static PetscErrorCode MatNorm_ADA(Mat mat,NormType type,PetscReal *norm)
363 {
364   PetscErrorCode ierr;
365   TaoMatADACtx   ctx;
366 
367   PetscFunctionBegin;
368   ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr);
369   if (type == NORM_FROBENIUS) {
370     *norm = 1.0;
371   } else if (type == NORM_1 || type == NORM_INFINITY) {
372     *norm = 1.0;
373   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
374   PetscFunctionReturn(0);
375 }
376 
377 /*@C
378    MatCreateADA - Creates a matrix M=A^T D1 A + D2 where D1, D2 are diagonal
379 
380    Collective on matrix
381 
382    Input Parameters:
383 +  mat - matrix of arbitrary type
384 .  d1 - A vector defining a diagonal matrix
385 -  d2 - A vector defining a diagonal matrix
386 
387    Output Parameters:
388 .  J - New matrix whose operations are defined in terms of mat, D1, and D2.
389 
390    Notes:
391    The user provides the input data and is responsible for destroying
392    this data after matrix J has been destroyed.
393 
394    Level: developer
395 
396 .seealso: MatCreate()
397 @*/
MatCreateADA(Mat mat,Vec d1,Vec d2,Mat * J)398 PetscErrorCode MatCreateADA(Mat mat,Vec d1, Vec d2, Mat *J)
399 {
400   MPI_Comm       comm = PetscObjectComm((PetscObject)mat);
401   TaoMatADACtx   ctx;
402   PetscErrorCode ierr;
403   PetscInt       nloc,n;
404 
405   PetscFunctionBegin;
406   ierr = PetscNew(&ctx);CHKERRQ(ierr);
407   ctx->A=mat;
408   ctx->D1=d1;
409   ctx->D2=d2;
410   if (d1){
411     ierr = VecDuplicate(d1,&ctx->W);CHKERRQ(ierr);
412     ierr = PetscObjectReference((PetscObject)d1);CHKERRQ(ierr);
413   } else {
414     ctx->W = NULL;
415   }
416   if (d2){
417     ierr = VecDuplicate(d2,&ctx->W2);CHKERRQ(ierr);
418     ierr = VecDuplicate(d2,&ctx->ADADiag);CHKERRQ(ierr);
419     ierr =  PetscObjectReference((PetscObject)d2);CHKERRQ(ierr);
420   } else {
421     ctx->W2      = NULL;
422     ctx->ADADiag = NULL;
423   }
424 
425   ctx->GotDiag = 0;
426   ierr =  PetscObjectReference((PetscObject)mat);CHKERRQ(ierr);
427 
428   ierr=VecGetLocalSize(d2,&nloc);CHKERRQ(ierr);
429   ierr=VecGetSize(d2,&n);CHKERRQ(ierr);
430 
431   ierr = MatCreateShell(comm,nloc,nloc,n,n,ctx,J);CHKERRQ(ierr);
432   ierr = MatShellSetManageScalingShifts(*J);CHKERRQ(ierr);
433   ierr = MatShellSetOperation(*J,MATOP_MULT,(void(*)(void))MatMult_ADA);CHKERRQ(ierr);
434   ierr = MatShellSetOperation(*J,MATOP_DESTROY,(void(*)(void))MatDestroy_ADA);CHKERRQ(ierr);
435   ierr = MatShellSetOperation(*J,MATOP_VIEW,(void(*)(void))MatView_ADA);CHKERRQ(ierr);
436   ierr = MatShellSetOperation(*J,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_ADA);CHKERRQ(ierr);
437   ierr = MatShellSetOperation(*J,MATOP_DIAGONAL_SET,(void(*)(void))MatDiagonalSet_ADA);CHKERRQ(ierr);
438   ierr = MatShellSetOperation(*J,MATOP_SHIFT,(void(*)(void))MatShift_ADA);CHKERRQ(ierr);
439   ierr = MatShellSetOperation(*J,MATOP_EQUAL,(void(*)(void))MatEqual_ADA);CHKERRQ(ierr);
440   ierr = MatShellSetOperation(*J,MATOP_SCALE,(void(*)(void))MatScale_ADA);CHKERRQ(ierr);
441   ierr = MatShellSetOperation(*J,MATOP_TRANSPOSE,(void(*)(void))MatTranspose_ADA);CHKERRQ(ierr);
442   ierr = MatShellSetOperation(*J,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_ADA);CHKERRQ(ierr);
443   ierr = MatShellSetOperation(*J,MATOP_CREATE_SUBMATRICES,(void(*)(void))MatCreateSubMatrices_ADA);CHKERRQ(ierr);
444   ierr = MatShellSetOperation(*J,MATOP_NORM,(void(*)(void))MatNorm_ADA);CHKERRQ(ierr);
445   ierr = MatShellSetOperation(*J,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_ADA);CHKERRQ(ierr);
446   ierr = MatShellSetOperation(*J,MATOP_CREATE_SUBMATRIX,(void(*)(void))MatCreateSubMatrix_ADA);CHKERRQ(ierr);
447 
448   ierr = PetscLogObjectParent((PetscObject)(*J),(PetscObject)ctx->W);CHKERRQ(ierr);
449   ierr = PetscLogObjectParent((PetscObject)mat,(PetscObject)(*J));CHKERRQ(ierr);
450 
451   ierr = MatSetOption(*J,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr);
452   PetscFunctionReturn(0);
453 }
454