Actual source code: taovec_petsc.c

  1: #include "tao_general.h"

  3: #ifdef TAO_USE_PETSC

  5: #include "taovec_petsc.h"
  6: #include "../indexset/taois_petsc.h"

  8: int VecMedian(Vec, Vec, Vec, Vec);
  9: int VecPointwiseMax(Vec, Vec, Vec);
 10: int VecPointwiseMin(Vec, Vec, Vec);
 11: int VecFischer(Vec, Vec, Vec, Vec, Vec);
 12: int VecSFischer(Vec, Vec, Vec, Vec, double, Vec);
 13: int VecBoundProjection(Vec, Vec, Vec, Vec, Vec);
 14: int VecCreateSubVec(Vec, IS, Vec *);
 15: int VecISAXPY(Vec, double, Vec, IS);
 16: int VecStepMax(Vec, Vec, double*);
 17: int VecStepMax2(Vec, Vec,Vec, Vec, double*);
 18: int VecCompare(Vec, Vec, PetscTruth *);
 19: int VecBoundStepInfo(Vec, Vec, Vec, Vec, PetscReal *, PetscReal *, PetscReal *);

 21: static int petscminmaxevent=0;

 25: /*@C
 26:    TaoWrapPetscVec - Creates a new TaoVec object using an existing
 27:    PETSc Vec.

 29:    Input Parameter:
 30: +  V -  a PETSc vector
 31: -  TV -  the address of a pointer to a TaoVecPetsc

 33:    Output Parameter:
 34: .  TV - pointer to a new TaoVecPetsc

 36:    Note:  A TaoVecPetsc is an object with the methods of an abstract
 37:    TaoVec object.  A TaoVecPetsc contains an implementation of the TaoVec
 38:    methods.  Routines using these vectors should declare a pointer to 
 39:    a TaoVec, assign this pointer to the address of a TaoVec object, 
 40:    use the pointer to invoke methods on the object, and use this pointer
 41:    as an argument when calling other routines.  This usage is different
 42:    from the usage of a PETSc Vec.  In PETSc, applications will typically
 43:    declare a Vec, and pass it as an argument into routines.  That is,
 44:    applications will typically declare a pointer to a TaoVec and use the
 45:    pointer, or declare a Vec and use it directly.
 46:   
 47:    Note:
 48:    The user us repsonsible for destroying the Vec V, in addition to
 49:    to TaoVecPetsc vector TV.  The Vec can be destroyed immediately
 50:    after this routine.

 52:    Level: developer

 54: .seealso TaoVecGetPetscVec(), TaoVecDestroy()

 56: @*/
 57: int TaoWrapPetscVec( Vec V, TaoVecPetsc* *TV){
 58:   int info;
 60:   *TV = new  TaoVecPetsc(V);
 61:   info=(*TV)->SetVec(V); CHKERRQ(info);
 62:   return(0);
 63: }

 67: /*@C
 68:    TaoVecGetPetscVec - Set the address of a Vec equal to the location
 69:    of the underlying Vec structure in this TaoVecPetsc object. 

 71:    Input Parameter:
 72: +  TV - the TaoVecPetsc 
 73: -  V -  the address of Vec

 75:    Output Parameter:
 76: .  V -  the address of a specific Vec

 78:    Note:  
 79:    This routine does not create a Vec.  It sets a pointer
 80:    to the location of an existing Vec.

 82:    Note:
 83:    The underlying PETSc Vec may also be obtained by
 84:    Vec VV = (Vec) TV->VecObject;

 86:    Note:
 87:    
 88:    VecObject is a void* pointer that can be cast to the underlying Vec object.

 90:    Level: advanced

 92: @*/
 93: int TaoVecGetPetscVec( TaoVec* TV, Vec *V){
 95:   if (TV && V){
 96:     *V=(Vec)(TV->VecObject);
 97:   } else if (V){
 98:     *V=0;
 99:   }
100:   return(0);
101: }

