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