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 PetscLogEvent petscisevents=0;
 13: static PetscLogEvent 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, PetscInt imax, TaoIndexSetPetsc ** SS){
 49:   //  if (0==0) return 1;
 50:   //  *SS = new  TaoIndexSetPetsc(imax, S);
 51:   PetscFunctionReturn(1);
 52: }

 56: /* @C
 57:    TaoIndexSetPetsc - Create a new TaoIndexSet object using PETSc.

 59:    Input Parameter:
 60: -  V = A vector from the space that this indexset will describe.
 61: +  SS - an Index Set

 63:    Level: beginner

 65:    Note: 
 66:    This index set will be deleted when the TaoIndexSet is deleted

 68:    Note: 
 69:    The constuctor can also be called without the IndexSet
 70:    
 71:    Note:
 72:    The method TaoIndexSetPetsc::GetIS() returns a pointer to the current PETSc Index Set.
 73: @ */
 74: TaoIndexSetPetsc::TaoIndexSetPetsc(Vec V, IS SS):TaoIndexSet(){
 75:   Vec V2;
 76:   isp=0;isp2=0;ISGathered=0;ispComplement=0;scatter=0;
 77:   VecDuplicate(V,&V2);
 78:   this->SetUp(V2,SS);
 79:   return;
 80: }

 84: /* @C
 85:    TaoIndexSetPetsc - Create a new TaoIndexSet object using PETSc.

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

 90:    Level: beginner

 92:    Note:
 93:    The method TaoIndexSetPetsc::GetIS() returns a pointer to the current PETSc Index Set.
 94: @ */
 95: TaoIndexSetPetsc::TaoIndexSetPetsc(Vec V):TaoIndexSet(){
 96:   int info;
 97:   PetscInt low,high;
 98:   IS is;
 99:   Vec V2;
100:   isp=0;isp2=0;ISGathered=0;ispComplement=0;scatter=0;
101:   VecDuplicate(V,&V2);
102:   info=VecGetOwnershipRange(V,&low,&high);
103:   info=ISCreateStride(PETSC_COMM_WORLD,1,low,1,&is);
104:   this->SetUp(V2,is);
105:   info = PetscObjectDereference((PetscObject)is);
106:   return;
107: }

111: TaoIndexSetPetsc::~TaoIndexSetPetsc(){
112:   clearit();
113:   delete[] iptr; 
114:   if (this->VSpace){
115:     VecDestroy(this->VSpace);
116:   }
117: };

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

123:   int info;
124:   int size;
125:   PetscInt imax;
126:   PetscTruth flg,ivalue;
127:   MPI_Comm comm;


132:   info=VecGetLocalSize(V2,&imax); CHKERRQ(info);
133:   VSpace=V2;

135:   iptr = new PetscInt[imax];

137:   isp=0;isp2=0;ISGathered=0;ispComplement=0;scatter=0;
138:   this->ispviewer=PETSC_VIEWER_STDOUT_WORLD;
139:   reducedtype=defaultistype;

141:   info=this->SetIS(SS); CHKERRQ(info);
142:   nlocal=imax;

144:   info = PetscObjectGetComm((PetscObject)V2,&comm);CHKERRQ(info);
145:   MPI_Comm_size(comm,&size);
146:   if (size==1 && defaultistype!=TaoMaskFullSpace &&defaultistype!=TaoMatrixFree ){ 
147:     reducedtype=TaoSingleProcessor;
148:   } else {
149:     PetscOptionsGetTruth(0,"-redistribute",&ivalue,&flg);
150:     if (flg && ivalue==PETSC_FALSE) reducedtype=TaoNoRedistributeSubset;
151:     else if (flg && ivalue==PETSC_TRUE)reducedtype=TaoRedistributeSubset;
152:     else reducedtype=defaultistype;
153:   }
154:   PetscOptionsHasName(0,"-taomask",&flg);
155:   if (flg) { reducedtype= TaoMaskFullSpace; }
156:   PetscOptionsHasName(0,"-taosubmatrixfree",&flg);
157:   if (flg) { reducedtype= TaoMatrixFree; }

159:   if (petscisevents==0){
160:     PetscLogEventRegister("Identify Indices",IS_COOKIE,&petscisevents);
161:   }
162:   if (tredistribute==0){
163:     PetscLogEventRegister("IS Redistribute ",IS_COOKIE,&tredistribute);
164:   }