105: TaoVecPetsc::TaoVecPetsc(Vec VV){
106:   this->pv=0;
107:   this->pvecviewer=0;
108:   if (VV){
109:     this->SetVec(VV);
110:   }
111:   if (petscminmaxevent==0){
112:     PetscLogEventRegister(&petscminmaxevent,"PointwiseMinMax",VEC_COOKIE);
113:   }
114:   return;
115: }

119: TaoVecPetsc::~TaoVecPetsc(){
120:   this->SetVec(0);
121:   return;
122: }

126: /*@C
127:    SetVec - Set or replace the underlying Vec object in this TaoVec.

129:    Input Parameter:
130: .  V -  a PETSc vector

132:    Note:
133:    The user us repsonsible for destroying the Vec V, in addition to
134:    to TaoVecPetsc vector TV.  The Vec can be destroyed immediately
135:    after this routine.

137:    Level: developer

139: .seealso TaoWrapPetscVec(), TaoVecDestroy()

141: @*/
142: int TaoVecPetsc::SetVec(Vec VV){
143:   int info;
145:   if (VV){
147:     PetscObjectReference((PetscObject)VV);
148:     //    info = PetscLogInfo((VV,"Set the PETSc Vec in a TaoVec .\n")); CHKERRQ(info);
149:   }
150:   if (this->pv){
151:     info=VecDestroy(this->pv); CHKERRQ(info);
152:   }
153:   this->pv=VV;
154:   this->VecObject=(void*)VV;
155:   return(0);
156: }

160: int TaoVecPetsc::Clone(TaoVec**ntv){
161:   int info;
162:   TaoVecPetsc *nptv;
163:   Vec npv;
165:   info=VecDuplicate(pv,&npv); CHKERRQ(info);
166:   info = TaoWrapPetscVec(npv,&nptv); CHKERRQ(info);
167:   *ntv=nptv;
168:   info = VecDestroy(npv); CHKERRQ(info);
169:   nptv->pvecviewer=this->pvecviewer;
170:   return(0);
171: }


176: int TaoVecPetsc::Compatible(TaoVec *tv, TaoTruth *flag){
177:   TaoVecPetsc *vv = (TaoVecPetsc *)(tv);
178:   PetscTruth flg=PETSC_TRUE;
179:   int info=0;

182:   if (vv==0){
183:     *flag=TAO_FALSE;
184:     return(0);
185:   }
186:   info = VecCompare(pv, vv->pv, &flg); CHKERRQ(info);
187:   if (flg==PETSC_FALSE){
188:     *flag=TAO_FALSE;
189:   } else {
190:     *flag=TAO_TRUE;
191:   }
192:   return(0);
193: }

197: int TaoVecPetsc::SetToConstant( double cc ){
198:   int info;
199:   PetscScalar c=cc;
201:   info=VecSet(pv, c);  CHKERRQ(info);
202:   return(0);
203: }

207: int TaoVecPetsc::SetToZero(){
208:   PetscScalar zero=0.0;
209:   int info;
211:   info=VecSet(pv, zero);  CHKERRQ(info);
212:   return(0);
213: }

217: int TaoVecPetsc::CopyFrom( TaoVec* tv ){
218:   Vec pv2 = (Vec)(tv->VecObject);
219:   int info;
221:   info=VecCopy(pv2,pv);  CHKERRQ(info);
222:   return(0);
223: }

227: int TaoVecPetsc::ScaleCopyFrom( double a, TaoVec* tv ){
228:   Vec v = (Vec) tv->VecObject;
229:   PetscScalar alpha=a, zero=0.0;
230:   int info;
232:   info=VecAXPBY(this->pv,alpha,zero,v);  CHKERRQ(info);
233:   return(0);
234: }

238: int TaoVecPetsc::NormInfinity(double *vnorm){
239:   int info;
240:   PetscReal vv;
242:   info=VecNorm(pv,NORM_INFINITY,&vv);  CHKERRQ(info);
243:   *vnorm=(double)vv;
244:   return(0);
245: }

