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