Actual source code: matadic.c
petsc-3.3-p7 2013-05-11
2: /*
3: ADIC matrix-free matrix implementation
4: */
6: #include <petsc-private/matimpl.h>
7: #include <petscdmda.h> /*I "petscdmda.h" I*/
8: #include <petscsnes.h> /*I "petscsnes.h" I*/
9: EXTERN_C_BEGIN
10: #include <adic/ad_utils.h>
11: EXTERN_C_END
13: typedef struct {
14: DM da;
15: Vec localu; /* point at which Jacobian is evaluated */
16: void *ctx;
17: SNES snes;
18: Vec diagonal; /* current matrix diagonal */
19: PetscBool diagonalvalid; /* indicates if diagonal matches current base vector */
20: } Mat_DAAD;
24: PetscErrorCode MatAssemblyEnd_DAAD(Mat A,MatAssemblyType atype)
25: {
26: Mat_DAAD *a = (Mat_DAAD*)A->data;
28: Vec u;
31: a->diagonalvalid = PETSC_FALSE;
32: if (a->snes) {
33: SNESGetSolution(a->snes,&u);
34: DMGlobalToLocalBegin(a->da,u,INSERT_VALUES,a->localu);
35: DMGlobalToLocalEnd(a->da,u,INSERT_VALUES,a->localu);
36: }
37: return(0);
38: }
42: PetscErrorCode MatMult_DAAD(Mat A,Vec xx,Vec yy)
43: {
44: Mat_DAAD *a = (Mat_DAAD*)A->data;
45: Vec localxx;
49: DMGetLocalVector(a->da,&localxx);
50: DMGlobalToLocalBegin(a->da,xx,INSERT_VALUES,localxx);
51: DMGlobalToLocalEnd(a->da,xx,INSERT_VALUES,localxx);
52: DMDAMultiplyByJacobian1WithAD(a->da,a->localu,localxx,yy,a->ctx);
53: DMRestoreLocalVector(a->da,&localxx);
54: return(0);
55: }
57: #include <../src/dm/da/daimpl.h>
61: PetscErrorCode MatGetDiagonal_DAAD(Mat A,Vec dd)
62: {
63: Mat_DAAD *a = (Mat_DAAD*)A->data;
65: int j,nI,gI,gtdof;
66: PetscScalar *avu,*ad_vustart,ad_f[2],*d;
67: DMDALocalInfo info;
68: MatStencil stencil;
69: void* *ad_vu;
73: /* get space for derivative object. */
74: DMDAGetAdicMFArray(a->da,PETSC_TRUE,(void **)&ad_vu,(void**)&ad_vustart,>dof);
76: /* copy input vector into derivative object */
77: VecGetArray(a->localu,&avu);
78: for (j=0; j<gtdof; j++) {
79: ad_vustart[2*j] = avu[j];
80: ad_vustart[2*j+1] = 0.0;
81: }
82: VecRestoreArray(a->localu,&avu);
84: PetscADResetIndep();
85: PetscADIncrementTotalGradSize(1);
86: PetscADSetIndepDone();
88: VecGetArray(dd,&d);
90: DMDAGetLocalInfo(a->da,&info);
91: nI = 0;
92: for (stencil.k = info.zs; stencil.k<info.zs+info.zm; stencil.k++) {
93: for (stencil.j = info.ys; stencil.j<info.ys+info.ym; stencil.j++) {
94: for (stencil.i = info.xs; stencil.i<info.xs+info.xm; stencil.i++) {
95: for (stencil.c = 0; stencil.c<info.dof; stencil.c++) {
96: 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;
97: ad_vustart[1+2*gI] = 1.0;
98: (*a->da->adicmf_lfi)(&info,&stencil,ad_vu,ad_f,a->ctx);
99: d[nI] = ad_f[1];
100: ad_vustart[1+2*gI] = 0.0;
101: nI++;
102: }
103: }
104: }
105: }
107: VecRestoreArray(dd,&d);
108: DMDARestoreAdicMFArray(a->da,PETSC_TRUE,(void **)&ad_vu,(void**)&ad_vustart,>dof);
109: return(0);
110: }
115: PetscErrorCode MatSOR_DAAD(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,int its,int lits,Vec xx)
116: {
117: Mat_DAAD *a = (Mat_DAAD*)A->data;
119: int j,gtdof,nI,gI;
120: PetscScalar *avu,*av,*ad_vustart,ad_f[2],*d,*b;
121: Vec localxx,dd;
122: DMDALocalInfo info;
123: MatStencil stencil;
124: void* *ad_vu;
127: if (omega != 1.0) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONG,"Currently only support omega of 1.0");
128: if (fshift) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONG,"Currently do not support fshift");
129: if (its <= 0 || lits <= 0) SETERRQ2(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits);
131: if (!a->diagonal) {
132: DMCreateGlobalVector(a->da,&a->diagonal);
133: }
134: if (!a->diagonalvalid) {
135: MatGetDiagonal(A,a->diagonal);
136: a->diagonalvalid = PETSC_TRUE;
137: }
138: dd = a->diagonal;
141: DMGetLocalVector(a->da,&localxx);
142: if (flag & SOR_ZERO_INITIAL_GUESS) {
143: VecSet(localxx,0.0);
144: } else {
145: DMGlobalToLocalBegin(a->da,xx,INSERT_VALUES,localxx);
146: DMGlobalToLocalEnd(a->da,xx,INSERT_VALUES,localxx);
147: }
149: /* get space for derivative object. */
150: DMDAGetAdicMFArray(a->da,PETSC_TRUE,(void **)&ad_vu,(void**)&ad_vustart,>dof);
152: /* copy input vector into derivative object */
153: VecGetArray(a->localu,&avu);
154: VecGetArray(localxx,&av);
155: for (j=0; j<gtdof; j++) {
156: ad_vustart[2*j] = avu[j];
157: ad_vustart[2*j+1] = av[j];
158: }
159: VecRestoreArray(a->localu,&avu);
160: VecRestoreArray(localxx,&av);
162: PetscADResetIndep();
163: PetscADIncrementTotalGradSize(1);
164: PetscADSetIndepDone();
166: VecGetArray(dd,&d);
167: VecGetArray(bb,&b);
169: DMDAGetLocalInfo(a->da,&info);
170: while (its--) {
171: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){
172: nI = 0;
173: for (stencil.k = info.zs; stencil.k<info.zs+info.zm; stencil.k++) {
174: for (stencil.j = info.ys; stencil.j<info.ys+info.ym; stencil.j++) {
175: for (stencil.i = info.xs; stencil.i<info.xs+info.xm; stencil.i++) {
176: for (stencil.c = 0; stencil.c<info.dof; stencil.c++) {
177: (*a->da->adicmf_lfi)(&info,&stencil,ad_vu,ad_f,a->ctx);
178: 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;
179: ad_vustart[1+2*gI] += (b[nI] - ad_f[1])/d[nI];
180: nI++;
181: }
182: }
183: }
184: }
185: }
186: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){
187: nI = info.dof*info.xm*info.ym*info.zm - 1;
188: for (stencil.k = info.zs+info.zm-1; stencil.k>=info.zs; stencil.k--) {
189: for (stencil.j = info.ys+info.ym-1; stencil.j>=info.ys; stencil.j--) {
190: for (stencil.i = info.xs+info.xm-1; stencil.i>=info.xs; stencil.i--) {
191: for (stencil.c = info.dof-1; stencil.c>=0; stencil.c--) {
192: (*a->da->adicmf_lfi)(&info,&stencil,ad_vu,ad_f,a->ctx);
193: 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;
194: ad_vustart[1+2*gI] += (b[nI] - ad_f[1])/d[nI];
195: nI--;
196: }
197: }
198: }
199: }
200: }
201: }
203: VecRestoreArray(dd,&d);
204: VecRestoreArray(bb,&b);
206: VecGetArray(localxx,&av);
207: for (j=0; j<gtdof; j++) {
208: av[j] = ad_vustart[2*j+1];
209: }
210: VecRestoreArray(localxx,&av);
211: DMLocalToGlobalBegin(a->da,localxx,INSERT_VALUES,xx);
212: DMLocalToGlobalEnd(a->da,localxx,INSERT_VALUES,xx);
214: DMRestoreLocalVector(a->da,&localxx);
215: DMDARestoreAdicMFArray(a->da,PETSC_TRUE,(void **)&ad_vu,(void**)&ad_vustart,>dof);
216: return(0);
217: }
223: PetscErrorCode MatDestroy_DAAD(Mat A)
224: {
225: Mat_DAAD *a = (Mat_DAAD*)A->data;
229: DMDestroy(&a->da);
230: VecDestroy(&a->localu);
231: VecDestroy(&a->diagonal);
232: PetscFree(A->data);
233: PetscObjectChangeTypeName((PetscObject)A,0);
234: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetBase_C","",PETSC_NULL);
235: PetscObjectComposeFunction((PetscObject)A,"MatDAADSetDA_C","",PETSC_NULL);
236: PetscObjectComposeFunction((PetscObject)A,"MatDAADSetSNES_C","",PETSC_NULL);
237: PetscObjectComposeFunction((PetscObject)A,"MatDAADSetCtx_C","",PETSC_NULL);
238: return(0);
239: }
241: /* -------------------------------------------------------------------*/
242: static struct _MatOps MatOps_Values = {0,
243: 0,
244: 0,
245: MatMult_DAAD,
246: /* 4*/ 0,
247: 0,
248: 0,
249: 0,
250: 0,
251: 0,
252: /*10*/ 0,
253: 0,
254: 0,
255: MatSOR_DAAD,
256: 0,
257: /*15*/ 0,
258: 0,
259: MatGetDiagonal_DAAD,
260: 0,
261: 0,
262: /*20*/ 0,
263: MatAssemblyEnd_DAAD,
264: 0,
265: 0,
266: /*24*/ 0,
267: 0,
268: 0,
269: 0,
270: 0,
271: /*29*/ 0,
272: 0,
273: 0,
274: 0,
275: 0,
276: /*34*/ 0,
277: 0,
278: 0,
279: 0,
280: 0,
281: /*39*/ 0,
282: 0,
283: 0,
284: 0,
285: 0,
286: /*44*/ 0,
287: 0,
288: 0,
289: 0,
290: 0,
291: /*49*/ 0,
292: 0,
293: 0,
294: 0,
295: 0,
296: /*54*/ 0,
297: 0,
298: 0,
299: 0,
300: 0,
301: /*59*/ 0,
302: MatDestroy_DAAD,
303: 0,
304: 0,
305: 0,
306: /*64*/ 0,
307: 0,
308: 0,
309: 0,
310: 0,
311: /*69*/ 0,
312: 0,
313: 0,
314: 0,
315: 0,
316: /*74*/ 0,
317: 0,
318: 0,
319: 0,
320: 0,
321: /*79*/ 0,
322: 0,
323: 0,
324: 0,
325: 0,
326: /*84*/ 0,
327: 0,
328: 0,
329: 0,
330: 0,
331: /*89*/ 0,
332: 0,
333: 0,
334: 0,
335: 0,
336: /*94*/ 0,
337: 0,
338: 0,
339: 0};
341: /* --------------------------------------------------------------------------------*/
343: EXTERN_C_BEGIN
346: PetscErrorCode MatMFFDSetBase_AD(Mat J,Vec U,Vec F)
347: {
349: Mat_DAAD *a = (Mat_DAAD*)J->data;
352: a->diagonalvalid = PETSC_FALSE;
353: DMGlobalToLocalBegin(a->da,U,INSERT_VALUES,a->localu);
354: DMGlobalToLocalEnd(a->da,U,INSERT_VALUES,a->localu);
355: return(0);
356: }
357: EXTERN_C_END
359: EXTERN_C_BEGIN
362: PetscErrorCode MatDAADSetDA_AD(Mat A,DM da)
363: {
364: Mat_DAAD *a = (Mat_DAAD*)A->data;
366: PetscInt nc,nx,ny,nz,Nx,Ny,Nz;
369: PetscObjectReference((PetscObject)da);
370: if (a->da) { DMDestroy(a->da); }
371: a->da = da;
372: DMDAGetInfo(da,0,&Nx,&Ny,&Nz,0,0,0,&nc,0,0,0);
373: DMDAGetCorners(da,0,0,0,&nx,&ny,&nz);
374: A->rmap->n = A->cmap->n = nc*nx*ny*nz;
375: A->rmap->N = A->cmap->N = nc*Nx*Ny*Nz;
376: DMCreateLocalVector(da,&a->localu);
377: return(0);
378: }
379: EXTERN_C_END
381: EXTERN_C_BEGIN
384: PetscErrorCode MatDAADSetSNES_AD(Mat A,SNES snes)
385: {
386: Mat_DAAD *a = (Mat_DAAD*)A->data;
389: a->snes = snes;
390: return(0);
391: }
392: EXTERN_C_END
394: EXTERN_C_BEGIN
397: PetscErrorCode MatDAADSetCtx_AD(Mat A,void *ctx)
398: {
399: Mat_DAAD *a = (Mat_DAAD*)A->data;
402: a->ctx = ctx;
403: return(0);
404: }
405: EXTERN_C_END
407: /*MC
408: MATDAAD - MATDAAD = "daad" - A matrix type that can do matrix-vector products using a local function that
409: is differentiated with ADIFOR or ADIC.
411: Level: intermediate
413: .seealso: MatCreateDAAD
414: M*/
416: EXTERN_C_BEGIN
419: PetscErrorCode MatCreate_DAAD(Mat B)
420: {
421: Mat_DAAD *b;
425: PetscNewLog(B,Mat_DAAD,&b);
426: B->data = (void*)b;
427: PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
428:
429: PetscLayoutSetUp(B->rmap);
430: PetscLayoutSetUp(B->cmap);
432: PetscObjectChangeTypeName((PetscObject)B,MATDAAD);
433: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMFFDSetBase_C","MatMFFDSetBase_AD",MatMFFDSetBase_AD);
434: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatDAADSetDA_C","MatDAADSetDA_AD",MatDAADSetDA_AD);
435: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatDAADSetSNES_C","MatDAADSetSNES_AD",MatDAADSetSNES_AD);
436: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatDAADSetCtx_C","MatDAADSetCtx_AD",MatDAADSetCtx_AD);
437: PetscObjectChangeTypeName((PetscObject)B,MATDAAD);
438: return(0);
439: }
440: EXTERN_C_END
445: /*@C
446: MatDAADSetDA - Tells the matrix what DMDA it is using for layout and Jacobian.
448: Logically Collective on Mat and DMDA
450: Input Parameters:
451: + mat - the matrix
452: - da - the DMDA
454: Level: intermediate
456: .seealso: MatCreate(), DMDASetLocalAdicMFFunction(), MatCreateDAAD()
458: @*/
459: PetscErrorCode MatDAADSetDA(Mat A,DM da)
460: {
466: PetscTryMethod(A,"MatDAADSetDA_C",(Mat,void*),(A,da));
467: return(0);
468: }
472: /*@C
473: MatDAADSetSNES - Tells the matrix what SNES it is using for the base U.
475: Logically Collective on Mat and SNES
477: Input Parameters:
478: + mat - the matrix
479: - snes - the SNES
481: Level: intermediate
483: Notes: this is currently turned off for Fortran usage
485: .seealso: MatCreate(), DMDASetLocalAdicMFFunction(), MatCreateDAAD(), MatDAADSetDA()
487: @*/
488: PetscErrorCode MatDAADSetSNES(Mat A,SNES snes)
489: {
495: PetscTryMethod(A,"MatDAADSetSNES_C",(Mat,void*),(A,snes));
496: return(0);
497: }
501: /*@C
502: MatDAADSetCtx - Sets the user context for a DMDAAD (ADIC matrix-free) matrix.
504: Logically Collective on Mat
506: Input Parameters:
507: + mat - the matrix
508: - ctx - the context
510: Level: intermediate
512: .seealso: MatCreate(), DMDASetLocalAdicMFFunction(), MatCreateDAAD(), MatDAADSetDA()
514: @*/
515: PetscErrorCode MatDAADSetCtx(Mat A,void *ctx)
516: {
521: PetscTryMethod(A,"MatDAADSetCtx_C",(Mat,void*),(A,ctx));
522: return(0);
523: }
527: /*@C
528: MatCreateDAAD - Creates a matrix that can do matrix-vector products using a local
529: function that is differentiated with ADIFOR or ADIC.
531: Collective on DMDA
533: Input Parameters:
534: . da - the DMDA that defines the distribution of the vectors
536: Output Parameter:
537: . A - the matrix
539: Level: intermediate
541: Notes: this is currently turned off for Fortran
543: .seealso: MatCreate(), DMDASetLocalAdicMFFunction()
545: @*/
546: PetscErrorCode MatCreateDAAD(DM da,Mat *A)
547: {
549: MPI_Comm comm;
552: PetscObjectGetComm((PetscObject)da,&comm);
553: MatCreate(comm,A);
554: MatSetType(*A,MATDAAD);
555: MatDAADSetDA(*A,da);
556: return(0);
557: }
562: /*@
563: MatRegisterDAAD - Registers DMDAAD matrix type
565: Level: advanced
567: .seealso: MatCreateDAAD(), DMDASetLocalAdicMFFunction()
569: @*/
570: PetscErrorCode MatRegisterDAAD(void)
571: {
574: MatRegisterDynamic(MATDAAD,PETSC_NULL,"MatCreate_DAAD",MatCreate_DAAD);
575: return(0);
576: }