249: int TaoVecPetsc::Norm1(double *vnorm){
250:   int info;
251:   PetscReal vv;
253:   info=VecNorm(pv,NORM_1,&vv);  CHKERRQ(info);
254:   *vnorm=(double)vv;
255:   return(0);
256: }

260: int TaoVecPetsc::Norm2(double *vnorm){
261:   int info;
262:   PetscReal vv;
264:   info=VecNorm(pv,NORM_2,&vv);  CHKERRQ(info);
265:   *vnorm=(double)vv;
266:   return(0);
267: }

271: int TaoVecPetsc::Norm2squared(double *vnorm){
272:   int info;
273:   PetscReal vv;
275:   info=VecNorm(pv,NORM_2,&vv);  CHKERRQ(info);
276:   *vnorm=(double)(vv*vv);
277:   return(0);
278: }

282: int TaoVecPetsc::Scale( double alpha ){
283:   int info;
284:   PetscScalar aalpha=alpha;
286:   info=VecScale(pv, aalpha);  CHKERRQ(info);
287:   return(0);
288: }


293: int TaoVecPetsc::Axpy( double alpha, TaoVec* tv ){
294:   Vec v = (Vec) tv->VecObject;
295:   PetscScalar aalpha=alpha;
296:   int info;
298:   info=VecAXPY(this->pv, aalpha, v);  CHKERRQ(info);
299:   return(0);
300: }


305: int TaoVecPetsc::Axpby( double alpha, TaoVec* tv, double beta ){
306:   Vec v = (Vec) tv->VecObject;
307:   PetscScalar aalpha = alpha, bbeta = beta;
308:   int info;
310:   info=VecAXPBY(this->pv,aalpha,bbeta,v);  CHKERRQ(info);
311:   return(0);
312: }

316: int TaoVecPetsc::Aypx( double alpha, TaoVec* x ){
317:   Vec X = (Vec) x->VecObject;
318:   PetscScalar aalpha = alpha;
319:   int info;
321:   info=VecAYPX(this->pv,aalpha,X);  CHKERRQ(info);
322:   return(0);
323: }
324: 
327: int TaoVecPetsc::AddConstant( double cc ){
328:   int info;
329:   PetscScalar c=cc;
331:   info=VecShift(pv, c);  CHKERRQ(info);
332:   return(0);
333: }

337: int TaoVecPetsc::Dot( TaoVec* tv, double *vDotv ){
338:   Vec v = (Vec) tv->VecObject;
339:   PetscScalar c;
340:   int info;
342:   info=VecDot(v,this->pv,&c);  CHKERRQ(info);
343:   *vDotv=c;
344:   return(0);
345: }

349: int TaoVecPetsc::Negate(){
350:   PetscScalar m1=-1.0;
351:   int info;
353:   info=VecScale(pv, m1);  CHKERRQ(info);
354:   return(0);
355: }

359: int TaoVecPetsc::Reciprocal(){
360:   int info;
362:   info=VecReciprocal(pv);  CHKERRQ(info);
363:   return(0);
364: }

368: int TaoVecPetsc::GetDimension(int *n){
369:   int info;
371:   info=VecGetSize(pv,n);  CHKERRQ(info);
372:   return(0);
373: }

377: int TaoVecPetsc::PointwiseMultiply( TaoVec* tv, TaoVec* tw ){
378:   Vec v = (Vec) tv->VecObject;
379:   Vec w = (Vec) tw->VecObject;
380:   int info;
382:   info=VecPointwiseMult(this->pv, w, v);  CHKERRQ(info);
383:   return(0);
384: }


389: int TaoVecPetsc::PointwiseDivide( TaoVec* tv , TaoVec* tw){
390:   Vec v = (Vec) tv->VecObject;
391:   Vec w = (Vec) tw->VecObject;
392:   int info;
394:   info=VecPointwiseDivide(this->pv, v, w);  CHKERRQ(info);
395:   return(0);
396: }

