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