Actual source code: matadic.c
1: #define PETSCMAT_DLL
3: /*
4: ADIC matrix-free matrix implementation
5: */
7: #include 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,"MatSNESMFSetBase_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 MatSNESMFSetBase_AD(Mat J,Vec U)
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: int nc,nx,ny,nz,Nx,Ny,Nz;
371: a->da = da;
372: PetscObjectReference((PetscObject)da);
373: DAGetInfo(da,0,&Nx,&Ny,&Nz,0,0,0,&nc,0,0,0);
374: DAGetCorners(da,0,0,0,&nx,&ny,&nz);
375: A->rmap.n = A->cmap.n = nc*nx*ny*nz;
376: A->rmap.N = A->cmap.N = nc*Nx*Ny*Nz;
377: DACreateLocalVector(da,&a->localu);
378: return(0);
379: }
385: PetscErrorCode MatDAADSetSNES_AD(Mat A,SNES snes)
386: {
387: Mat_DAAD *a = (Mat_DAAD*)A->data;
390: a->snes = snes;
391: return(0);
392: }
398: PetscErrorCode MatDAADSetCtx_AD(Mat A,void *ctx)
399: {
400: Mat_DAAD *a = (Mat_DAAD*)A->data;
403: a->ctx = ctx;
404: return(0);
405: }
408: /*MC
409: MATDAAD - MATDAAD = "daad" - A matrix type that can do matrix-vector products using a local function that
410: is differentiated with ADIFOR or ADIC.
412: Level: intermediate
414: .seealso: MatCreateDAAD
415: M*/
420: PetscErrorCode MatCreate_DAAD(Mat B)
421: {
422: Mat_DAAD *b;
426: PetscNew(Mat_DAAD,&b);
427: B->data = (void*)b;
428: PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
429:
430: PetscMapInitialize(B->comm,&B->rmap);
431: PetscMapInitialize(B->comm,&B->cmap);
433: PetscObjectChangeTypeName((PetscObject)B,MATDAAD);
434: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSNESMFSetBase_C","MatSNESMFSetBase_AD",MatSNESMFSetBase_AD);
435: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatDAADSetDA_C","MatDAADSetDA_AD",MatDAADSetDA_AD);
436: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatDAADSetSNES_C","MatDAADSetSNES_AD",MatDAADSetSNES_AD);
437: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatDAADSetCtx_C","MatDAADSetCtx_AD",MatDAADSetCtx_AD);
438: PetscObjectChangeTypeName((PetscObject)B,MATDAAD);
439: return(0);
440: }
446: /*@C
447: MatDAADSetDA - Tells the matrix what DA it is using for layout and Jacobian.
449: Collective on Mat and DA
451: Input Parameters:
452: + mat - the matrix
453: - da - the DA
455: Level: intermediate
457: .seealso: MatCreate(), DASetLocalAdicMFFunction(), MatCreateDAAD()
459: @*/
460: PetscErrorCode MatDAADSetDA(Mat A,DA da)
461: {
462: PetscErrorCode ierr,(*f)(Mat,void*);
467: PetscObjectQueryFunction((PetscObject)A,"MatDAADSetDA_C",(void (**)(void))&f);
468: if (f) {
469: (*f)(A,da);
470: }
471: return(0);
472: }
476: /*@
477: MatDAADSetSNES - Tells the matrix what SNES it is using for the base U.
479: Collective on Mat and SNES
481: Input Parameters:
482: + mat - the matrix
483: - snes - the SNES
485: Level: intermediate
487: .seealso: MatCreate(), DASetLocalAdicMFFunction(), MatCreateDAAD(), MatDAADSetDA()
489: @*/
490: PetscErrorCode MatDAADSetSNES(Mat A,SNES snes)
491: {
492: PetscErrorCode ierr,(*f)(Mat,void*);
497: PetscObjectQueryFunction((PetscObject)A,"MatDAADSetSNES_C",(void (**)(void))&f);
498: if (f) {
499: (*f)(A,snes);
500: }
501: return(0);
502: }
506: /*@C
507: MatDAADSetCtx - Sets the user context for a DAAD (ADIC matrix-free) matrix.
509: Collective on Mat
511: Input Parameters:
512: + mat - the matrix
513: - ctx - the context
515: Level: intermediate
517: .seealso: MatCreate(), DASetLocalAdicMFFunction(), MatCreateDAAD(), MatDAADSetDA()
519: @*/
520: PetscErrorCode MatDAADSetCtx(Mat A,void *ctx)
521: {
522: PetscErrorCode ierr,(*f)(Mat,void*);
526: PetscObjectQueryFunction((PetscObject)A,"MatDAADSetCtx_C",(void (**)(void))&f);
527: if (f) {
528: (*f)(A,ctx);
529: }
530: return(0);
531: }
535: /*@
536: MatCreateDAAD - Creates a matrix that can do matrix-vector products using a local
537: function that is differentiated with ADIFOR or ADIC.
539: Collective on DA
541: Input Parameters:
542: . da - the DA that defines the distribution of the vectors
544: Output Parameter:
545: . A - the matrix
547: Level: intermediate
549: .seealso: MatCreate(), DASetLocalAdicMFFunction()
551: @*/
552: PetscErrorCode MatCreateDAAD(DA da,Mat *A)
553: {
555: MPI_Comm comm;
558: PetscObjectGetComm((PetscObject)da,&comm);
559: MatCreate(comm,A);
560: MatSetType(*A,MATDAAD);
561: MatDAADSetDA(*A,da);
562: return(0);
563: }
568: /*@
569: MatRegisterDAAD - Registers DAAD matrix type
571: Level: advanced
573: .seealso: MatCreateDAAD(), DASetLocalAdicMFFunction()
575: @*/
576: PetscErrorCode MatRegisterDAAD(void)
577: {
580: MatRegisterDynamic(MATDAAD,PETSC_NULL,"MatCreate_DAAD",MatCreate_DAAD);
581: return(0);
582: }