Actual source code: matadic.c
1: /*$Id: matadic.c,v 1.13 2001/08/07 03:03:06 balay Exp $*/
2: /*
3: ADIC matrix-free matrix implementation
4: */
6: #include src/mat/matimpl.h
7: #include petscda.h
8: #include petscsnes.h
9: #include petscsys.h
10: EXTERN_C_BEGIN
11: #include "adic/ad_utils.h"
12: EXTERN_C_END
14: typedef struct {
15: DA da;
16: Vec localu; /* point at which Jacobian is evaluated */
17: void *ctx;
18: SNES snes;
19: Vec diagonal; /* current matrix diagonal */
20: PetscTruth diagonalvalid; /* indicates if diagonal matches current base vector */
21: } Mat_DAAD;
23: #undef __FUNCT__
25: int MatAssemblyEnd_DAAD(Mat A,MatAssemblyType atype)
26: {
27: Mat_DAAD *a = (Mat_DAAD*)A->data;
28: int ierr;
29: Vec u;
32: a->diagonalvalid = PETSC_FALSE;
33: if (a->snes) {
34: SNESGetSolution(a->snes,&u);
35: DAGlobalToLocalBegin(a->da,u,INSERT_VALUES,a->localu);
36: DAGlobalToLocalEnd(a->da,u,INSERT_VALUES,a->localu);
37: }
38: return(0);
39: }
41: #undef __FUNCT__
43: int MatMult_DAAD(Mat A,Vec xx,Vec yy)
44: {
45: Mat_DAAD *a = (Mat_DAAD*)A->data;
46: Vec localxx;
47: int ierr;
50: DAGetLocalVector(a->da,&localxx);
51: DAGlobalToLocalBegin(a->da,xx,INSERT_VALUES,localxx);
52: DAGlobalToLocalEnd(a->da,xx,INSERT_VALUES,localxx);
53: DAMultiplyByJacobian1WithAD(a->da,a->localu,localxx,yy,a->ctx);
54: DARestoreLocalVector(a->da,&localxx);
55: return(0);
56: }
58: #include src/dm/da/daimpl.h
60: #undef __FUNCT__
62: int MatGetDiagonal_DAAD(Mat A,Vec dd)
63: {
64: Mat_DAAD *a = (Mat_DAAD*)A->data;
65: int ierr,j,nI,gI,gtdof;
66: PetscScalar *avu,*ad_vustart,ad_f[2],*d;
67: DALocalInfo info;
68: MatStencil stencil;
69: void* *ad_vu;
73: /* get space for derivative object. */
74: DAGetAdicMFArray(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: DAGetLocalInfo(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.xm + (stencil.k - info.gzs)*info.dof*info.xm*info.ym;
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: DARestoreAdicMFArray(a->da,PETSC_TRUE,(void **)&ad_vu,(void**)&ad_vustart,>dof);
109: return(0);
110: }
113: #undef __FUNCT__
115: int MatRelax_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;
118: int ierr,j,gtdof,nI,gI;
119: PetscScalar *avu,*av,*ad_vustart,ad_f[2],zero = 0.0,*d,*b;
120: Vec localxx,dd;
121: DALocalInfo info;
122: MatStencil stencil;
123: void* *ad_vu;
126: if (its <= 0 || lits <= 0) SETERRQ2(PETSC_ERR_ARG_WRONG,"Relaxation requires global its %d and local its %d both positive",its,lits);
127: if (!a->diagonal) {
128: DACreateGlobalVector(a->da,&a->diagonal);
129: }
130: if (!a->diagonalvalid) {
131: ierr = MatGetDiagonal(A,a->diagonal);
132: a->diagonalvalid = PETSC_TRUE;
133: }
134: dd = a->diagonal;
137: DAGetLocalVector(a->da,&localxx);
138: if (flag & SOR_ZERO_INITIAL_GUESS) {
139: VecSet(&zero,localxx);
140: } else {
141: DAGlobalToLocalBegin(a->da,xx,INSERT_VALUES,localxx);
142: DAGlobalToLocalEnd(a->da,xx,INSERT_VALUES,localxx);
143: }
145: /* get space for derivative object. */
146: DAGetAdicMFArray(a->da,PETSC_TRUE,(void **)&ad_vu,(void**)&ad_vustart,>dof);
148: /* copy input vector into derivative object */
149: VecGetArray(a->localu,&avu);
150: VecGetArray(localxx,&av);
151: for (j=0; j<gtdof; j++) {
152: ad_vustart[2*j] = avu[j];
153: ad_vustart[2*j+1] = av[j];
154: }
155: VecRestoreArray(a->localu,&avu);
156: VecRestoreArray(localxx,&av);
158: PetscADResetIndep();
159: PetscADIncrementTotalGradSize(1);
160: PetscADSetIndepDone();
162: VecGetArray(dd,&d);
163: VecGetArray(bb,&b);
165: DAGetLocalInfo(a->da,&info);
166: while (its--) {
167: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){
168: nI = 0;
169: for (stencil.k = info.zs; stencil.k<info.zs+info.zm; stencil.k++) {
170: for (stencil.j = info.ys; stencil.j<info.ys+info.ym; stencil.j++) {
171: for (stencil.i = info.xs; stencil.i<info.xs+info.xm; stencil.i++) {
172: for (stencil.c = 0; stencil.c<info.dof; stencil.c++) {
173: (*a->da->adicmf_lfi)(&info,&stencil,ad_vu,ad_f,a->ctx);
174: gI = stencil.c + (stencil.i - info.gxs)*info.dof + (stencil.j - info.gys)*info.dof*info.xm + (stencil.k - info.gzs)*info.dof*info.xm*info.ym;
175: ad_vustart[1+2*gI] += (b[nI] - ad_f[1])/d[nI];
176: nI++;
177: }
178: }
179: }
180: }
181: }
182: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){
183: nI = info.dof*info.xm*info.ym*info.zm - 1;
184: for (stencil.k = info.zs+info.zm-1; stencil.k>=info.zs; stencil.k--) {
185: for (stencil.j = info.ys+info.ym-1; stencil.j>=info.ys; stencil.j--) {
186: for (stencil.i = info.xs+info.xm-1; stencil.i>=info.xs; stencil.i--) {
187: for (stencil.c = info.dof-1; stencil.c>=0; stencil.c--) {
188: (*a->da->adicmf_lfi)(&info,&stencil,ad_vu,ad_f,a->ctx);
189: gI = stencil.c + (stencil.i - info.gxs)*info.dof + (stencil.j - info.gys)*info.dof*info.xm + (stencil.k - info.gzs)*info.dof*info.xm*info.ym;
190: ad_vustart[1+2*gI] += (b[nI] - ad_f[1])/d[nI];
191: nI--;
192: }
193: }
194: }
195: }
196: }
197: }
199: VecRestoreArray(dd,&d);
200: VecRestoreArray(bb,&b);
202: VecGetArray(localxx,&av);
203: for (j=0; j<gtdof; j++) {
204: av[j] = ad_vustart[2*j+1];
205: }
206: VecRestoreArray(localxx,&av);
207: DALocalToGlobal(a->da,localxx,INSERT_VALUES,xx);
209: DARestoreLocalVector(a->da,&localxx);
210: DARestoreAdicMFArray(a->da,PETSC_TRUE,(void **)&ad_vu,(void**)&ad_vustart,>dof);
211: return(0);
212: }
216: #undef __FUNCT__
218: int MatDestroy_DAAD(Mat A)
219: {
220: Mat_DAAD *a = (Mat_DAAD*)A->data;
221: int ierr;
224: DADestroy(a->da);
225: VecDestroy(a->localu);
226: if (a->diagonal) {VecDestroy(a->diagonal);}
227: PetscFree(a);
228: return(0);
229: }
231: /* -------------------------------------------------------------------*/
232: static struct _MatOps MatOps_Values = {0,
233: 0,
234: 0,
235: MatMult_DAAD,
236: 0,
237: 0,
238: 0,
239: 0,
240: 0,
241: 0,
242: 0,
243: 0,
244: 0,
245: MatRelax_DAAD,
246: 0,
247: 0,
248: 0,
249: MatGetDiagonal_DAAD,
250: 0,
251: 0,
252: 0,
253: MatAssemblyEnd_DAAD,
254: 0,
255: 0,
256: 0,
257: 0,
258: 0,
259: 0,
260: 0,
261: 0,
262: 0,
263: 0,
264: 0,
265: 0,
266: 0,
267: 0,
268: 0,
269: 0,
270: 0,
271: 0,
272: 0,
273: 0,
274: 0,
275: 0,
276: 0,
277: 0,
278: 0,
279: 0,
280: 0,
281: 0,
282: 0,
283: 0,
284: 0,
285: 0,
286: 0,
287: 0,
288: 0,
289: 0,
290: 0,
291: 0,
292: 0,
293: MatDestroy_DAAD,
294: 0,
295: 0,
296: 0,
297: 0,
298: 0,
299: 0,
300: 0,
301: 0,
302: 0,
303: 0,
304: 0,
305: 0,
306: 0};
308: /* --------------------------------------------------------------------------------*/
310: EXTERN_C_BEGIN
311: #undef __FUNCT__
313: int MatSNESMFSetBase_AD(Mat J,Vec U)
314: {
315: int ierr;
316: Mat_DAAD *a = (Mat_DAAD*)J->data;
319: a->diagonalvalid = PETSC_FALSE;
320: DAGlobalToLocalBegin(a->da,U,INSERT_VALUES,a->localu);
321: DAGlobalToLocalEnd(a->da,U,INSERT_VALUES,a->localu);
322: return(0);
323: }
324: EXTERN_C_END
326: EXTERN_C_BEGIN
327: #undef __FUNCT__
329: int MatDAADSetDA_AD(Mat A,DA da)
330: {
331: Mat_DAAD *a = (Mat_DAAD*)A->data;
332: int ierr,nc,nx,ny,nz,Nx,Ny,Nz;
335: a->da = da;
336: ierr = PetscObjectReference((PetscObject)da);
337: ierr = DAGetInfo(da,0,&Nx,&Ny,&Nz,0,0,0,&nc,0,0,0);
338: ierr = DAGetCorners(da,0,0,0,&nx,&ny,&nz);
339: A->m = A->n = nc*nx*ny*nz;
340: A->M = A->N = nc*Nx*Ny*Nz;
341: ierr = DACreateLocalVector(da,&a->localu);
342: return(0);
343: }
344: EXTERN_C_END
346: EXTERN_C_BEGIN
347: #undef __FUNCT__
349: int MatDAADSetSNES_AD(Mat A,SNES snes)
350: {
351: Mat_DAAD *a = (Mat_DAAD*)A->data;
354: a->snes = snes;
355: return(0);
356: }
357: EXTERN_C_END
359: EXTERN_C_BEGIN
360: #undef __FUNCT__
362: int MatDAADSetCtx_AD(Mat A,void *ctx)
363: {
364: Mat_DAAD *a = (Mat_DAAD*)A->data;
367: a->ctx = ctx;
368: return(0);
369: }
370: EXTERN_C_END
372: EXTERN_C_BEGIN
373: #undef __FUNCT__
375: int MatCreate_DAAD(Mat B)
376: {
377: Mat_DAAD *b;
378: int ierr;
381: ierr = PetscNew(Mat_DAAD,&b);
382: B->data = (void*)b;
383: PetscMemzero(b,sizeof(Mat_DAAD));
384: PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));
385:
386: PetscMapCreateMPI(B->comm,B->m,B->m,&B->rmap);
387: PetscMapCreateMPI(B->comm,B->n,B->n,&B->cmap);
389: PetscObjectChangeTypeName((PetscObject)B,MATDAAD);
390: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSNESMFSetBase_C","MatSNESMFSetBase_AD",MatSNESMFSetBase_AD);
391: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatDAADSetDA_C","MatDAADSetDA_AD",MatDAADSetDA_AD);
392: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatDAADSetSNES_C","MatDAADSetSNES_AD",MatDAADSetSNES_AD);
393: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatDAADSetCtx_C","MatDAADSetCtx_AD",MatDAADSetCtx_AD);
394: return(0);
395: }
396: EXTERN_C_END
399: #undef __FUNCT__
401: /*@C
402: MatDAADSetDA - Tells the matrix what DA it is using for layout and Jacobian.
404: Collective on Mat and DA
406: Input Parameters:
407: + mat - the matrix
408: - da - the DA
410: Level: intermediate
412: .seealso: MatCreate(), DASetLocalAdicMFFunction(), MatCreateDAAD()
414: @*/
415: int MatDAADSetDA(Mat A,DA da)
416: {
417: int ierr,(*f)(Mat,void*);
422: PetscObjectQueryFunction((PetscObject)A,"MatDAADSetDA_C",(void (**)(void))&f);
423: if (f) {
424: (*f)(A,da);
425: }
426: return(0);
427: }
429: #undef __FUNCT__
431: /*@C
432: MatDAADSetSNES - Tells the matrix what SNES it is using for the base U.
434: Collective on Mat and SNES
436: Input Parameters:
437: + mat - the matrix
438: - snes - the SNES
440: Level: intermediate
442: .seealso: MatCreate(), DASetLocalAdicMFFunction(), MatCreateDAAD(), MatDAADSetDA()
444: @*/
445: int MatDAADSetSNES(Mat A,SNES snes)
446: {
447: int ierr,(*f)(Mat,void*);
452: PetscObjectQueryFunction((PetscObject)A,"MatDAADSetSNES_C",(void (**)(void))&f);
453: if (f) {
454: (*f)(A,snes);
455: }
456: return(0);
457: }
459: #undef __FUNCT__
461: /*@C
462: MatDAADSetCtx - Sets the user context for a DAAD (ADIC matrix-free) matrix.
464: Collective on Mat
466: Input Parameters:
467: + mat - the matrix
468: - ctx - the context
470: Level: intermediate
472: .seealso: MatCreate(), DASetLocalAdicMFFunction(), MatCreateDAAD(), MatDAADSetDA()
474: @*/
475: int MatDAADSetCtx(Mat A,void *ctx)
476: {
477: int ierr,(*f)(Mat,void*);
481: PetscObjectQueryFunction((PetscObject)A,"MatDAADSetCtx_C",(void (**)(void))&f);
482: if (f) {
483: (*f)(A,ctx);
484: }
485: return(0);
486: }
488: #undef __FUNCT__
490: /*@C
491: MatCreateDAAD - Creates a matrix that can do matrix-vector products using a local
492: function that is differentiated with ADIFOR or ADIC.
494: Collective on DA
496: Input Parameters:
497: . da - the DA that defines the distribution of the vectors
499: Output Parameter:
500: . A - the matrix
502: Level: intermediate
504: .seealso: MatCreate(), DASetLocalAdicMFFunction()
506: @*/
507: int MatCreateDAAD(DA da,Mat *A)
508: {
509: int ierr;
510: MPI_Comm comm;
513: PetscObjectGetComm((PetscObject)da,&comm);
514: MatCreate(comm,0,0,0,0,A);
515: MatSetType(*A,MATDAAD);
516: MatDAADSetDA(*A,da);
517: return(0);
518: }
520: #undef __FUNCT__
522: int MatRegisterDAAD(void)
523: {
526: MatRegisterDynamic(MATDAAD,PETSC_NULL,"MatCreate_DAAD",MatCreate_DAAD);
527: return(0);
528: }