166:   return(0);
167: }

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

173:   int info;
175:   if (SS){
177:     PetscObjectReference((PetscObject)SS); 
178:   }
179:   info=this->clearit(); CHKERRQ(info);
180:   this->isp=SS;
181:   return(0);
182: }

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

189:    Input Parameter:
190: .  type - the manner is which subsets are chosen.  

192:    Choice of types:
193: +  TaoRedistributeSubset - When using multiple processors, redistribute the subset of variables
194: over the processors.
195: .  TaoNoRedistributeSubset - When using multiple processors, keep variables on current processors
196: .  TaoMaskFullSpace - Use the full set of variables, and mask the unneeded ones.
197: -  TaoMatrixFree - Use the full set of variables, and mask the unneeded ones.

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

202:    Level: advanced
203: @*/
204: int TaoSelectSubset( TaoPetscISType type){
206:   defaultistype=type;
207:   return(0);
208: }

212: int TaoIndexSetPetsc::GetMask(Vec *VMask){
213:   int info;
214:   PetscScalar zero=0.0;

217:   info=VecSet(this->VSpace, zero); CHKERRQ(info);
218:   info=VecISSetToConstant(this->isp,1.0,this->VSpace); CHKERRQ(info);
219:   *VMask=this->VSpace;
220:   return(0);
221: }

225: int TaoIndexSetPetsc::UnionOf(TaoIndexSet* S1, TaoIndexSet* S2){
226:   IS SS1=((TaoIndexSetPetsc*)S1)->GetIS();
227:   IS SS2=((TaoIndexSetPetsc*)S2)->GetIS();
228:   int info;
230:   info=PetscLogEventBegin(petscisevents,0,0,0,0); CHKERRQ(info);
231:   info=ISSum(SS1,SS2,&SS1); CHKERRQ(info);
232:   info=PetscLogEventEnd(petscisevents,0,0,0,0); CHKERRQ(info);
233:   return(0);
234: }

238: int TaoIndexSetPetsc::IntersectionOf(TaoIndexSet* S1, TaoIndexSet* S2){
239:   IS SS1=((TaoIndexSetPetsc*)S1)->GetIS();
240:   IS SS2=((TaoIndexSetPetsc*)S2)->GetIS();
241:   IS SS3;
242:   int info;
244:   info=PetscLogEventBegin(petscisevents,0,0,0,0); CHKERRQ(info);
245:   info=ISDifference(SS1,SS2,&SS3); CHKERRQ(info);
246:   info=this->SetIS(SS3); CHKERRQ(info);
247:   info=ISDestroy(SS3); CHKERRQ(info);
248:   info=PetscLogEventEnd(petscisevents,0,0,0,0); CHKERRQ(info);
249:   return(0);
250: }

254: int TaoIndexSetPetsc::ComplementOf(TaoIndexSet* SS){
255:   IS SS1=((TaoIndexSetPetsc*)SS)->GetIS();
256:   IS SS2;
257:   int info;
259:   info=PetscLogEventBegin(petscisevents,0,0,0,0); CHKERRQ(info);
260:   info=ISCreateComplement(SS1,this->VSpace,&SS2); CHKERRQ(info);
261:   info=this->SetIS(SS2); CHKERRQ(info);
262:   info=ISDestroy(SS2); CHKERRQ(info);
263:   info=PetscLogEventEnd(petscisevents,0,0,0,0); CHKERRQ(info);
264:   return(0);
265: }

269: int TaoIndexSetPetsc::Duplicate(TaoIndexSet**SS){
270:   int info;
271:   TaoIndexSetPetsc *S;
272:   IS is;

275:   if (isp){
276:     info = ISDuplicate(isp,&is); CHKERRQ(info);
277:     S = new TaoIndexSetPetsc(this->VSpace, is);
278:     info=ISDestroy(is); CHKERRQ(info);
279:   } else {
280:     S = new TaoIndexSetPetsc(this->VSpace);
281:   }
282:   *SS=S;

284:   return(0);
285: }

