Actual source code: taomat_ga.c

      1: //File: taomat_ga.c

      3: /**************************************************************

      5: Author: Limin Zhang, Ph.D.
      6:         Mathematics Department
      7:         Columbia Basin College
      8:         Pasco, WA 99301
      9:         Limin.Zhang@cbc2.org

     11: Mentor: Jarek Naplocha, Ph.D.
     12:         Environmental Molecular Science Laboratory
     13:         Pacific Northwest National Laboratory
     14:         Richland, WA 99352

     16: Date: 7/11/2002

     18: Purpose:
     19:       to design and implement some interfaces between TAO and
     20:       global arrays.

     22: Revise History:

     24: 8/8/02
     25:         To clean Petsc function calls and marcos.
     26: **************************************************************/

     29: //#include "tao_general.h"

     31: #include "taomat_ga.h"
     32: #include "taovec_ga.h"

     34: /*
     35: int MatD_Fischer(Mat, Vec, Vec, Vec, Vec, Vec, Vec, Vec, Vec);
     36: int MatCreateADA(Mat,Vec,Vec,Mat*);
     37: int MatSMFResetRowColumn(Mat,IS,IS);
     38: */

     40: /*@C
     41:    TaoWrapGaMat - Creates a new TaoMat object using a GA matrix.

     43:    Input Parameter:
     44: .  M -  a GA matrix

     46:    Output Parameter:
     47: .  MM - new TaoMat

     49:    Level: advanced

     51: .seealso TaoMatGetGaMat(), TaoMatDestroy()
     52: @*/
     53: int TaoWrapGaMat (GAMat M, TaoMatGa ** MM)
     54: {
     55:   TaoFunctionBegin;
     56:   *MM = new TaoMatGa (M);
     57:   TaoFunctionReturn (0);
     58: }

     60: /*@C
     61:    TaoMatGetGaMat - If the TaoMat is of the TaoMatGa type, this routine
     62:    gets the underlying GA matrix.

     64:    Input Parameter:
     65: .  MM - the TaoMat

     67:    Output Parameter:
     68: .  M -  the GA mat

     70:    Note:
     71:    The function TaoMatGa::GetMat() will also return the Mat M.

     73:    Level: advanced
     74: @*/
     75: int TaoMatGetGaMat (TaoMat * MM, GAMat * M)
     76: {
     77:   TaoFunctionBegin;
     78:   if (MM)
     79:     {
     80:       *M = ((TaoMatGa *) MM)->GetMat ();
     81:     }
     82:   else
     83:     {
     84:       *M = 0;
     85:     }
     86:   TaoFunctionReturn (0);
     87: }

     89: TaoMatGa::TaoMatGa (GAMat MM):TaoMat ()
     90: {
     91:   this->pm = MM;
     92:   this->MatObject = (void *) (&MM);
     93:   this->pm_pre = MM;
     94:   return;
     95: }

     97: int TaoMatGa::Compatible (TaoVec * xx, TaoVec * yy, TaoTruth *flag)
     98: {
     99:   int
    100:     nx,
    101:     ny,
    102:     m,
    103:     n;
    104:   GAVec
    105:     x,
    106:     y;
    107:   if (xx == 0 || yy == 0)
    108:     { *flag=TAO_FALSE; return 0; }

    110:   x = ((TaoVecGa *) xx)->GetVec ();
    111:   y = ((TaoVecGa *) yy)->GetVec ();
    112:   if (x == 0 || y == 0 || this->pm == 0)
    113:     { *flag=TAO_FALSE; return 0; }

    115:   int
    116:     xtype,
    117:     ytype,
    118:     type,
    119:     ndim,
    120:     dims[2];

    122:   NGA_Inquire (x, &xtype, &ndim, &nx);
    123:   if (ndim != 1)
    124:     GA_Error ("TaoMatGa::Compatible: Wrong dimension", ndim);

    126:   NGA_Inquire (y, &ytype, &ndim, &ny);
    127:   if (ndim != 1)
    128:     GA_Error ("TaoMatGa::Compatible:Wrong dimension", ndim);

    130:   NGA_Inquire (this->pm, &type, &ndim, dims);
    131:   if (ndim != 2)
    132:     GA_Error ("TaoMatGa::Compatible:Wrong dimension", ndim);
    133:   n = dims[0];
    134:   m = dims[1];

    136:   if (xtype != type || ytype != type)
    137:     GA_Error ("TaoMatGa::Compatible:Wrong data type", type);

    139:   if (n != nx || m != ny)
    140:     { *flag=TAO_FALSE; return 0; }

    142:   *flag=TAO_TRUE;
    143:   return 0;

    145: }

    147: int
    148: TaoMatGa::Clone (TaoMat ** ntm)
    149: {
    150:   TaoMatGa *nptm;
    151:   GAMat npm;
    152:   TaoFunctionBegin;
    153:   npm = GA_Duplicate (this->pm, "Cloned");
    154:   if (!npm)
    155:     GA_Error ("TaoMatGa::Clone:duplicate failed: ", npm);
    156:   TaoWrapGaMat (npm, &nptm);
    157:   *ntm = nptm;
    158:   TaoFunctionReturn (0);
    159: }

    161: int
    162: TaoMatGa::CopyFrom (TaoMat * tm)
    163: {
    164:   TaoFunctionBegin;
    165: #if 0 //try a different way to get GAMat
    166:   TaoMatGa *temp2 = dynamic_cast < TaoMatGa * >(tm);
    167:   GAMat pm2 = temp2->GetMat ();
    168: #endif
    169:   GAMat pm2 = ((TaoMatGa *)tm)->GetMat ();
    170:   GA_Copy (pm2, this->pm);
    171:   TaoFunctionReturn (0);
    172: }

    174: int
    175: TaoMatGa::GetDimensions (int *m, int *n)
    176: {
    177:   TaoFunctionBegin;
    178:   int type, ndim, dims[2];
    179:   NGA_Inquire (this->pm, &type, &ndim, dims);
    180:   if (ndim != 2)
    181:     GA_Error ("TaoMatGa::GetDimension: wrong dimension", ndim);
    182:   else
    183:     {
    184:       *n = dims[0];
    185:       *m = dims[1];
    186:     }
    187:   TaoFunctionReturn (0);
    188: }

    190: int
    191: TaoMatGa::Multiply (TaoVec * tv, TaoVec * tw)
    192: {
    193:   double alpha = 1.0;
    194:   double beta = 0.0;
    195:   TaoTruth flag;
    196:   TaoFunctionBegin;
    197:   GAVec vv = ((TaoVecGa *) tv)->GetVec ();
    198:   GAVec ww = ((TaoVecGa *) tw)->GetVec ();
    199:   int m, n = 1, k;                //n = 1 due to vecotor ww that is k by 1 array
    200:   this->GetDimensions (&m, &k);
    201:   this->Compatible (tv, tw,&flag);
    202:   if (flag)
    203:     GA_Dgemm ('N', 'N', m, n, k, alpha, this->pm, vv, beta, ww);
    204:   else
    205:     TaoFunctionReturn (1);
    206:   TaoFunctionReturn (0);
    207: }

    209: int
    210: TaoMatGa::MultiplyAdd (TaoVec * tv, TaoVec * tw, TaoVec * ty)
    211: {
    212:   TaoTruth flag1,flag2;
    213:   TaoFunctionBegin;
    214:   GAVec xx = ((TaoVecGa *) tv)->GetVec ();
    215:   GAVec ww = ((TaoVecGa *) tw)->GetVec ();
    216:   GAVec yy = ((TaoVecGa *) ty)->GetVec ();
    217:   GAVec temp_ww = GA_Duplicate (ww, "temp_WW");
    218:   if (!temp_ww)
    219:     GA_Error ("TaoMatGa::MultiplyAdd:Duplicate:", temp_ww);
    220:   GA_Copy (ww, temp_ww);
    221:   double alpha = 1.0;
    222:   double beta = 1.0;
    223:   int m, n = 1, k;                //n = 1 due to vecotor ww that is k by 1 array
    224:   this->GetDimensions (&m, &k);
    225:   this->Compatible (tv, ty, &flag1);
    226:   this->Compatible (tv, tw, &flag2);
    227:   if (flag1 && flag2)
    228:     {
    229:       GA_Dgemm ('N', 'N', m, n, k, alpha, this->pm, xx, beta, temp_ww);
    230:       //copy from temp_ww to yy. Doing this way, the data in ww is not dirty!!!!!
    231:       GA_Copy (temp_ww, yy);
    232:       //free the temporary memory
    233:       GA_Destroy(temp_ww);
    234:     }
    235:   else
    236:     TaoFunctionReturn (1);
    237:   TaoFunctionReturn (0);
    238: }

    240: int
    241: TaoMatGa::MultiplyTranspose (TaoVec * tv, TaoVec * tw)
    242: {
    243:   double alpha = 1.0;
    244:   double beta = 0.0;
    245:   int m, n = 1, k;                /*n = 1 due to vecotor ww that is k */
    246:   TaoFunctionBegin;
    247:   GAVec vv = ((TaoVecGa *) tv)->GetVec ();
    248:   GAVec ww = ((TaoVecGa *) tw)->GetVec ();
    249:   this->GetDimensions (&m, &k);
    250:   GA_Dgemm ('T', 'N', k, n, m, alpha, this->pm, vv, beta, ww);
    251:   TaoFunctionReturn (0);
    252: }

    254: int
    255: TaoMatGa::MultiplyTransposeAdd (TaoVec * tv, TaoVec * tw, TaoVec * ty)
    256: {
    257:   GAVec xx = ((TaoVecGa *) tv)->GetVec ();
    258:   GAVec ww = ((TaoVecGa *) tw)->GetVec ();
    259:   GAVec yy = ((TaoVecGa *) ty)->GetVec ();
    260:   TaoFunctionBegin;
    261:   TaoFunctionBegin;
    262:   GAVec temp_ww = GA_Duplicate (ww, "temp_WW");
    263:   if (!temp_ww)
    264:     GA_Error ("TaoMatGa::MultiplyAdd:Dupli cate:", temp_ww);
    265:   GA_Copy (ww, temp_ww);
    266:   double alpha = 1.0;
    267:   double beta = 1.0;
    268:   int m, n = 1, k;
    269:   //n = 1 due to vecotor ww that is k by 1 array
    270:   this->GetDimensions (&m, &k);
    271:   GA_Dgemm ('T', 'N', k, n, m, alpha, this->pm, xx, beta, temp_ww);
    272:   //copy from temp_ww to yy. Doing this way, the data in ww is not dirty!!!!!
    273:   GA_Copy (temp_ww, yy);
    274:   //free the memory held by the temporary GA array temp_ww
    275:   GA_Destroy(temp_ww);
    276:   TaoFunctionReturn (0);
    277: }

    279: int
    280: TaoMatGa::SetDiagonal (TaoVec * tv)
    281: {
    282:   GAVec vv = ((TaoVecGa *) tv)->GetVec ();
    283:   TaoFunctionBegin;
    284:   GA_Set_diagonal (this->pm, vv);
    285:   TaoFunctionReturn (0);
    286: }

    288: int
    289: TaoMatGa::AddDiagonal (TaoVec * tv)
    290: {
    291:   GAVec vv = ((TaoVecGa *) tv)->GetVec ();
    292:   TaoFunctionBegin;
    293:   GA_Add_diagonal (this->pm, vv);
    294:   TaoFunctionReturn (0);
    295: }

    297: int
    298: TaoMatGa::GetDiagonal (TaoVec * tv)
    299: {
    300:   GAVec vv = ((TaoVecGa *) tv)->GetVec ();
    301:   TaoFunctionBegin;
    302:   GA_Get_diag (this->pm, vv);
    303:   TaoFunctionReturn (0);
    304: }

    306: int
    307: TaoMatGa::ShiftDiagonal (TaoScalar c)
    308: {
    309:   TaoFunctionBegin;
    310:   GA_Shift_diagonal (this->pm, &c);
    311:   TaoFunctionReturn (0);
    312: }

    314: int
    315: TaoMatGa::View ()
    316: {
    317:   TaoFunctionBegin;
    318:   GA_Print (this->pm);
    319:   TaoFunctionReturn (0);
    320: }

    323: int
    324: TaoMatGa::Presolve ()
    325: {
    326:   TaoFunctionBegin;
    327:   TaoFunctionReturn (0);
    328: }

    330: int
    331: TaoMatGa::Solve (TaoVec * tv, TaoVec * tw, TaoTruth * tt)
    332: {
    333:   TaoFunctionBegin;
    334:   GA_Error ("not yet implemented.", 0);
    335:   TaoFunctionReturn (1);
    336: }

    338: int
    339: TaoMatGa::RowScale (TaoVec * tv)
    340: {
    341:   GAVec vv = ((TaoVecGa *) tv)->GetVec ();
    342:   TaoFunctionBegin;
    343:   GA_Scale_rows (this->pm, vv);
    344:   TaoFunctionReturn (0);
    345: }

    347: int
    348: TaoMatGa::ColScale (TaoVec * tv)
    349: {
    350:   GAVec vv = ((TaoVecGa *) tv)->GetVec ();
    351:   TaoFunctionBegin;
    352:   GA_Scale_cols (this->pm, vv);
    353:   TaoFunctionReturn (0);
    354: }

    356: int
    357: TaoMatGa::CreateReducedMatrix (TaoIndexSet * S1,
    358:                                TaoIndexSet * S2, TaoMat ** MM)
    359: {
    360:   int info;
    361:   TaoMatGa *M;
    362:   GAMat B = 0;
    363:   TaoFunctionBegin;
    364:   info = TaoWrapGaMat (B, &M);
    365:   CHKERRQ (info);
    366:   *MM = M;
    367:   TaoFunctionReturn (0);
    368: }

    370: int
    371: TaoMatGa::SetReducedMatrix (TaoMat * M, TaoIndexSet * S1, TaoIndexSet * S2)
    372: {
    373:   TaoFunctionBegin;
    374:   GA_Error ("Not yet implemented.", 0);
    375:   TaoFunctionReturn (1);
    376: }

    380: int
    381: TaoMatGa::D_Fischer (TaoVec * tx, TaoVec * tf, TaoVec * tl,
    382:                      TaoVec * tu, TaoVec * tda,
    383:                      TaoVec * tdb, TaoVec * tt1, TaoVec * tt2)
    384: {
    385:   TaoFunctionBegin;
    386:   GA_Error ("Not yet implemented.", 0);
    387:   TaoFunctionReturn (1);
    388: }

    390: int
    391: TaoMatGa::Norm1 (double *nm)
    392: {
    393:   TaoFunctionBegin;
    394:   GA_Norm1 (this->pm, nm);
    395:   TaoFunctionReturn (0);
    396: }

    399: int
    400: TaoMatGa::Fill (double c)
    401: {

    403:   TaoFunctionBegin;
    404:   GA_Fill (this->pm, &c);
    405:   TaoFunctionReturn (0);
    406: }

    409: int
    410: TaoMatGa::Zero ()
    411: {

    413:   TaoFunctionBegin;
    414:   GA_Zero (this->pm);
    415:   TaoFunctionReturn (0);

    417: }