400: int TaoVecPetsc::Median( TaoVec* tv, TaoVec* tw, TaoVec* tx){
401:   Vec v = (Vec) tv->VecObject;
402:   Vec w = (Vec) tw->VecObject;
403:   Vec x = (Vec) tx->VecObject;
404:   int info;
406:   info=PetscLogEventBegin(petscminmaxevent,0,0,0,0); CHKERRQ(info);
407:   info=VecMedian(v,w,x,this->pv); CHKERRQ(info);
408:   info=PetscLogEventEnd(petscminmaxevent,0,0,0,0); CHKERRQ(info);
409:   return(0);
410: }

414: int TaoVecPetsc::PointwiseMinimum( TaoVec* tv, TaoVec* tw){
415:   Vec v = (Vec) tv->VecObject;
416:   Vec w = (Vec) tw->VecObject;
417:   int info;
419:   info=PetscLogEventBegin(petscminmaxevent,0,0,0,0); CHKERRQ(info);
420:   info=VecPointwiseMin(this->pv, w, v); CHKERRQ(info);
421:   info=PetscLogEventEnd(petscminmaxevent,0,0,0,0); CHKERRQ(info);
422:   return(0);
423: }

427: int TaoVecPetsc::Fischer(TaoVec *tx, TaoVec *tf, TaoVec *tl, TaoVec *tu)
428: {
429:   TaoVecPetsc *xx = (TaoVecPetsc *)(tx);
430:   TaoVecPetsc *ff = (TaoVecPetsc *)(tf);
431:   TaoVecPetsc *ll = (TaoVecPetsc *)(tl);
432:   TaoVecPetsc *uu = (TaoVecPetsc *)(tu);
433:   int info;

436:   info=VecFischer(xx->pv, ff->pv, ll->pv, uu->pv, pv); CHKERRQ(info);
437:   return(0);
438: }

442: int TaoVecPetsc::SFischer(TaoVec *tx, TaoVec *tf,
443:                           TaoVec *tl, TaoVec *tu, double mu)
444: {
445:   TaoVecPetsc *xx = (TaoVecPetsc *)(tx);
446:   TaoVecPetsc *ff = (TaoVecPetsc *)(tf);
447:   TaoVecPetsc *ll = (TaoVecPetsc *)(tl);
448:   TaoVecPetsc *uu = (TaoVecPetsc *)(tu);
449:   int info;

452:   if ((mu >= -TAO_EPSILON) && (mu <= TAO_EPSILON)) {
453:     Fischer(tx, tf, tl, tu);
454:   }
455:   else {
456:     info=VecSFischer(xx->pv, ff->pv, ll->pv, uu->pv, mu, pv); CHKERRQ(info);
457:   }
458:   return(0);
459: }

463: int TaoVecPetsc::PointwiseMaximum( TaoVec* tv, TaoVec* tw){
464:   Vec v = (Vec) tv->VecObject;
465:   Vec w = (Vec) tw->VecObject;
466:   int info;
468:   info=PetscLogEventBegin(petscminmaxevent,0,0,0,0); CHKERRQ(info);
469:   info=VecPointwiseMax(this->pv, w, v);  CHKERRQ(info);
470:   info=PetscLogEventEnd(petscminmaxevent,0,0,0,0); CHKERRQ(info);
471:   return(0);
472: }


