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: }