289: int TaoIndexSetPetsc::IsSame(TaoIndexSet*SS, TaoTruth*flg){
290:   int info;
291:   PetscTruth pflg;
292:   IS SS2=((TaoIndexSetPetsc*)SS)->GetIS();
294:   info=PetscLogEventBegin(petscisevents,0,0,0,0); CHKERRQ(info);
295:   info=ISEqual(isp,SS2,&pflg); CHKERRQ(info);
296:   if (pflg==PETSC_TRUE) *flg=TAO_TRUE; else *flg=TAO_FALSE;
297:   info=PetscLogEventEnd(petscisevents,0,0,0,0); CHKERRQ(info);
298:   return(0);
299: }

303: int TaoIndexSetPetsc::View(){
304:   int info;
306:   info=ISView(this->isp,this->ispviewer);CHKERRQ(info); CHKERRQ(info);
307:   return(0);
308: }


313: int TaoIndexSetPetsc::WhichEqual(TaoVec*v1,TaoVec*v2){
314:   int info;
315:   TaoVecPetsc *tv1, *tv2;
316:   IS ISnew;
318:   info=PetscLogEventBegin(petscisevents,0,0,0,0); CHKERRQ(info);
319:   tv1 = dynamic_cast<TaoVecPetsc*>(v1);
320:   tv2 = dynamic_cast<TaoVecPetsc*>(v2);
321:   if (tv1 == NULL) {
322:     SETERRQ(1,"Could not cast argument 1 to TaoVecPetsc*");
323:   }
324:   if (tv2 == NULL) {
325:     SETERRQ(1,"Could not cast argument 2 to TaoVecPetsc*");
326:   }
327:   info=VecWhichEqual(tv1->GetVec(),tv2->GetVec(),&ISnew); CHKERRQ(info);
328:   info=this->SetIS(ISnew); CHKERRQ(info);
329:   info=ISDestroy(ISnew); CHKERRQ(info);
330:   info=PetscLogEventEnd(petscisevents,0,0,0,0); CHKERRQ(info);
331:   return(0);
332: }

336: int TaoIndexSetPetsc::WhichLessThan(TaoVec*v1,TaoVec*v2){
337:   int info;
338:   Vec V1=((TaoVecPetsc *)v1)->GetVec();
339:   Vec V2=((TaoVecPetsc *)v2)->GetVec();
340:   IS ISnew;
342:   info=PetscLogEventBegin(petscisevents,0,0,0,0); CHKERRQ(info);
343:   info=VecWhichLessThan(V1,V2,&ISnew); CHKERRQ(info);
344:   info=this->SetIS(ISnew); CHKERRQ(info);
345:   info=ISDestroy(ISnew); CHKERRQ(info);
346:   info=PetscLogEventEnd(petscisevents,0,0,0,0); CHKERRQ(info);
347:   return(0);
348: }

352: int TaoIndexSetPetsc::WhichGreaterThan(TaoVec*v1,TaoVec*v2){
353:   int info;
354:   Vec V1=((TaoVecPetsc *)v1)->GetVec();
355:   Vec V2=((TaoVecPetsc *)v2)->GetVec();
356:   IS ISnew;
358:   info=PetscLogEventBegin(petscisevents,0,0,0,0); CHKERRQ(info);
359:   info=VecWhichGreaterThan(V1,V2,&ISnew); CHKERRQ(info);
360:   info=this->SetIS(ISnew); CHKERRQ(info);
361:   info=ISDestroy(ISnew); CHKERRQ(info);
362:   info=PetscLogEventEnd(petscisevents,0,0,0,0); CHKERRQ(info);
363:   return(0);
364: }

368: int TaoIndexSetPetsc::WhichBetween(TaoVec*V1,TaoVec*V2,TaoVec*V3){
369:   int info;
370:   Vec VecLow=((TaoVecPetsc *)V1)->GetVec();
371:   Vec VecHigh=((TaoVecPetsc *)V3)->GetVec();
372:   Vec V=((TaoVecPetsc *)V2)->GetVec();
373:   IS ISnew;

376:   info=PetscLogEventBegin(petscisevents,0,0,0,0); CHKERRQ(info);
377:   info=VecWhichBetween(VecLow,V,VecHigh,&ISnew); CHKERRQ(info);
378:   info=this->SetIS(ISnew); CHKERRQ(info);
379:   info=ISDestroy(ISnew); CHKERRQ(info);
380:   info=PetscLogEventEnd(petscisevents,0,0,0,0);CHKERRQ(info);          

382:   return(0);
383: }