477: int TaoVecPetsc::Waxpby  ( double aa, TaoVec* tv, double bb, TaoVec* tw){
478:   Vec v = (Vec) tv->VecObject;
479:   Vec w = (Vec) tw->VecObject;
480:   int info;
481:   PetscScalar a=aa, b=bb, zero=0.0;

484:   if (a==1){
485:     info = VecWAXPY(this->pv,b,w,v); CHKERRQ(info);
486:   } else if (b==1) {
487:     info = VecWAXPY(this->pv,a,v,w); CHKERRQ(info);
488:   } else if (a==0) {
489:     info = VecSet(pv, zero); CHKERRQ(info);
490:     info = VecAXPY(this->pv, b, w); CHKERRQ(info);
491:   } else if (b==0) {
492:     info = VecSet(pv, zero); CHKERRQ(info);
493:     info = VecAXPY(this->pv, a, v); CHKERRQ(info);
494:   } else {
495:     info = VecCopy(w,pv); CHKERRQ(info);
496:     info = VecAXPBY(this->pv,a,b,v); CHKERRQ(info);
497:   }
498:   return(0);
499: }

503: int TaoVecPetsc::AbsoluteValue(){
504:   int info;
506:   info=VecAbs(pv);  CHKERRQ(info);
507:   return(0);
508: }

512: int TaoVecPetsc::MinElement(double*val){
513:   int info, p;
514:   PetscScalar a;
516:   info = VecMin(this->pv,&p, &a);CHKERRQ(info);
517:   *val=a;
518:   return(0);
519: }


524: int TaoVecPetsc::GetArray(TaoScalar **dptr, int *n){
525:   int info;
526:   PetscScalar *pptr;
528:   if (sizeof(TaoScalar)==sizeof(PetscScalar)){
529:     info = VecGetLocalSize(pv,n);CHKERRQ(info);
530:     info = VecGetArray(pv,&pptr);CHKERRQ(info);
531:     *dptr=(TaoScalar*)pptr;
532:   } else {
533:     SETERRQ(1,"Incompatible data types");
534:   }
535:   return(0);
536: }

540: int TaoVecPetsc::RestoreArray(TaoScalar **dptr,int *n){
541:   int info;
542:   PetscScalar *pptr;
544:   if (sizeof(TaoScalar)==sizeof(PetscScalar)){
545:     pptr=(PetscScalar*)(*dptr);
546:     info = VecRestoreArray(pv,&pptr); CHKERRQ(info);
547:     *dptr=PETSC_NULL;
548:     *n=0;
549:   } else {
550:     SETERRQ(1,"Incompatible data types");
551:   }
552:   return(0);
553: }


558: int TaoVecPetsc::SetReducedVec(TaoVec*VV,TaoIndexSet*SS){
559:   int info,nn,nnn,nred;
560:   TaoPetscISType type;
561:   Vec R, V = (Vec) VV->VecObject;
562:   TaoIndexSetPetsc *TSS = (TaoIndexSetPetsc*) SS;
563:   IS S2;
564:   MPI_Comm comm;

567:   info = SS->GetSize(&nn); CHKERRQ(info);
568:   info = VV->GetDimension(&nnn); CHKERRQ(info);
569:   info = TSS->GetReducedType(&type); CHKERRQ(info);
570:   if (type==TaoMaskFullSpace){

572:   } else if (nn==nnn){
573:     info = VecDuplicate(V,&R); CHKERRQ(info);
574:     info = this->SetVec(R); CHKERRQ(info);
575:     PetscObjectDereference((PetscObject)R);
576:   } else {
577:     info = TSS->RedistributeIS(&S2); CHKERRQ(info);
578:     info = ISGetLocalSize(S2,&nred); CHKERRQ(info);
579:     info = PetscObjectGetComm((PetscObject)V,&comm);CHKERRQ(info);
580:     info = VecCreateMPI(comm,nred,nn,&R);CHKERRQ(info);
581:     info = this->SetVec(R); CHKERRQ(info);
582:     PetscObjectDereference((PetscObject)R);
583:   }
584:   info=this->ReducedCopyFromFull(VV,SS);CHKERRQ(info);

586:   return(0);
587: }

