Actual source code: taovec_ga.c
1: /**************************************************************
2: File: taovec_ga.c
4: TAO/GLobal Array Project
6: Author: Limin Zhang, Ph.D.
7: Mathematics Department
8: Columbia Basin College
9: Pasco, WA 99301
10: Limin.Zhang@cbc2.org
12: Mentor: Jarek Naplocha, Ph.D.
13: Environmental Molecular Science Laboratory
14: Pacific Northwest National Laboratory
15: Richland, WA 99352
17: Purpose:
18: to design and implement some interfaces between TAO and
19: global arrays.
21: Date: 4/22/2002
24: Revise history:
26: 7/15/02
27: replace pv with this->pv
28: introduce GAVec as int
30: 8/8/02
31: clean Petsc function calls and marcos.
33: **************************************************************/
35: #include "math.h"
38: #include "tao_general.h"
39: #include "taovec_ga.h" /*I "taovec_ga.h" I*/
44: /*@C
45: TaoWrapGaVec - Creates a new TaoVec object using an existing
46: GlobalArray vector.
48: Input Parameter:
49: + V - a GlobalArray vector
50: - TV - the address of a pointer to a TaoGaVec
52: Output Parameter:
53: . TV - pointer to a new TaoGaVec
55: Level: advanced
57: .seealso TaoVecGetGaVec(), TaoWrapPetscVec()
58: @*/
59: int TaoWrapGaVec (GAVec V, TaoVecGa ** TV)
60: {
61: TaoFunctionBegin;
62: *TV = new TaoVecGa (V);
63: TaoFunctionReturn(0);
64: }
68: /*@C
69: TaoVecGetGaVec - If the TaoVec is of the TaoVecGa type, this routine gets the underlying Ga Vec.
71: Input Parameter:
72: . TV - the TaoVecGa
74: Output Parameter:
75: . V - the Ga vector
78: Level: advanced
80: .seealso TaoWrapGaVec()
82: @*/
84: int TaoVecGetGaVec (TaoVec * TV, GAVec * V)
85: {
86: TaoFunctionBegin;
87: if (TV)
88: {
89: TaoVecGa *tmp = dynamic_cast<TaoVecGa *>(TV);
90: if (!tmp) {
91: SETERRQ(1,"Could not cast from TaoVec* to TaoVecGa*.");
92: }
93: *V = tmp->GetVec ();
94: }
95: else
96: {
97: *V = 0;
98: }
99: TaoFunctionReturn(0);
100: }
104: TaoVecGa::TaoVecGa (GAVec VV)
105: {
106: this->pv = VV;
107: this->VecObject = (void *) (&VV);
108: return;
109: }
113: int
114: TaoVecGa::Clone (TaoVec ** ntv)
115: {
116: TaoVecGa *nptv;
117: GAVec npv;
119: TaoFunctionBegin;
120: npv = GA_Duplicate (this->pv, "Cloned"); //pv is the ga handle of the current object.
121: if (!npv)
122: GA_Error ("TaoVecGa::Clone:duplicate failed: ", npv);
123: TaoWrapGaVec (npv, &nptv);
124: *ntv = nptv;
125: npv = 0;
126: TaoFunctionReturn(0);
127: }
132: int TaoVecGa::Compatible (TaoVec * tv, TaoTruth *flag)
133: {
134: int info;
135: int pv2 = 0;
137: TaoFunctionBegin;
138: info = TaoVecGetGaVec (tv, &pv2); CHKERRQ(info);
139: info = GA_Compare_distr (this->pv, pv2);
141: if (info == 0)
142: *flag = TAO_TRUE;
143: else
144: *flag = TAO_FALSE;
146: TaoFunctionReturn(0);
147: }
152: int
153: TaoVecGa::SetToConstant (TaoScalar c)
154: {
155: TaoFunctionBegin;
156: GA_Fill (this->pv, &c);
157: TaoFunctionReturn(0);
158: }
163: int
164: TaoVecGa::SetToZero ()
165: {
166: TaoFunctionBegin;
167: GA_Zero (this->pv);
168: TaoFunctionReturn(0);
169: }
173: int
174: TaoVecGa::CopyFrom (TaoVec * tv)
175: {
176: int pv2 = 0;
177: TaoFunctionBegin;
178: TaoVecGetGaVec (tv, &pv2);
179: GA_Copy (pv2, this->pv);
180: TaoFunctionReturn(0);
181: }
186: int
187: TaoVecGa::ScaleCopyFrom (TaoScalar a, TaoVec * tv)
188: {
189: int info, pv2 = 0;
190: TaoFunctionBegin;
191: info = TaoVecGetGaVec (tv, &pv2);
192: GA_Copy (pv2, this->pv);
193: GA_Scale (this->pv, &a);
194: TaoFunctionReturn(0);
195: }
201: int
202: TaoVecGa::NormInfinity (TaoScalar * vnorm)
203: {
204: TaoFunctionBegin;
205: GA_Norm_infinity (this->pv, vnorm);
206: TaoFunctionReturn(0);
207: }
211: int
212: TaoVecGa::Norm1 (TaoScalar * vnorm)
213: {
214: TaoFunctionBegin;
215: GA_Norm1 (this->pv, vnorm);
216: TaoFunctionReturn(0);
217: }
222: int
223: TaoVecGa::Norm2 (TaoScalar * vnorm)
224: {
225: TaoFunctionBegin;
226: *vnorm = sqrt (GA_Ddot (this->pv, this->pv));
227: TaoFunctionReturn(0);
228: }
232: int
233: TaoVecGa::Norm2squared (TaoScalar * vnorm)
234: {
235: TaoFunctionBegin;
236: *vnorm = GA_Ddot (this->pv, this->pv);
237: TaoFunctionReturn(0);
238: }
242: int
243: TaoVecGa::Scale (TaoScalar alpha)
244: {
245: TaoFunctionBegin;
246: GA_Scale (this->pv, &alpha);
247: TaoFunctionReturn(0);
248: }
253: int
254: TaoVecGa::Axpy (TaoScalar alpha, TaoVec * tv)
255: {
256: int info, pv2 = 0;
257: TaoScalar a = 1.0, b = alpha;
258: TaoFunctionBegin;
259: info = TaoVecGetGaVec (tv, &pv2); CHKERRQ(info);
260: GA_Add (&a, this->pv, &b, pv2, this->pv);
261: TaoFunctionReturn(0);
262: }
267: int
268: TaoVecGa::Axpby (TaoScalar alpha, TaoVec * tv, TaoScalar beta)
269: {
270: int info, pv2 = 0;
271: TaoScalar a = alpha, b = beta;
272: TaoFunctionBegin;
274: info = TaoVecGetGaVec (tv, &pv2); CHKERRQ(info);
275: GA_Add (&a, pv2, &b, this->pv, this->pv);
276: TaoFunctionReturn(0);
277: }
282: int
283: TaoVecGa::Aypx (TaoScalar alpha, TaoVec * x)
284: {
285: int info, pv2 = 0;
286: TaoScalar a = alpha, b = 1.0;
287: TaoFunctionBegin;
289: info = TaoVecGetGaVec (x, &pv2); CHKERRQ(info);
290: GA_Add (&a, this->pv, &b, pv2, this->pv);
291: TaoFunctionReturn(0);
292: }
297: int
298: TaoVecGa::AddConstant (TaoScalar c)
299: {
300: TaoFunctionBegin;
301: GA_Add_constant (this->pv, &c);
302: TaoFunctionReturn(0);
303: }
307: int
308: TaoVecGa::Dot (TaoVec * tv, TaoScalar * vDotv)
309: {
310: int info, pv2 = 0;
311: TaoFunctionBegin;
312: info = TaoVecGetGaVec (tv, &pv2); CHKERRQ(info);
313: *vDotv = GA_Ddot (this->pv, pv2);
314: TaoFunctionReturn(0);
315: }
319: int
320: TaoVecGa::Negate ()
321: {
322: TaoScalar m1 = -1.0;
323: TaoFunctionBegin;
324: GA_Scale (this->pv, &m1);
325: TaoFunctionReturn(0);
326: }
330: int
331: TaoVecGa::Reciprocal ()
332: {
333: TaoFunctionBegin;
334: GA_Recip (this->pv);
335: TaoFunctionReturn(0);
336: }
341: int
342: TaoVecGa::GetDimension (int *n)
343: {
344: int type, ndim, dims;
345: TaoFunctionBegin;
346: NGA_Inquire (this->pv, &type, &ndim, &dims);
347: if (ndim != 1)
348: SETERRQ(1,"GA Vector has bad dimension");
349: *n = dims;
350: TaoFunctionReturn(0);
351: }
356: int
357: TaoVecGa::PointwiseMultiply (TaoVec * tv, TaoVec * tw)
358: {
359: int info, pv1 = 0, pv2 = 0;
360: TaoFunctionBegin;
361: info = TaoVecGetGaVec (tv, &pv1); CHKERRQ(info);
362: info = TaoVecGetGaVec (tw, &pv2); CHKERRQ(info);
363: GA_Elem_multiply (pv1, pv2, this->pv);
364: TaoFunctionReturn(0);
365: }
370: int
371: TaoVecGa::PointwiseDivide (TaoVec * tv, TaoVec * tw)
372: {
373: int info, pv1 = 0, pv2 = 0;
374: TaoFunctionBegin;
375: info = TaoVecGetGaVec (tv, &pv1); CHKERRQ(info);
376: info = TaoVecGetGaVec (tw, &pv2); CHKERRQ(info);
377: GA_Elem_divide (pv1, pv2, this->pv);
378: TaoFunctionReturn(0);
379: }
383: int
384: TaoVecGa::Median (TaoVec * tv, TaoVec * tw, TaoVec * tx)
385: {
386: int info, pv1 = 0, pv2 = 0, pv3 = 0;
387: TaoFunctionBegin;
388: info = TaoVecGetGaVec (tv, &pv1); CHKERRQ(info);
389: info = TaoVecGetGaVec (tw, &pv2); CHKERRQ(info);
390: info = TaoVecGetGaVec (tx, &pv3); CHKERRQ(info);
392: GA_Median (pv1, pv2, pv3, this->pv);
394: TaoFunctionReturn(0);
395: }
399: int
400: TaoVecGa::PointwiseMinimum (TaoVec * tv, TaoVec * tw)
401: {
402: int info, pv1 = 0, pv2 = 0;
403: info = TaoVecGetGaVec (tv, &pv1); CHKERRQ(info);
404: info = TaoVecGetGaVec (tw, &pv2); CHKERRQ(info);
406: GA_Elem_minimum (pv1, pv2, this->pv);
407: TaoFunctionReturn(0);
408: }
412: int
413: TaoVecGa::Fischer (TaoVec * tx, TaoVec * tf, TaoVec * tl, TaoVec * tu)
414: {
415: TaoFunctionBegin;
416: SETERRQ(1," TaoVecGa::Fischer not implemented.");
417: TaoFunctionReturn(0);
418: }
422: int
423: TaoVecGa::PointwiseMaximum (TaoVec * tv, TaoVec * tw)
424: {
425: int info, pv1 = 0, pv2 = 0;
426: TaoFunctionBegin;
427: info = TaoVecGetGaVec (tv, &pv1); CHKERRQ(info);
428: info = TaoVecGetGaVec (tw, &pv2); CHKERRQ(info);
429: GA_Elem_maximum (pv1, pv2, this->pv);
430: TaoFunctionReturn(0);
431: }
435: //Sums two scaled vectors and stores the result in this vector. (this=alpha*xx+beta*yy)
436: int
437: TaoVecGa::Waxpby (TaoScalar a, TaoVec * tv, TaoScalar b, TaoVec * tw)
438: {
439: int info, pv1 = 0, pv2 = 0;
440: TaoFunctionBegin;
441: info = TaoVecGetGaVec (tv, &pv1); CHKERRQ(info);
442: info = TaoVecGetGaVec (tw, &pv2); CHKERRQ(info);
444: if ((a == 0) && (b == 0))
445: GA_Zero (this->pv);
446: else if (a == 0)
447: {
448: GA_Scale (pv2, &b);
449: GA_Copy (pv2, this->pv);
450: }
451: else if (b == 0)
452: {
453: GA_Scale (pv1, &a);
454: GA_Copy (pv1, this->pv);
455: }
456: else
457: GA_Add (&a, pv1, &b, pv2, this->pv);
459: TaoFunctionReturn(0);
460: }
465: int
466: TaoVecGa::AbsoluteValue ()
467: {
468: TaoFunctionBegin;
469: GA_Abs_value (this->pv);
470: TaoFunctionReturn(0);
471: }
476: int
477: TaoVecGa::MinElement (TaoScalar * val)
478: {
479: int index;
480: TaoFunctionBegin;
481: NGA_Select_elem (this->pv, "min", val, &index);
482: TaoFunctionReturn(0);
483: }
485: /*
487: int TaoVecGa::GetArray (TaoScalar ** dptr, int *n)
488: dptr - a pointer to the pointer to the first memory location on the local patch (output)
489: n - the number of elements of the GA vector at the location pointed by dptr. (output)
490: */
493: int
494: TaoVecGa::GetArray (TaoScalar ** dptr, int *n)
495: {
496: int type, ndim, dims;
497: int lo, hi, ld;
498: TaoScalar *ptr = 0;
500: TaoFunctionBegin;
501: int me = GA_Nodeid ();
503: GA_Sync ();
505: NGA_Inquire (this->pv, &type, &ndim, &dims);
506: if (type != MT_C_DBL) {
507: //the only datatype for taovec is TaoScalar
508: SETERRQ(1,"Global Array is wrong data type");
509: }
511: if (ndim != 1) {
512: //we only deal with one dimensional global array as a vector
513: SETERRQ(1,"Global Array has wrong dimension");
514: }
516: NGA_Distribution (this->pv, me, &lo, &hi);
518: //Get the size of the patch on this processor
519: if (lo >= 0 && hi >= 0) {
520: NGA_Access (this->pv, &lo, &hi, &ptr, &ld);
521: *n = hi - lo + 1;
522: } else
523: *n = 0;
525: *dptr = ptr;
526: TaoFunctionReturn(0);
527: }
531: int
532: TaoVecGa::RestoreArray (TaoScalar ** dptr, int *n)
533: {
534: int lo, hi;
535: int me = GA_Nodeid ();
537: TaoFunctionBegin;
538: GA_Sync ();
539: NGA_Distribution (this->pv, me, &lo, &hi);
541: if (lo != 0 && hi != -1)
542: NGA_Release_update (this->pv, &lo, &hi);
544: GA_Sync ();
546: *dptr = 0;
547: *n = 0;
548: TaoFunctionReturn(0);
549: }
553: int
554: TaoVecGa::SetReducedVec (TaoVec * VV, TaoIndexSet * SS)
555: {
556: TaoFunctionBegin;
557: SETERRQ(1,"TaoVecGa::SetReducedVec not yet implemented");
558: TaoFunctionReturn(0);
559: }
564: int
565: TaoVecGa::ReducedXPY (TaoVec * VV, TaoIndexSet * SS)
566: {
567: TaoFunctionBegin;
568: SETERRQ(1,"TaoVecGa::ReducedXPY not yet implemented.");
569: TaoFunctionReturn(0);
570: }
574: int
575: TaoVecGa::ReducedCopyFromFull (TaoVec * VV, TaoIndexSet * SS)
576: {
577: TaoFunctionBegin;
578: SETERRQ(1,"TaoVecGa::ReducedCopyFromFull not yet implemented.");
579: TaoFunctionReturn(0);
580: }
584: int
585: TaoVecGa::StepMax (TaoVec * tv, double * step)
586: {
587: int info;
588: TaoFunctionBegin;
589: GAVec X = this->GetVec ();
590: GAVec DX = 0;
591: info = TaoVecGetGaVec (tv, &DX); CHKERRQ(info);
592: GA_Step_max (X, DX, step);
593: TaoFunctionReturn(0);
594: }
599: int
600: TaoVecGa::StepMax2 (TaoVec * tv, TaoVec * xl, TaoVec * xu, TaoScalar * step)
601: {
602: int info;
603: TaoFunctionBegin;
604: SETERRQ(1,"TaoVecGa::CreatIndexSet not yet implemented.");
605: TaoFunctionReturn(0);
607: }
613: int
614: TaoVecGa::View ()
615: {
616: TaoFunctionBegin;
617: GA_Print (this->pv);
618: TaoFunctionReturn(0);
619: }
623: int
624: TaoVecGa::CreateIndexSet (TaoIndexSet ** SSS)
625: {
626: TaoFunctionBegin;
627: SETERRQ(1,"TaoVecGa::CreatIndexSet not yet implemented.");
628: TaoFunctionReturn(0);
629: }
633: int
634: TaoVecGa::StepBoundInfo(TaoVec *TXL, TaoVec *TXU, TaoVec *S,
635: double *bmin1, double *bmin2, double *bmax) {
636: int i, n, n2, n3, n4, info;
637: double tt,t1=1.0e+20, t2=1.0e+20, t3=0;
638: GAScalar *x,*xl,*xu,*dx;
640: TaoVecGa *X=this;
642: TaoFunctionBegin;
643: info = X->GetArray((TaoScalar **)&x, &n); CHKERRQ(info);
644: info = TXL->GetArray((TaoScalar **)&xl, &n2); CHKERRQ(info);
645: info = TXU->GetArray((TaoScalar **)&xu, &n3); CHKERRQ(info);
646: info = S->GetArray((TaoScalar **)&dx, &n4); CHKERRQ(info);
648: if ( (n != n2) || (n != n3) || (n != n4) ) {
649: SETERRQ(1,"Vectors must all be same size");
650: }
651:
652: for (i=0;i<n;i++){
653: if (dx[i]>0){
654: tt=(xu[i]-x[i])/dx[i];
655: t1=TaoMin(t1,tt);
656: if (t1>0){
657: t2=TaoMin(t2,tt);
658: }
659: t3=TaoMax(t3,tt);
660: } else if (dx[i]<0){
661: tt=(xl[i]-x[i])/dx[i];
662: t1=TaoMin(t1,tt);
663: if (t1>0){
664: t2=TaoMin(t2,tt);
665: }
666: t3=TaoMax(t3,tt);
667: }
668: }
669:
670: info = X->RestoreArray((TaoScalar **)&x, &n); CHKERRQ(info);
671: info = TXL->RestoreArray((TaoScalar **)&xl, &n2); CHKERRQ(info);
672: info = TXU->RestoreArray((TaoScalar **)&xu, &n3); CHKERRQ(info);
673: info = S->RestoreArray((TaoScalar **)&dx, &n4); CHKERRQ(info);
674:
675: if (bmin1) GA_Dgop(&t1, 1, "min");
676: if (bmin2) GA_Dgop(&t2, 1, "min");
677: if (bmax) GA_Dgop(&t3, 1, "max");
678:
679: *bmin1 = t1; /* boundmin */
680: *bmin2 = t2; /* wolfemin */
681: *bmax = t3; /* boundmax */
682:
683: TaoFunctionReturn(0);
684: }
690: int
691: TaoVecGa::BoundGradientProjection(TaoVec *GG, TaoVec *XXL, TaoVec *XX,
692: TaoVec *XXU) {
693:
694: TaoVecGa *GP = this;
695:
696: int i, n, n2, n3, n4, n5, info;
697: GAScalar *xptr,*xlptr,*xuptr,*gptr,*gpptr;
698: GAScalar xval,gpval;
699:
700: TaoFunctionBegin;
701: /* Project variables at the lower and upper bound */
703: info = XX->GetArray((TaoScalar **)&xptr, &n); CHKERRQ(info);
704: info = XXL->GetArray((TaoScalar **)&xlptr, &n2); CHKERRQ(info);
705: info = XXU->GetArray((TaoScalar **)&xuptr, &n3); CHKERRQ(info);
706: info = GG->GetArray((TaoScalar **)&gptr, &n4); CHKERRQ(info);
708: if ( (n != n2) || (n != n3) || (n != n4) ) {
709: SETERRQ(1,"Vectors must all be same size");
710: }
711:
713: if (GG!=GP){
714: info = GP->GetArray((TaoScalar **)&gpptr,&n5); CHKERRQ(info);
715: if (n != n5) {SETERRQ(1, "Vectors must all be same size");}
716: } else { gpptr=gptr; }
718: for (i=0; i<n; ++i){
719: gpval = gptr[i]; xval = xptr[i];
720:
721: if (gpval>0 && xval<=xlptr[i]){
722: gpval = 0;
723: } else if (gpval<0 && xval>=xuptr[i]){
724: gpval = 0;
725: }
726: gpptr[i] = gpval;
727: }
728:
729: info = XX->RestoreArray((TaoScalar **)&xptr, &n); CHKERRQ(info);
730: info = XXL->RestoreArray((TaoScalar **)&xlptr, &n2); CHKERRQ(info);
731: info = XXU->RestoreArray((TaoScalar **)&xuptr, &n3); CHKERRQ(info);
732: info = GG->RestoreArray((TaoScalar **)&gptr, &n4); CHKERRQ(info);
733: if (GG!=GP) {
734: info = GP->RestoreArray((TaoScalar **)&gpptr, &n5); CHKERRQ(info);
735: }
736: TaoFunctionReturn(0);
737: }