Actual source code: veccomp.c
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
11: #include <slepc/private/vecimplslepc.h>
13: /* Private MPI datatypes and operators */
14: static MPI_Datatype MPIU_NORM2=0, MPIU_NORM1_AND_2=0;
15: static PetscBool VecCompInitialized = PETSC_FALSE;
16: MPI_Op MPIU_NORM2_SUM=0;
18: /* Private functions */
19: static inline void SumNorm2(PetscReal*,PetscReal*,PetscReal*,PetscReal*);
20: static inline PetscReal GetNorm2(PetscReal,PetscReal);
21: static inline void AddNorm2(PetscReal*,PetscReal*,PetscReal);
22: static PetscErrorCode VecCompSetSubVecs_Comp(Vec,PetscInt,Vec*);
23: static PetscErrorCode VecCompGetSubVecs_Comp(Vec,PetscInt*,const Vec**);
25: #include "veccomp0.h"
28: #include "veccomp0.h"
30: static inline void SumNorm2(PetscReal *ssq0,PetscReal *scale0,PetscReal *ssq1,PetscReal *scale1)
31: {
32: PetscReal q;
33: if (*scale0 > *scale1) {
34: q = *scale1/(*scale0);
35: *ssq1 = *ssq0 + q*q*(*ssq1);
36: *scale1 = *scale0;
37: } else {
38: q = *scale0/(*scale1);
39: *ssq1 += q*q*(*ssq0);
40: }
41: }
43: static inline PetscReal GetNorm2(PetscReal ssq,PetscReal scale)
44: {
45: return scale*PetscSqrtReal(ssq);
46: }
48: static inline void AddNorm2(PetscReal *ssq,PetscReal *scale,PetscReal x)
49: {
50: PetscReal absx,q;
51: if (x != 0.0) {
52: absx = PetscAbs(x);
53: if (*scale < absx) {
54: q = *scale/absx;
55: *ssq = 1.0 + *ssq*q*q;
56: *scale = absx;
57: } else {
58: q = absx/(*scale);
59: *ssq += q*q;
60: }
61: }
62: }
64: SLEPC_EXTERN void MPIAPI SlepcSumNorm2_Local(void *in,void *out,PetscMPIInt *cnt,MPI_Datatype *datatype)
65: {
66: PetscInt i,count = *cnt;
67: PetscReal *xin = (PetscReal*)in,*xout = (PetscReal*)out;
69: PetscFunctionBegin;
70: if (*datatype == MPIU_NORM2) {
71: for (i=0;i<count;i++) {
72: SumNorm2(&xin[i*2],&xin[i*2+1],&xout[i*2],&xout[i*2+1]);
73: }
74: } else if (*datatype == MPIU_NORM1_AND_2) {
75: for (i=0;i<count;i++) {
76: xout[i*3] += xin[i*3];
77: SumNorm2(&xin[i*3+1],&xin[i*3+2],&xout[i*3+1],&xout[i*3+2]);
78: }
79: } else {
80: (void)(*PetscErrorPrintf)("Can only handle MPIU_NORM* data types");
81: MPI_Abort(PETSC_COMM_WORLD,1);
82: }
83: PetscFunctionReturnVoid();
84: }
86: static PetscErrorCode VecCompNormEnd(void)
87: {
88: PetscFunctionBegin;
89: PetscCallMPI(MPI_Type_free(&MPIU_NORM2));
90: PetscCallMPI(MPI_Type_free(&MPIU_NORM1_AND_2));
91: PetscCallMPI(MPI_Op_free(&MPIU_NORM2_SUM));
92: VecCompInitialized = PETSC_FALSE;
93: PetscFunctionReturn(PETSC_SUCCESS);
94: }
96: static PetscErrorCode VecCompNormInit(void)
97: {
98: PetscFunctionBegin;
99: PetscCallMPI(MPI_Type_contiguous(2,MPIU_REAL,&MPIU_NORM2));
100: PetscCallMPI(MPI_Type_commit(&MPIU_NORM2));
101: PetscCallMPI(MPI_Type_contiguous(3,MPIU_REAL,&MPIU_NORM1_AND_2));
102: PetscCallMPI(MPI_Type_commit(&MPIU_NORM1_AND_2));
103: PetscCallMPI(MPI_Op_create(SlepcSumNorm2_Local,PETSC_TRUE,&MPIU_NORM2_SUM));
104: PetscCall(PetscRegisterFinalize(VecCompNormEnd));
105: PetscFunctionReturn(PETSC_SUCCESS);
106: }
108: PetscErrorCode VecDestroy_Comp(Vec v)
109: {
110: Vec_Comp *vs = (Vec_Comp*)v->data;
111: PetscInt i;
113: PetscFunctionBegin;
114: #if defined(PETSC_USE_LOG)
115: PetscCall(PetscLogObjectState((PetscObject)v,"Length=%" PetscInt_FMT,v->map->n));
116: #endif
117: for (i=0;i<vs->nx;i++) PetscCall(VecDestroy(&vs->x[i]));
118: if (--vs->n->friends <= 0) PetscCall(PetscFree(vs->n));
119: PetscCall(PetscFree(vs->x));
120: PetscCall(PetscFree(vs));
121: PetscCall(PetscObjectComposeFunction((PetscObject)v,"VecCompSetSubVecs_C",NULL));
122: PetscCall(PetscObjectComposeFunction((PetscObject)v,"VecCompGetSubVecs_C",NULL));
123: PetscFunctionReturn(PETSC_SUCCESS);
124: }
126: static struct _VecOps DvOps = {
127: PetscDesignatedInitializer(duplicate,VecDuplicate_Comp),
128: PetscDesignatedInitializer(duplicatevecs,VecDuplicateVecs_Comp),
129: PetscDesignatedInitializer(destroyvecs,VecDestroyVecs_Comp),
130: PetscDesignatedInitializer(dot,VecDot_Comp_MPI),
131: PetscDesignatedInitializer(mdot,VecMDot_Comp_MPI),
132: PetscDesignatedInitializer(norm,VecNorm_Comp_MPI),
133: PetscDesignatedInitializer(tdot,VecTDot_Comp_MPI),
134: PetscDesignatedInitializer(mtdot,VecMTDot_Comp_MPI),
135: PetscDesignatedInitializer(scale,VecScale_Comp),
136: PetscDesignatedInitializer(copy,VecCopy_Comp),
137: PetscDesignatedInitializer(set,VecSet_Comp),
138: PetscDesignatedInitializer(swap,VecSwap_Comp),
139: PetscDesignatedInitializer(axpy,VecAXPY_Comp),
140: PetscDesignatedInitializer(axpby,VecAXPBY_Comp),
141: PetscDesignatedInitializer(maxpy,VecMAXPY_Comp),
142: PetscDesignatedInitializer(aypx,VecAYPX_Comp),
143: PetscDesignatedInitializer(waxpy,VecWAXPY_Comp),
144: PetscDesignatedInitializer(axpbypcz,VecAXPBYPCZ_Comp),
145: PetscDesignatedInitializer(pointwisemult,VecPointwiseMult_Comp),
146: PetscDesignatedInitializer(pointwisedivide,VecPointwiseDivide_Comp),
147: PetscDesignatedInitializer(setvalues,NULL),
148: PetscDesignatedInitializer(assemblybegin,NULL),
149: PetscDesignatedInitializer(assemblyend,NULL),
150: PetscDesignatedInitializer(getarray,NULL),
151: PetscDesignatedInitializer(getsize,VecGetSize_Comp),
152: PetscDesignatedInitializer(getlocalsize,VecGetLocalSize_Comp),
153: PetscDesignatedInitializer(restorearray,NULL),
154: PetscDesignatedInitializer(max,VecMax_Comp),
155: PetscDesignatedInitializer(min,VecMin_Comp),
156: PetscDesignatedInitializer(setrandom,VecSetRandom_Comp),
157: PetscDesignatedInitializer(setoption,NULL),
158: PetscDesignatedInitializer(setvaluesblocked,NULL),
159: PetscDesignatedInitializer(destroy,VecDestroy_Comp),
160: PetscDesignatedInitializer(view,VecView_Comp),
161: PetscDesignatedInitializer(placearray,NULL),
162: PetscDesignatedInitializer(replacearray,NULL),
163: PetscDesignatedInitializer(dot_local,VecDot_Comp_Seq),
164: PetscDesignatedInitializer(tdot_local,VecTDot_Comp_Seq),
165: PetscDesignatedInitializer(norm_local,VecNorm_Comp_Seq),
166: PetscDesignatedInitializer(mdot_local,VecMDot_Comp_Seq),
167: PetscDesignatedInitializer(mtdot_local,VecMTDot_Comp_Seq),
168: PetscDesignatedInitializer(load,NULL),
169: PetscDesignatedInitializer(reciprocal,VecReciprocal_Comp),
170: PetscDesignatedInitializer(conjugate,VecConjugate_Comp),
171: PetscDesignatedInitializer(setlocaltoglobalmapping,NULL),
172: PetscDesignatedInitializer(getlocaltoglobalmapping,NULL),
173: PetscDesignatedInitializer(resetarray,NULL),
174: PetscDesignatedInitializer(setfromoptions,NULL),
175: PetscDesignatedInitializer(maxpointwisedivide,VecMaxPointwiseDivide_Comp),
176: PetscDesignatedInitializer(pointwisemax,VecPointwiseMax_Comp),
177: PetscDesignatedInitializer(pointwisemaxabs,VecPointwiseMaxAbs_Comp),
178: PetscDesignatedInitializer(pointwisemin,VecPointwiseMin_Comp),
179: PetscDesignatedInitializer(getvalues,NULL),
180: PetscDesignatedInitializer(sqrt,VecSqrtAbs_Comp),
181: PetscDesignatedInitializer(abs,VecAbs_Comp),
182: PetscDesignatedInitializer(exp,VecExp_Comp),
183: PetscDesignatedInitializer(log,VecLog_Comp),
184: PetscDesignatedInitializer(shift,VecShift_Comp),
185: PetscDesignatedInitializer(create,NULL),
186: PetscDesignatedInitializer(stridegather,NULL),
187: PetscDesignatedInitializer(stridescatter,NULL),
188: PetscDesignatedInitializer(dotnorm2,VecDotNorm2_Comp_MPI),
189: PetscDesignatedInitializer(getsubvector,NULL),
190: PetscDesignatedInitializer(restoresubvector,NULL),
191: PetscDesignatedInitializer(getarrayread,NULL),
192: PetscDesignatedInitializer(restorearrayread,NULL),
193: PetscDesignatedInitializer(stridesubsetgather,NULL),
194: PetscDesignatedInitializer(stridesubsetscatter,NULL),
195: PetscDesignatedInitializer(viewnative,NULL),
196: PetscDesignatedInitializer(loadnative,NULL),
197: PetscDesignatedInitializer(getlocalvector,NULL)
198: };
200: PetscErrorCode VecDuplicateVecs_Comp(Vec w,PetscInt m,Vec *V[])
201: {
202: PetscInt i;
204: PetscFunctionBegin;
206: PetscAssertPointer(V,3);
207: PetscCheck(m>0,PetscObjectComm((PetscObject)w),PETSC_ERR_ARG_OUTOFRANGE,"m must be > 0: m = %" PetscInt_FMT,m);
208: PetscCall(PetscMalloc1(m,V));
209: for (i=0;i<m;i++) PetscCall(VecDuplicate(w,*V+i));
210: PetscFunctionReturn(PETSC_SUCCESS);
211: }
213: PetscErrorCode VecDestroyVecs_Comp(PetscInt m,Vec v[])
214: {
215: PetscInt i;
217: PetscFunctionBegin;
218: PetscAssertPointer(v,2);
219: PetscCheck(m>0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"m must be > 0: m = %" PetscInt_FMT,m);
220: for (i=0;i<m;i++) PetscCall(VecDestroy(&v[i]));
221: PetscCall(PetscFree(v));
222: PetscFunctionReturn(PETSC_SUCCESS);
223: }
225: static PetscErrorCode VecCreate_Comp_Private(Vec v,Vec *x,PetscInt nx,PetscBool x_to_me,Vec_Comp_N *n)
226: {
227: Vec_Comp *s;
228: PetscInt N=0,lN=0,i,k;
230: PetscFunctionBegin;
231: if (!VecCompInitialized) {
232: VecCompInitialized = PETSC_TRUE;
233: PetscCall(VecRegister(VECCOMP,VecCreate_Comp));
234: PetscCall(VecCompNormInit());
235: }
237: /* Allocate a new Vec_Comp */
238: if (v->data) PetscCall(PetscFree(v->data));
239: PetscCall(PetscNew(&s));
240: PetscCall(PetscMemcpy(v->ops,&DvOps,sizeof(DvOps)));
241: v->data = (void*)s;
242: v->petscnative = PETSC_FALSE;
244: /* Allocate the array of Vec, if it is needed to be done */
245: if (!x_to_me) {
246: if (nx) PetscCall(PetscMalloc1(nx,&s->x));
247: if (x) PetscCall(PetscArraycpy(s->x,x,nx));
248: } else s->x = x;
250: s->nx = nx;
252: if (nx && x) {
253: /* Allocate the shared structure, if it is not given */
254: if (!n) {
255: for (i=0;i<nx;i++) {
256: PetscCall(VecGetSize(x[i],&k));
257: N+= k;
258: PetscCall(VecGetLocalSize(x[i],&k));
259: lN+= k;
260: }
261: PetscCall(PetscNew(&n));
262: s->n = n;
263: n->n = nx;
264: n->N = N;
265: n->lN = lN;
266: n->friends = 1;
267: } else { /* If not, check in the vector in the shared structure */
268: s->n = n;
269: s->n->friends++;
270: }
272: /* Set the virtual sizes as the real sizes of the vector */
273: PetscCall(VecSetSizes(v,s->n->lN,s->n->N));
274: }
276: PetscCall(PetscObjectChangeTypeName((PetscObject)v,VECCOMP));
277: PetscCall(PetscObjectComposeFunction((PetscObject)v,"VecCompSetSubVecs_C",VecCompSetSubVecs_Comp));
278: PetscCall(PetscObjectComposeFunction((PetscObject)v,"VecCompGetSubVecs_C",VecCompGetSubVecs_Comp));
279: PetscFunctionReturn(PETSC_SUCCESS);
280: }
282: SLEPC_EXTERN PetscErrorCode VecCreate_Comp(Vec V)
283: {
284: PetscFunctionBegin;
285: PetscCall(VecCreate_Comp_Private(V,NULL,0,PETSC_FALSE,NULL));
286: PetscFunctionReturn(PETSC_SUCCESS);
287: }
289: /*MC
290: VECCOMP - VECCOMP = "comp" - Vector type consisting of several subvectors,
291: each stored separately.
293: Level: developer
295: Notes:
296: This is similar to PETSc's `VECNEST` but customized for SLEPc's needs. In particular,
297: the number of child vectors can be modified dynamically, with `VecCompSetSubVecs()`.
299: .seealso: `Vec`, `VecType`, `VecCreateComp()`, `VecCreateCompWithVecs()`
300: M*/
302: /*@
303: VecCreateComp - Creates a new vector containing several subvectors,
304: each stored separately.
306: Collective
308: Input Parameters:
309: + comm - communicator for the new `Vec`
310: . Nx - array of (initial) global sizes of child vectors
311: . n - number of child vectors
312: . t - type of the child vectors
313: - Vparent - (optional) template vector
315: Output Parameter:
316: . V - new vector
318: Notes:
319: This is similar to PETSc's `VECNEST` but customized for SLEPc's needs. In particular,
320: the number of child vectors can be modified dynamically, with `VecCompSetSubVecs()`.
322: Level: developer
324: .seealso: `VecCreateCompWithVecs()`, `VecCompSetSubVecs()`
325: @*/
326: PetscErrorCode VecCreateComp(MPI_Comm comm,PetscInt Nx[],PetscInt n,VecType t,Vec Vparent,Vec *V)
327: {
328: Vec *x;
329: PetscInt i;
331: PetscFunctionBegin;
332: PetscCall(VecCreate(comm,V));
333: PetscCall(PetscMalloc1(n,&x));
334: for (i=0;i<n;i++) {
335: PetscCall(VecCreate(comm,&x[i]));
336: PetscCall(VecSetSizes(x[i],PETSC_DECIDE,Nx[i]));
337: PetscCall(VecSetType(x[i],t));
338: }
339: PetscCall(VecCreate_Comp_Private(*V,x,n,PETSC_TRUE,Vparent?((Vec_Comp*)Vparent->data)->n:NULL));
340: PetscFunctionReturn(PETSC_SUCCESS);
341: }
343: /*@
344: VecCreateCompWithVecs - Creates a new vector containing several subvectors,
345: each stored separately, from an array of `Vec`s.
347: Collective
349: Input Parameters:
350: + x - array of `Vec`s
351: . n - number of child vectors
352: - Vparent - (optional) template vector
354: Output Parameter:
355: . V - new vector
357: Level: developer
359: .seealso: `VecCreateComp()`
360: @*/
361: PetscErrorCode VecCreateCompWithVecs(Vec x[],PetscInt n,Vec Vparent,Vec *V)
362: {
363: PetscInt i;
365: PetscFunctionBegin;
366: PetscAssertPointer(x,1);
369: PetscCall(VecCreate(PetscObjectComm((PetscObject)x[0]),V));
370: for (i=0;i<n;i++) PetscCall(PetscObjectReference((PetscObject)x[i]));
371: PetscCall(VecCreate_Comp_Private(*V,x,n,PETSC_FALSE,Vparent?((Vec_Comp*)Vparent->data)->n:NULL));
372: PetscFunctionReturn(PETSC_SUCCESS);
373: }
375: PetscErrorCode VecDuplicate_Comp(Vec win,Vec *V)
376: {
377: Vec *x;
378: PetscInt i;
379: Vec_Comp *s = (Vec_Comp*)win->data;
381: PetscFunctionBegin;
382: SlepcValidVecComp(win,1);
383: PetscCall(VecCreate(PetscObjectComm((PetscObject)win),V));
384: PetscCall(PetscMalloc1(s->nx,&x));
385: for (i=0;i<s->nx;i++) {
386: if (s->x[i]) PetscCall(VecDuplicate(s->x[i],&x[i]));
387: else x[i] = NULL;
388: }
389: PetscCall(VecCreate_Comp_Private(*V,x,s->nx,PETSC_TRUE,s->n));
390: PetscFunctionReturn(PETSC_SUCCESS);
391: }
393: static PetscErrorCode VecCompGetSubVecs_Comp(Vec win,PetscInt *n,const Vec **x)
394: {
395: Vec_Comp *s = (Vec_Comp*)win->data;
397: PetscFunctionBegin;
398: if (x) *x = s->x;
399: if (n) *n = s->n->n;
400: PetscFunctionReturn(PETSC_SUCCESS);
401: }
403: /*@C
404: VecCompGetSubVecs - Returns the entire array of vectors defining a
405: compound vector.
407: Collective
409: Input Parameter:
410: . win - compound vector
412: Output Parameters:
413: + n - number of child vectors
414: - x - array of child vectors
416: Level: developer
418: .seealso: `VecCreateComp()`
419: @*/
420: PetscErrorCode VecCompGetSubVecs(Vec win,PetscInt *n,const Vec *x[])
421: {
422: PetscFunctionBegin;
424: PetscUseMethod(win,"VecCompGetSubVecs_C",(Vec,PetscInt*,const Vec**),(win,n,x));
425: PetscFunctionReturn(PETSC_SUCCESS);
426: }
428: static PetscErrorCode VecCompSetSubVecs_Comp(Vec win,PetscInt n,Vec *x)
429: {
430: Vec_Comp *s = (Vec_Comp*)win->data;
431: PetscInt i,N,nlocal;
432: Vec_Comp_N *nn;
434: PetscFunctionBegin;
435: PetscCheck(s,PetscObjectComm((PetscObject)win),PETSC_ERR_ORDER,"Must call VecSetSizes first");
436: if (!s->nx) {
437: /* vector has been created via VecCreate+VecSetType+VecSetSizes, so allocate data structures */
438: PetscCall(PetscMalloc1(n,&s->x));
439: PetscCall(VecGetSize(win,&N));
440: PetscCheck(N%n==0,PetscObjectComm((PetscObject)win),PETSC_ERR_SUP,"Global dimension %" PetscInt_FMT " is not divisible by %" PetscInt_FMT,N,n);
441: PetscCall(VecGetLocalSize(win,&nlocal));
442: PetscCheck(nlocal%n==0,PetscObjectComm((PetscObject)win),PETSC_ERR_SUP,"Local dimension %" PetscInt_FMT " is not divisible by %" PetscInt_FMT,nlocal,n);
443: s->nx = n;
444: for (i=0;i<n;i++) {
445: PetscCall(VecCreate(PetscObjectComm((PetscObject)win),&s->x[i]));
446: PetscCall(VecSetSizes(s->x[i],nlocal/n,N/n));
447: PetscCall(VecSetFromOptions(s->x[i]));
448: }
449: if (!s->n) {
450: PetscCall(PetscNew(&nn));
451: s->n = nn;
452: nn->N = N;
453: nn->lN = nlocal;
454: nn->friends = 1;
455: }
456: } else PetscCheck(n<=s->nx,PetscObjectComm((PetscObject)win),PETSC_ERR_SUP,"Number of child vectors cannot be larger than %" PetscInt_FMT,s->nx);
457: if (x) PetscCall(PetscArraycpy(s->x,x,n));
458: s->n->n = n;
459: PetscFunctionReturn(PETSC_SUCCESS);
460: }
462: /*@
463: VecCompSetSubVecs - Resets the number of subvectors defining a compound vector,
464: or replaces the subvectors.
466: Collective
468: Input Parameters:
469: + win - compound vector
470: . n - number of child vectors
471: - x - array of child vectors
473: Note:
474: It is not possible to increase the number of subvectors with respect to the
475: number set at its creation.
477: Level: developer
479: .seealso: `VecCreateComp()`, `VecCompGetSubVecs()`
480: @*/
481: PetscErrorCode VecCompSetSubVecs(Vec win,PetscInt n,Vec x[])
482: {
483: PetscFunctionBegin;
486: PetscTryMethod(win,"VecCompSetSubVecs_C",(Vec,PetscInt,Vec*),(win,n,x));
487: PetscFunctionReturn(PETSC_SUCCESS);
488: }
490: PetscErrorCode VecAXPY_Comp(Vec v,PetscScalar alpha,Vec w)
491: {
492: Vec_Comp *vs = (Vec_Comp*)v->data,*ws = (Vec_Comp*)w->data;
493: PetscInt i;
495: PetscFunctionBegin;
496: SlepcValidVecComp(v,1);
497: SlepcValidVecComp(w,3);
498: for (i=0;i<vs->n->n;i++) PetscCall(VecAXPY(vs->x[i],alpha,ws->x[i]));
499: PetscFunctionReturn(PETSC_SUCCESS);
500: }
502: PetscErrorCode VecAYPX_Comp(Vec v,PetscScalar alpha,Vec w)
503: {
504: Vec_Comp *vs = (Vec_Comp*)v->data,*ws = (Vec_Comp*)w->data;
505: PetscInt i;
507: PetscFunctionBegin;
508: SlepcValidVecComp(v,1);
509: SlepcValidVecComp(w,3);
510: for (i=0;i<vs->n->n;i++) PetscCall(VecAYPX(vs->x[i],alpha,ws->x[i]));
511: PetscFunctionReturn(PETSC_SUCCESS);
512: }
514: PetscErrorCode VecAXPBY_Comp(Vec v,PetscScalar alpha,PetscScalar beta,Vec w)
515: {
516: Vec_Comp *vs = (Vec_Comp*)v->data,*ws = (Vec_Comp*)w->data;
517: PetscInt i;
519: PetscFunctionBegin;
520: SlepcValidVecComp(v,1);
521: SlepcValidVecComp(w,4);
522: for (i=0;i<vs->n->n;i++) PetscCall(VecAXPBY(vs->x[i],alpha,beta,ws->x[i]));
523: PetscFunctionReturn(PETSC_SUCCESS);
524: }
526: PetscErrorCode VecMAXPY_Comp(Vec v,PetscInt n,const PetscScalar *alpha,Vec *w)
527: {
528: Vec_Comp *vs = (Vec_Comp*)v->data;
529: Vec *wx;
530: PetscInt i,j;
532: PetscFunctionBegin;
533: SlepcValidVecComp(v,1);
534: for (i=0;i<n;i++) SlepcValidVecComp(w[i],4);
536: PetscCall(PetscMalloc1(n,&wx));
538: for (j=0;j<vs->n->n;j++) {
539: for (i=0;i<n;i++) wx[i] = ((Vec_Comp*)w[i]->data)->x[j];
540: PetscCall(VecMAXPY(vs->x[j],n,alpha,wx));
541: }
543: PetscCall(PetscFree(wx));
544: PetscFunctionReturn(PETSC_SUCCESS);
545: }
547: PetscErrorCode VecWAXPY_Comp(Vec v,PetscScalar alpha,Vec w,Vec z)
548: {
549: Vec_Comp *vs = (Vec_Comp*)v->data,*ws = (Vec_Comp*)w->data,*zs = (Vec_Comp*)z->data;
550: PetscInt i;
552: PetscFunctionBegin;
553: SlepcValidVecComp(v,1);
554: SlepcValidVecComp(w,3);
555: SlepcValidVecComp(z,4);
556: for (i=0;i<vs->n->n;i++) PetscCall(VecWAXPY(vs->x[i],alpha,ws->x[i],zs->x[i]));
557: PetscFunctionReturn(PETSC_SUCCESS);
558: }
560: PetscErrorCode VecAXPBYPCZ_Comp(Vec v,PetscScalar alpha,PetscScalar beta,PetscScalar gamma,Vec w,Vec z)
561: {
562: Vec_Comp *vs = (Vec_Comp*)v->data,*ws = (Vec_Comp*)w->data,*zs = (Vec_Comp*)z->data;
563: PetscInt i;
565: PetscFunctionBegin;
566: SlepcValidVecComp(v,1);
567: SlepcValidVecComp(w,5);
568: SlepcValidVecComp(z,6);
569: for (i=0;i<vs->n->n;i++) PetscCall(VecAXPBYPCZ(vs->x[i],alpha,beta,gamma,ws->x[i],zs->x[i]));
570: PetscFunctionReturn(PETSC_SUCCESS);
571: }
573: PetscErrorCode VecGetSize_Comp(Vec v,PetscInt *size)
574: {
575: Vec_Comp *vs = (Vec_Comp*)v->data;
577: PetscFunctionBegin;
578: PetscAssertPointer(size,2);
579: if (vs->n) {
580: SlepcValidVecComp(v,1);
581: *size = vs->n->N;
582: } else *size = v->map->N;
583: PetscFunctionReturn(PETSC_SUCCESS);
584: }
586: PetscErrorCode VecGetLocalSize_Comp(Vec v,PetscInt *size)
587: {
588: Vec_Comp *vs = (Vec_Comp*)v->data;
590: PetscFunctionBegin;
591: PetscAssertPointer(size,2);
592: if (vs->n) {
593: SlepcValidVecComp(v,1);
594: *size = vs->n->lN;
595: } else *size = v->map->n;
596: PetscFunctionReturn(PETSC_SUCCESS);
597: }
599: PetscErrorCode VecMax_Comp(Vec v,PetscInt *idx,PetscReal *z)
600: {
601: Vec_Comp *vs = (Vec_Comp*)v->data;
602: PetscInt idxp,s=0,s0;
603: PetscReal zp,z0;
604: PetscInt i;
606: PetscFunctionBegin;
607: SlepcValidVecComp(v,1);
608: if (!idx && !z) PetscFunctionReturn(PETSC_SUCCESS);
610: if (vs->n->n > 0) PetscCall(VecMax(vs->x[0],idx?&idxp:NULL,&zp));
611: else {
612: zp = PETSC_MIN_REAL;
613: if (idx) idxp = -1;
614: }
615: for (i=1;i<vs->n->n;i++) {
616: PetscCall(VecGetSize(vs->x[i-1],&s0));
617: s += s0;
618: PetscCall(VecMax(vs->x[i],idx?&idxp:NULL,&z0));
619: if (zp < z0) {
620: if (idx) *idx = s+idxp;
621: zp = z0;
622: }
623: }
624: if (z) *z = zp;
625: PetscFunctionReturn(PETSC_SUCCESS);
626: }
628: PetscErrorCode VecMin_Comp(Vec v,PetscInt *idx,PetscReal *z)
629: {
630: Vec_Comp *vs = (Vec_Comp*)v->data;
631: PetscInt idxp,s=0,s0;
632: PetscReal zp,z0;
633: PetscInt i;
635: PetscFunctionBegin;
636: SlepcValidVecComp(v,1);
637: if (!idx && !z) PetscFunctionReturn(PETSC_SUCCESS);
639: if (vs->n->n > 0) PetscCall(VecMin(vs->x[0],idx?&idxp:NULL,&zp));
640: else {
641: zp = PETSC_MAX_REAL;
642: if (idx) idxp = -1;
643: }
644: for (i=1;i<vs->n->n;i++) {
645: PetscCall(VecGetSize(vs->x[i-1],&s0));
646: s += s0;
647: PetscCall(VecMin(vs->x[i],idx?&idxp:NULL,&z0));
648: if (zp > z0) {
649: if (idx) *idx = s+idxp;
650: zp = z0;
651: }
652: }
653: if (z) *z = zp;
654: PetscFunctionReturn(PETSC_SUCCESS);
655: }
657: PetscErrorCode VecMaxPointwiseDivide_Comp(Vec v,Vec w,PetscReal *m)
658: {
659: Vec_Comp *vs = (Vec_Comp*)v->data,*ws = (Vec_Comp*)w->data;
660: PetscReal work;
661: PetscInt i;
663: PetscFunctionBegin;
664: SlepcValidVecComp(v,1);
665: SlepcValidVecComp(w,2);
666: if (!m || vs->n->n == 0) PetscFunctionReturn(PETSC_SUCCESS);
667: PetscCall(VecMaxPointwiseDivide(vs->x[0],ws->x[0],m));
668: for (i=1;i<vs->n->n;i++) {
669: PetscCall(VecMaxPointwiseDivide(vs->x[i],ws->x[i],&work));
670: *m = PetscMax(*m,work);
671: }
672: PetscFunctionReturn(PETSC_SUCCESS);
673: }
680: PetscErrorCode __COMPOSE3__(Vec,NAME,_Comp)(Vec v) \
681: { \
682: Vec_Comp *vs = (Vec_Comp*)v->data; \
683: PetscInt i; \
684: \
685: PetscFunctionBegin; \
686: SlepcValidVecComp(v,1); \
687: for (i=0;i<vs->n->n;i++) { \
688: PetscCall(__COMPOSE2__(Vec,NAME)(vs->x[i])); \
689: } \
690: PetscFunctionReturn(PETSC_SUCCESS);\
691: }
693: __FUNC_TEMPLATE1__(Conjugate)
694: __FUNC_TEMPLATE1__(Reciprocal)
695: __FUNC_TEMPLATE1__(SqrtAbs)
696: __FUNC_TEMPLATE1__(Abs)
697: __FUNC_TEMPLATE1__(Exp)
698: __FUNC_TEMPLATE1__(Log)
701: PetscErrorCode __COMPOSE3__(Vec,NAME,_Comp)(Vec v,T0 __a) \
702: { \
703: Vec_Comp *vs = (Vec_Comp*)v->data; \
704: PetscInt i; \
705: \
706: PetscFunctionBegin; \
707: SlepcValidVecComp(v,1); \
708: for (i=0;i<vs->n->n;i++) { \
709: PetscCall(__COMPOSE2__(Vec,NAME)(vs->x[i],__a)); \
710: } \
711: PetscFunctionReturn(PETSC_SUCCESS);\
712: }
714: __FUNC_TEMPLATE2__(Set,PetscScalar)
715: __FUNC_TEMPLATE2__(View,PetscViewer)
716: __FUNC_TEMPLATE2__(Scale,PetscScalar)
717: __FUNC_TEMPLATE2__(SetRandom,PetscRandom)
718: __FUNC_TEMPLATE2__(Shift,PetscScalar)
721: PetscErrorCode __COMPOSE3__(Vec,NAME,_Comp)(Vec v,Vec w) \
722: { \
723: Vec_Comp *vs = (Vec_Comp*)v->data,\
724: *ws = (Vec_Comp*)w->data; \
725: PetscInt i; \
726: \
727: PetscFunctionBegin; \
728: SlepcValidVecComp(v,1); \
729: SlepcValidVecComp(w,2); \
730: for (i=0;i<vs->n->n;i++) { \
731: PetscCall(__COMPOSE2__(Vec,NAME)(vs->x[i],ws->x[i])); \
732: } \
733: PetscFunctionReturn(PETSC_SUCCESS);\
734: }
736: __FUNC_TEMPLATE3__(Copy)
737: __FUNC_TEMPLATE3__(Swap)
740: PetscErrorCode __COMPOSE3__(Vec,NAME,_Comp)(Vec v,Vec w,Vec z) \
741: { \
742: Vec_Comp *vs = (Vec_Comp*)v->data, \
743: *ws = (Vec_Comp*)w->data, \
744: *zs = (Vec_Comp*)z->data; \
745: PetscInt i; \
746: \
747: PetscFunctionBegin; \
748: SlepcValidVecComp(v,1); \
749: SlepcValidVecComp(w,2); \
750: SlepcValidVecComp(z,3); \
751: for (i=0;i<vs->n->n;i++) { \
752: PetscCall(__COMPOSE2__(Vec,NAME)(vs->x[i],ws->x[i],zs->x[i])); \
753: } \
754: PetscFunctionReturn(PETSC_SUCCESS);\
755: }
757: __FUNC_TEMPLATE4__(PointwiseMax)
758: __FUNC_TEMPLATE4__(PointwiseMaxAbs)
759: __FUNC_TEMPLATE4__(PointwiseMin)
760: __FUNC_TEMPLATE4__(PointwiseMult)
761: __FUNC_TEMPLATE4__(PointwiseDivide)