591: int TaoVecPetsc::ReducedXPY(TaoVec*VV,TaoIndexSet*SS){
592:   VecScatter scatterit;
593:   TaoIndexSetPetsc *TSS = (TaoIndexSetPetsc*) SS;
594:   Vec R = (Vec) VV->VecObject;
595:   Vec V=this->pv;
596:   int info;


600:   info = TSS->GetReducedVecScatter(R,V, &scatterit);CHKERRQ(info);
601: 
602:   info = VecScatterBegin(R,V,ADD_VALUES,SCATTER_REVERSE,scatterit);
603:   CHKERRQ(info);
604:   info = VecScatterEnd(R,V,ADD_VALUES,SCATTER_REVERSE, scatterit);
605:   CHKERRQ(info);

607:   return(0);
608: }

612: int TaoVecPetsc::ReducedCopyFromFull(TaoVec*VV,TaoIndexSet*SS){
613:   int info;
614:   PetscScalar zero=0.0;
615:   TaoIndexSetPetsc *TSS = (TaoIndexSetPetsc*) SS;
616:   VecScatter scatterit;
617:   Vec V = (Vec) VV->VecObject;
618:   Vec R=this->pv;


622:   info = VecSet(R, zero); CHKERRQ(info);
623:   info = TSS->GetReducedVecScatter(V,R, &scatterit);CHKERRQ(info);
624:   info = VecScatterBegin(V,R,INSERT_VALUES,SCATTER_FORWARD,scatterit);
625:   CHKERRQ(info);
626:   info = VecScatterEnd(V,R,INSERT_VALUES,SCATTER_FORWARD, scatterit);
627:   CHKERRQ(info);

629:   return(0);
630: }



636: int TaoVecPetsc::StepMax(TaoVec*tv,double*step){
637:   Vec dx = (Vec) tv->VecObject;
638:   PetscReal pstep;
639:   int info;
641:   info = VecStepMax(this->pv, dx,&pstep);CHKERRQ(info);
642:   *step=pstep;
643:   return(0);
644: }

648: int TaoVecPetsc::StepBoundInfo(TaoVec *TXL ,TaoVec *TXU, TaoVec*S, double* bmin1,double* bmin2, double* bmax){
649:   int info;
650:   PetscReal p1,p2,p3;
651:   Vec X=(Vec)this->VecObject, DX=(Vec)S->VecObject;
652:   Vec XL=(Vec)TXL->VecObject, XU=(Vec)TXU->VecObject;
654:   info=VecBoundStepInfo(X,XL,XU,DX,&p1,&p2,&p3); CHKERRQ(info);
655:   *bmin1=p1; *bmin2=p2; *bmax=p3;
656:   return(0);
657: }

661: /*@C
662:    SetPetscViewer

664:    Input Parameter:
665: .  viewer - a viewer object to be used with this->View() and VecView()

667:    Level: advanced

669: @*/
670: int TaoVecPetsc::SetPetscViewer(PetscViewer viewer){
672:   this->pvecviewer= viewer;
673:   return(0);
674: }


679: int TaoVecPetsc::View(){
680:   int info;
682:   info=VecView(this->pv,this->pvecviewer);CHKERRQ(info);
683:   return(0);
684: }

688: int TaoVecPetsc::BoundGradientProjection(TaoVec*GG,TaoVec*XXL,TaoVec*XX, TaoVec*XXU){
689:   Vec v = (Vec) XXL->VecObject;
690:   Vec w = (Vec) XXU->VecObject;
691:   Vec x = (Vec) XX->VecObject;
692:   Vec g = (Vec) GG->VecObject;
693:   int info;
695:   info = PetscLogEventBegin(petscminmaxevent,0,0,0,0);CHKERRQ(info);
696:   info = VecBoundProjection(g,x,v,w,this->pv);CHKERRQ(info);
697:   info = PetscLogEventEnd(petscminmaxevent,0,0,0,0);CHKERRQ(info);
698:   return(0);
699: }

703: int TaoVecPetsc::CreateIndexSet(TaoIndexSet**SSS){
704:   TaoIndexSetPetsc *SS;

707:   SS = new TaoIndexSetPetsc(this->pv);
708:   *SSS=SS;
709:   return(0);
710: }


713: #endif