Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
6 : This file is part of SLEPc.
7 : SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9 : */
10 :
11 : #include <slepc/private/vecimplslepc.h> /*I "slepcvec.h" I*/
12 :
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;
17 :
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**);
24 :
25 : #include "veccomp0.h"
26 :
27 : #define __WITH_MPI__
28 : #include "veccomp0.h"
29 :
30 1608 : static inline void SumNorm2(PetscReal *ssq0,PetscReal *scale0,PetscReal *ssq1,PetscReal *scale1)
31 : {
32 1608 : PetscReal q;
33 1608 : if (*scale0 > *scale1) {
34 804 : q = *scale1/(*scale0);
35 804 : *ssq1 = *ssq0 + q*q*(*ssq1);
36 804 : *scale1 = *scale0;
37 : } else {
38 804 : q = *scale0/(*scale1);
39 804 : *ssq1 += q*q*(*ssq0);
40 : }
41 1608 : }
42 :
43 17860 : static inline PetscReal GetNorm2(PetscReal ssq,PetscReal scale)
44 : {
45 17860 : return scale*PetscSqrtReal(ssq);
46 : }
47 :
48 18171 : static inline void AddNorm2(PetscReal *ssq,PetscReal *scale,PetscReal x)
49 : {
50 18171 : PetscReal absx,q;
51 18171 : if (x != 0.0) {
52 18168 : absx = PetscAbs(x);
53 18168 : if (*scale < absx) {
54 17980 : q = *scale/absx;
55 17980 : *ssq = 1.0 + *ssq*q*q;
56 17980 : *scale = absx;
57 : } else {
58 188 : q = absx/(*scale);
59 188 : *ssq += q*q;
60 : }
61 : }
62 18171 : }
63 :
64 1608 : SLEPC_EXTERN void MPIAPI SlepcSumNorm2_Local(void *in,void *out,PetscMPIInt *cnt,MPI_Datatype *datatype)
65 : {
66 1608 : PetscInt i,count = *cnt;
67 1608 : PetscReal *xin = (PetscReal*)in,*xout = (PetscReal*)out;
68 :
69 1608 : PetscFunctionBegin;
70 1608 : if (*datatype == MPIU_NORM2) {
71 3212 : for (i=0;i<count;i++) {
72 1606 : SumNorm2(&xin[i*2],&xin[i*2+1],&xout[i*2],&xout[i*2+1]);
73 : }
74 2 : } else if (*datatype == MPIU_NORM1_AND_2) {
75 4 : for (i=0;i<count;i++) {
76 2 : xout[i*3] += xin[i*3];
77 2 : SumNorm2(&xin[i*3+1],&xin[i*3+2],&xout[i*3+1],&xout[i*3+2]);
78 : }
79 : } else {
80 0 : (void)(*PetscErrorPrintf)("Can only handle MPIU_NORM* data types");
81 0 : MPI_Abort(PETSC_COMM_WORLD,1);
82 : }
83 1608 : PetscFunctionReturnVoid();
84 : }
85 :
86 34 : static PetscErrorCode VecCompNormEnd(void)
87 : {
88 34 : PetscFunctionBegin;
89 34 : PetscCallMPI(MPI_Type_free(&MPIU_NORM2));
90 34 : PetscCallMPI(MPI_Type_free(&MPIU_NORM1_AND_2));
91 34 : PetscCallMPI(MPI_Op_free(&MPIU_NORM2_SUM));
92 34 : VecCompInitialized = PETSC_FALSE;
93 34 : PetscFunctionReturn(PETSC_SUCCESS);
94 : }
95 :
96 34 : static PetscErrorCode VecCompNormInit(void)
97 : {
98 34 : PetscFunctionBegin;
99 34 : PetscCallMPI(MPI_Type_contiguous(2,MPIU_REAL,&MPIU_NORM2));
100 34 : PetscCallMPI(MPI_Type_commit(&MPIU_NORM2));
101 34 : PetscCallMPI(MPI_Type_contiguous(3,MPIU_REAL,&MPIU_NORM1_AND_2));
102 34 : PetscCallMPI(MPI_Type_commit(&MPIU_NORM1_AND_2));
103 34 : PetscCallMPI(MPI_Op_create(SlepcSumNorm2_Local,PETSC_TRUE,&MPIU_NORM2_SUM));
104 34 : PetscCall(PetscRegisterFinalize(VecCompNormEnd));
105 34 : PetscFunctionReturn(PETSC_SUCCESS);
106 : }
107 :
108 1885 : PetscErrorCode VecDestroy_Comp(Vec v)
109 : {
110 1885 : Vec_Comp *vs = (Vec_Comp*)v->data;
111 1885 : PetscInt i;
112 :
113 1885 : PetscFunctionBegin;
114 : #if defined(PETSC_USE_LOG)
115 1885 : PetscCall(PetscLogObjectState((PetscObject)v,"Length=%" PetscInt_FMT,v->map->n));
116 : #endif
117 4754 : for (i=0;i<vs->nx;i++) PetscCall(VecDestroy(&vs->x[i]));
118 1885 : if (--vs->n->friends <= 0) PetscCall(PetscFree(vs->n));
119 1885 : PetscCall(PetscFree(vs->x));
120 1885 : PetscCall(PetscFree(vs));
121 1885 : PetscCall(PetscObjectComposeFunction((PetscObject)v,"VecCompSetSubVecs_C",NULL));
122 1885 : PetscCall(PetscObjectComposeFunction((PetscObject)v,"VecCompGetSubVecs_C",NULL));
123 1885 : PetscFunctionReturn(PETSC_SUCCESS);
124 : }
125 :
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 : };
200 :
201 38 : PetscErrorCode VecDuplicateVecs_Comp(Vec w,PetscInt m,Vec *V[])
202 : {
203 38 : PetscInt i;
204 :
205 38 : PetscFunctionBegin;
206 38 : PetscValidHeaderSpecific(w,VEC_CLASSID,1);
207 38 : PetscAssertPointer(V,3);
208 38 : PetscCheck(m>0,PetscObjectComm((PetscObject)w),PETSC_ERR_ARG_OUTOFRANGE,"m must be > 0: m = %" PetscInt_FMT,m);
209 38 : PetscCall(PetscMalloc1(m,V));
210 404 : for (i=0;i<m;i++) PetscCall(VecDuplicate(w,*V+i));
211 38 : PetscFunctionReturn(PETSC_SUCCESS);
212 : }
213 :
214 38 : PetscErrorCode VecDestroyVecs_Comp(PetscInt m,Vec v[])
215 : {
216 38 : PetscInt i;
217 :
218 38 : PetscFunctionBegin;
219 38 : PetscAssertPointer(v,2);
220 38 : PetscCheck(m>0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"m must be > 0: m = %" PetscInt_FMT,m);
221 404 : for (i=0;i<m;i++) PetscCall(VecDestroy(&v[i]));
222 38 : PetscCall(PetscFree(v));
223 38 : PetscFunctionReturn(PETSC_SUCCESS);
224 : }
225 :
226 1885 : static PetscErrorCode VecCreate_Comp_Private(Vec v,Vec *x,PetscInt nx,PetscBool x_to_me,Vec_Comp_N *n)
227 : {
228 1885 : Vec_Comp *s;
229 1885 : PetscInt N=0,lN=0,i,k;
230 :
231 1885 : PetscFunctionBegin;
232 1885 : if (!VecCompInitialized) {
233 34 : VecCompInitialized = PETSC_TRUE;
234 34 : PetscCall(VecRegister(VECCOMP,VecCreate_Comp));
235 34 : PetscCall(VecCompNormInit());
236 : }
237 :
238 : /* Allocate a new Vec_Comp */
239 1885 : if (v->data) PetscCall(PetscFree(v->data));
240 1885 : PetscCall(PetscNew(&s));
241 1885 : PetscCall(PetscMemcpy(v->ops,&DvOps,sizeof(DvOps)));
242 1885 : v->data = (void*)s;
243 1885 : v->petscnative = PETSC_FALSE;
244 :
245 : /* Allocate the array of Vec, if it is needed to be done */
246 1885 : if (!x_to_me) {
247 1513 : if (nx) PetscCall(PetscMalloc1(nx,&s->x));
248 1513 : if (x) PetscCall(PetscArraycpy(s->x,x,nx));
249 372 : } else s->x = x;
250 :
251 1885 : s->nx = nx;
252 :
253 1885 : if (nx && x) {
254 : /* Allocate the shared structure, if it is not given */
255 1882 : if (!n) {
256 121 : for (i=0;i<nx;i++) {
257 73 : PetscCall(VecGetSize(x[i],&k));
258 73 : N+= k;
259 73 : PetscCall(VecGetLocalSize(x[i],&k));
260 73 : lN+= k;
261 : }
262 48 : PetscCall(PetscNew(&n));
263 48 : s->n = n;
264 48 : n->n = nx;
265 48 : n->N = N;
266 48 : n->lN = lN;
267 48 : n->friends = 1;
268 : } else { /* If not, check in the vector in the shared structure */
269 1834 : s->n = n;
270 1834 : s->n->friends++;
271 : }
272 :
273 : /* Set the virtual sizes as the real sizes of the vector */
274 1882 : PetscCall(VecSetSizes(v,s->n->lN,s->n->N));
275 : }
276 :
277 1885 : PetscCall(PetscObjectChangeTypeName((PetscObject)v,VECCOMP));
278 1885 : PetscCall(PetscObjectComposeFunction((PetscObject)v,"VecCompSetSubVecs_C",VecCompSetSubVecs_Comp));
279 1885 : PetscCall(PetscObjectComposeFunction((PetscObject)v,"VecCompGetSubVecs_C",VecCompGetSubVecs_Comp));
280 1885 : PetscFunctionReturn(PETSC_SUCCESS);
281 : }
282 :
283 3 : SLEPC_EXTERN PetscErrorCode VecCreate_Comp(Vec V)
284 : {
285 3 : PetscFunctionBegin;
286 3 : PetscCall(VecCreate_Comp_Private(V,NULL,0,PETSC_FALSE,NULL));
287 3 : PetscFunctionReturn(PETSC_SUCCESS);
288 : }
289 :
290 : /*@
291 : VecCreateComp - Creates a new vector containing several subvectors,
292 : each stored separately.
293 :
294 : Collective
295 :
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
302 :
303 : Output Parameter:
304 : . V - new vector
305 :
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().
309 :
310 : Level: developer
311 :
312 : .seealso: VecCreateCompWithVecs(), VecCompSetSubVecs()
313 : @*/
314 3 : PetscErrorCode VecCreateComp(MPI_Comm comm,PetscInt Nx[],PetscInt n,VecType t,Vec Vparent,Vec *V)
315 : {
316 3 : Vec *x;
317 3 : PetscInt i;
318 :
319 3 : PetscFunctionBegin;
320 3 : PetscCall(VecCreate(comm,V));
321 3 : PetscCall(PetscMalloc1(n,&x));
322 9 : for (i=0;i<n;i++) {
323 6 : PetscCall(VecCreate(comm,&x[i]));
324 6 : PetscCall(VecSetSizes(x[i],PETSC_DECIDE,Nx[i]));
325 6 : PetscCall(VecSetType(x[i],t));
326 : }
327 3 : PetscCall(VecCreate_Comp_Private(*V,x,n,PETSC_TRUE,Vparent?((Vec_Comp*)Vparent->data)->n:NULL));
328 3 : PetscFunctionReturn(PETSC_SUCCESS);
329 : }
330 :
331 : /*@
332 : VecCreateCompWithVecs - Creates a new vector containing several subvectors,
333 : each stored separately, from an array of Vecs.
334 :
335 : Collective
336 :
337 : Input Parameters:
338 : + x - array of Vecs
339 : . n - number of child vectors
340 : - Vparent - (optional) template vector
341 :
342 : Output Parameter:
343 : . V - new vector
344 :
345 : Level: developer
346 :
347 : .seealso: VecCreateComp()
348 : @*/
349 1510 : PetscErrorCode VecCreateCompWithVecs(Vec x[],PetscInt n,Vec Vparent,Vec *V)
350 : {
351 1510 : PetscInt i;
352 :
353 1510 : PetscFunctionBegin;
354 1510 : PetscAssertPointer(x,1);
355 1510 : PetscValidHeaderSpecific(*x,VEC_CLASSID,1);
356 4530 : PetscValidLogicalCollectiveInt(*x,n,2);
357 1510 : PetscCall(VecCreate(PetscObjectComm((PetscObject)x[0]),V));
358 3813 : for (i=0;i<n;i++) PetscCall(PetscObjectReference((PetscObject)x[i]));
359 1510 : PetscCall(VecCreate_Comp_Private(*V,x,n,PETSC_FALSE,Vparent?((Vec_Comp*)Vparent->data)->n:NULL));
360 1510 : PetscFunctionReturn(PETSC_SUCCESS);
361 : }
362 :
363 369 : PetscErrorCode VecDuplicate_Comp(Vec win,Vec *V)
364 : {
365 369 : Vec *x;
366 369 : PetscInt i;
367 369 : Vec_Comp *s = (Vec_Comp*)win->data;
368 :
369 369 : PetscFunctionBegin;
370 369 : SlepcValidVecComp(win,1);
371 369 : PetscCall(VecCreate(PetscObjectComm((PetscObject)win),V));
372 369 : PetscCall(PetscMalloc1(s->nx,&x));
373 923 : for (i=0;i<s->nx;i++) {
374 554 : if (s->x[i]) PetscCall(VecDuplicate(s->x[i],&x[i]));
375 0 : else x[i] = NULL;
376 : }
377 369 : PetscCall(VecCreate_Comp_Private(*V,x,s->nx,PETSC_TRUE,s->n));
378 369 : PetscFunctionReturn(PETSC_SUCCESS);
379 : }
380 :
381 64660 : static PetscErrorCode VecCompGetSubVecs_Comp(Vec win,PetscInt *n,const Vec **x)
382 : {
383 64660 : Vec_Comp *s = (Vec_Comp*)win->data;
384 :
385 64660 : PetscFunctionBegin;
386 64660 : if (x) *x = s->x;
387 64660 : if (n) *n = s->n->n;
388 64660 : PetscFunctionReturn(PETSC_SUCCESS);
389 : }
390 :
391 : /*@C
392 : VecCompGetSubVecs - Returns the entire array of vectors defining a
393 : compound vector.
394 :
395 : Collective
396 :
397 : Input Parameter:
398 : . win - compound vector
399 :
400 : Output Parameters:
401 : + n - number of child vectors
402 : - x - array of child vectors
403 :
404 : Level: developer
405 :
406 : .seealso: VecCreateComp()
407 : @*/
408 64660 : PetscErrorCode VecCompGetSubVecs(Vec win,PetscInt *n,const Vec **x)
409 : {
410 64660 : PetscFunctionBegin;
411 64660 : PetscValidHeaderSpecific(win,VEC_CLASSID,1);
412 64660 : PetscUseMethod(win,"VecCompGetSubVecs_C",(Vec,PetscInt*,const Vec**),(win,n,x));
413 64660 : PetscFunctionReturn(PETSC_SUCCESS);
414 : }
415 :
416 718 : static PetscErrorCode VecCompSetSubVecs_Comp(Vec win,PetscInt n,Vec *x)
417 : {
418 718 : Vec_Comp *s = (Vec_Comp*)win->data;
419 718 : PetscInt i,N,nlocal;
420 718 : Vec_Comp_N *nn;
421 :
422 718 : PetscFunctionBegin;
423 718 : PetscCheck(s,PetscObjectComm((PetscObject)win),PETSC_ERR_ORDER,"Must call VecSetSizes first");
424 718 : if (!s->nx) {
425 : /* vector has been created via VecCreate+VecSetType+VecSetSizes, so allocate data structures */
426 3 : PetscCall(PetscMalloc1(n,&s->x));
427 3 : PetscCall(VecGetSize(win,&N));
428 3 : PetscCheck(N%n==0,PetscObjectComm((PetscObject)win),PETSC_ERR_SUP,"Global dimension %" PetscInt_FMT " is not divisible by %" PetscInt_FMT,N,n);
429 3 : PetscCall(VecGetLocalSize(win,&nlocal));
430 3 : PetscCheck(nlocal%n==0,PetscObjectComm((PetscObject)win),PETSC_ERR_SUP,"Local dimension %" PetscInt_FMT " is not divisible by %" PetscInt_FMT,nlocal,n);
431 3 : s->nx = n;
432 9 : for (i=0;i<n;i++) {
433 6 : PetscCall(VecCreate(PetscObjectComm((PetscObject)win),&s->x[i]));
434 6 : PetscCall(VecSetSizes(s->x[i],nlocal/n,N/n));
435 6 : PetscCall(VecSetFromOptions(s->x[i]));
436 : }
437 3 : if (!s->n) {
438 3 : PetscCall(PetscNew(&nn));
439 3 : s->n = nn;
440 3 : nn->N = N;
441 3 : nn->lN = nlocal;
442 3 : nn->friends = 1;
443 : }
444 715 : } else PetscCheck(n<=s->nx,PetscObjectComm((PetscObject)win),PETSC_ERR_SUP,"Number of child vectors cannot be larger than %" PetscInt_FMT,s->nx);
445 718 : if (x) PetscCall(PetscArraycpy(s->x,x,n));
446 718 : s->n->n = n;
447 718 : PetscFunctionReturn(PETSC_SUCCESS);
448 : }
449 :
450 : /*@
451 : VecCompSetSubVecs - Resets the number of subvectors defining a compound vector,
452 : or replaces the subvectors.
453 :
454 : Collective
455 :
456 : Input Parameters:
457 : + win - compound vector
458 : . n - number of child vectors
459 : - x - array of child vectors
460 :
461 : Note:
462 : It is not possible to increase the number of subvectors with respect to the
463 : number set at its creation.
464 :
465 : Level: developer
466 :
467 : .seealso: VecCreateComp(), VecCompGetSubVecs()
468 : @*/
469 718 : PetscErrorCode VecCompSetSubVecs(Vec win,PetscInt n,Vec x[])
470 : {
471 718 : PetscFunctionBegin;
472 718 : PetscValidHeaderSpecific(win,VEC_CLASSID,1);
473 2154 : PetscValidLogicalCollectiveInt(win,n,2);
474 718 : PetscTryMethod(win,"VecCompSetSubVecs_C",(Vec,PetscInt,Vec*),(win,n,x));
475 718 : PetscFunctionReturn(PETSC_SUCCESS);
476 : }
477 :
478 27046 : PetscErrorCode VecAXPY_Comp(Vec v,PetscScalar alpha,Vec w)
479 : {
480 27046 : Vec_Comp *vs = (Vec_Comp*)v->data,*ws = (Vec_Comp*)w->data;
481 27046 : PetscInt i;
482 :
483 27046 : PetscFunctionBegin;
484 27046 : SlepcValidVecComp(v,1);
485 27046 : SlepcValidVecComp(w,3);
486 54559 : for (i=0;i<vs->n->n;i++) PetscCall(VecAXPY(vs->x[i],alpha,ws->x[i]));
487 27046 : PetscFunctionReturn(PETSC_SUCCESS);
488 : }
489 :
490 15542 : PetscErrorCode VecAYPX_Comp(Vec v,PetscScalar alpha,Vec w)
491 : {
492 15542 : Vec_Comp *vs = (Vec_Comp*)v->data,*ws = (Vec_Comp*)w->data;
493 15542 : PetscInt i;
494 :
495 15542 : PetscFunctionBegin;
496 15542 : SlepcValidVecComp(v,1);
497 15542 : SlepcValidVecComp(w,3);
498 31339 : for (i=0;i<vs->n->n;i++) PetscCall(VecAYPX(vs->x[i],alpha,ws->x[i]));
499 15542 : PetscFunctionReturn(PETSC_SUCCESS);
500 : }
501 :
502 3 : PetscErrorCode VecAXPBY_Comp(Vec v,PetscScalar alpha,PetscScalar beta,Vec w)
503 : {
504 3 : Vec_Comp *vs = (Vec_Comp*)v->data,*ws = (Vec_Comp*)w->data;
505 3 : PetscInt i;
506 :
507 3 : PetscFunctionBegin;
508 3 : SlepcValidVecComp(v,1);
509 3 : SlepcValidVecComp(w,4);
510 9 : for (i=0;i<vs->n->n;i++) PetscCall(VecAXPBY(vs->x[i],alpha,beta,ws->x[i]));
511 3 : PetscFunctionReturn(PETSC_SUCCESS);
512 : }
513 :
514 16734 : PetscErrorCode VecMAXPY_Comp(Vec v,PetscInt n,const PetscScalar *alpha,Vec *w)
515 : {
516 16734 : Vec_Comp *vs = (Vec_Comp*)v->data;
517 16734 : Vec *wx;
518 16734 : PetscInt i,j;
519 :
520 16734 : PetscFunctionBegin;
521 16734 : SlepcValidVecComp(v,1);
522 64504 : for (i=0;i<n;i++) SlepcValidVecComp(w[i],4);
523 :
524 16734 : PetscCall(PetscMalloc1(n,&wx));
525 :
526 33723 : for (j=0;j<vs->n->n;j++) {
527 65269 : for (i=0;i<n;i++) wx[i] = ((Vec_Comp*)w[i]->data)->x[j];
528 16989 : PetscCall(VecMAXPY(vs->x[j],n,alpha,wx));
529 : }
530 :
531 16734 : PetscCall(PetscFree(wx));
532 16734 : PetscFunctionReturn(PETSC_SUCCESS);
533 : }
534 :
535 177 : PetscErrorCode VecWAXPY_Comp(Vec v,PetscScalar alpha,Vec w,Vec z)
536 : {
537 177 : Vec_Comp *vs = (Vec_Comp*)v->data,*ws = (Vec_Comp*)w->data,*zs = (Vec_Comp*)z->data;
538 177 : PetscInt i;
539 :
540 177 : PetscFunctionBegin;
541 177 : SlepcValidVecComp(v,1);
542 177 : SlepcValidVecComp(w,3);
543 177 : SlepcValidVecComp(z,4);
544 357 : for (i=0;i<vs->n->n;i++) PetscCall(VecWAXPY(vs->x[i],alpha,ws->x[i],zs->x[i]));
545 177 : PetscFunctionReturn(PETSC_SUCCESS);
546 : }
547 :
548 177 : PetscErrorCode VecAXPBYPCZ_Comp(Vec v,PetscScalar alpha,PetscScalar beta,PetscScalar gamma,Vec w,Vec z)
549 : {
550 177 : Vec_Comp *vs = (Vec_Comp*)v->data,*ws = (Vec_Comp*)w->data,*zs = (Vec_Comp*)z->data;
551 177 : PetscInt i;
552 :
553 177 : PetscFunctionBegin;
554 177 : SlepcValidVecComp(v,1);
555 177 : SlepcValidVecComp(w,5);
556 177 : SlepcValidVecComp(z,6);
557 357 : for (i=0;i<vs->n->n;i++) PetscCall(VecAXPBYPCZ(vs->x[i],alpha,beta,gamma,ws->x[i],zs->x[i]));
558 177 : PetscFunctionReturn(PETSC_SUCCESS);
559 : }
560 :
561 9 : PetscErrorCode VecGetSize_Comp(Vec v,PetscInt *size)
562 : {
563 9 : Vec_Comp *vs = (Vec_Comp*)v->data;
564 :
565 9 : PetscFunctionBegin;
566 9 : PetscAssertPointer(size,2);
567 9 : if (vs->n) {
568 6 : SlepcValidVecComp(v,1);
569 6 : *size = vs->n->N;
570 3 : } else *size = v->map->N;
571 9 : PetscFunctionReturn(PETSC_SUCCESS);
572 : }
573 :
574 28353 : PetscErrorCode VecGetLocalSize_Comp(Vec v,PetscInt *size)
575 : {
576 28353 : Vec_Comp *vs = (Vec_Comp*)v->data;
577 :
578 28353 : PetscFunctionBegin;
579 28353 : PetscAssertPointer(size,2);
580 28353 : if (vs->n) {
581 28350 : SlepcValidVecComp(v,1);
582 28350 : *size = vs->n->lN;
583 3 : } else *size = v->map->n;
584 28353 : PetscFunctionReturn(PETSC_SUCCESS);
585 : }
586 :
587 3 : PetscErrorCode VecMax_Comp(Vec v,PetscInt *idx,PetscReal *z)
588 : {
589 3 : Vec_Comp *vs = (Vec_Comp*)v->data;
590 3 : PetscInt idxp,s=0,s0;
591 3 : PetscReal zp,z0;
592 3 : PetscInt i;
593 :
594 3 : PetscFunctionBegin;
595 3 : SlepcValidVecComp(v,1);
596 3 : if (!idx && !z) PetscFunctionReturn(PETSC_SUCCESS);
597 :
598 6 : if (vs->n->n > 0) PetscCall(VecMax(vs->x[0],idx?&idxp:NULL,&zp));
599 : else {
600 0 : zp = PETSC_MIN_REAL;
601 0 : if (idx) idxp = -1;
602 : }
603 6 : for (i=1;i<vs->n->n;i++) {
604 3 : PetscCall(VecGetSize(vs->x[i-1],&s0));
605 3 : s += s0;
606 6 : PetscCall(VecMax(vs->x[i],idx?&idxp:NULL,&z0));
607 3 : if (zp < z0) {
608 0 : if (idx) *idx = s+idxp;
609 0 : zp = z0;
610 : }
611 : }
612 3 : if (z) *z = zp;
613 3 : PetscFunctionReturn(PETSC_SUCCESS);
614 : }
615 :
616 3 : PetscErrorCode VecMin_Comp(Vec v,PetscInt *idx,PetscReal *z)
617 : {
618 3 : Vec_Comp *vs = (Vec_Comp*)v->data;
619 3 : PetscInt idxp,s=0,s0;
620 3 : PetscReal zp,z0;
621 3 : PetscInt i;
622 :
623 3 : PetscFunctionBegin;
624 3 : SlepcValidVecComp(v,1);
625 3 : if (!idx && !z) PetscFunctionReturn(PETSC_SUCCESS);
626 :
627 6 : if (vs->n->n > 0) PetscCall(VecMin(vs->x[0],idx?&idxp:NULL,&zp));
628 : else {
629 0 : zp = PETSC_MAX_REAL;
630 0 : if (idx) idxp = -1;
631 : }
632 6 : for (i=1;i<vs->n->n;i++) {
633 3 : PetscCall(VecGetSize(vs->x[i-1],&s0));
634 3 : s += s0;
635 6 : PetscCall(VecMin(vs->x[i],idx?&idxp:NULL,&z0));
636 3 : if (zp > z0) {
637 3 : if (idx) *idx = s+idxp;
638 3 : zp = z0;
639 : }
640 : }
641 3 : if (z) *z = zp;
642 3 : PetscFunctionReturn(PETSC_SUCCESS);
643 : }
644 :
645 3 : PetscErrorCode VecMaxPointwiseDivide_Comp(Vec v,Vec w,PetscReal *m)
646 : {
647 3 : Vec_Comp *vs = (Vec_Comp*)v->data,*ws = (Vec_Comp*)w->data;
648 3 : PetscReal work;
649 3 : PetscInt i;
650 :
651 3 : PetscFunctionBegin;
652 3 : SlepcValidVecComp(v,1);
653 3 : SlepcValidVecComp(w,2);
654 3 : if (!m || vs->n->n == 0) PetscFunctionReturn(PETSC_SUCCESS);
655 3 : PetscCall(VecMaxPointwiseDivide(vs->x[0],ws->x[0],m));
656 6 : for (i=1;i<vs->n->n;i++) {
657 3 : PetscCall(VecMaxPointwiseDivide(vs->x[i],ws->x[i],&work));
658 4 : *m = PetscMax(*m,work);
659 : }
660 3 : PetscFunctionReturn(PETSC_SUCCESS);
661 : }
662 :
663 : #define __QUOTEME__(x) #x
664 : #define __COMPOSE2__(A,B) A##B
665 : #define __COMPOSE3__(A,B,C) A##B##C
666 :
667 : #define __FUNC_TEMPLATE1__(NAME) \
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 : }
680 :
681 0 : __FUNC_TEMPLATE1__(Conjugate)
682 9 : __FUNC_TEMPLATE1__(Reciprocal)
683 9 : __FUNC_TEMPLATE1__(SqrtAbs)
684 9 : __FUNC_TEMPLATE1__(Abs)
685 9 : __FUNC_TEMPLATE1__(Exp)
686 9 : __FUNC_TEMPLATE1__(Log)
687 :
688 : #define __FUNC_TEMPLATE2__(NAME,T0) \
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 : }
701 :
702 2884 : __FUNC_TEMPLATE2__(Set,PetscScalar)
703 9 : __FUNC_TEMPLATE2__(View,PetscViewer)
704 2151 : __FUNC_TEMPLATE2__(Scale,PetscScalar)
705 9 : __FUNC_TEMPLATE2__(SetRandom,PetscRandom)
706 9 : __FUNC_TEMPLATE2__(Shift,PetscScalar)
707 :
708 : #define __FUNC_TEMPLATE3__(NAME) \
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 : }
723 :
724 24294 : __FUNC_TEMPLATE3__(Copy)
725 9 : __FUNC_TEMPLATE3__(Swap)
726 :
727 : #define __FUNC_TEMPLATE4__(NAME) \
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 : }
744 :
745 9 : __FUNC_TEMPLATE4__(PointwiseMax)
746 9 : __FUNC_TEMPLATE4__(PointwiseMaxAbs)
747 9 : __FUNC_TEMPLATE4__(PointwiseMin)
748 9 : __FUNC_TEMPLATE4__(PointwiseMult)
749 9 : __FUNC_TEMPLATE4__(PointwiseDivide)
|