387: int TaoIndexSetPetsc::WhichBetweenOrEqual(TaoVec *V1, TaoVec *V2, TaoVec *V3)
388: {
389:   int info;
390:   Vec VecLow  = ((TaoVecPetsc *)V1)->GetVec();
391:   Vec VecHigh = ((TaoVecPetsc *)V3)->GetVec();
392:   Vec V       = ((TaoVecPetsc *)V2)->GetVec();
393:   IS ISnew;

396:   info=PetscLogEventBegin(petscisevents,0,0,0,0); CHKERRQ(info);
397:   info=VecWhichBetweenOrEqual(VecLow,V,VecHigh,&ISnew); CHKERRQ(info);
398:   info=this->SetIS(ISnew); CHKERRQ(info);
399:   info=ISDestroy(ISnew); CHKERRQ(info);
400:   info=PetscLogEventEnd(petscisevents,0,0,0,0);CHKERRQ(info);
401:   return(0);
402: }

406: int TaoIndexSetPetsc::GetSize(TaoInt *nn){
408:   int info=ISGetSize(isp,nn); CHKERRQ(info);
409:   return(0);
410: }

414: int TaoIndexSetPetsc::clearit(){
415:   int info;
417:   if (scatter){
418:     info=VecScatterDestroy(scatter); CHKERRQ(info);
419:     scatter=0;
420:   }
421:   if (ISGathered){
422:     info=ISDestroy(ISGathered); CHKERRQ(info);
423:     ISGathered=0;
424:   }
425:   if (isp){
426:     info=ISDestroy(isp); CHKERRQ(info);
427:     isp=0;
428:   }
429:   if (isp2){
430:     info=ISDestroy(isp2); CHKERRQ(info);
431:     isp2=0;
432:   }
433:   if (ispComplement){
434:     info=ISDestroy(ispComplement); CHKERRQ(info);
435:     ispComplement=0;
436:   }
437:   isp=0;isp2=0;ISGathered=0;ispComplement=0;scatter=0;

439:   return(0);
440: }


445: int TaoIndexSetPetsc::RedistributeIS(IS*ISNew){
446:   int info;
447:   PetscInt n,nn,low,high,i;
448:   const PetscInt *is;
449:   PetscInt *gl;
450:   Vec R;
451:   IS ISAll;
452:   MPI_Comm comm;


456:   info=PetscLogEventBegin(petscisevents,0,0,0,0);CHKERRQ(info);
457:   info = VecGetSize(VSpace,&n); CHKERRQ(info);
458:   info = this->GetSize(&nn); CHKERRQ(info);
459:   if (reducedtype==TaoMaskFullSpace){
460:     *ISNew = isp;
461:   } else if (reducedtype==TaoSingleProcessor){
462:     *ISNew = isp;
463:   } else if (reducedtype==TaoNoRedistributeSubset || n==nn){
464:     *ISNew = isp;
465:   } else if (reducedtype==TaoRedistributeSubset){

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

469:     if (isp2==0){
470:       info = this->GetWholeIS(&ISAll); CHKERRQ(info);
471:       
472:       info = PetscObjectGetComm((PetscObject)isp,&comm);CHKERRQ(info);
473:       info = VecCreateMPI(comm,PETSC_DECIDE,nn,&R);CHKERRQ(info);
474:       info = VecGetLocalSize(R,&n); CHKERRQ(info);
475:       info = VecGetOwnershipRange(R, &low, &high); CHKERRQ(info);
476:       info = ISGetIndices(ISAll, &is); CHKERRQ(info);
477:       if (n>0){
478:         info = PetscMalloc( n*sizeof(PetscInt),&gl ); CHKERRQ(info);
479:         for (i=0; i<n; i++){
480:           gl[i]= is[low+i];
481:         } 
482:       } else gl=NULL;
483:       
484:       info = PetscObjectGetComm((PetscObject)isp,&comm);CHKERRQ(info);
485:       info = ISCreateGeneral(comm,n,gl,ISNew); CHKERRQ(info);
486:       info = ISRestoreIndices(ISGathered, &is); CHKERRQ(info);
487:       if (n>0) {
488:         info = PetscFree(gl); CHKERRQ(info);
489:       }

491:       isp2=*ISNew;
492:       info = VecDestroy(R); CHKERRQ(info);

494:     }
495:     info=PetscLogEventEnd(tredistribute,0,0,0,0);CHKERRQ(info);
496:     *ISNew = isp2;
497:     
498:   }

500:   info=PetscLogEventEnd(petscisevents,0,0,0,0);CHKERRQ(info);
501:   return(0);
502: }


