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