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