Actual source code: taovec_petsc.c
1: #include "tao_general.h"
3: #ifdef TAO_USE_PETSC
5: #include "taovec_petsc.h"
6: #include "../indexset/taois_petsc.h"
8: int VecMedian(Vec, Vec, Vec, Vec);
9: int VecPointwiseMax(Vec, Vec, Vec);
10: int VecPointwiseMin(Vec, Vec, Vec);
11: int VecFischer(Vec, Vec, Vec, Vec, Vec);
12: int VecSFischer(Vec, Vec, Vec, Vec, double, Vec);
13: int VecBoundProjection(Vec, Vec, Vec, Vec, Vec);
14: int VecCreateSubVec(Vec, IS, Vec *);
15: int VecISAXPY(Vec, double, Vec, IS);
16: int VecStepMax(Vec, Vec, double*);
17: int VecStepMax2(Vec, Vec,Vec, Vec, double*);
18: int VecCompare(Vec, Vec, PetscTruth *);
19: int VecBoundStepInfo(Vec, Vec, Vec, Vec, PetscReal *, PetscReal *, PetscReal *);
21: static int petscminmaxevent=0;
25: /*@C
26: TaoWrapPetscVec - Creates a new TaoVec object using an existing
27: PETSc Vec.
29: Input Parameter:
30: + V - a PETSc vector
31: - TV - the address of a pointer to a TaoVecPetsc
33: Output Parameter:
34: . TV - pointer to a new TaoVecPetsc
36: Note: A TaoVecPetsc is an object with the methods of an abstract
37: TaoVec object. A TaoVecPetsc contains an implementation of the TaoVec
38: methods. Routines using these vectors should declare a pointer to
39: a TaoVec, assign this pointer to the address of a TaoVec object,
40: use the pointer to invoke methods on the object, and use this pointer
41: as an argument when calling other routines. This usage is different
42: from the usage of a PETSc Vec. In PETSc, applications will typically
43: declare a Vec, and pass it as an argument into routines. That is,
44: applications will typically declare a pointer to a TaoVec and use the
45: pointer, or declare a Vec and use it directly.
46:
47: Note:
48: The user us repsonsible for destroying the Vec V, in addition to
49: to TaoVecPetsc vector TV. The Vec can be destroyed immediately
50: after this routine.
52: Level: developer
54: .seealso TaoVecGetPetscVec(), TaoVecDestroy()
56: @*/
57: int TaoWrapPetscVec( Vec V, TaoVecPetsc* *TV){
58: int info;
60: *TV = new TaoVecPetsc(V);
61: info=(*TV)->SetVec(V); CHKERRQ(info);
62: return(0);
63: }
67: /*@C
68: TaoVecGetPetscVec - Set the address of a Vec equal to the location
69: of the underlying Vec structure in this TaoVecPetsc object.
71: Input Parameter:
72: + TV - the TaoVecPetsc
73: - V - the address of Vec
75: Output Parameter:
76: . V - the address of a specific Vec
78: Note:
79: This routine does not create a Vec. It sets a pointer
80: to the location of an existing Vec.
82: Note:
83: The underlying PETSc Vec may also be obtained by
84: Vec VV = (Vec) TV->VecObject;
86: Note:
87:
88: VecObject is a void* pointer that can be cast to the underlying Vec object.
90: Level: advanced
92: @*/
93: int TaoVecGetPetscVec( TaoVec* TV, Vec *V){
95: if (TV && V){
96: *V=(Vec)(TV->VecObject);
97: } else if (V){
98: *V=0;
99: }
100: return(0);
101: }
105: TaoVecPetsc::TaoVecPetsc(Vec VV){
106: this->pv=0;
107: this->pvecviewer=0;
108: if (VV){
109: this->SetVec(VV);
110: }
111: if (petscminmaxevent==0){
112: PetscLogEventRegister(&petscminmaxevent,"PointwiseMinMax",VEC_COOKIE);
113: }
114: return;
115: }
119: TaoVecPetsc::~TaoVecPetsc(){
120: this->SetVec(0);
121: return;
122: }
126: /*@C
127: SetVec - Set or replace the underlying Vec object in this TaoVec.
129: Input Parameter:
130: . V - a PETSc vector
132: Note:
133: The user us repsonsible for destroying the Vec V, in addition to
134: to TaoVecPetsc vector TV. The Vec can be destroyed immediately
135: after this routine.
137: Level: developer
139: .seealso TaoWrapPetscVec(), TaoVecDestroy()
141: @*/
142: int TaoVecPetsc::SetVec(Vec VV){
143: int info;
145: if (VV){
147: PetscObjectReference((PetscObject)VV);
148: // info = PetscLogInfo((VV,"Set the PETSc Vec in a TaoVec .\n")); CHKERRQ(info);
149: }
150: if (this->pv){
151: info=VecDestroy(this->pv); CHKERRQ(info);
152: }
153: this->pv=VV;
154: this->VecObject=(void*)VV;
155: return(0);
156: }
160: int TaoVecPetsc::Clone(TaoVec**ntv){
161: int info;
162: TaoVecPetsc *nptv;
163: Vec npv;
165: info=VecDuplicate(pv,&npv); CHKERRQ(info);
166: info = TaoWrapPetscVec(npv,&nptv); CHKERRQ(info);
167: *ntv=nptv;
168: info = VecDestroy(npv); CHKERRQ(info);
169: nptv->pvecviewer=this->pvecviewer;
170: return(0);
171: }
176: int TaoVecPetsc::Compatible(TaoVec *tv, TaoTruth *flag){
177: TaoVecPetsc *vv = (TaoVecPetsc *)(tv);
178: PetscTruth flg=PETSC_TRUE;
179: int info=0;
182: if (vv==0){
183: *flag=TAO_FALSE;
184: return(0);
185: }
186: info = VecCompare(pv, vv->pv, &flg); CHKERRQ(info);
187: if (flg==PETSC_FALSE){
188: *flag=TAO_FALSE;
189: } else {
190: *flag=TAO_TRUE;
191: }
192: return(0);
193: }
197: int TaoVecPetsc::SetToConstant( double cc ){
198: int info;
199: PetscScalar c=cc;
201: info=VecSet(pv, c); CHKERRQ(info);
202: return(0);
203: }
207: int TaoVecPetsc::SetToZero(){
208: PetscScalar zero=0.0;
209: int info;
211: info=VecSet(pv, zero); CHKERRQ(info);
212: return(0);
213: }
217: int TaoVecPetsc::CopyFrom( TaoVec* tv ){
218: Vec pv2 = (Vec)(tv->VecObject);
219: int info;
221: info=VecCopy(pv2,pv); CHKERRQ(info);
222: return(0);
223: }
227: int TaoVecPetsc::ScaleCopyFrom( double a, TaoVec* tv ){
228: Vec v = (Vec) tv->VecObject;
229: PetscScalar alpha=a, zero=0.0;
230: int info;
232: info=VecAXPBY(this->pv,alpha,zero,v); CHKERRQ(info);
233: return(0);
234: }
238: int TaoVecPetsc::NormInfinity(double *vnorm){
239: int info;
240: PetscReal vv;
242: info=VecNorm(pv,NORM_INFINITY,&vv); CHKERRQ(info);
243: *vnorm=(double)vv;
244: return(0);
245: }
249: int TaoVecPetsc::Norm1(double *vnorm){
250: int info;
251: PetscReal vv;
253: info=VecNorm(pv,NORM_1,&vv); CHKERRQ(info);
254: *vnorm=(double)vv;
255: return(0);
256: }
260: int TaoVecPetsc::Norm2(double *vnorm){
261: int info;
262: PetscReal vv;
264: info=VecNorm(pv,NORM_2,&vv); CHKERRQ(info);
265: *vnorm=(double)vv;
266: return(0);
267: }
271: int TaoVecPetsc::Norm2squared(double *vnorm){
272: int info;
273: PetscReal vv;
275: info=VecNorm(pv,NORM_2,&vv); CHKERRQ(info);
276: *vnorm=(double)(vv*vv);
277: return(0);
278: }
282: int TaoVecPetsc::Scale( double alpha ){
283: int info;
284: PetscScalar aalpha=alpha;
286: info=VecScale(pv, aalpha); CHKERRQ(info);
287: return(0);
288: }
293: int TaoVecPetsc::Axpy( double alpha, TaoVec* tv ){
294: Vec v = (Vec) tv->VecObject;
295: PetscScalar aalpha=alpha;
296: int info;
298: info=VecAXPY(this->pv, aalpha, v); CHKERRQ(info);
299: return(0);
300: }
305: int TaoVecPetsc::Axpby( double alpha, TaoVec* tv, double beta ){
306: Vec v = (Vec) tv->VecObject;
307: PetscScalar aalpha = alpha, bbeta = beta;
308: int info;
310: info=VecAXPBY(this->pv,aalpha,bbeta,v); CHKERRQ(info);
311: return(0);
312: }
316: int TaoVecPetsc::Aypx( double alpha, TaoVec* x ){
317: Vec X = (Vec) x->VecObject;
318: PetscScalar aalpha = alpha;
319: int info;
321: info=VecAYPX(this->pv,aalpha,X); CHKERRQ(info);
322: return(0);
323: }
324:
327: int TaoVecPetsc::AddConstant( double cc ){
328: int info;
329: PetscScalar c=cc;
331: info=VecShift(pv, c); CHKERRQ(info);
332: return(0);
333: }
337: int TaoVecPetsc::Dot( TaoVec* tv, double *vDotv ){
338: Vec v = (Vec) tv->VecObject;
339: PetscScalar c;
340: int info;
342: info=VecDot(v,this->pv,&c); CHKERRQ(info);
343: *vDotv=c;
344: return(0);
345: }
349: int TaoVecPetsc::Negate(){
350: PetscScalar m1=-1.0;
351: int info;
353: info=VecScale(pv, m1); CHKERRQ(info);
354: return(0);
355: }
359: int TaoVecPetsc::Reciprocal(){
360: int info;
362: info=VecReciprocal(pv); CHKERRQ(info);
363: return(0);
364: }
368: int TaoVecPetsc::GetDimension(int *n){
369: int info;
371: info=VecGetSize(pv,n); CHKERRQ(info);
372: return(0);
373: }
377: int TaoVecPetsc::PointwiseMultiply( TaoVec* tv, TaoVec* tw ){
378: Vec v = (Vec) tv->VecObject;
379: Vec w = (Vec) tw->VecObject;
380: int info;
382: info=VecPointwiseMult(this->pv, w, v); CHKERRQ(info);
383: return(0);
384: }
389: int TaoVecPetsc::PointwiseDivide( TaoVec* tv , TaoVec* tw){
390: Vec v = (Vec) tv->VecObject;
391: Vec w = (Vec) tw->VecObject;
392: int info;
394: info=VecPointwiseDivide(this->pv, v, w); CHKERRQ(info);
395: return(0);
396: }
400: int TaoVecPetsc::Median( TaoVec* tv, TaoVec* tw, TaoVec* tx){
401: Vec v = (Vec) tv->VecObject;
402: Vec w = (Vec) tw->VecObject;
403: Vec x = (Vec) tx->VecObject;
404: int info;
406: info=PetscLogEventBegin(petscminmaxevent,0,0,0,0); CHKERRQ(info);
407: info=VecMedian(v,w,x,this->pv); CHKERRQ(info);
408: info=PetscLogEventEnd(petscminmaxevent,0,0,0,0); CHKERRQ(info);
409: return(0);
410: }
414: int TaoVecPetsc::PointwiseMinimum( TaoVec* tv, TaoVec* tw){
415: Vec v = (Vec) tv->VecObject;
416: Vec w = (Vec) tw->VecObject;
417: int info;
419: info=PetscLogEventBegin(petscminmaxevent,0,0,0,0); CHKERRQ(info);
420: info=VecPointwiseMin(this->pv, w, v); CHKERRQ(info);
421: info=PetscLogEventEnd(petscminmaxevent,0,0,0,0); CHKERRQ(info);
422: return(0);
423: }
427: int TaoVecPetsc::Fischer(TaoVec *tx, TaoVec *tf, TaoVec *tl, TaoVec *tu)
428: {
429: TaoVecPetsc *xx = (TaoVecPetsc *)(tx);
430: TaoVecPetsc *ff = (TaoVecPetsc *)(tf);
431: TaoVecPetsc *ll = (TaoVecPetsc *)(tl);
432: TaoVecPetsc *uu = (TaoVecPetsc *)(tu);
433: int info;
436: info=VecFischer(xx->pv, ff->pv, ll->pv, uu->pv, pv); CHKERRQ(info);
437: return(0);
438: }
442: int TaoVecPetsc::SFischer(TaoVec *tx, TaoVec *tf,
443: TaoVec *tl, TaoVec *tu, double mu)
444: {
445: TaoVecPetsc *xx = (TaoVecPetsc *)(tx);
446: TaoVecPetsc *ff = (TaoVecPetsc *)(tf);
447: TaoVecPetsc *ll = (TaoVecPetsc *)(tl);
448: TaoVecPetsc *uu = (TaoVecPetsc *)(tu);
449: int info;
452: if ((mu >= -TAO_EPSILON) && (mu <= TAO_EPSILON)) {
453: Fischer(tx, tf, tl, tu);
454: }
455: else {
456: info=VecSFischer(xx->pv, ff->pv, ll->pv, uu->pv, mu, pv); CHKERRQ(info);
457: }
458: return(0);
459: }
463: int TaoVecPetsc::PointwiseMaximum( TaoVec* tv, TaoVec* tw){
464: Vec v = (Vec) tv->VecObject;
465: Vec w = (Vec) tw->VecObject;
466: int info;
468: info=PetscLogEventBegin(petscminmaxevent,0,0,0,0); CHKERRQ(info);
469: info=VecPointwiseMax(this->pv, w, v); CHKERRQ(info);
470: info=PetscLogEventEnd(petscminmaxevent,0,0,0,0); CHKERRQ(info);
471: return(0);
472: }
477: int TaoVecPetsc::Waxpby ( double aa, TaoVec* tv, double bb, TaoVec* tw){
478: Vec v = (Vec) tv->VecObject;
479: Vec w = (Vec) tw->VecObject;
480: int info;
481: PetscScalar a=aa, b=bb, zero=0.0;
484: if (a==1){
485: info = VecWAXPY(this->pv,b,w,v); CHKERRQ(info);
486: } else if (b==1) {
487: info = VecWAXPY(this->pv,a,v,w); CHKERRQ(info);
488: } else if (a==0) {
489: info = VecSet(pv, zero); CHKERRQ(info);
490: info = VecAXPY(this->pv, b, w); CHKERRQ(info);
491: } else if (b==0) {
492: info = VecSet(pv, zero); CHKERRQ(info);
493: info = VecAXPY(this->pv, a, v); CHKERRQ(info);
494: } else {
495: info = VecCopy(w,pv); CHKERRQ(info);
496: info = VecAXPBY(this->pv,a,b,v); CHKERRQ(info);
497: }
498: return(0);
499: }
503: int TaoVecPetsc::AbsoluteValue(){
504: int info;
506: info=VecAbs(pv); CHKERRQ(info);
507: return(0);
508: }
512: int TaoVecPetsc::MinElement(double*val){
513: int info, p;
514: PetscScalar a;
516: info = VecMin(this->pv,&p, &a);CHKERRQ(info);
517: *val=a;
518: return(0);
519: }
524: int TaoVecPetsc::GetArray(TaoScalar **dptr, int *n){
525: int info;
526: PetscScalar *pptr;
528: if (sizeof(TaoScalar)==sizeof(PetscScalar)){
529: info = VecGetLocalSize(pv,n);CHKERRQ(info);
530: info = VecGetArray(pv,&pptr);CHKERRQ(info);
531: *dptr=(TaoScalar*)pptr;
532: } else {
533: SETERRQ(1,"Incompatible data types");
534: }
535: return(0);
536: }
540: int TaoVecPetsc::RestoreArray(TaoScalar **dptr,int *n){
541: int info;
542: PetscScalar *pptr;
544: if (sizeof(TaoScalar)==sizeof(PetscScalar)){
545: pptr=(PetscScalar*)(*dptr);
546: info = VecRestoreArray(pv,&pptr); CHKERRQ(info);
547: *dptr=PETSC_NULL;
548: *n=0;
549: } else {
550: SETERRQ(1,"Incompatible data types");
551: }
552: return(0);
553: }
558: int TaoVecPetsc::SetReducedVec(TaoVec*VV,TaoIndexSet*SS){
559: int info,nn,nnn,nred;
560: TaoPetscISType type;
561: Vec R, V = (Vec) VV->VecObject;
562: TaoIndexSetPetsc *TSS = (TaoIndexSetPetsc*) SS;
563: IS S2;
564: MPI_Comm comm;
567: info = SS->GetSize(&nn); CHKERRQ(info);
568: info = VV->GetDimension(&nnn); CHKERRQ(info);
569: info = TSS->GetReducedType(&type); CHKERRQ(info);
570: if (type==TaoMaskFullSpace){
572: } else if (nn==nnn){
573: info = VecDuplicate(V,&R); CHKERRQ(info);
574: info = this->SetVec(R); CHKERRQ(info);
575: PetscObjectDereference((PetscObject)R);
576: } else {
577: info = TSS->RedistributeIS(&S2); CHKERRQ(info);
578: info = ISGetLocalSize(S2,&nred); CHKERRQ(info);
579: info = PetscObjectGetComm((PetscObject)V,&comm);CHKERRQ(info);
580: info = VecCreateMPI(comm,nred,nn,&R);CHKERRQ(info);
581: info = this->SetVec(R); CHKERRQ(info);
582: PetscObjectDereference((PetscObject)R);
583: }
584: info=this->ReducedCopyFromFull(VV,SS);CHKERRQ(info);
586: return(0);
587: }
591: int TaoVecPetsc::ReducedXPY(TaoVec*VV,TaoIndexSet*SS){
592: VecScatter scatterit;
593: TaoIndexSetPetsc *TSS = (TaoIndexSetPetsc*) SS;
594: Vec R = (Vec) VV->VecObject;
595: Vec V=this->pv;
596: int info;
600: info = TSS->GetReducedVecScatter(R,V, &scatterit);CHKERRQ(info);
601:
602: info = VecScatterBegin(R,V,ADD_VALUES,SCATTER_REVERSE,scatterit);
603: CHKERRQ(info);
604: info = VecScatterEnd(R,V,ADD_VALUES,SCATTER_REVERSE, scatterit);
605: CHKERRQ(info);
607: return(0);
608: }
612: int TaoVecPetsc::ReducedCopyFromFull(TaoVec*VV,TaoIndexSet*SS){
613: int info;
614: PetscScalar zero=0.0;
615: TaoIndexSetPetsc *TSS = (TaoIndexSetPetsc*) SS;
616: VecScatter scatterit;
617: Vec V = (Vec) VV->VecObject;
618: Vec R=this->pv;
622: info = VecSet(R, zero); CHKERRQ(info);
623: info = TSS->GetReducedVecScatter(V,R, &scatterit);CHKERRQ(info);
624: info = VecScatterBegin(V,R,INSERT_VALUES,SCATTER_FORWARD,scatterit);
625: CHKERRQ(info);
626: info = VecScatterEnd(V,R,INSERT_VALUES,SCATTER_FORWARD, scatterit);
627: CHKERRQ(info);
629: return(0);
630: }
636: int TaoVecPetsc::StepMax(TaoVec*tv,double*step){
637: Vec dx = (Vec) tv->VecObject;
638: PetscReal pstep;
639: int info;
641: info = VecStepMax(this->pv, dx,&pstep);CHKERRQ(info);
642: *step=pstep;
643: return(0);
644: }
648: int TaoVecPetsc::StepBoundInfo(TaoVec *TXL ,TaoVec *TXU, TaoVec*S, double* bmin1,double* bmin2, double* bmax){
649: int info;
650: PetscReal p1,p2,p3;
651: Vec X=(Vec)this->VecObject, DX=(Vec)S->VecObject;
652: Vec XL=(Vec)TXL->VecObject, XU=(Vec)TXU->VecObject;
654: info=VecBoundStepInfo(X,XL,XU,DX,&p1,&p2,&p3); CHKERRQ(info);
655: *bmin1=p1; *bmin2=p2; *bmax=p3;
656: return(0);
657: }
661: /*@C
662: SetPetscViewer
664: Input Parameter:
665: . viewer - a viewer object to be used with this->View() and VecView()
667: Level: advanced
669: @*/
670: int TaoVecPetsc::SetPetscViewer(PetscViewer viewer){
672: this->pvecviewer= viewer;
673: return(0);
674: }
679: int TaoVecPetsc::View(){
680: int info;
682: info=VecView(this->pv,this->pvecviewer);CHKERRQ(info);
683: return(0);
684: }
688: int TaoVecPetsc::BoundGradientProjection(TaoVec*GG,TaoVec*XXL,TaoVec*XX, TaoVec*XXU){
689: Vec v = (Vec) XXL->VecObject;
690: Vec w = (Vec) XXU->VecObject;
691: Vec x = (Vec) XX->VecObject;
692: Vec g = (Vec) GG->VecObject;
693: int info;
695: info = PetscLogEventBegin(petscminmaxevent,0,0,0,0);CHKERRQ(info);
696: info = VecBoundProjection(g,x,v,w,this->pv);CHKERRQ(info);
697: info = PetscLogEventEnd(petscminmaxevent,0,0,0,0);CHKERRQ(info);
698: return(0);
699: }
703: int TaoVecPetsc::CreateIndexSet(TaoIndexSet**SSS){
704: TaoIndexSetPetsc *SS;
707: SS = new TaoIndexSetPetsc(this->pv);
708: *SSS=SS;
709: return(0);
710: }
713: #endif