Actual source code: taois_petsc.c

      1: #include "tao_general.h"            /*I "tao.h" I*/
      2: #ifdef TAO_USE_PETSC
      3: #include "taois_petsc.h"
      4: #include "../vector/taovec_petsc.h"

      6: int VecWhichBetween(Vec, Vec, Vec, IS *);
      7: int VecWhichBetweenOrEqual(Vec, Vec, Vec, IS *);
      8: int VecWhichGreaterThan(Vec, Vec, IS * );
      9: int VecWhichLessThan(Vec, Vec, IS *);
     10: int VecWhichEqual(Vec, Vec, IS *);

     12: static int petscisevents=0;
     13: static int tredistribute=0;
     14: //static TaoPetscISType defaultistype=TaoRedistributeSubset;
     15: static TaoPetscISType defaultistype=TaoNoRedistributeSubset;
     16: //static TaoPetscISType defaultistype=TaoMaskFullSpace;

     20: /* @C
     21:    TaoWrapPetscIS - Create a new TaoIndexSet object using PETSc.

     23:    Input Parameter:
     24: +  S -  a PETSc IS
     25: -  imax - the maximum local length of the index set (Should be equal to or greater than the local length of the vector)

     27:    Output Parameter:
     28: .  SS - address of a new TaoIndexSet

     30:    Note:
     31:    A TaoIndexSetPetsc is an object with the methods of an abstract
     32:    TaoIndexSet object.  A TaoIndexSetPetsc contains an implementation of the
     33:    TaoIndexSet methods.  Routines using these vectors should declare a
     34:    pointer to a TaoIndexSet, assign this pointer to the address of a
     35:    TaoIndexSet object,use the pointer to invoke methods on the object, and use
     36:    this pointer as an argument when calling other routines.  This usage is
     37:    different from the usage of a PETSc IS.  In PETSc, applications will
     38:    typically declare a IS, and pass it as an argument into routines.  That is,
     39:    applications will typically declare a pointer to a TaoIndexSet and use the
     40:    pointer, or declare a IS and use it directly.

     42: .seealso TaoIndexSetGetPetscIS(), TaoIndexSetDestroy()

     44:    Level: developer
     45: @ */
     46: int TaoWrapPetscIS( IS S, int imax, TaoIndexSetPetsc ** SS){
     49:   //  if (0==0) return 1;
     50:   //  *SS = new  TaoIndexSetPetsc(imax, S);
     51:   PetscFunctionReturn(1);
     52: }

     56: /*@C
     57:    TaoIndexSetGetPetscIS - If the TaoIndexSet is of the TaoIndexSetPetsc type, it gets the underlying PETSc IS

     59:    Input Parameter:
     60: +  SS - the TaoIndexSetPetsc
     61: -  S -  the address of an IS

     63:    Output Parameter:
     64: .  S -  the PETSc IS

     66:    Note:
     67:    This routine does not create a IS.  It sets a pointer
     68:    to the location of an existing IS.

     70:    Level: advanced

     72: .seealso TaoVecGetPetscVec(), TaoMatGetPetscMat()
     73: @*/
     74: int TaoIndexSetGetPetscIS( TaoIndexSet *SS, IS * S){
     76:   if (SS){
     77:     *S=((TaoIndexSetPetsc *)SS)->GetIS();
     79:   } else {
     80:     *S=0;
     81:   }
     82:   return(0);
     83: }

     87: /* @C
     88:    TaoIndexSetPetsc - Create a new TaoIndexSet object using PETSc.

     90:    Input Parameter:
     91: -  V = A vector from the space that this indexset will describe.
     92: +  SS - an Index Set

     94:    Level: beginner

     96:    Note:
     97:    This index set will be deleted when the TaoIndexSet is deleted

     99:    Note:
    100:    The constuctor can also be called without the IndexSet
    101:
    102:    Note:
    103:    The method TaoIndexSetPetsc::GetIS() returns a pointer to the current PETSc Index Set.
    104: @ */
    105: TaoIndexSetPetsc::TaoIndexSetPetsc(Vec V, IS SS):TaoIndexSet(){
    106:   Vec V2;
    107:   VecDuplicate(V,&V2);
    108:   this->SetUp(V2,SS);
    109:   return;
    110: }

    114: /* @C
    115:    TaoIndexSetPetsc - Create a new TaoIndexSet object using PETSc.

    117:    Input Parameter:
    118: .  V = A vector from the space that this indexset will describe.

    120:    Level: beginner

    122:    Note:
    123:    The method TaoIndexSetPetsc::GetIS() returns a pointer to the current PETSc Index Set.
    124: @ */
    125: TaoIndexSetPetsc::TaoIndexSetPetsc(Vec V):TaoIndexSet(){
    126:   Vec V2;
    127:   VecDuplicate(V,&V2);
    128:   this->SetUp(V2,0);
    129:   return;
    130: }

    134: TaoIndexSetPetsc::~TaoIndexSetPetsc(){
    135:     clearit();
    136:     delete[] iptr;
    137:     if (this->VSpace){
    138:       VecDestroy(this->VSpace);
    139:     }
    140: };

    144: int TaoIndexSetPetsc::SetUp(Vec V2, IS SS){

    146:   int info,size,imax;
    147:   PetscTruth flg,ivalue;
    148:   MPI_Comm comm;

    153:   info=VecGetLocalSize(V2,&imax); CHKERRQ(info);
    154:   VSpace=V2;

    156:   iptr = new int[imax];

    158:   isp=0;isp2=0;ISGathered=0;ispComplement=0;scatter=0;
    159:   this->ispviewer=PETSC_VIEWER_STDOUT_WORLD;
    160:   reducedtype=defaultistype;

    162:   info=this->SetIS(SS); CHKERRQ(info);
    163:   nlocal=imax;

    165:   info = PetscObjectGetComm((PetscObject)V2,&comm);CHKERRQ(info);
    166:   MPI_Comm_size(comm,&size);
    167:   if (size==1 && defaultistype!=TaoMaskFullSpace ){
    168:     reducedtype=TaoSingleProcessor;
    169:   } else {
    170:     PetscOptionsGetTruth(0,"-redistribute",&ivalue,&flg);
    171:     if (flg && ivalue==PETSC_FALSE) reducedtype=TaoNoRedistributeSubset;
    172:     else if (flg && ivalue==PETSC_TRUE)reducedtype=TaoRedistributeSubset;
    173:     else reducedtype=defaultistype;
    174:   }
    175:   PetscOptionsHasName(0,"-submatrixfree",&flg);
    176:   if (flg) { reducedtype= TaoMaskFullSpace; }

    178:   if (petscisevents==0){
    179:     PetscLogEventRegister(&petscisevents,"Identify Indices",IS_COOKIE);
    180:   }
    181:   if (tredistribute==0){
    182:     PetscLogEventRegister(&tredistribute,"IS Redistribute ",IS_COOKIE);
    183:   }

    185:   return(0);
    186: }

    190: int TaoIndexSetPetsc::SetIS(IS SS){

    192:   int info;
    194:   if (SS){
    196:     PetscObjectReference((PetscObject)SS);
    197:   }
    198:   info=this->clearit(); CHKERRQ(info);
    199:   this->isp=SS;
    200:   this->ISObject=(void*)SS;
    201:   return(0);
    202: }

    206: /*@C
    207:    TaoSelectSubset - Select the manner in which subsets of variables in a matrix are selected

    209:    Input Parameter:
    210: .  type - the manner is which subsets are chosen.

    212:    Choice of types:
    213: +  TaoRedistributeSubset - When using multiple processors, redistribute the subset of variables
    214: over the processors.
    215: .  TaoNoRedistributeSubset - When using multiple processors, keep variables on current processors
    216: -  TaoMaskSubset - Use the full set of variables, and mask the unneeeded ones.

    218:    Note:
    219:    For use with active set methods in bound constrained optimization such as TRON and GPCG

    221:    Level: advanced
    222: @*/
    223: int TaoSelectSubset( TaoPetscISType type){
    225:   defaultistype=type;
    226:   return(0);
    227: }

    231: int TaoIndexSetPetsc::GetMask(Vec *VMask){
    232:   int info;
    233:   PetscScalar zero=0.0;

    236:   info=VecSet(this->VSpace, zero); CHKERRQ(info);
    237:   info=VecISSetToConstant(this->isp,1.0,this->VSpace); CHKERRQ(info);
    238:   *VMask=this->VSpace;
    239:   return(0);
    240: }

    244: int TaoIndexSetPetsc::UnionOf(TaoIndexSet* S1, TaoIndexSet* S2){
    245:   IS SS1=((TaoIndexSetPetsc*)S1)->GetIS();
    246:   IS SS2=((TaoIndexSetPetsc*)S2)->GetIS();
    247:   int info;
    249:   info=PetscLogEventBegin(petscisevents,0,0,0,0); CHKERRQ(info);
    250:   info=ISSum(&SS1,SS2); CHKERRQ(info);
    251:   info=PetscLogEventEnd(petscisevents,0,0,0,0); CHKERRQ(info);
    252:   return(0);
    253: }

    257: int TaoIndexSetPetsc::IntersectionOf(TaoIndexSet* S1, TaoIndexSet* S2){
    258:   IS SS1=((TaoIndexSetPetsc*)S1)->GetIS();
    259:   IS SS2=((TaoIndexSetPetsc*)S2)->GetIS();
    260:   IS SS3;
    261:   int info;
    263:   info=PetscLogEventBegin(petscisevents,0,0,0,0); CHKERRQ(info);
    264:   info=ISDifference(SS1,SS2,&SS3); CHKERRQ(info);
    265:   info=this->SetIS(SS3); CHKERRQ(info);
    266:   info=ISDestroy(SS3); CHKERRQ(info);
    267:   info=PetscLogEventEnd(petscisevents,0,0,0,0); CHKERRQ(info);
    268:   return(0);
    269: }

    273: int TaoIndexSetPetsc::ComplementOf(TaoIndexSet* SS){
    274:   IS SS1=((TaoIndexSetPetsc*)SS)->GetIS();
    275:   IS SS2;
    276:   int info;
    278:   info=PetscLogEventBegin(petscisevents,0,0,0,0); CHKERRQ(info);
    279:   info=ISCreateComplement(SS1,this->VSpace,&SS2); CHKERRQ(info);
    280:   info=this->SetIS(SS2); CHKERRQ(info);
    281:   info=ISDestroy(SS2); CHKERRQ(info);
    282:   info=PetscLogEventEnd(petscisevents,0,0,0,0); CHKERRQ(info);
    283:   return(0);
    284: }

    288: int TaoIndexSetPetsc::Duplicate(TaoIndexSet**SS){
    289:   int info;
    290:   TaoIndexSetPetsc *S;
    291:   IS is;

    294:   if (isp){
    295:     info = ISDuplicate(isp,&is); CHKERRQ(info);
    296:     S = new TaoIndexSetPetsc(this->VSpace, is);
    297:     info=ISDestroy(is); CHKERRQ(info);
    298:   } else {
    299:     S = new TaoIndexSetPetsc(this->VSpace);
    300:   }
    301:   *SS=S;

    303:   return(0);
    304: }

    308: int TaoIndexSetPetsc::IsSame(TaoIndexSet*SS, TaoTruth*flg){
    309:   int info;
    310:   PetscTruth pflg;
    311:   IS SS2=((TaoIndexSetPetsc*)SS)->GetIS();
    313:   info=PetscLogEventBegin(petscisevents,0,0,0,0); CHKERRQ(info);
    314:   info=ISEqual(isp,SS2,&pflg); CHKERRQ(info);
    315:   if (pflg==PETSC_TRUE) *flg=TAO_TRUE; else *flg=TAO_FALSE;
    316:   info=PetscLogEventEnd(petscisevents,0,0,0,0); CHKERRQ(info);
    317:   return(0);
    318: }

    322: int TaoIndexSetPetsc::View(){
    323:   int info;
    325:   info=ISView(this->isp,this->ispviewer);CHKERRQ(info); CHKERRQ(info);
    326:   return(0);
    327: }

    332: int TaoIndexSetPetsc::WhichEqual(TaoVec*v1,TaoVec*v2){
    333:   int info;
    335:   Vec V1=((TaoVecPetsc *)v1)->GetVec();
    336:   Vec V2=((TaoVecPetsc *)v2)->GetVec();
    337:   IS ISnew;
    339:   info=PetscLogEventBegin(petscisevents,0,0,0,0); CHKERRQ(info);
    340:   info=VecWhichEqual(V1,V2,&ISnew); CHKERRQ(info);
    341:   info=this->SetIS(ISnew); CHKERRQ(info);
    342:   info=ISDestroy(ISnew); CHKERRQ(info);
    343:   info=PetscLogEventEnd(petscisevents,0,0,0,0); CHKERRQ(info);
    344:   return(0);
    345: }

    349: int TaoIndexSetPetsc::WhichLessThan(TaoVec*v1,TaoVec*v2){
    350:   int info;
    351:   Vec V1=((TaoVecPetsc *)v1)->GetVec();
    352:   Vec V2=((TaoVecPetsc *)v2)->GetVec();
    353:   IS ISnew;
    355:   info=PetscLogEventBegin(petscisevents,0,0,0,0); CHKERRQ(info);
    356:   info=VecWhichLessThan(V1,V2,&ISnew); CHKERRQ(info);
    357:   info=this->SetIS(ISnew); CHKERRQ(info);
    358:   info=ISDestroy(ISnew); CHKERRQ(info);
    359:   info=PetscLogEventEnd(petscisevents,0,0,0,0); CHKERRQ(info);
    360:   return(0);
    361: }

    365: int TaoIndexSetPetsc::WhichGreaterThan(TaoVec*v1,TaoVec*v2){
    366:   int info;
    367:   Vec V1=((TaoVecPetsc *)v1)->GetVec();
    368:   Vec V2=((TaoVecPetsc *)v2)->GetVec();
    369:   IS ISnew;
    371:   info=PetscLogEventBegin(petscisevents,0,0,0,0); CHKERRQ(info);
    372:   info=VecWhichGreaterThan(V1,V2,&ISnew); CHKERRQ(info);
    373:   info=this->SetIS(ISnew); CHKERRQ(info);
    374:   info=ISDestroy(ISnew); CHKERRQ(info);
    375:   info=PetscLogEventEnd(petscisevents,0,0,0,0); CHKERRQ(info);
    376:   return(0);
    377: }

    381: int TaoIndexSetPetsc::WhichBetween(TaoVec*V1,TaoVec*V2,TaoVec*V3){
    382:   int info;
    383:   Vec VecLow=((TaoVecPetsc *)V1)->GetVec();
    384:   Vec VecHigh=((TaoVecPetsc *)V3)->GetVec();
    385:   Vec V=((TaoVecPetsc *)V2)->GetVec();
    386:   IS ISnew;

    389:   info=PetscLogEventBegin(petscisevents,0,0,0,0); CHKERRQ(info);
    390:   info=VecWhichBetween(VecLow,V,VecHigh,&ISnew); CHKERRQ(info);
    391:   info=this->SetIS(ISnew); CHKERRQ(info);
    392:   info=ISDestroy(ISnew); CHKERRQ(info);
    393:   info=PetscLogEventEnd(petscisevents,0,0,0,0);CHKERRQ(info);

    395:   return(0);
    396: }

    400: int TaoIndexSetPetsc::WhichBetweenOrEqual(TaoVec *V1, TaoVec *V2, TaoVec *V3)
    401: {
    402:   int info;
    403:   Vec VecLow  = ((TaoVecPetsc *)V1)->GetVec();
    404:   Vec VecHigh = ((TaoVecPetsc *)V3)->GetVec();
    405:   Vec V       = ((TaoVecPetsc *)V2)->GetVec();
    406:   IS ISnew;

    409:   info=PetscLogEventBegin(petscisevents,0,0,0,0); CHKERRQ(info);
    410:   info=VecWhichBetweenOrEqual(VecLow,V,VecHigh,&ISnew); CHKERRQ(info);
    411:   info=this->SetIS(ISnew); CHKERRQ(info);
    412:   info=ISDestroy(ISnew); CHKERRQ(info);
    413:   info=PetscLogEventEnd(petscisevents,0,0,0,0);CHKERRQ(info);
    414:   return(0);
    415: }

    419: int TaoIndexSetPetsc::GetSize(int *nn){
    421:   int info=ISGetSize(isp,nn); CHKERRQ(info);
    422:   return(0);
    423: }

    427: int TaoIndexSetPetsc::clearit(){
    428:   int info;
    430:   if (scatter){
    431:     info=VecScatterDestroy(scatter); CHKERRQ(info);
    432:     scatter=0;
    433:   }
    434:   if (ISGathered){
    435:     info=ISDestroy(ISGathered); CHKERRQ(info);
    436:     ISGathered=0;
    437:   }
    438:   if (isp){
    439:     info=ISDestroy(isp); CHKERRQ(info);
    440:     isp=0;
    441:   }
    442:   if (isp2){
    443:     info=ISDestroy(isp2); CHKERRQ(info);
    444:     isp2=0;
    445:   }
    446:   if (ispComplement){
    447:     info=ISDestroy(ispComplement); CHKERRQ(info);
    448:     ispComplement=0;
    449:   }
    450:   isp=0;isp2=0;ISGathered=0;ispComplement=0;scatter=0;

    452:   return(0);
    453: }

    458: int TaoIndexSetPetsc::RedistributeIS(IS*ISNew){
    459:   int info,i,n,nn,low,high;
    460:   int *gl,*is;
    461:   Vec R;
    462:   IS ISAll;
    463:   MPI_Comm comm;

    467:   info=PetscLogEventBegin(petscisevents,0,0,0,0);CHKERRQ(info);
    468:   info = VecGetSize(VSpace,&n); CHKERRQ(info);
    469:   info = this->GetSize(&nn); CHKERRQ(info);
    470:   if (reducedtype==TaoMaskFullSpace){
    471:     *ISNew = isp;
    472:   } else if (reducedtype==TaoSingleProcessor){
    473:     *ISNew = isp;
    474:   } else if (reducedtype==TaoNoRedistributeSubset || n==nn){
    475:     *ISNew = isp;
    476:   } else if (reducedtype==TaoRedistributeSubset){

    478:     info=PetscLogEventBegin(tredistribute,0,0,0,0);CHKERRQ(info);

    480:     if (isp2==0){
    481:       info = this->GetWholeIS(&ISAll); CHKERRQ(info);
    482:
    483:       info = PetscObjectGetComm((PetscObject)isp,&comm);CHKERRQ(info);
    484:       info = VecCreateMPI(comm,PETSC_DECIDE,nn,&R);CHKERRQ(info);
    485:       info = VecGetLocalSize(R,&n); CHKERRQ(info);
    486:       info = VecGetOwnershipRange(R, &low, &high); CHKERRQ(info);
    487:       info = ISGetIndices(ISAll, &is); CHKERRQ(info);
    488:       if (n>0){
    489:         info = PetscMalloc( n*sizeof(int),&gl ); CHKERRQ(info);
    490:         for (i=0; i<n; i++){
    491:           gl[i]= is[low+i];
    492:         }
    493:       } else gl=NULL;
    494:
    495:       info = PetscObjectGetComm((PetscObject)isp,&comm);CHKERRQ(info);
    496:       info = ISCreateGeneral(comm,n,gl,ISNew); CHKERRQ(info);
    497:       info = ISRestoreIndices(ISGathered, &is); CHKERRQ(info);
    498:       if (n>0) PetscFree(gl);

    500:       isp2=*ISNew;
    501:       info = VecDestroy(R); CHKERRQ(info);

    503:     }
    504:     info=PetscLogEventEnd(tredistribute,0,0,0,0);CHKERRQ(info);
    505:     *ISNew = isp2;
    506:
    507:   }

    509:   info=PetscLogEventEnd(petscisevents,0,0,0,0);CHKERRQ(info);
    510:   return(0);
    511: }

    516: int TaoIndexSetPetsc::GetWholeIS(IS*ISAll){
    517:   int info;
    519:   info=PetscLogEventBegin(petscisevents,0,0,0,0);CHKERRQ(info);
    520:   if (reducedtype==TaoSingleProcessor){
    521:     *ISAll = isp;
    522:   } else {
    523:     info=PetscLogEventBegin(tredistribute,0,0,0,0);CHKERRQ(info);
    524:     if (ISGathered==0){
    525:       info = ISAllGather(isp,&ISGathered); CHKERRQ(info);
    526:       info = ISSort(ISGathered); CHKERRQ(info);
    527:     } else if (ISGathered==0){
    528:       int      *indices,*sizes,size,*offsets,n,*lindices,i,N,ierr;
    529:       MPI_Comm comm;

    531:       ierr=PetscObjectGetComm((PetscObject)isp,&comm);
    532:       MPI_Comm_size(comm,&size);
    533:       PetscMalloc(2*size*sizeof(int),&sizes);
    534:       offsets = sizes + size;

    536:       ierr=ISGetLocalSize(isp,&n);
    537:       MPI_Allgather(&n,1,MPI_INT,sizes,1,MPI_INT,comm);
    538:       offsets[0] = 0;
    539:       for (i=1;i<size; i++) offsets[i] = offsets[i-1] + sizes[i-1];
    540:       N = offsets[size-1] + sizes[size-1];

    542:       PetscMalloc((N+1)*sizeof(int),&indices);
    543:       ierr=ISGetIndices(isp,&lindices);
    544:       MPI_Allgatherv(lindices,n,MPI_INT,indices,sizes,offsets,MPI_INT,comm);
    545:       ierr=ISRestoreIndices(isp,&lindices);

    547:       ierr=ISCreateGeneral(comm,N,indices,&ISGathered);
    548:       PetscFree(indices);
    549:       PetscFree(sizes);
    550:     }
    551:     *ISAll = ISGathered;
    552:     info=PetscLogEventEnd(tredistribute,0,0,0,0);CHKERRQ(info);
    553:   }

    555:   info=PetscLogEventEnd(petscisevents,0,0,0,0);CHKERRQ(info);
    556:   return(0);
    557: }

    562: int TaoIndexSetPetsc::GetReducedVecScatter(Vec VFull, Vec VSub, VecScatter*scatterit){
    563:   int info,n,nlow,nhigh;
    564:   IS S1, Ident;
    565:   MPI_Comm comm;

    569:   if (scatter==0){
    570:     info = VecGetSize(VFull,&n); CHKERRQ(info);
    571:     info = VecGetSize(VSub,&nlow); CHKERRQ(info);
    572:     if (reducedtype==TaoMaskFullSpace || n==nlow){

    574:       info = VecScatterCreate(VFull,isp,VSub,isp,&scatter); CHKERRQ(info);

    576:     } else if (reducedtype==TaoSingleProcessor || reducedtype==TaoNoRedistributeSubset){

    578:       S1=isp;
    579:       info = VecGetOwnershipRange(VSub,&nlow,&nhigh); CHKERRQ(info);
    580:       n=nhigh-nlow;

    582:       info = PetscObjectGetComm((PetscObject)VFull,&comm);CHKERRQ(info);
    583:       info = ISCreateStride(comm,n,nlow,1,&Ident); CHKERRQ(info);
    584:       info = VecScatterCreate(VFull,S1,VSub,Ident,&scatter); CHKERRQ(info);
    585:       info = ISDestroy(Ident); CHKERRQ(info);

    587:     } else {

    589:       info = this->GetWholeIS(&S1); CHKERRQ(info);
    590:       nlow = 0;
    591:       info = VecGetSize(VSub,&n);CHKERRQ(info);

    593:       info = PetscObjectGetComm((PetscObject)VFull,&comm);CHKERRQ(info);
    594:       info = ISCreateStride(comm,n,nlow,1,&Ident); CHKERRQ(info);
    595:       info = VecScatterCreate(VFull,S1,VSub,Ident,&scatter); CHKERRQ(info);
    596:       info = ISDestroy(Ident); CHKERRQ(info);

    598:     }
    599:   }
    600:
    601:   *scatterit=scatter;

    603:   return(0);
    604: }

    608: int TaoIndexSetPetsc::GetReducedType(TaoPetscISType* rtype){
    610:   *rtype=this->reducedtype;
    611:   return(0);
    612: }

    616: int TaoIndexSetPetsc::GetComplementIS(IS*isppp){
    617:   int info;
    619:   if (this->ispComplement==0){
    620:     info=PetscLogEventBegin(petscisevents,0,0,0,0);CHKERRQ(info);
    621:     info = ISCreateComplement(this->isp, this->VSpace, &this->ispComplement);CHKERRQ(info);
    622:     info=PetscLogEventEnd(petscisevents,0,0,0,0);CHKERRQ(info);
    623:   }
    624:   *isppp = this->ispComplement;
    625:   return(0);
    626: }

    630: int VecWhichEqual(Vec Vec1, Vec Vec2, IS * S)
    631: {
    632:   /*
    633:      Create an index set containing the indices of
    634:      the vectors Vec1 and Vec2 with identical elements.
    635:   */
    636:   int    i,n,low,high,low2,high2,n_same=0,info;
    637:   int    *same;
    638:   PetscScalar *v1,*v2;
    639:   MPI_Comm comm;

    647:   info = VecGetOwnershipRange(Vec1, &low, &high); CHKERRQ(info);
    648:   info = VecGetOwnershipRange(Vec2, &low2, &high2); CHKERRQ(info);

    650:   if ( low != low2 || high != high2 )
    651:     SETERRQ(1,"Vectors must be identically loaded over processors");

    653:   info = VecGetLocalSize(Vec1,&n); CHKERRQ(info);

    655:   if (n>0){
    656:
    657:     if (Vec1 == Vec2){
    658:       info = VecGetArray(Vec1,&v1); CHKERRQ(info);
    659:       v2=v1;
    660:     } else {
    661:       info = VecGetArray(Vec1,&v1); CHKERRQ(info);
    662:       info = VecGetArray(Vec2,&v2); CHKERRQ(info);
    663:     }

    665:     info = PetscMalloc( n*sizeof(int),&same ); CHKERRQ(info);
    666:
    667:     for (i=0; i<n; i++){
    668:       if (v1[i] == v2[i]) {same[n_same]=low+i; n_same++;}
    669:     }
    670:
    671:     if (Vec1 == Vec2){
    672:       info = VecRestoreArray(Vec1,&v1); CHKERRQ(info);
    673:     } else {
    674:       info = VecRestoreArray(Vec1,&v1); CHKERRQ(info);
    675:       info = VecRestoreArray(Vec2,&v2); CHKERRQ(info);
    676:     }

    678:   } else {

    680:     n_same = 0; same=NULL;

    682:   }

    684:   info = PetscObjectGetComm((PetscObject)Vec1,&comm);CHKERRQ(info);
    685:   info = ISCreateGeneral(comm,n_same,same,S);CHKERRQ(info);

    687:   if (same) PetscFree(same);

    689:   return(0);
    690: }

    694: int VecWhichLessThan(Vec Vec1, Vec Vec2, IS * S)
    695: {
    696:   int i,n,low,high,low2,high2,n_lt=0,info;
    697:   int *lt;
    698:   PetscScalar *v1,*v2;
    699:   MPI_Comm comm;

    706:   info = VecGetOwnershipRange(Vec1, &low, &high); CHKERRQ(info);
    707:   info = VecGetOwnershipRange(Vec2, &low2, &high2); CHKERRQ(info);

    709:   if ( low != low2 || high != high2 )
    710:     SETERRQ(1,"Vectors must be identically loaded over processors");

    712:   info = VecGetLocalSize(Vec1,&n); CHKERRQ(info);

    714:   if (n>0){

    716:     if (Vec1 == Vec2){
    717:       info = VecGetArray(Vec1,&v1); CHKERRQ(info);
    718:       v2=v1;
    719:     } else {
    720:       info = VecGetArray(Vec1,&v1); CHKERRQ(info);
    721:       info = VecGetArray(Vec2,&v2); CHKERRQ(info);
    722:     }
    723:     info = PetscMalloc( n*sizeof(int),&lt ); CHKERRQ(info);
    724:
    725:     for (i=0; i<n; i++){
    726:       if (v1[i] < v2[i]) {lt[n_lt]=high+i; n_lt++;}
    727:     }

    729:     if (Vec1 == Vec2){
    730:       info = VecRestoreArray(Vec1,&v1); CHKERRQ(info);
    731:     } else {
    732:       info = VecRestoreArray(Vec1,&v1); CHKERRQ(info);
    733:       info = VecRestoreArray(Vec2,&v2); CHKERRQ(info);
    734:     }
    735:
    736:   } else {
    737:     n_lt=0; lt=NULL;
    738:   }

    740:   info = PetscObjectGetComm((PetscObject)Vec1,&comm);CHKERRQ(info);
    741:   info = ISCreateGeneral(comm,n_lt,lt,S);CHKERRQ(info);

    743:   if (lt) PetscFree(lt);

    745:   return(0);
    746: }

    750: int VecWhichGreaterThan(Vec Vec1, Vec Vec2, IS * S)
    751: {
    752:   int    i,n,low,high,low2,high2,n_gt=0,info;
    753:   int    *gt=NULL;
    754:   PetscScalar *v1,*v2;
    755:   MPI_Comm comm;

    762:   info = VecGetOwnershipRange(Vec1, &low, &high); CHKERRQ(info);
    763:   info = VecGetOwnershipRange(Vec2, &low2, &high2); CHKERRQ(info);

    765:   if ( low != low2 || high != high2 )
    766:     SETERRQ(1,"Vectors must be identically loaded over processors");

    768:   info = VecGetLocalSize(Vec1,&n); CHKERRQ(info);

    770:   if (n>0){

    772:     if (Vec1 == Vec2){
    773:       info = VecGetArray(Vec1,&v1); CHKERRQ(info);
    774:       v2=v1;
    775:     } else {
    776:       info = VecGetArray(Vec1,&v1); CHKERRQ(info);
    777:       info = VecGetArray(Vec2,&v2); CHKERRQ(info);
    778:     }

    780:     info = PetscMalloc( n*sizeof(int), &gt ); CHKERRQ(info);
    781:
    782:     for (i=0; i<n; i++){
    783:       if (v1[i] > v2[i]) {gt[n_gt]=low+i; n_gt++;}
    784:     }

    786:     if (Vec1 == Vec2){
    787:       info = VecRestoreArray(Vec1,&v1); CHKERRQ(info);
    788:     } else {
    789:       info = VecRestoreArray(Vec1,&v1); CHKERRQ(info);
    790:       info = VecRestoreArray(Vec2,&v2); CHKERRQ(info);
    791:     }
    792:
    793:   } else{
    794:
    795:     n_gt=0; gt=NULL;

    797:   }

    799:   info = PetscObjectGetComm((PetscObject)Vec1,&comm);CHKERRQ(info);
    800:   info = ISCreateGeneral(comm,n_gt,gt,S);CHKERRQ(info);

    802:   if (gt) PetscFree(gt);

    804:   return(0);
    805: }

    809: int VecWhichBetween(Vec VecLow, Vec V, Vec VecHigh, IS * S)
    810: {
    811:   /*
    812:      Creates an index set with the indices of V whose
    813:      elements are stictly between the corresponding elements
    814:      of the vector VecLow and the Vector VecHigh
    815:   */
    816:   int i,n,low,high,low2,high2,low3,high3,n_vm=0,info;
    817:   int *vm;
    818:   PetscScalar *v1,*v2,*vmiddle;
    819:   MPI_Comm comm;

    824:   info = VecGetOwnershipRange(VecLow, &low, &high); CHKERRQ(info);
    825:   info = VecGetOwnershipRange(VecHigh, &low2, &high2); CHKERRQ(info);
    826:   info = VecGetOwnershipRange(V, &low3, &high3); CHKERRQ(info);

    828:   if ( low!=low2 || high!=high2 || low!=low3 || high!=high3 )
    829:     SETERRQ(1,"Vectors must be identically loaded over processors");

    831:   info = VecGetLocalSize(VecLow,&n); CHKERRQ(info);

    833:   if (n>0){

    835:     info = VecGetArray(VecLow,&v1); CHKERRQ(info);
    836:     if (VecLow != VecHigh){
    837:       info = VecGetArray(VecHigh,&v2); CHKERRQ(info);
    838:     } else {
    839:       v2=v1;
    840:     }
    841:     if ( V != VecLow && V != VecHigh){
    842:       info = VecGetArray(V,&vmiddle); CHKERRQ(info);
    843:     } else if ( V==VecLow ){
    844:       vmiddle=v1;
    845:     } else {
    846:       vmiddle =v2;
    847:     }

    849:     info = PetscMalloc( n*sizeof(int), &vm ); CHKERRQ(info);
    850:
    851:     for (i=0; i<n; i++){
    852:       if (v1[i] < vmiddle[i] && vmiddle[i] < v2[i]) {vm[n_vm]=low+i; n_vm++;}
    853:     }

    855:     info = VecRestoreArray(VecLow,&v1); CHKERRQ(info);
    856:     if (VecLow != VecHigh){
    857:       info = VecRestoreArray(VecHigh,&v2); CHKERRQ(info);
    858:     }
    859:     if ( V != VecLow && V != VecHigh){
    860:       info = VecRestoreArray(V,&vmiddle); CHKERRQ(info);
    861:     }

    863:   } else {

    865:     n_vm=0; vm=NULL;

    867:   }

    869:   info = PetscObjectGetComm((PetscObject)V,&comm);CHKERRQ(info);
    870:   info = ISCreateGeneral(comm,n_vm,vm,S);CHKERRQ(info);

    872:   if (vm) PetscFree(vm);

    874:   return(0);
    875: }

    880: int VecWhichBetweenOrEqual(Vec VecLow, Vec V, Vec VecHigh, IS * S)
    881: {
    882:   /*
    883:      Creates an index set with the indices of V whose
    884:      elements are stictly between the corresponding elements
    885:      of the vector VecLow and the Vector VecHigh
    886:   */
    887:   int i,n,low,high,low2,high2,low3,high3,n_vm=0,info;
    888:   int *vm;
    889:   PetscScalar *v1,*v2,*vmiddle;
    890:   MPI_Comm comm;

    895:   info = VecGetOwnershipRange(VecLow, &low, &high); CHKERRQ(info);
    896:   info = VecGetOwnershipRange(VecHigh, &low2, &high2); CHKERRQ(info);
    897:   info = VecGetOwnershipRange(V, &low3, &high3); CHKERRQ(info);

    899:   if ( low!=low2 || high!=high2 || low!=low3 || high!=high3 )
    900:     SETERRQ(1,"Vectors must be identically loaded over processors");

    902:   info = VecGetLocalSize(VecLow,&n); CHKERRQ(info);

    904:   if (n>0){

    906:     info = VecGetArray(VecLow,&v1); CHKERRQ(info);
    907:     if (VecLow != VecHigh){
    908:       info = VecGetArray(VecHigh,&v2); CHKERRQ(info);
    909:     } else {
    910:       v2=v1;
    911:     }
    912:     if ( V != VecLow && V != VecHigh){
    913:       info = VecGetArray(V,&vmiddle); CHKERRQ(info);
    914:     } else if ( V==VecLow ){
    915:       vmiddle=v1;
    916:     } else {
    917:       vmiddle =v2;
    918:     }

    920:     info = PetscMalloc( n*sizeof(int), &vm ); CHKERRQ(info);
    921:
    922:     for (i=0; i<n; i++){
    923:       if (v1[i] <= vmiddle[i] && vmiddle[i] <= v2[i]) {vm[n_vm]=low+i; n_vm++;}
    924:     }

    926:     info = VecRestoreArray(VecLow,&v1); CHKERRQ(info);
    927:     if (VecLow != VecHigh){
    928:       info = VecRestoreArray(VecHigh,&v2); CHKERRQ(info);
    929:     }
    930:     if ( V != VecLow && V != VecHigh){
    931:       info = VecRestoreArray(V,&vmiddle); CHKERRQ(info);
    932:     }

    934:   } else {

    936:     n_vm=0; vm=NULL;

    938:   }

    940:   info = PetscObjectGetComm((PetscObject)V,&comm);CHKERRQ(info);
    941:   info = ISCreateGeneral(comm,n_vm,vm,S);CHKERRQ(info);

    943:   if (vm) PetscFree(vm);

    945:   return(0);
    946: }

    951: /*@C
    952:    ISCreateComplement - Creates the complement of the the index set

    954:    Input Parameter:
    955: +  S -  a PETSc IS
    956: -  V - the reference vector space

    958:    Output Parameter:
    959: .  T -  the complement of S

    962: .seealso ISCreateGeneral()

    964:    Level: advanced
    965: @*/
    966: int ISCreateComplement(IS S, Vec V, IS *T){
    967:   int info,i,nis,nloc,high,low,n=0;
    968:   int *s,*tt,*ss;
    969:   MPI_Comm comm;

    975:   info = VecGetOwnershipRange(V,&low,&high); CHKERRQ(info);
    976:   info = VecGetLocalSize(V,&nloc); CHKERRQ(info);
    977:   info = ISGetLocalSize(S,&nis); CHKERRQ(info);
    978:   info = ISGetIndices(S, &s); CHKERRQ(info);
    979:   info = PetscMalloc( nloc*sizeof(int),&tt ); CHKERRQ(info);
    980:   info = PetscMalloc( nloc*sizeof(int),&ss ); CHKERRQ(info);

    982:   for (i=low; i<high; i++){ tt[i-low]=i; }

    984:   for (i=0; i<nis; i++){ tt[s[i]-low] = -2; }
    985:
    986:   for (i=0; i<nloc; i++){
    987:     if (tt[i]>-1){ ss[n]=tt[i]; n++; }
    988:   }

    990:   info = ISRestoreIndices(S, &s); CHKERRQ(info);
    991:
    992:   info = PetscObjectGetComm((PetscObject)S,&comm);CHKERRQ(info);
    993:   info = ISCreateGeneral(comm,n,ss,T);CHKERRQ(info);
    994:
    995:   if (tt) PetscFree(tt);
    996:   if (ss) PetscFree(ss);

    998:   return(0);
    999: }
    1000: #endif