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,&gtdof);

 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,&gtdof);
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,&gtdof);

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,&gtdof);
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: }