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 *);
20: int VecPow(Vec, double);
22: static int petscminmaxevent=0;
26: /*@C
27: TaoWrapPetscVec - Creates a new TaoVec object using an existing
28: PETSc Vec.
30: Input Parameter:
31: + V - a PETSc vector
32: - TV - the address of a pointer to a TaoVecPetsc
34: Output Parameter:
35: . TV - pointer to a new TaoVecPetsc
37: Note: A TaoVecPetsc is an object with the methods of an abstract
38: TaoVec object. A TaoVecPetsc contains an implementation of the TaoVec
39: methods. Routines using these vectors should declare a pointer to
40: a TaoVec, assign this pointer to the address of a TaoVec object,
41: use the pointer to invoke methods on the object, and use this pointer
42: as an argument when calling other routines. This usage is different
43: from the usage of a PETSc Vec. In PETSc, applications will typically
44: declare a Vec, and pass it as an argument into routines. That is,
45: applications will typically declare a pointer to a TaoVec and use the
46: pointer, or declare a Vec and use it directly.
47:
48: Note:
49: The user is responsible for destroying the Vec V, in addition to
50: to TaoVecPetsc vector TV. The Vec can be destroyed immediately
51: after this routine.
53: Level: developer
55: .seealso TaoVecGetPetscVec(), TaoVecDestroy()
57: @*/
58: int TaoWrapPetscVec( Vec V, TaoVecPetsc* *TV){
59: int info;
61: *TV = new TaoVecPetsc(V);
62: info=(*TV)->SetVec(V); CHKERRQ(info);
63: return(0);
64: }
68: TaoVecPetsc::TaoVecPetsc(Vec VV)
69: {
70: pv=0;
71: pvecviewer=0;
72: if (VV) {
73: SetVec(VV);
74: }
75: if (petscminmaxevent==0){
76: PetscLogEventRegister("PointwiseMinMax",VEC_COOKIE,&petscminmaxevent);
77: }
78: return;
79: }
83: TaoVecPetsc::~TaoVecPetsc()
84: {
85: SetVec(0);
86: return;
87: }
91: /*@C
92: SetVec - Set or replace the underlying Vec object in this TaoVec.
94: Input Parameter:
95: . V - a PETSc vector
97: Note:
98: The user us repsonsible for destroying the Vec V, in addition to
99: to TaoVecPetsc vector TV. The Vec can be destroyed immediately
100: after this routine.
102: Level: developer
104: .seealso TaoWrapPetscVec(), TaoVecDestroy()
106: @*/
107: int TaoVecPetsc::SetVec(Vec VV){
108: int info;
110: if (VV){
112: PetscObjectReference((PetscObject)VV);
113: // info = PetscInfo(VV,"Set the PETSc Vec in a TaoVec .\n"); CHKERRQ(info);
114: }
115: if (pv) {
116: info=VecDestroy(pv); CHKERRQ(info);
117: }
118: pv=VV;
119: return(0);
120: }
124: int TaoVecPetsc::Clone(TaoVec**ntv){
125: int info;
126: TaoVecPetsc *nptv;
127: Vec npv;
129: info=VecDuplicate(pv,&npv); CHKERRQ(info);
130: info = TaoWrapPetscVec(npv,&nptv); CHKERRQ(info);
131: *ntv=nptv;
132: info = VecDestroy(npv); CHKERRQ(info);
133: nptv->pvecviewer=pvecviewer;
134: return(0);
135: }
140: int TaoVecPetsc::Compatible(TaoVec *tv, TaoTruth *flag){
141: TaoVecPetsc *vv = (TaoVecPetsc *)(tv);
142: PetscTruth flg=PETSC_TRUE;
143: int info=0;
146: if (vv==0){
147: *flag=TAO_FALSE;
148: return(0);
149: }
150: info = VecCompare(pv, vv->pv, &flg); CHKERRQ(info);
151: if (flg==PETSC_FALSE){
152: *flag=TAO_FALSE;
153: } else {
154: *flag=TAO_TRUE;
155: }
156: return(0);
157: }
161: int TaoVecPetsc::SetToConstant( double cc ){
162: int info;
163: PetscScalar c=cc;
165: info=VecSet(pv, c); CHKERRQ(info);
166: return(0);
167: }
171: int TaoVecPetsc::SetToZero(){
172: PetscScalar zero=0.0;
173: int info;
175: info=VecSet(pv, zero); CHKERRQ(info);
176: return(0);
177: }
181: int TaoVecPetsc::CopyFrom(TaoVec* tv)
182: {
183: TaoVecPetsc *mv = dynamic_cast <TaoVecPetsc *> (tv);
184: int info;
187: info = VecCopy(mv->pv, pv); CHKERRQ(info);
188: return(0);
189: }
193: int TaoVecPetsc::ScaleCopyFrom(double a, TaoVec *tv)
194: {
195: TaoVecPetsc *mv = dynamic_cast <TaoVecPetsc *> (tv);
196: PetscScalar alpha=a, zero=0.0;
197: int info;
199: info=VecAXPBY(pv,alpha,zero,mv->pv); CHKERRQ(info);
200: return(0);
201: }
205: int TaoVecPetsc::NormInfinity(double *vnorm){
206: int info;
207: PetscReal vv;
209: info=VecNorm(pv,NORM_INFINITY,&vv); CHKERRQ(info);
210: *vnorm=(double)vv;
211: return(0);
212: }
216: int TaoVecPetsc::Norm1(double *vnorm){
217: int info;
218: PetscReal vv;
220: info=VecNorm(pv,NORM_1,&vv); CHKERRQ(info);
221: *vnorm=(double)vv;
222: return(0);
223: }
227: int TaoVecPetsc::Norm2(double *vnorm){
228: int info;
229: PetscReal vv;
231: info=VecNorm(pv,NORM_2,&vv); CHKERRQ(info);
232: *vnorm=(double)vv;
233: return(0);
234: }
238: int TaoVecPetsc::Norm2squared(double *vnorm){
239: int info;
240: PetscReal vv;
242: info=VecNorm(pv,NORM_2,&vv); CHKERRQ(info);
243: *vnorm=(double)(vv*vv);
244: return(0);
245: }
249: int TaoVecPetsc::Scale( double alpha ){
250: int info;
251: PetscScalar aalpha=alpha;
253: info=VecScale(pv, aalpha); CHKERRQ(info);
254: return(0);
255: }
260: int TaoVecPetsc::Axpy(double alpha, TaoVec *tv)
261: {
262: TaoVecPetsc *mv = dynamic_cast <TaoVecPetsc *> (tv);
263: PetscScalar aalpha=alpha;
264: int info;
266: info=VecAXPY(pv, aalpha, mv->pv); CHKERRQ(info);
267: return(0);
268: }
273: int TaoVecPetsc::Axpby(double alpha, TaoVec *tv, double beta)
274: {
275: TaoVecPetsc *mv = dynamic_cast <TaoVecPetsc *> (tv);
276: PetscScalar aalpha = alpha, bbeta = beta;
277: int info;
279: info=VecAXPBY(pv,aalpha,bbeta,mv->pv); CHKERRQ(info);
280: return(0);
281: }
285: int TaoVecPetsc::Aypx(double alpha, TaoVec *tv)
286: {
287: TaoVecPetsc *mv = dynamic_cast <TaoVecPetsc *> (tv);
288: PetscScalar aalpha = alpha;
289: int info;
291: info=VecAYPX(pv,aalpha,mv->pv); CHKERRQ(info);
292: return(0);
293: }
294:
297: int TaoVecPetsc::AddConstant( double cc )
298: {
299: int info;
300: PetscScalar c=cc;
302: info=VecShift(pv, c); CHKERRQ(info);
303: return(0);
304: }
308: int TaoVecPetsc::Dot( TaoVec* tv, double *vDotv )
309: {
310: TaoVecPetsc *mv = dynamic_cast <TaoVecPetsc *> (tv);
311: PetscScalar c;
312: int info;
314: info=VecDot(pv, mv->pv, &c); CHKERRQ(info);
315: *vDotv=c;
316: return(0);
317: }
321: int TaoVecPetsc::Negate(){
322: PetscScalar m1=-1.0;
323: int info;
325: info=VecScale(pv, m1); CHKERRQ(info);
326: return(0);
327: }
331: int TaoVecPetsc::Reciprocal(){
332: int info;
334: info=VecReciprocal(pv); CHKERRQ(info);
335: return(0);
336: }
341: int TaoVecPetsc::Sqrt(){
342: int info;
344: info=VecSqrt(pv); CHKERRQ(info);
345: return(0);
346: }
351: int TaoVecPetsc::Pow(double p){
352: int info;
354: info=VecPow(pv, p); CHKERRQ(info);
355: return(0);
356: }
360: int TaoVecPetsc::GetDimension(TaoInt *n){
361: int info;
363: info=VecGetSize(pv,n); CHKERRQ(info);
364: return(0);
365: }
369: int TaoVecPetsc::PointwiseMultiply(TaoVec *tv, TaoVec *tw)
370: {
371: TaoVecPetsc *mv = dynamic_cast <TaoVecPetsc *> (tv);
372: TaoVecPetsc *mw = dynamic_cast <TaoVecPetsc *> (tw);
373: int info;
375: info=VecPointwiseMult(pv, mv->pv, mw->pv); CHKERRQ(info);
376: return(0);
377: }
382: int TaoVecPetsc::PointwiseDivide( TaoVec* tv , TaoVec* tw)
383: {
384: TaoVecPetsc *mv = dynamic_cast <TaoVecPetsc *> (tv);
385: TaoVecPetsc *mw = dynamic_cast <TaoVecPetsc *> (tw);
386: int info;
388: info=VecPointwiseDivide(pv, mv->pv, mw->pv); CHKERRQ(info);
389: return(0);
390: }
394: int TaoVecPetsc::Median( TaoVec* tv, TaoVec* tw, TaoVec* tx)
395: {
396: TaoVecPetsc *mv = dynamic_cast <TaoVecPetsc *> (tv);
397: TaoVecPetsc *mw = dynamic_cast <TaoVecPetsc *> (tw);
398: TaoVecPetsc *mx = dynamic_cast <TaoVecPetsc *> (tx);
399: int info;
401: info=PetscLogEventBegin(petscminmaxevent,0,0,0,0); CHKERRQ(info);
402: info=VecMedian(mv->pv,mw->pv,mx->pv,pv); CHKERRQ(info);
403: info=PetscLogEventEnd(petscminmaxevent,0,0,0,0); CHKERRQ(info);
404: return(0);
405: }
409: int TaoVecPetsc::PointwiseMinimum( TaoVec* tv, TaoVec* tw)
410: {
411: TaoVecPetsc *mv = dynamic_cast <TaoVecPetsc *> (tv);
412: TaoVecPetsc *mw = dynamic_cast <TaoVecPetsc *> (tw);
413: int info;
415: info=PetscLogEventBegin(petscminmaxevent,0,0,0,0); CHKERRQ(info);
416: info=VecPointwiseMin(pv, mv->pv, mw->pv); CHKERRQ(info);
417: info=PetscLogEventEnd(petscminmaxevent,0,0,0,0); CHKERRQ(info);
418: return(0);
419: }
423: int TaoVecPetsc::Fischer(TaoVec *tx, TaoVec *tf, TaoVec *tl, TaoVec *tu)
424: {
425: TaoVecPetsc *xx = dynamic_cast <TaoVecPetsc *> (tx);
426: TaoVecPetsc *ff = dynamic_cast <TaoVecPetsc *> (tf);
427: TaoVecPetsc *ll = dynamic_cast <TaoVecPetsc *> (tl);
428: TaoVecPetsc *uu = dynamic_cast <TaoVecPetsc *> (tu);
429: int info;
432: info=VecFischer(xx->pv, ff->pv, ll->pv, uu->pv, pv); CHKERRQ(info);
433: return(0);
434: }
438: int TaoVecPetsc::SFischer(TaoVec *tx, TaoVec *tf,
439: TaoVec *tl, TaoVec *tu, double mu)
440: {
441: TaoVecPetsc *xx = dynamic_cast <TaoVecPetsc *> (tx);
442: TaoVecPetsc *ff = dynamic_cast <TaoVecPetsc *> (tf);
443: TaoVecPetsc *ll = dynamic_cast <TaoVecPetsc *> (tl);
444: TaoVecPetsc *uu = dynamic_cast <TaoVecPetsc *> (tu);
445: int info;
448: if ((mu >= -TAO_EPSILON) && (mu <= TAO_EPSILON)) {
449: Fischer(tx, tf, tl, tu);
450: }
451: else {
452: info=VecSFischer(xx->pv, ff->pv, ll->pv, uu->pv, mu, pv); CHKERRQ(info);
453: }
454: return(0);
455: }
459: int TaoVecPetsc::PointwiseMaximum( TaoVec* tv, TaoVec* tw)
460: {
461: TaoVecPetsc *mv = dynamic_cast <TaoVecPetsc *> (tv);
462: TaoVecPetsc *mw = dynamic_cast <TaoVecPetsc *> (tw);
463: int info;
465: info=PetscLogEventBegin(petscminmaxevent,0,0,0,0); CHKERRQ(info);
466: info=VecPointwiseMax(pv, mv->pv, mw->pv); CHKERRQ(info);
467: info=PetscLogEventEnd(petscminmaxevent,0,0,0,0); CHKERRQ(info);
468: return(0);
469: }
474: int TaoVecPetsc::Waxpby ( double aa, TaoVec* tv, double bb, TaoVec* tw)
475: {
476: TaoVecPetsc *mv = dynamic_cast <TaoVecPetsc *> (tv);
477: TaoVecPetsc *mw = dynamic_cast <TaoVecPetsc *> (tw);
478: int info;
479: PetscScalar a=aa, b=bb, zero=0.0;
482: if (a==1){
483: info = VecWAXPY(pv,b,mw->pv,mv->pv); CHKERRQ(info);
484: }
485: else if (b==1) {
486: info = VecWAXPY(pv,a,mv->pv,mw->pv); CHKERRQ(info);
487: }
488: else if (a==0) {
489: info = VecSet(pv, zero); CHKERRQ(info);
490: info = VecAXPY(pv, b, mw->pv); CHKERRQ(info);
491: }
492: else if (b==0) {
493: info = VecSet(pv, zero); CHKERRQ(info);
494: info = VecAXPY(pv, a, mv->pv); CHKERRQ(info);
495: }
496: else {
497: info = VecCopy(mw->pv,pv); CHKERRQ(info);
498: info = VecAXPBY(pv,a,b,mv->pv); CHKERRQ(info);
499: }
500: return(0);
501: }
505: int TaoVecPetsc::AbsoluteValue(){
506: int info;
508: info=VecAbs(pv); CHKERRQ(info);
509: return(0);
510: }
514: int TaoVecPetsc::MinElement(double*val){
515: int info;
516: PetscInt p;
517: PetscScalar a;
519: info = VecMin(pv,&p, &a);CHKERRQ(info);
520: *val=a;
521: return(0);
522: }
527: int TaoVecPetsc::GetArray(TaoScalar **dptr, TaoInt *n){
528: int info;
529: PetscScalar *pptr;
531: if (sizeof(TaoScalar)==sizeof(PetscScalar)){
532: info = VecGetLocalSize(pv,n);CHKERRQ(info);
533: info = VecGetArray(pv,&pptr);CHKERRQ(info);
534: *dptr=(TaoScalar*)pptr;
535: } else {
536: SETERRQ(1,"Incompatible data types");
537: }
538: return(0);
539: }
543: int TaoVecPetsc::RestoreArray(TaoScalar **dptr,TaoInt *n){
544: int info;
545: PetscScalar *pptr;
547: if (sizeof(TaoScalar)==sizeof(PetscScalar)){
548: pptr=(PetscScalar*)(*dptr);
549: info = VecRestoreArray(pv,&pptr); CHKERRQ(info);
550: *dptr=PETSC_NULL;
551: *n=0;
552: } else {
553: SETERRQ(1,"Incompatible data types");
554: }
555: return(0);
556: }
561: int TaoVecPetsc::SetReducedVec(TaoVec *tv, TaoIndexSet *ti)
562: {
563: TaoVecPetsc *mv = dynamic_cast <TaoVecPetsc *> (tv);
564: TaoIndexSetPetsc *TSS = dynamic_cast <TaoIndexSetPetsc *> (ti);
566: int info;
567: PetscInt nn,nnn;
568: TaoPetscISType type;
571: info = TSS->GetSize(&nn); CHKERRQ(info);
572: info = TSS->GetReducedType(&type); CHKERRQ(info);
573: info = mv->GetDimension(&nnn); CHKERRQ(info);
574: if (type==TaoMaskFullSpace){
575: }
576: else if (type==TaoMatrixFree){
577: info = VecCopy(mv->pv,pv); CHKERRQ(info);
578: info = VecISSetToConstant(TSS->GetIS(),0.0,pv);CHKERRQ(info);
579: }
581: info=ReducedCopyFromFull(mv,TSS);CHKERRQ(info);
583: return(0);
584: }
588: int TaoVecPetsc::ReducedXPY(TaoVec *tv, TaoIndexSet *ti)
589: {
590: TaoVecPetsc *mv = dynamic_cast <TaoVecPetsc *> (tv);
591: TaoIndexSetPetsc *TSS = dynamic_cast <TaoIndexSetPetsc *> (ti);
593: VecScatter scatterit;
594: int info;
598: info = TSS->GetReducedVecScatter(mv->pv,pv, &scatterit);CHKERRQ(info);
599:
600: info = VecScatterBegin(scatterit,mv->pv,pv,ADD_VALUES,SCATTER_REVERSE);
601: CHKERRQ(info);
602: info = VecScatterEnd(scatterit,mv->pv,pv,ADD_VALUES,SCATTER_REVERSE);
603: CHKERRQ(info);
605: return(0);
606: }
610: int TaoVecPetsc::ReducedCopyFromFull(TaoVec *tv, TaoIndexSet *ti)
611: {
612: TaoVecPetsc *mv = dynamic_cast <TaoVecPetsc *> (tv);
613: TaoIndexSetPetsc *TSS = dynamic_cast <TaoIndexSetPetsc *> (ti);
615: int info;
616: PetscScalar zero=0.0;
617: VecScatter scatterit;
621: info = VecSet(pv, zero); CHKERRQ(info);
622: info = TSS->GetReducedVecScatter(mv->pv,pv, &scatterit);CHKERRQ(info);
623: info = VecScatterBegin(scatterit,mv->pv,pv,INSERT_VALUES,SCATTER_FORWARD);
624: CHKERRQ(info);
625: info = VecScatterEnd(scatterit,mv->pv,pv,INSERT_VALUES,SCATTER_FORWARD);
626: CHKERRQ(info);
628: return(0);
629: }
635: int TaoVecPetsc::StepMax(TaoVec *tv, double *step)
636: {
637: TaoVecPetsc *mv = dynamic_cast <TaoVecPetsc *> (tv);
638: PetscReal pstep;
639: int info;
641: info = VecStepMax(pv, mv->pv, &pstep);CHKERRQ(info);
642: *step=pstep;
643: return(0);
644: }
648: int TaoVecPetsc::StepBoundInfo(TaoVec *txl ,TaoVec *txu, TaoVec *ts,
649: double *bmin1, double *bmin2, double *bmax)
650: {
651: TaoVecPetsc *mxl = dynamic_cast <TaoVecPetsc *> (txl);
652: TaoVecPetsc *mxu = dynamic_cast <TaoVecPetsc *> (txu);
653: TaoVecPetsc *ms = dynamic_cast <TaoVecPetsc *> (ts);
655: int info;
656: PetscReal p1,p2,p3;
658: info=VecBoundStepInfo(pv,mxl->pv,mxu->pv,ms->pv,&p1,&p2,&p3); CHKERRQ(info);
659: *bmin1=p1; *bmin2=p2; *bmax=p3;
660: return(0);
661: }
665: /*@C
666: SetPetscViewer
668: Input Parameter:
669: . viewer - a viewer object to be used with View() and VecView()
671: Level: advanced
673: @*/
674: int TaoVecPetsc::SetPetscViewer(PetscViewer viewer)
675: {
677: pvecviewer= viewer;
678: return(0);
679: }
684: int TaoVecPetsc::View()
685: {
686: int info;
688: info=VecView(pv,pvecviewer);CHKERRQ(info);
689: return(0);
690: }
694: int TaoVecPetsc::BoundGradientProjection(TaoVec *tg, TaoVec *txl,
695: TaoVec *tx, TaoVec *txu)
696: {
697: TaoVecPetsc *mg = dynamic_cast <TaoVecPetsc *> (tg);
698: TaoVecPetsc *mx = dynamic_cast <TaoVecPetsc *> (tx);
699: TaoVecPetsc *mxl = dynamic_cast <TaoVecPetsc *> (txl);
700: TaoVecPetsc *mxu = dynamic_cast <TaoVecPetsc *> (txu);
701: int info;
703: info = PetscLogEventBegin(petscminmaxevent,0,0,0,0);CHKERRQ(info);
704: info = VecBoundProjection(mg->pv,mx->pv,mxl->pv,mxu->pv,pv);CHKERRQ(info);
705: info = PetscLogEventEnd(petscminmaxevent,0,0,0,0);CHKERRQ(info);
706: return(0);
707: }
711: int TaoVecPetsc::CreateIndexSet(TaoIndexSet**SSS)
712: {
713: TaoIndexSetPetsc *SS;
716: SS = new TaoIndexSetPetsc(pv);
717: *SSS=SS;
718: return(0);
719: }
722: #endif