Actual source code: matadic.c
1: #define PETSCMAT_DLL
3: /*
4: ADIC matrix-free matrix implementation
5: */
7: #include private/matimpl.h
8: #include petscda.h
9: #include petscsnes.h
10: #include petscsys.h
12: #include "adic/ad_utils.h"
15: typedef struct {
16: DA da;
17: Vec localu; /* point at which Jacobian is evaluated */
18: void *ctx;
19: SNES snes;
20: Vec diagonal; /* current matrix diagonal */
21: PetscTruth diagonalvalid; /* indicates if diagonal matches current base vector */
22: } Mat_DAAD;
26: PetscErrorCode MatAssemblyEnd_DAAD(Mat A,MatAssemblyType atype)
27: {
28: Mat_DAAD *a = (Mat_DAAD*)A->data;
30: Vec u;
33: a->diagonalvalid = PETSC_FALSE;
34: if (a->snes) {
35: SNESGetSolution(a->snes,&u);
36: DAGlobalToLocalBegin(a->da,u,INSERT_VALUES,a->localu);
37: DAGlobalToLocalEnd(a->da,u,INSERT_VALUES,a->localu);
38: }
39: return(0);
40: }
44: PetscErrorCode MatMult_DAAD(Mat A,Vec xx,Vec yy)
45: {
46: Mat_DAAD *a = (Mat_DAAD*)A->data;
47: Vec localxx;
51: DAGetLocalVector(a->da,&localxx);
52: DAGlobalToLocalBegin(a->da,xx,INSERT_VALUES,localxx);
53: DAGlobalToLocalEnd(a->da,xx,INSERT_VALUES,localxx);
54: DAMultiplyByJacobian1WithAD(a->da,a->localu,localxx,yy,a->ctx);
55: DARestoreLocalVector(a->da,&localxx);
56: return(0);
57: }
59: #include ../src/dm/da/daimpl.h
63: PetscErrorCode MatGetDiagonal_DAAD(Mat A,Vec dd)
64: {
65: Mat_DAAD *a = (Mat_DAAD*)A->data;
67: int j,nI,gI,gtdof;
68: PetscScalar *avu,*ad_vustart,ad_f[2],*d;
69: DALocalInfo info;
70: MatStencil stencil;
71: void* *ad_vu;
75: /* get space for derivative object. */
76: DAGetAdicMFArray(a->da,PETSC_TRUE,(void **)&ad_vu,(void**)&ad_vustart,>dof);
78: /* copy input vector into derivative object */
79: VecGetArray(a->localu,&avu);
80: for (j=0; j<gtdof; j++) {
81: ad_vustart[2*j] = avu[j];
82: ad_vustart[2*j+1] = 0.0;
83: }
84: VecRestoreArray(a->localu,&avu);
86: PetscADResetIndep();
87: PetscADIncrementTotalGradSize(1);
88: PetscADSetIndepDone();
90: VecGetArray(dd,&d);
92: DAGetLocalInfo(a->da,&info);
93: nI = 0;
94: for (stencil.k = info.zs; stencil.k<info.zs+info.zm; stencil.k++) {
95: for (stencil.j = info.ys; stencil.j<info.ys+info.ym; stencil.j++) {
96: for (stencil.i = info.xs; stencil.i<info.xs+info.xm; stencil.i++) {
97: for (stencil.c = 0; stencil.c<info.dof; stencil.c++) {
98: gI = stencil.c + (stencil.i - info.gxs)*info.dof + (stencil.j - info.gys)*info.dof*info.gxm + (stencil.k - info.gzs)*info.dof*info.gxm*info.gym;
99: ad_vustart[1+2*gI] = 1.0;
100: (*a->da->adicmf_lfi)(&info,&stencil,ad_vu,ad_f,a->ctx);
101: d[nI] = ad_f[1];
102: ad_vustart[1+2*gI] = 0.0;
103: nI++;
104: }
105: }
106: }
107: }
109: VecRestoreArray(dd,&d);
110: DARestoreAdicMFArray(a->da,PETSC_TRUE,(void **)&ad_vu,(void**)&ad_vustart,>dof);
111: return(0);
112: }
117: PetscErrorCode MatRelax_DAAD(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,int its,int lits,Vec xx)
118: {
119: Mat_DAAD *a = (Mat_DAAD*)A->data;
121: int j,gtdof,nI,gI;
122: PetscScalar *avu,*av,*ad_vustart,ad_f[2],*d,*b;
123: Vec localxx,dd;
124: DALocalInfo info;
125: MatStencil stencil;
126: void* *ad_vu;
129: if (omega != 1.0) SETERRQ(PETSC_ERR_ARG_WRONG,"Currently only support omega of 1.0");
130: if (fshift) SETERRQ(PETSC_ERR_ARG_WRONG,"Currently do not support fshift");
131: if (its <= 0 || lits <= 0) SETERRQ2(PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits);
133: if (!a->diagonal) {
134: DACreateGlobalVector(a->da,&a->diagonal);
135: }
136: if (!a->diagonalvalid) {
137: MatGetDiagonal(A,a->diagonal);
138: a->diagonalvalid = PETSC_TRUE;
139: }
140: dd = a->diagonal;
143: DAGetLocalVector(a->da,&localxx);
144: if (flag & SOR_ZERO_INITIAL_GUESS) {
145: VecSet(localxx,0.0);
146: } else {
147: DAGlobalToLocalBegin(a->da,xx,INSERT_VALUES,localxx);
148: DAGlobalToLocalEnd(a->da,xx,INSERT_VALUES,localxx);
149: }
151: /* get space for derivative object. */
152: DAGetAdicMFArray(a->da,PETSC_TRUE,(void **)&ad_vu,(void**)&ad_vustart,>dof);
154: /* copy input vector into derivative object */
155: VecGetArray(a->localu,&avu);
156: VecGetArray(localxx,&av);
157: for (j=0; j<gtdof; j++) {
158: ad_vustart[2*j] = avu[j];
159: ad_vustart[2*j+1] = av[j];
160: }
161: VecRestoreArray(a->localu,&avu);
162: VecRestoreArray(localxx,&av);
164: PetscADResetIndep();
165: PetscADIncrementTotalGradSize(1);
166: PetscADSetIndepDone();
168: VecGetArray(dd,&d);
169: VecGetArray(bb,&b);
171: DAGetLocalInfo(a->da,&info);
172: while (its--) {
173: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){
174: nI = 0;
175: for (stencil.k = info.zs; stencil.k<info.zs+info.zm; stencil.k++) {
176: for (stencil.j = info.ys; stencil.j<info.ys+info.ym; stencil.j++) {
177: for (stencil.i = info.xs; stencil.i<info.xs+info.xm; stencil.i++) {
178: for (stencil.c = 0; stencil.c<info.dof; stencil.c++) {
179: (*a->da->adicmf_lfi)(&info,&stencil,ad_vu,ad_f,a->ctx);
180: gI = stencil.c + (stencil.i - info.gxs)*info.dof + (stencil.j - info.gys)*info.dof*info.gxm + (stencil.k - info.gzs)*info.dof*info.gxm*info.gym;
181: ad_vustart[1+2*gI] += (b[nI] - ad_f[1])/d[nI];
182: nI++;
183: }
184: }
185: }
186: }
187: }
188: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){
189: nI = info.dof*info.xm*info.ym*info.zm - 1;
190: for (stencil.k = info.zs+info.zm-1; stencil.k>=info.zs; stencil.k--) {
191: for (stencil.j = info.ys+info.ym-1; stencil.j>=info.ys; stencil.j--) {
192: for (stencil.i = info.xs+info.xm-1; stencil.i>=info.xs; stencil.i--) {
193: for (stencil.c = info.dof-1; stencil.c>=0; stencil.c--) {
194: (*a->da->adicmf_lfi)(&info,&stencil,ad_vu,ad_f,a->ctx);
195: gI = stencil.c + (stencil.i - info.gxs)*info.dof + (stencil.j - info.gys)*info.dof*info.gxm + (stencil.k - info.gzs)*info.dof*info.gxm*info.gym;
196: ad_vustart[1+2*gI] += (b[nI] - ad_f[1])/d[nI];
197: nI--;
198: }
199: }
200: }
201: }
202: }
203: }
205: VecRestoreArray(dd,&d);
206: VecRestoreArray(bb,&b);
208: VecGetArray(localxx,&av);
209: for (j=0; j<gtdof; j++) {
210: av[j] = ad_vustart[2*j+1];
211: }
212: VecRestoreArray(localxx,&av);
213: DALocalToGlobal(a->da,localxx,INSERT_VALUES,xx);
215: DARestoreLocalVector(a->da,&localxx);
216: DARestoreAdicMFArray(a->da,PETSC_TRUE,(void **)&ad_vu,(void**)&ad_vustart,>dof);
217: return(0);
218: }
224: PetscErrorCode MatDestroy_DAAD(Mat A)
225: {
226: Mat_DAAD *a = (Mat_DAAD*)A->data;
230: DADestroy(a->da);
231: VecDestroy(a->localu);
232: if (a->diagonal) {VecDestroy(a->diagonal);}
233: PetscFree(a);
234: PetscObjectChangeTypeName((PetscObject)A,0);
235: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetBase_C","",PETSC_NULL);
236: PetscObjectComposeFunction((PetscObject)A,"MatDAADSetDA_C","",PETSC_NULL);
237: PetscObjectComposeFunction((PetscObject)A,"MatDAADSetSNES_C","",PETSC_NULL);
238: PetscObjectComposeFunction((PetscObject)A,"MatDAADSetCtx_C","",PETSC_NULL);
239: return(0);
240: }
242: /* -------------------------------------------------------------------*/
243: static struct _MatOps MatOps_Values = {0,
244: 0,
245: 0,
246: MatMult_DAAD,
247: /* 4*/ 0,
248: 0,
249: 0,
250: 0,
251: 0,
252: 0,
253: /*10*/ 0,
254: 0,
255: 0,
256: MatRelax_DAAD,
257: 0,
258: /*15*/ 0,
259: 0,
260: MatGetDiagonal_DAAD,
261: 0,
262: 0,
263: /*20*/ 0,
264: MatAssemblyEnd_DAAD,
265: 0,
266: 0,
267: 0,
268: /*25*/ 0,
269: 0,
270: 0,
271: 0,
272: 0,
273: /*30*/ 0,
274: 0,
275: 0,
276: 0,
277: 0,
278: /*35*/ 0,
279: 0,
280: 0,
281: 0,
282: 0,
283: /*40*/ 0,
284: 0,
285: 0,
286: 0,
287: 0,
288: /*45*/ 0,
289: 0,
290: 0,
291: 0,
292: 0,
293: /*50*/ 0,
294: 0,
295: 0,
296: 0,
297: 0,
298: /*55*/ 0,
299: 0,
300: 0,
301: 0,
302: 0,
303: /*60*/ 0,
304: MatDestroy_DAAD,
305: 0,
306: 0,
307: 0,
308: /*65*/ 0,
309: 0,
310: 0,
311: 0,
312: 0,
313: /*70*/ 0,
314: 0,
315: 0,
316: 0,
317: 0,
318: /*75*/ 0,
319: 0,
320: 0,
321: 0,
322: 0,
323: /*80*/ 0,
324: 0,
325: 0,
326: 0,
327: 0,
328: /*85*/ 0,
329: 0,
330: 0,
331: 0,
332: 0,
333: /*90*/ 0,
334: 0,
335: 0,
336: 0,
337: 0,
338: /*95*/ 0,
339: 0,
340: 0,
341: 0};
343: /* --------------------------------------------------------------------------------*/
348: PetscErrorCode MatMFFDSetBase_AD(Mat J,Vec U,Vec F)
349: {
351: Mat_DAAD *a = (Mat_DAAD*)J->data;
354: a->diagonalvalid = PETSC_FALSE;
355: DAGlobalToLocalBegin(a->da,U,INSERT_VALUES,a->localu);
356: DAGlobalToLocalEnd(a->da,U,INSERT_VALUES,a->localu);
357: return(0);
358: }
364: PetscErrorCode MatDAADSetDA_AD(Mat A,DA da)
365: {
366: Mat_DAAD *a = (Mat_DAAD*)A->data;
368: PetscInt nc,nx,ny,nz,Nx,Ny,Nz;
371: PetscObjectReference((PetscObject)da);
372: if (a->da) { DADestroy(a->da); }
373: a->da = da;
374: DAGetInfo(da,0,&Nx,&Ny,&Nz,0,0,0,&nc,0,0,0);
375: DAGetCorners(da,0,0,0,&nx,&ny,&nz);
376: A->rmap->n = A->cmap->n = nc*nx*ny*nz;
377: A->rmap->N = A->cmap->N = nc*Nx*Ny*Nz;
378: DACreateLocalVector(da,&a->localu);
379: return(0);
380: }
386: PetscErrorCode MatDAADSetSNES_AD(Mat A,SNES snes)
387: {
388: Mat_DAAD *a = (Mat_DAAD*)A->data;
391: a->snes = snes;
392: return(0);
393: }
399: PetscErrorCode MatDAADSetCtx_AD(Mat A,void *ctx)
400: {
401: Mat_DAAD *a = (Mat_DAAD*)A->data;
404: a->ctx = ctx;
405: return(0);
406: }
409: /*MC
410: MATDAAD - MATDAAD = "daad" - A matrix type that can do matrix-vector products using a local function that
411: is differentiated with ADIFOR or ADIC.
413: Level: intermediate
415: .seealso: MatCreateDAAD
416: M*/
421: PetscErrorCode MatCreate_DAAD(Mat B)
422: {
423: Mat_DAAD *b;
427: PetscNewLog(B,Mat_DAAD,&b);
428: B->data = (void*)b;
429: PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
430:
431: PetscMapSetBlockSize(B->rmap,1);
432: PetscMapSetBlockSize(B->cmap,1);
433: PetscMapSetUp(B->rmap);
434: PetscMapSetUp(B->cmap);
436: PetscObjectChangeTypeName((PetscObject)B,MATDAAD);
437: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMFFDSetBase_C","MatMFFDSetBase_AD",MatMFFDSetBase_AD);
438: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatDAADSetDA_C","MatDAADSetDA_AD",MatDAADSetDA_AD);
439: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatDAADSetSNES_C","MatDAADSetSNES_AD",MatDAADSetSNES_AD);
440: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatDAADSetCtx_C","MatDAADSetCtx_AD",MatDAADSetCtx_AD);
441: PetscObjectChangeTypeName((PetscObject)B,MATDAAD);
442: return(0);
443: }
449: /*@C
450: MatDAADSetDA - Tells the matrix what DA it is using for layout and Jacobian.
452: Collective on Mat and DA
454: Input Parameters:
455: + mat - the matrix
456: - da - the DA
458: Level: intermediate
460: .seealso: MatCreate(), DASetLocalAdicMFFunction(), MatCreateDAAD()
462: @*/
463: PetscErrorCode MatDAADSetDA(Mat A,DA da)
464: {
465: PetscErrorCode ierr,(*f)(Mat,void*);
470: PetscObjectQueryFunction((PetscObject)A,"MatDAADSetDA_C",(void (**)(void))&f);
471: if (f) {
472: (*f)(A,da);
473: }
474: return(0);
475: }
479: /*@C
480: MatDAADSetSNES - Tells the matrix what SNES it is using for the base U.
482: Collective on Mat and SNES
484: Input Parameters:
485: + mat - the matrix
486: - snes - the SNES
488: Level: intermediate
490: Notes: this is currently turned off for Fortran usage
492: .seealso: MatCreate(), DASetLocalAdicMFFunction(), MatCreateDAAD(), MatDAADSetDA()
494: @*/
495: PetscErrorCode MatDAADSetSNES(Mat A,SNES snes)
496: {
497: PetscErrorCode ierr,(*f)(Mat,void*);
502: PetscObjectQueryFunction((PetscObject)A,"MatDAADSetSNES_C",(void (**)(void))&f);
503: if (f) {
504: (*f)(A,snes);
505: }
506: return(0);
507: }
511: /*@C
512: MatDAADSetCtx - Sets the user context for a DAAD (ADIC matrix-free) matrix.
514: Collective on Mat
516: Input Parameters:
517: + mat - the matrix
518: - ctx - the context
520: Level: intermediate
522: .seealso: MatCreate(), DASetLocalAdicMFFunction(), MatCreateDAAD(), MatDAADSetDA()
524: @*/
525: PetscErrorCode MatDAADSetCtx(Mat A,void *ctx)
526: {
527: PetscErrorCode ierr,(*f)(Mat,void*);
531: PetscObjectQueryFunction((PetscObject)A,"MatDAADSetCtx_C",(void (**)(void))&f);
532: if (f) {
533: (*f)(A,ctx);
534: }
535: return(0);
536: }
540: /*@C
541: MatCreateDAAD - Creates a matrix that can do matrix-vector products using a local
542: function that is differentiated with ADIFOR or ADIC.
544: Collective on DA
546: Input Parameters:
547: . da - the DA that defines the distribution of the vectors
549: Output Parameter:
550: . A - the matrix
552: Level: intermediate
554: Notes: this is currently turned off for Fortran
556: .seealso: MatCreate(), DASetLocalAdicMFFunction()
558: @*/
559: PetscErrorCode MatCreateDAAD(DA da,Mat *A)
560: {
562: MPI_Comm comm;
565: PetscObjectGetComm((PetscObject)da,&comm);
566: MatCreate(comm,A);
567: MatSetType(*A,MATDAAD);
568: MatDAADSetDA(*A,da);
569: return(0);
570: }
575: /*@
576: MatRegisterDAAD - Registers DAAD matrix type
578: Level: advanced
580: .seealso: MatCreateDAAD(), DASetLocalAdicMFFunction()
582: @*/
583: PetscErrorCode MatRegisterDAAD(void)
584: {
587: MatRegisterDynamic(MATDAAD,PETSC_NULL,"MatCreate_DAAD",MatCreate_DAAD);
588: return(0);
589: }