Actual source code: ebvec1.c
4: /*$Id: ebvec1.c,v 1.9 2001/09/26 17:10:29 balay Exp $*/
7: #include src/vec/vecimpl.h
8: #include src/vec/impls/dvecimpl.h
9: #include esi/petsc/vector.h
11: typedef struct {
12: ::esi::Vector<double,int> *evec;
13: } Vec_ESI;
15: /*
16: Wraps a PETSc vector to look like an ESI vector and stashes the wrapper inside the
17: PETSc vector. If PETSc vector already had wrapper uses that instead.
18: */
19: #undef __FUNCT__
21: int VecESIWrap(Vec xin,::esi::Vector<double,int> **v)
22: {
23: esi::petsc::Vector<double,int> *t;
24: int ierr;
27: if (!xin->esivec) {
28: t = new esi::petsc::Vector<double,int>(xin);
29: t->getInterface("esi::Vector",xin->esivec);
30: }
31: *v = reinterpret_cast<esi::Vector<double,int>* >(xin->esivec);
32: return(0);
33: }
35: #undef __FUNCT__
37: /*@C
38: VecESISetVector - Takes a PETSc vector sets it to type ESI and
39: provides the ESI vector that it wraps to look like a PETSc vector.
41: Level: intermediate
43: @*/
44: int VecESISetVector(Vec xin,::esi::Vector<double,int> *v)
45: {
46: Vec_ESI *x;
47: PetscTruth tesi;
48: int ierr;
51: PetscTypeCompare((PetscObject)xin,0,&tesi);
52: if (tesi) {
53: VecSetType(xin,VECESI);
54: }
55: PetscTypeCompare((PetscObject)xin,VECESI,&tesi);
56: if (tesi) {
57: int n,N;
58: ::esi::IndexSpace<int> *map;
60: v->getGlobalSize(N);
61: if (xin->N == -1) xin->N = N;
62: else if (xin->N != N) SETERRQ2(1,"Global size of Vec %d not equal size of esi::Vector %d",xin->N,N);
64: v->getIndexSpace(map);
65: map->getLocalSize(n);
66: if (xin->n == -1) xin->n = n;
67: else if (xin->n != n) SETERRQ2(1,"Local size of Vec %d not equal size of esi::Vector %d",xin->n,n);
69: x = (Vec_ESI*)xin->data;
70: x->evec = v;
71: v->addReference();
72: if (!xin->map){
73: PetscMapCreateMPI(xin->comm,n,N,&xin->map);
74: }
75: VecStashCreate_Private(xin->comm,1,&xin->stash);
76: VecStashCreate_Private(xin->comm,xin->bs,&xin->bstash);
77: (v)->getInterface("esi::Vector",xin->esivec);
78: }
79: return(0);
80: }
82: /* ---------------------------------------------------------------------------------*/
84: #undef __FUNCT__
86: int VecPlaceArray_ESI(Vec vin,const PetscScalar *a)
87: {
88: Vec_ESI *v = (Vec_ESI *)vin->data;
89: ::esi::VectorReplaceAccess<double,int> *vr;
90: int ierr;
93: v->evec->getInterface("esi::VectorReplaceAccess",reinterpret_cast<void *&>(vr));
94: vr->setArrayPointer((PetscScalar*)a,vin->n);
95: return(0);
96: }
98: #undef __FUNCT__
100: int VecSet_ESI(const PetscScalar *alpha,Vec xin)
101: {
102: Vec_ESI *x = (Vec_ESI*)xin->data;
103: int ierr;
106: x->evec->put(*alpha);
107: return(0);
108: }
110: #undef __FUNCT__
112: int VecDuplicate_ESI(Vec xin,Vec *xout)
113: {
114: Vec_ESI *x = (Vec_ESI*)xin->data;
115: int ierr;
116: ::esi::Vector<double,int> *nevec;
119: VecCreate(xin->comm,xout);
120: VecSetSizes(*xout,xin->n,xin->N);
121: x->evec->clone(nevec);
122: VecESISetVector(*xout,nevec);
123: nevec->deleteReference();
124: return(0);
125: }
127: #undef __FUNCT__
129: int VecDot_ESI(Vec xin,Vec yin,PetscScalar *z)
130: {
131: Vec_ESI *x = (Vec_ESI*)xin->data;
132: int ierr;
133: ::esi::Vector<double,int> *ytmp;
136: /* Make yin look like an esi:Vector */
137: VecESIWrap(yin,&ytmp);
138: x->evec->dot(*ytmp,*z);
139: return(0);
140: }
142: #undef __FUNCT__
144: int VecAXPY_ESI(const PetscScalar *a,Vec xin,Vec yin)
145: {
146: Vec_ESI *x = (Vec_ESI*)xin->data;
147: int ierr;
148: ::esi::Vector<double,int> *ytmp;
151: /* Make yin look like an esi:Vector */
152: VecESIWrap(yin,&ytmp);
153: ytmp->axpy(*x->evec,*a);
154: return(0);
155: }
157: #undef __FUNCT__
159: int VecAYPX_ESI(const PetscScalar *a,Vec xin,Vec yin)
160: {
161: Vec_ESI *x = (Vec_ESI*)xin->data;
162: int ierr;
163: ::esi::Vector<double,int> *ytmp;
166: /* Make yin look like an esi:Vector */
167: VecESIWrap(yin,&ytmp);
168: ytmp->aypx(*a,*x->evec);
169: return(0);
170: }
172: #undef __FUNCT__
174: int VecWAXPY_ESI(const PetscScalar *a,Vec xin,Vec yin,Vec win)
175: {
176: Vec_ESI *x = (Vec_ESI*)xin->data;
177: int ierr;
178: ::esi::Vector<double,int> *ytmp,*wtmp;
181: /* Make yin look like an esi:Vector */
182: VecESIWrap(yin,&ytmp);
183: VecESIWrap(win,&wtmp);
184: wtmp->axpby(*a,*x->evec,1.0,*ytmp);
185: return(0);
186: }
188: #undef __FUNCT__
190: int VecCopy_ESI(Vec xin,Vec yin)
191: {
192: Vec_ESI *x = (Vec_ESI*)xin->data;
193: int ierr;
194: ::esi::Vector<double,int> *ytmp;
197: if (xin != yin) {
198: /* Make yin look like an esi:Vector */
199: VecESIWrap(yin,&ytmp);
200: ytmp->copy(*x->evec);
201: }
202: return(0);
203: }
205: #undef __FUNCT__
207: int VecPointwiseMult_ESI(Vec xin,Vec yin,Vec zin)
208: {
209: Vec_ESI *x = (Vec_ESI*)xin->data;
210: int ierr;
211: ::esi::Vector<double,int> *ztmp;
214: if (zin != yin) {
215: VecCopy(yin,zin);
216: }
218: /* Make zin look like an esi:Vector */
219: VecESIWrap(zin,&ztmp);
220: ztmp->scaleDiagonal(*x->evec);
221: return(0);
222: }
224: #undef __FUNCT__
226: int VecPointwiseDivide_ESI(Vec xin,Vec yin,Vec win)
227: {
228: int n = xin->n,i,ierr;
229: PetscScalar *xx,*yy,*ww;
232: VecGetArrayFast(yin,&yy);
233: if (yin != xin) {VecGetArrayFast(xin,&xx);}
234: else xx = yy;
235: if (yin != win) {VecGetArrayFast(win,&ww);}
236: else ww = yy;
237: for (i=0; i<n; i++) ww[i] = xx[i] / yy[i];
238: VecRestoreArrayFast(yin,&yy);
239: if (yin != win) {VecRestoreArrayFast(win,&ww);}
240: if (xin != win) {VecRestoreArrayFast(xin,&xx);}
241: return(0);
242: }
244: #include petscblaslapack.h
245: /*
246: ESI does not provide a method for this
247: */
248: #undef __FUNCT__
250: int VecSwap_ESI(Vec xin,Vec yin)
251: {
255: if (xin != yin) {
256: PetscScalar *ya,*xa;
257: int one = 1;
259: VecGetArrayFast(yin,&ya);
260: VecGetArrayFast(xin,&xa);
261: BLswap_(&xin->n,xa,&one,ya,&one);
262: VecRestoreArrayFast(xin,&xa);
263: VecRestoreArrayFast(yin,&ya);
264: }
265: return(0);
266: }
268: #undef __FUNCT__
270: int VecMDot_ESI(int nv,Vec xin,const Vec yin[],PetscScalar *z)
271: {
272: Vec_ESI *x = (Vec_ESI *)xin->data;
273: int ierr,i;
274: ::esi::Vector<double,int> *ytmp;
277: for (i=0; i<nv; i++) {
278: /* Make yin look like an esi:Vector */
279: VecESIWrap(yin[i],&ytmp);
280: x->evec->dot(*ytmp,z[i]);
281: }
282: return(0);
283: }
286: #undef __FUNCT__
288: int VecMAXPY_ESI(int nv,const PetscScalar *a,Vec xin, Vec yin[])
289: {
290: Vec_ESI *x = (Vec_ESI *)xin->data;
291: int ierr,i;
292: ::esi::Vector<double,int> *ytmp;
295: for (i=0; i<nv; i++) {
296: /* Make yin look like an esi:Vector */
297: VecESIWrap(yin[i],&ytmp);
298: x->evec->axpy(*ytmp,a[i]);
299: }
300: return(0);
301: }
304: #undef __FUNCT__
306: int VecGetSize_ESI(Vec vin,int *size)
307: {
308: Vec_ESI *x = (Vec_ESI*)vin->data;
309: int ierr;
312: x->evec->getGlobalSize(*size);
313: return(0);
314: }
316: #undef __FUNCT__
318: int VecGetLocalSize_ESI(Vec vin,int *size)
319: {
320: Vec_ESI *x = (Vec_ESI*)vin->data;
321: int ierr;
322: ::esi::IndexSpace<int> *map;
325: x->evec->getIndexSpace(map);
326: map->getLocalSize(*size);
327: return(0);
328: }
330: #undef __FUNCT__
332: int VecGetArray_ESI(Vec vin,PetscScalar **array)
333: {
334: Vec_ESI *x = (Vec_ESI*)vin->data;
335: int ierr;
338: x->evec->getCoefPtrReadWriteLock(*array);
339: return(0);
340: }
342: #undef __FUNCT__
344: int VecRestoreArray_ESI(Vec vin,PetscScalar **array)
345: {
346: Vec_ESI *x = (Vec_ESI*)vin->data;
347: int ierr;
350: x->evec->releaseCoefPtrLock(*array);
351: return(0);
352: }
354: #undef __FUNCT__
356: int VecScale_ESI(const PetscScalar *array,Vec vin)
357: {
358: Vec_ESI *x = (Vec_ESI*)vin->data;
359: int ierr;
362: x->evec->scale(*array);
363: return(0);
364: }
366: #undef __FUNCT__
368: int VecNorm_ESI(Vec vin,NormType ntype,PetscReal *norm)
369: {
370: Vec_ESI *x = (Vec_ESI*)vin->data;
371: int ierr;
374: if (ntype == NORM_2) {
375: x->evec->norm2(*norm);
376: } else if (ntype == NORM_1) {
377: x->evec->norm1(*norm);
378: } else if (ntype == NORM_INFINITY) {
379: x->evec->normInfinity(*norm);
380: } else SETERRQ1(1,"Unknown NormType %d",ntype);
381: return(0);
382: }
384: extern int VecSetValues_MPI(Vec,int,const int[],const PetscScalar[],InsertMode);
385: extern int VecAssemblyBegin_MPI(Vec);
386: extern int VecAssemblyEnd_MPI(Vec);
387: extern int VecView_MPI(Vec,PetscViewer);
388: extern int VecReciprocal_Default(Vec);
389: extern int VecSetRandom_Seq(PetscRandom,Vec);
390: extern int VecSetValuesBlocked_MPI(Vec,int,const int[],const PetscScalar[],InsertMode);
391: extern int VecMax_MPI(Vec,int*,PetscReal*);
392: extern int VecMin_MPI(Vec,int*,PetscReal*);
394: /* ---------------------------------------------------------------------------------*/
396: #undef __FUNCT__
398: int VecDestroy_ESI(Vec v)
399: {
400: Vec_ESI *vs = (Vec_ESI*)v->data;
401: int ierr;
404: if (vs->evec) {
405: vs->evec->deleteReference();
406: }
407: VecStashDestroy_Private(&v->bstash);
408: VecStashDestroy_Private(&v->stash);
409: PetscFree(vs);
410: return(0);
411: }
413: EXTERN_C_BEGIN
414: #undef __FUNCT__
416: int VecCreate_PetscESI(Vec V)
417: {
418: int ierr;
419: Vec v;
420: esi::petsc::Vector<double,int> *ve;
423: V->ops->destroy = 0; /* since this is called from VecSetType() we have to make sure it doesn't get destroyed twice */
424: VecSetType(V,VECESI);
425: VecCreate(V->comm,&v);
426: VecSetSizes(v,V->n,V->N);
427: if (V->bs > 1) {VecSetBlockSize(v,V->bs);}
428: PetscObjectSetOptionsPrefix((PetscObject)v,"esi_");
429: VecSetFromOptions(v);
430: ve = new esi::petsc::Vector<double,int>(v);
431: VecESISetVector(V,ve);
432: ve->deleteReference();
433: PetscObjectDereference((PetscObject)v);
434: return(0);
435: }
436: EXTERN_C_END
438: static struct _VecOps EvOps = {VecDuplicate_ESI,
439: VecDuplicateVecs_Default,
440: VecDestroyVecs_Default,
441: VecDot_ESI,
442: VecMDot_ESI,
443: VecNorm_ESI,
444: 0,
445: 0,
446: VecScale_ESI,
447: VecCopy_ESI,
448: VecSet_ESI,
449: VecSwap_ESI,
450: VecAXPY_ESI,
451: 0,
452: VecMAXPY_ESI,
453: VecAYPX_ESI,
454: VecWAXPY_ESI,
455: VecPointwiseMult_ESI,
456: VecPointwiseDivide_ESI,
457: VecSetValues_MPI,
458: VecAssemblyBegin_MPI,
459: VecAssemblyEnd_MPI,
460: VecGetArray_ESI,
461: VecGetSize_ESI,
462: VecGetLocalSize_ESI,
463: VecRestoreArray_ESI,
464: VecMax_MPI,
465: VecMin_MPI,
466: VecSetRandom_Seq,
467: 0,
468: VecSetValuesBlocked_MPI,
469: VecDestroy_ESI,
470: VecView_MPI,
471: VecPlaceArray_ESI,
472: 0,
473: VecDot_Seq,
474: VecTDot_Seq,
475: VecNorm_Seq,
476: 0,
477: VecReciprocal_Default};
479: #undef __FUNCT__
481: int VecESISetFromOptions(Vec V)
482: {
483: char string[1024];
484: PetscTruth flg;
485: int ierr;
486:
488: PetscTypeCompare((PetscObject)V,VECESI,&flg);
489: if (flg) {
490: PetscOptionsGetString(V->prefix,"-vec_esi_type",string,1024,&flg);
491: if (flg) {
492: VecESISetType(V,string);
493: }
494: }
495: return(0);
496: }
499: EXTERN_C_BEGIN
500: #undef __FUNCT__
502: int VecCreate_ESI(Vec V)
503: {
504: Vec_ESI *s;
505: int ierr;
506:
508: ierr = PetscNew(Vec_ESI,&s);
509: ierr = PetscMemzero(s,sizeof(Vec_ESI));
511: s->evec = 0;
512: V->data = (void*)s;
513: V->petscnative = PETSC_FALSE;
514: V->esivec = 0;
515: ierr = PetscMemcpy(V->ops,&EvOps,sizeof(EvOps));
516: return(0);
517: }
518: EXTERN_C_END
520: extern PetscFList CCAList;
522: #undef __FUNCT__
524: /*@C
525: ESICreateIndexSpace - Creates an esi::IndexSpace using the -is_esi_type type
527: Level: beginner
528:
529: @*/
530: int ESICreateIndexSpace(const char * commname,void *comm,int m,::esi::IndexSpace<int>*&v)
531: {
532: int ierr;
533: ::esi::IndexSpace<int>::Factory *f;
534: ::esi::IndexSpace<int>::Factory *(*r)(void);
535: char name[PETSC_MAX_PATH_LEN];
536: PetscTruth found;
539: PetscOptionsGetString(PETSC_NULL,"-is_esi_type",name,1024,&found);
540: if (!found) {
541: PetscStrcpy(name,"esi::petsc::IndexSpace");
542: }
543: PetscFListFind(*(MPI_Comm*)comm,CCAList,name,(void(**)(void))&r);
544: if (!r) SETERRQ1(1,"Unable to load esi::IndexSpace Factory constructor %s",name);
545: f = (*r)();
546: f->create(commname,comm,m,PETSC_DECIDE,PETSC_DECIDE,v);
547: delete f;
548: return(0);
549: }
551: #undef __FUNCT__
553: /*@C
554: VecESISetType - Given a PETSc vector of type ESI loads the ESI constructor
555: by name and wraps the ESI vector to look like a PETSc vector.
557: Level: intermediate
558: @*/
559: int VecESISetType(Vec V,char *name)
560: {
561: int ierr;
562: ::esi::Vector<double,int> *ve;
563: ::esi::Vector<double,int>::Factory *f,*(*r)(void);
564: ::esi::IndexSpace<int> *map;
567: PetscFListFind(V->comm,CCAList,name,(void(**)(void))&r);
568: if (!r) SETERRQ1(1,"Unable to load esi::VectorFactory constructor %s",name);
569: f = (*r)();
570: if (V->n == PETSC_DECIDE) {
571: PetscSplitOwnership(V->comm,&V->n,&V->N);
572: }
573: ESICreateIndexSpace("MPI",&V->comm,V->n,map);
574: f->create(*map,ve);
575: map->deleteReference();
576: delete f;
577: VecESISetVector(V,ve);
578: ve->deleteReference();
579: return(0);
580: }
582: #undef __FUNCT__
584: /*@C
585: ESILoadFactory - Creates an object for any ESI factory class. Generally does
586: this by dynamicly loading the constructor.
588: Collective on MPI_Comm
590: Input Parameters:
591: + commname - name of parallel context; currently only "MPI" is supported
592: . comm - communicator for parallel computing model, currently only MPI_Comm's are supported
593: - name - name of class the factory constructs, for example "esi::petsc::Vector"
595: Output Parameter:
596: . f - the factory object
598: Notes: The name of the class is the name of the class the factory CONSTRUCTS, not the name
599: of the factory class
601: Level: intermediate
602: @*/
603: int ESILoadFactory(char *commname,void *comm,char *classname,void *&f)
604: {
605: int ierr;
606: void *(*r)();
607: PetscTruth flag;
610: PetscStrcmp(commname,"MPI",&flag);
611: if (!flag) SETERRQ1(1,"Parallel computing model %s not supported",commname);
612: PetscFListFind(*(MPI_Comm*)comm,CCAList,classname,(void(**)(void))&r);
613: if (!r) SETERRQ1(1,"Unable to load constructor %s",classname);
614: f = (*r)();
615: return(0);
616: }