507: int TaoIndexSetPetsc::GetWholeIS(IS*ISAll){
508:   int info;
510:   info=PetscLogEventBegin(petscisevents,0,0,0,0);CHKERRQ(info);
511:   if (reducedtype==TaoSingleProcessor){
512:     *ISAll = isp;
513:   } else {
514:     info=PetscLogEventBegin(tredistribute,0,0,0,0);CHKERRQ(info);
515:     if (ISGathered==0){
516:       info = ISAllGather(isp,&ISGathered); CHKERRQ(info);
517:       info = ISSort(ISGathered); CHKERRQ(info);
518:     } else if (ISGathered==0){
519:       PetscInt  N;
520:       PetscMPIInt *sizes,*offsets,size,mpin;
521:       const PetscInt *lindices;
522:       PetscInt *indices;
523:       PetscInt i,n;
524:       MPI_Comm comm;

526:       info=PetscObjectGetComm((PetscObject)isp,&comm);CHKERRQ(info);
527:       MPI_Comm_size(comm,&size);
528:       PetscMalloc(2*size*sizeof(PetscInt),&sizes);
529:       offsets = sizes + size;

531:       info=ISGetLocalSize(isp,&n);CHKERRQ(info);
532:       mpin = (PetscMPIInt) n;
533:       if ((PetscInt)mpin != n) {
534:         SETERRQ(1,"Did not handle 64-bit-index correctly");
535:       }
536:       MPI_Allgather(&mpin,1,MPI_INT,sizes,1,MPI_INT,comm);
537:       offsets[0] = 0;
538:       for (i=1;i<size; i++) offsets[i] = offsets[i-1] + sizes[i-1];
539:       N = offsets[size-1] + sizes[size-1];

541:       PetscMalloc((N+1)*sizeof(PetscInt),&indices);
542:       info=ISGetIndices(isp,&lindices);CHKERRQ(info);
543:       MPI_Allgatherv((void*)lindices,n,MPIU_INT,indices,sizes,offsets,MPIU_INT,comm);
544:       info=ISRestoreIndices(isp,&lindices);CHKERRQ(info);

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

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


561: int TaoIndexSetPetsc::GetReducedVecScatter(Vec VFull, Vec VSub, VecScatter*scatterit){
562:   int info;
563:   PetscInt 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    info;
637:   PetscInt i,n_same=0;
638:   PetscInt n,low,high,low2,high2;
639:   PetscInt    *same;
640:   PetscScalar *v1,*v2;
641:   MPI_Comm comm;



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

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

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

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

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

680:   } else {

682:     n_same = 0; same=NULL;

684:   }

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

689:   if (same) {
690:     info = PetscFree(same); CHKERRQ(info);
691:   }

693:   return(0);
694: }

698: int VecWhichLessThan(Vec Vec1, Vec Vec2, IS * S)
699: {
700:   int info;
701:   PetscInt i;
702:   PetscInt n,low,high,low2,high2,n_lt=0;
703:   PetscInt *lt;
704:   PetscScalar *v1,*v2;
705:   MPI_Comm comm;


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

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

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

720:   if (n>0){

722:     if (Vec1 == Vec2){
723:       info = VecGetArray(Vec1,&v1); CHKERRQ(info);
724:       v2=v1;
725:     } else {
726:       info = VecGetArray(Vec1,&v1); CHKERRQ(info);
727:       info = VecGetArray(Vec2,&v2); CHKERRQ(info);
728:     }
729:     info = PetscMalloc( n*sizeof(PetscInt),&lt ); CHKERRQ(info);
730:     
731:     for (i=0; i<n; i++){
732:       if (v1[i] < v2[i]) {lt[n_lt]=high+i; n_lt++;}
733:     }

735:     if (Vec1 == Vec2){
736:       info = VecRestoreArray(Vec1,&v1); CHKERRQ(info);
737:     } else {
738:       info = VecRestoreArray(Vec1,&v1); CHKERRQ(info);
739:       info = VecRestoreArray(Vec2,&v2); CHKERRQ(info);
740:     }
741:       
742:   } else {
743:     n_lt=0; lt=NULL;
744:   }

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

749:   if (lt) {
750:     info = PetscFree(lt); CHKERRQ(info);
751:   }

753:   return(0);
754: }

758: int VecWhichGreaterThan(Vec Vec1, Vec Vec2, IS * S)
759: {
760:   int    info;
761:   PetscInt n,low,high,low2,high2,n_gt=0,i;
762:   PetscInt    *gt=NULL;
763:   PetscScalar *v1,*v2;
764:   MPI_Comm comm;


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

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

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

779:   if (n>0){

781:     if (Vec1 == Vec2){
782:       info = VecGetArray(Vec1,&v1); CHKERRQ(info);
783:       v2=v1;
784:     } else {
785:       info = VecGetArray(Vec1,&v1); CHKERRQ(info);
786:       info = VecGetArray(Vec2,&v2); CHKERRQ(info);
787:     }    

789:     info = PetscMalloc( n*sizeof(PetscInt), &gt ); CHKERRQ(info);
790:     
791:     for (i=0; i<n; i++){
792:       if (v1[i] > v2[i]) {gt[n_gt]=low+i; n_gt++;}
793:     }

795:     if (Vec1 == Vec2){
796:       info = VecRestoreArray(Vec1,&v1); CHKERRQ(info);
797:     } else {
798:       info = VecRestoreArray(Vec1,&v1); CHKERRQ(info);
799:       info = VecRestoreArray(Vec2,&v2); CHKERRQ(info);
800:     }
801:     
802:   } else{
803:     
804:     n_gt=0; gt=NULL;

806:   }

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

811:   if (gt) {
812:     info = PetscFree(gt); CHKERRQ(info);
813:   }
814:   return(0);
815: }

819: int VecWhichBetween(Vec VecLow, Vec V, Vec VecHigh, IS * S)
820: {
821:   /* 
822:      Creates an index set with the indices of V whose 
823:      elements are stictly between the corresponding elements 
824:      of the vector VecLow and the Vector VecHigh
825:   */
826:   int info;
827:   PetscInt n,low,high,low2,high2,low3,high3,n_vm=0;
828:   PetscInt *vm,i;
829:   PetscScalar *v1,*v2,*vmiddle;
830:   MPI_Comm comm;


835:   info = VecGetOwnershipRange(VecLow, &low, &high); CHKERRQ(info);
836:   info = VecGetOwnershipRange(VecHigh, &low2, &high2); CHKERRQ(info);
837:   info = VecGetOwnershipRange(V, &low3, &high3); CHKERRQ(info);

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

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

844:   if (n>0){

846:     info = VecGetArray(VecLow,&v1); CHKERRQ(info);
847:     if (VecLow != VecHigh){
848:       info = VecGetArray(VecHigh,&v2); CHKERRQ(info);
849:     } else {
850:       v2=v1;
851:     }
852:     if ( V != VecLow && V != VecHigh){
853:       info = VecGetArray(V,&vmiddle); CHKERRQ(info);
854:     } else if ( V==VecLow ){
855:       vmiddle=v1;
856:     } else {
857:       vmiddle =v2;
858:     }

860:     info = PetscMalloc( n*sizeof(PetscInt), &vm ); CHKERRQ(info);
861:     
862:     for (i=0; i<n; i++){
863:       if (v1[i] < vmiddle[i] && vmiddle[i] < v2[i]) {vm[n_vm]=low+i; n_vm++;}
864:     }

866:     info = VecRestoreArray(VecLow,&v1); CHKERRQ(info);
867:     if (VecLow != VecHigh){
868:       info = VecRestoreArray(VecHigh,&v2); CHKERRQ(info);
869:     }
870:     if ( V != VecLow && V != VecHigh){
871:       info = VecRestoreArray(V,&vmiddle); CHKERRQ(info);
872:     }

874:   } else {

876:     n_vm=0; vm=NULL;

878:   }

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

883:   if (vm) {
884:     info = PetscFree(vm); CHKERRQ(info);
885:   }

887:   return(0);
888: }


893: int VecWhichBetweenOrEqual(Vec VecLow, Vec V, Vec VecHigh, IS * S)
894: {
895:   /* 
896:      Creates an index set with the indices of V whose 
897:      elements are stictly between the corresponding elements 
898:      of the vector VecLow and the Vector VecHigh
899:   */
900:   int info;
901:   PetscInt n,low,high,low2,high2,low3,high3,n_vm=0,i;
902:   PetscInt *vm;
903:   PetscScalar *v1,*v2,*vmiddle;
904:   MPI_Comm comm;


909:   info = VecGetOwnershipRange(VecLow, &low, &high); CHKERRQ(info);
910:   info = VecGetOwnershipRange(VecHigh, &low2, &high2); CHKERRQ(info);
911:   info = VecGetOwnershipRange(V, &low3, &high3); CHKERRQ(info);

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

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

918:   if (n>0){

920:     info = VecGetArray(VecLow,&v1); CHKERRQ(info);
921:     if (VecLow != VecHigh){
922:       info = VecGetArray(VecHigh,&v2); CHKERRQ(info);
923:     } else {
924:       v2=v1;
925:     }
926:     if ( V != VecLow && V != VecHigh){
927:       info = VecGetArray(V,&vmiddle); CHKERRQ(info);
928:     } else if ( V==VecLow ){
929:       vmiddle=v1;
930:     } else {
931:       vmiddle =v2;
932:     }

934:     info = PetscMalloc( n*sizeof(PetscInt), &vm ); CHKERRQ(info);
935:     
936:     for (i=0; i<n; i++){
937:       if (v1[i] <= vmiddle[i] && vmiddle[i] <= v2[i]) {vm[n_vm]=low+i; n_vm++;}
938:     }

940:     info = VecRestoreArray(VecLow,&v1); CHKERRQ(info);
941:     if (VecLow != VecHigh){
942:       info = VecRestoreArray(VecHigh,&v2); CHKERRQ(info);
943:     }
944:     if ( V != VecLow && V != VecHigh){
945:       info = VecRestoreArray(V,&vmiddle); CHKERRQ(info);
946:     }

948:   } else {

950:     n_vm=0; vm=NULL;

952:   }

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

957:   if (vm) {
958:     info = PetscFree(vm); CHKERRQ(info);
959:   }

961:   return(0);
962: }


967: /*@C
968:    ISCreateComplement - Creates the complement of the the index set

970:    Input Parameter:
971: +  S -  a PETSc IS
972: -  V - the reference vector space

974:    Output Parameter:
975: .  T -  the complement of S


978: .seealso ISCreateGeneral()

980:    Level: advanced
981: @*/
982: int ISCreateComplement(IS S, Vec V, IS *T){
983:   int info;
984:   PetscInt i,nis,nloc,high,low,n=0;
985:   const PetscInt *s;
986:   PetscInt *tt,*ss;
987:   MPI_Comm comm;


993:   info = VecGetOwnershipRange(V,&low,&high); CHKERRQ(info);
994:   info = VecGetLocalSize(V,&nloc); CHKERRQ(info);
995:   info = ISGetLocalSize(S,&nis); CHKERRQ(info);
996:   info = ISGetIndices(S, &s); CHKERRQ(info);
997:   info = PetscMalloc( nloc*sizeof(PetscInt),&tt ); CHKERRQ(info);
998:   info = PetscMalloc( nloc*sizeof(PetscInt),&ss ); CHKERRQ(info);

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

1002:   for (i=0; i<nis; i++){ tt[s[i]-low] = -2; }
1003:   
1004:   for (i=0; i<nloc; i++){
1005:     if (tt[i]>-1){ ss[n]=tt[i]; n++; }
1006:   }

1008:   info = ISRestoreIndices(S, &s); CHKERRQ(info);
1009:   
1010:   info = PetscObjectGetComm((PetscObject)S,&comm);CHKERRQ(info);
1011:   info = ISCreateGeneral(comm,n,ss,T);CHKERRQ(info);
1012:   
1013:   if (tt) {
1014:     info = PetscFree(tt); CHKERRQ(info);
1015:   }
1016:   if (ss) {
1017:     info = PetscFree(ss); CHKERRQ(info);
1018:   }

1020:   return(0);
1021: }
1022: #endif