Actual source code: veccomp.c
slepc-3.22.1 2024-10-28
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(setvalueslocal,NULL),
174: PetscDesignatedInitializer(resetarray,NULL),
175: PetscDesignatedInitializer(setfromoptions,NULL),
176: PetscDesignatedInitializer(maxpointwisedivide,VecMaxPointwiseDivide_Comp),
177: PetscDesignatedInitializer(pointwisemax,VecPointwiseMax_Comp),
178: PetscDesignatedInitializer(pointwisemaxabs,VecPointwiseMaxAbs_Comp),
179: PetscDesignatedInitializer(pointwisemin,VecPointwiseMin_Comp),
180: PetscDesignatedInitializer(getvalues,NULL),
181: PetscDesignatedInitializer(sqrt,VecSqrtAbs_Comp),
182: PetscDesignatedInitializer(abs,VecAbs_Comp),
183: PetscDesignatedInitializer(exp,VecExp_Comp),
184: PetscDesignatedInitializer(log,VecLog_Comp),
185: PetscDesignatedInitializer(shift,VecShift_Comp),
186: PetscDesignatedInitializer(create,NULL),
187: PetscDesignatedInitializer(stridegather,NULL),
188: PetscDesignatedInitializer(stridescatter,NULL),
189: PetscDesignatedInitializer(dotnorm2,VecDotNorm2_Comp_MPI),
190: PetscDesignatedInitializer(getsubvector,NULL),
191: PetscDesignatedInitializer(restoresubvector,NULL),
192: PetscDesignatedInitializer(getarrayread,NULL),
193: PetscDesignatedInitializer(restorearrayread,NULL),
194: PetscDesignatedInitializer(stridesubsetgather,NULL),
195: PetscDesignatedInitializer(stridesubsetscatter,NULL),
196: PetscDesignatedInitializer(viewnative,NULL),
197: PetscDesignatedInitializer(loadnative,NULL),
198: PetscDesignatedInitializer(getlocalvector,NULL)
199: };
201: PetscErrorCode VecDuplicateVecs_Comp(Vec w,PetscInt m,Vec *V[])
202: {
203: PetscInt i;
205: PetscFunctionBegin;
207: PetscAssertPointer(V,3);
208: PetscCheck(m>0,PetscObjectComm((PetscObject)w),PETSC_ERR_ARG_OUTOFRANGE,"m must be > 0: m = %" PetscInt_FMT,m);
209: PetscCall(PetscMalloc1(m,V));
210: for (i=0;i<m;i++) PetscCall(VecDuplicate(w,*V+i));
211: PetscFunctionReturn(PETSC_SUCCESS);
212: }
214: PetscErrorCode VecDestroyVecs_Comp(PetscInt m,Vec v[])
215: {
216: PetscInt i;
218: PetscFunctionBegin;
219: PetscAssertPointer(v,2);
220: PetscCheck(m>0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"m must be > 0: m = %" PetscInt_FMT,m);
221: for (i=0;i<m;i++) PetscCall(VecDestroy(&v[i]));
222: PetscCall(PetscFree(v));
223: PetscFunctionReturn(PETSC_SUCCESS);
224: }
226: static PetscErrorCode VecCreate_Comp_Private(Vec v,Vec *x,PetscInt nx,PetscBool x_to_me,Vec_Comp_N *n)
227: {
228: Vec_Comp *s;
229: PetscInt N=0,lN=0,i,k;
231: PetscFunctionBegin;
232: if (!VecCompInitialized) {
233: VecCompInitialized = PETSC_TRUE;
234: PetscCall(VecRegister(VECCOMP,VecCreate_Comp));
235: PetscCall(VecCompNormInit());
236: }
238: /* Allocate a new Vec_Comp */
239: if (v->data) PetscCall(PetscFree(v->data));
240: PetscCall(PetscNew(&s));
241: PetscCall(PetscMemcpy(v->ops,&DvOps,sizeof(DvOps)));
242: v->data = (void*)s;
243: v->petscnative = PETSC_FALSE;
245: /* Allocate the array of Vec, if it is needed to be done */
246: if (!x_to_me) {
247: if (nx) PetscCall(PetscMalloc1(nx,&s->x));
248: if (x) PetscCall(PetscArraycpy(s->x,x,nx));
249: } else s->x = x;
251: s->nx = nx;
253: if (nx && x) {
254: /* Allocate the shared structure, if it is not given */
255: if (!n) {
256: for (i=0;i<nx;i++) {
257: PetscCall(VecGetSize(x[i],&k));
258: N+= k;
259: PetscCall(VecGetLocalSize(x[i],&k));
260: lN+= k;
261: }
262: PetscCall(PetscNew(&n));
263: s->n = n;
264: n->n = nx;
265: n->N = N;
266: n->lN = lN;
267: n->friends = 1;
268: } else { /* If not, check in the vector in the shared structure */
269: s->n = n;
270: s->n->friends++;
271: }
273: /* Set the virtual sizes as the real sizes of the vector */
274: PetscCall(VecSetSizes(v,s->n->lN,s->n->N));
275: }
277: PetscCall(PetscObjectChangeTypeName((PetscObject)v,VECCOMP));
278: PetscCall(PetscObjectComposeFunction((PetscObject)v,"VecCompSetSubVecs_C",VecCompSetSubVecs_Comp));
279: PetscCall(PetscObjectComposeFunction((PetscObject)v,"VecCompGetSubVecs_C",VecCompGetSubVecs_Comp));
280: PetscFunctionReturn(PETSC_SUCCESS);
281: }
283: SLEPC_EXTERN PetscErrorCode VecCreate_Comp(Vec V)
284: {
285: PetscFunctionBegin;
286: PetscCall(VecCreate_Comp_Private(V,NULL,0,PETSC_FALSE,NULL));
287: PetscFunctionReturn(PETSC_SUCCESS);
288: }
290: /*@
291: VecCreateComp - Creates a new vector containing several subvectors,
292: each stored separately.
294: Collective
296: Input Parameters:
297: + comm - communicator for the new Vec
298: . Nx - array of (initial) global sizes of child vectors
299: . n - number of child vectors
300: . t - type of the child vectors
301: - Vparent - (optional) template vector
303: Output Parameter:
304: . V - new vector
306: Notes:
307: This is similar to PETSc's VecNest but customized for SLEPc's needs. In particular,
308: the number of child vectors can be modified dynamically, with VecCompSetSubVecs().
310: Level: developer
312: .seealso: VecCreateCompWithVecs(), VecCompSetSubVecs()
313: @*/
314: PetscErrorCode VecCreateComp(MPI_Comm comm,PetscInt Nx[],PetscInt n,VecType t,Vec Vparent,Vec *V)
315: {
316: Vec *x;
317: PetscInt i;
319: PetscFunctionBegin;
320: PetscCall(VecCreate(comm,V));
321: PetscCall(PetscMalloc1(n,&x));
322: for (i=0;i<n;i++) {
323: PetscCall(VecCreate(comm,&x[i]));
324: PetscCall(VecSetSizes(x[i],PETSC_DECIDE,Nx[i]));
325: PetscCall(VecSetType(x[i],t));
326: }
327: PetscCall(VecCreate_Comp_Private(*V,x,n,PETSC_TRUE,Vparent?((Vec_Comp*)Vparent->data)->n:NULL));
328: PetscFunctionReturn(PETSC_SUCCESS);
329: }
331: /*@
332: VecCreateCompWithVecs - Creates a new vector containing several subvectors,
333: each stored separately, from an array of Vecs.
335: Collective
337: Input Parameters:
338: + x - array of Vecs
339: . n - number of child vectors
340: - Vparent - (optional) template vector
342: Output Parameter:
343: . V - new vector
345: Level: developer
347: .seealso: VecCreateComp()
348: @*/
349: PetscErrorCode VecCreateCompWithVecs(Vec x[],PetscInt n,Vec Vparent,Vec *V)
350: {
351: PetscInt i;
353: PetscFunctionBegin;
354: PetscAssertPointer(x,1);
357: PetscCall(VecCreate(PetscObjectComm((PetscObject)x[0]),V));
358: for (i=0;i<n;i++) PetscCall(PetscObjectReference((PetscObject)x[i]));
359: PetscCall(VecCreate_Comp_Private(*V,x,n,PETSC_FALSE,Vparent?((Vec_Comp*)Vparent->data)->n:NULL));
360: PetscFunctionReturn(PETSC_SUCCESS);
361: }
363: PetscErrorCode VecDuplicate_Comp(Vec win,Vec *V)
364: {
365: Vec *x;
366: PetscInt i;
367: Vec_Comp *s = (Vec_Comp*)win->data;
369: PetscFunctionBegin;
370: SlepcValidVecComp(win,1);
371: PetscCall(VecCreate(PetscObjectComm((PetscObject)win),V));
372: PetscCall(PetscMalloc1(s->nx,&x));
373: for (i=0;i<s->nx;i++) {
374: if (s->x[i]) PetscCall(VecDuplicate(s->x[i],&x[i]));
375: else x[i] = NULL;
376: }
377: PetscCall(VecCreate_Comp_Private(*V,x,s->nx,PETSC_TRUE,s->n));
378: PetscFunctionReturn(PETSC_SUCCESS);
379: }
381: static PetscErrorCode VecCompGetSubVecs_Comp(Vec win,PetscInt *n,const Vec **x)
382: {
383: Vec_Comp *s = (Vec_Comp*)win->data;
385: PetscFunctionBegin;
386: if (x) *x = s->x;
387: if (n) *n = s->n->n;
388: PetscFunctionReturn(PETSC_SUCCESS);
389: }
391: /*@C
392: VecCompGetSubVecs - Returns the entire array of vectors defining a
393: compound vector.
395: Collective
397: Input Parameter:
398: . win - compound vector
400: Output Parameters:
401: + n - number of child vectors
402: - x - array of child vectors
404: Level: developer
406: .seealso: VecCreateComp()
407: @*/
408: PetscErrorCode VecCompGetSubVecs(Vec win,PetscInt *n,const Vec **x)
409: {
410: PetscFunctionBegin;
412: PetscUseMethod(win,"VecCompGetSubVecs_C",(Vec,PetscInt*,const Vec**),(win,n,x));
413: PetscFunctionReturn(PETSC_SUCCESS);
414: }
416: static PetscErrorCode VecCompSetSubVecs_Comp(Vec win,PetscInt n,Vec *x)
417: {
418: Vec_Comp *s = (Vec_Comp*)win->data;
419: PetscInt i,N,nlocal;
420: Vec_Comp_N *nn;
422: PetscFunctionBegin;
423: PetscCheck(s,PetscObjectComm((PetscObject)win),PETSC_ERR_ORDER,"Must call VecSetSizes first");
424: if (!s->nx) {
425: /* vector has been created via VecCreate+VecSetType+VecSetSizes, so allocate data structures */
426: PetscCall(PetscMalloc1(n,&s->x));
427: PetscCall(VecGetSize(win,&N));
428: PetscCheck(N%n==0,PetscObjectComm((PetscObject)win),PETSC_ERR_SUP,"Global dimension %" PetscInt_FMT " is not divisible by %" PetscInt_FMT,N,n);
429: PetscCall(VecGetLocalSize(win,&nlocal));
430: PetscCheck(nlocal%n==0,PetscObjectComm((PetscObject)win),PETSC_ERR_SUP,"Local dimension %" PetscInt_FMT " is not divisible by %" PetscInt_FMT,nlocal,n);
431: s->nx = n;
432: for (i=0;i<n;i++) {
433: PetscCall(VecCreate(PetscObjectComm((PetscObject)win),&s->x[i]));
434: PetscCall(VecSetSizes(s->x[i],nlocal/n,N/n));
435: PetscCall(VecSetFromOptions(s->x[i]));
436: }
437: if (!s->n) {
438: PetscCall(PetscNew(&nn));
439: s->n = nn;
440: nn->N = N;
441: nn->lN = nlocal;
442: nn->friends = 1;
443: }
444: } else PetscCheck(n<=s->nx,PetscObjectComm((PetscObject)win),PETSC_ERR_SUP,"Number of child vectors cannot be larger than %" PetscInt_FMT,s->nx);
445: if (x) PetscCall(PetscArraycpy(s->x,x,n));
446: s->n->n = n;
447: PetscFunctionReturn(PETSC_SUCCESS);
448: }
450: /*@
451: VecCompSetSubVecs - Resets the number of subvectors defining a compound vector,
452: or replaces the subvectors.
454: Collective
456: Input Parameters:
457: + win - compound vector
458: . n - number of child vectors
459: - x - array of child vectors
461: Note:
462: It is not possible to increase the number of subvectors with respect to the
463: number set at its creation.
465: Level: developer
467: .seealso: VecCreateComp(), VecCompGetSubVecs()
468: @*/
469: PetscErrorCode VecCompSetSubVecs(Vec win,PetscInt n,Vec x[])
470: {
471: PetscFunctionBegin;
474: PetscTryMethod(win,"VecCompSetSubVecs_C",(Vec,PetscInt,Vec*),(win,n,x));
475: PetscFunctionReturn(PETSC_SUCCESS);
476: }
478: PetscErrorCode VecAXPY_Comp(Vec v,PetscScalar alpha,Vec w)
479: {
480: Vec_Comp *vs = (Vec_Comp*)v->data,*ws = (Vec_Comp*)w->data;
481: PetscInt i;
483: PetscFunctionBegin;
484: SlepcValidVecComp(v,1);
485: SlepcValidVecComp(w,3);
486: for (i=0;i<vs->n->n;i++) PetscCall(VecAXPY(vs->x[i],alpha,ws->x[i]));
487: PetscFunctionReturn(PETSC_SUCCESS);
488: }
490: PetscErrorCode VecAYPX_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(VecAYPX(vs->x[i],alpha,ws->x[i]));
499: PetscFunctionReturn(PETSC_SUCCESS);
500: }
502: PetscErrorCode VecAXPBY_Comp(Vec v,PetscScalar alpha,PetscScalar beta,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,4);
510: for (i=0;i<vs->n->n;i++) PetscCall(VecAXPBY(vs->x[i],alpha,beta,ws->x[i]));
511: PetscFunctionReturn(PETSC_SUCCESS);
512: }
514: PetscErrorCode VecMAXPY_Comp(Vec v,PetscInt n,const PetscScalar *alpha,Vec *w)
515: {
516: Vec_Comp *vs = (Vec_Comp*)v->data;
517: Vec *wx;
518: PetscInt i,j;
520: PetscFunctionBegin;
521: SlepcValidVecComp(v,1);
522: for (i=0;i<n;i++) SlepcValidVecComp(w[i],4);
524: PetscCall(PetscMalloc1(n,&wx));
526: for (j=0;j<vs->n->n;j++) {
527: for (i=0;i<n;i++) wx[i] = ((Vec_Comp*)w[i]->data)->x[j];
528: PetscCall(VecMAXPY(vs->x[j],n,alpha,wx));
529: }
531: PetscCall(PetscFree(wx));
532: PetscFunctionReturn(PETSC_SUCCESS);
533: }
535: PetscErrorCode VecWAXPY_Comp(Vec v,PetscScalar alpha,Vec w,Vec z)
536: {
537: Vec_Comp *vs = (Vec_Comp*)v->data,*ws = (Vec_Comp*)w->data,*zs = (Vec_Comp*)z->data;
538: PetscInt i;
540: PetscFunctionBegin;
541: SlepcValidVecComp(v,1);
542: SlepcValidVecComp(w,3);
543: SlepcValidVecComp(z,4);
544: for (i=0;i<vs->n->n;i++) PetscCall(VecWAXPY(vs->x[i],alpha,ws->x[i],zs->x[i]));
545: PetscFunctionReturn(PETSC_SUCCESS);
546: }
548: PetscErrorCode VecAXPBYPCZ_Comp(Vec v,PetscScalar alpha,PetscScalar beta,PetscScalar gamma,Vec w,Vec z)
549: {
550: Vec_Comp *vs = (Vec_Comp*)v->data,*ws = (Vec_Comp*)w->data,*zs = (Vec_Comp*)z->data;
551: PetscInt i;
553: PetscFunctionBegin;
554: SlepcValidVecComp(v,1);
555: SlepcValidVecComp(w,5);
556: SlepcValidVecComp(z,6);
557: for (i=0;i<vs->n->n;i++) PetscCall(VecAXPBYPCZ(vs->x[i],alpha,beta,gamma,ws->x[i],zs->x[i]));
558: PetscFunctionReturn(PETSC_SUCCESS);
559: }
561: PetscErrorCode VecGetSize_Comp(Vec v,PetscInt *size)
562: {
563: Vec_Comp *vs = (Vec_Comp*)v->data;
565: PetscFunctionBegin;
566: PetscAssertPointer(size,2);
567: if (vs->n) {
568: SlepcValidVecComp(v,1);
569: *size = vs->n->N;
570: } else *size = v->map->N;
571: PetscFunctionReturn(PETSC_SUCCESS);
572: }
574: PetscErrorCode VecGetLocalSize_Comp(Vec v,PetscInt *size)
575: {
576: Vec_Comp *vs = (Vec_Comp*)v->data;
578: PetscFunctionBegin;
579: PetscAssertPointer(size,2);
580: if (vs->n) {
581: SlepcValidVecComp(v,1);
582: *size = vs->n->lN;
583: } else *size = v->map->n;
584: PetscFunctionReturn(PETSC_SUCCESS);
585: }
587: PetscErrorCode VecMax_Comp(Vec v,PetscInt *idx,PetscReal *z)
588: {
589: Vec_Comp *vs = (Vec_Comp*)v->data;
590: PetscInt idxp,s=0,s0;
591: PetscReal zp,z0;
592: PetscInt i;
594: PetscFunctionBegin;
595: SlepcValidVecComp(v,1);
596: if (!idx && !z) PetscFunctionReturn(PETSC_SUCCESS);
598: if (vs->n->n > 0) PetscCall(VecMax(vs->x[0],idx?&idxp:NULL,&zp));
599: else {
600: zp = PETSC_MIN_REAL;
601: if (idx) idxp = -1;
602: }
603: for (i=1;i<vs->n->n;i++) {
604: PetscCall(VecGetSize(vs->x[i-1],&s0));
605: s += s0;
606: PetscCall(VecMax(vs->x[i],idx?&idxp:NULL,&z0));
607: if (zp < z0) {
608: if (idx) *idx = s+idxp;
609: zp = z0;
610: }
611: }
612: if (z) *z = zp;
613: PetscFunctionReturn(PETSC_SUCCESS);
614: }
616: PetscErrorCode VecMin_Comp(Vec v,PetscInt *idx,PetscReal *z)
617: {
618: Vec_Comp *vs = (Vec_Comp*)v->data;
619: PetscInt idxp,s=0,s0;
620: PetscReal zp,z0;
621: PetscInt i;
623: PetscFunctionBegin;
624: SlepcValidVecComp(v,1);
625: if (!idx && !z) PetscFunctionReturn(PETSC_SUCCESS);
627: if (vs->n->n > 0) PetscCall(VecMin(vs->x[0],idx?&idxp:NULL,&zp));
628: else {
629: zp = PETSC_MAX_REAL;
630: if (idx) idxp = -1;
631: }
632: for (i=1;i<vs->n->n;i++) {
633: PetscCall(VecGetSize(vs->x[i-1],&s0));
634: s += s0;
635: PetscCall(VecMin(vs->x[i],idx?&idxp:NULL,&z0));
636: if (zp > z0) {
637: if (idx) *idx = s+idxp;
638: zp = z0;
639: }
640: }
641: if (z) *z = zp;
642: PetscFunctionReturn(PETSC_SUCCESS);
643: }
645: PetscErrorCode VecMaxPointwiseDivide_Comp(Vec v,Vec w,PetscReal *m)
646: {
647: Vec_Comp *vs = (Vec_Comp*)v->data,*ws = (Vec_Comp*)w->data;
648: PetscReal work;
649: PetscInt i;
651: PetscFunctionBegin;
652: SlepcValidVecComp(v,1);
653: SlepcValidVecComp(w,2);
654: if (!m || vs->n->n == 0) PetscFunctionReturn(PETSC_SUCCESS);
655: PetscCall(VecMaxPointwiseDivide(vs->x[0],ws->x[0],m));
656: for (i=1;i<vs->n->n;i++) {
657: PetscCall(VecMaxPointwiseDivide(vs->x[i],ws->x[i],&work));
658: *m = PetscMax(*m,work);
659: }
660: PetscFunctionReturn(PETSC_SUCCESS);
661: }
668: PetscErrorCode __COMPOSE3__(Vec,NAME,_Comp)(Vec v) \
669: { \
670: Vec_Comp *vs = (Vec_Comp*)v->data; \
671: PetscInt i; \
672: \
673: PetscFunctionBegin; \
674: SlepcValidVecComp(v,1); \
675: for (i=0;i<vs->n->n;i++) { \
676: PetscCall(__COMPOSE2__(Vec,NAME)(vs->x[i])); \
677: } \
678: PetscFunctionReturn(PETSC_SUCCESS);\
679: }
681: __FUNC_TEMPLATE1__(Conjugate)
682: __FUNC_TEMPLATE1__(Reciprocal)
683: __FUNC_TEMPLATE1__(SqrtAbs)
684: __FUNC_TEMPLATE1__(Abs)
685: __FUNC_TEMPLATE1__(Exp)
686: __FUNC_TEMPLATE1__(Log)
689: PetscErrorCode __COMPOSE3__(Vec,NAME,_Comp)(Vec v,T0 __a) \
690: { \
691: Vec_Comp *vs = (Vec_Comp*)v->data; \
692: PetscInt i; \
693: \
694: PetscFunctionBegin; \
695: SlepcValidVecComp(v,1); \
696: for (i=0;i<vs->n->n;i++) { \
697: PetscCall(__COMPOSE2__(Vec,NAME)(vs->x[i],__a)); \
698: } \
699: PetscFunctionReturn(PETSC_SUCCESS);\
700: }
702: __FUNC_TEMPLATE2__(Set,PetscScalar)
703: __FUNC_TEMPLATE2__(View,PetscViewer)
704: __FUNC_TEMPLATE2__(Scale,PetscScalar)
705: __FUNC_TEMPLATE2__(SetRandom,PetscRandom)
706: __FUNC_TEMPLATE2__(Shift,PetscScalar)
709: PetscErrorCode __COMPOSE3__(Vec,NAME,_Comp)(Vec v,Vec w) \
710: { \
711: Vec_Comp *vs = (Vec_Comp*)v->data,\
712: *ws = (Vec_Comp*)w->data; \
713: PetscInt i; \
714: \
715: PetscFunctionBegin; \
716: SlepcValidVecComp(v,1); \
717: SlepcValidVecComp(w,2); \
718: for (i=0;i<vs->n->n;i++) { \
719: PetscCall(__COMPOSE2__(Vec,NAME)(vs->x[i],ws->x[i])); \
720: } \
721: PetscFunctionReturn(PETSC_SUCCESS);\
722: }
724: __FUNC_TEMPLATE3__(Copy)
725: __FUNC_TEMPLATE3__(Swap)
728: PetscErrorCode __COMPOSE3__(Vec,NAME,_Comp)(Vec v,Vec w,Vec z) \
729: { \
730: Vec_Comp *vs = (Vec_Comp*)v->data, \
731: *ws = (Vec_Comp*)w->data, \
732: *zs = (Vec_Comp*)z->data; \
733: PetscInt i; \
734: \
735: PetscFunctionBegin; \
736: SlepcValidVecComp(v,1); \
737: SlepcValidVecComp(w,2); \
738: SlepcValidVecComp(z,3); \
739: for (i=0;i<vs->n->n;i++) { \
740: PetscCall(__COMPOSE2__(Vec,NAME)(vs->x[i],ws->x[i],zs->x[i])); \
741: } \
742: PetscFunctionReturn(PETSC_SUCCESS);\
743: }
745: __FUNC_TEMPLATE4__(PointwiseMax)
746: __FUNC_TEMPLATE4__(PointwiseMaxAbs)
747: __FUNC_TEMPLATE4__(PointwiseMin)
748: __FUNC_TEMPLATE4__(PointwiseMult)
749: __FUNC_TEMPLATE4__(PointwiseDivide)