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),< ); 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), > ); 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