Actual source code: test1.c
slepc-3.21.2 2024-09-25
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: static char help[] = "Test VecComp.\n\n";
13: #include <slepcsys.h>
15: int main(int argc,char **argv)
16: {
17: Vec v,w,x,y,vc,wc,xc,yc,vparent,vchild[2],vecs[2];
18: const Vec *varray;
19: PetscMPIInt size,rank;
20: PetscInt i,n,k,Nx[2];
21: PetscReal norm,normc,norm12[2],norm12c[2],vmax,vmin;
22: PetscScalar dot[2],dotc[2];
24: PetscFunctionBeginUser;
25: PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
26: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size));
27: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank));
28: PetscCheck(size<=2,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This test needs one or two processes");
29: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"VecComp test\n"));
31: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
32: Create standard vectors
33: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
35: PetscCall(VecCreate(PETSC_COMM_WORLD,&v));
36: PetscCall(VecSetSizes(v,8/size,8));
37: PetscCall(VecSetFromOptions(v));
39: if (!rank) {
40: PetscCall(VecSetValue(v,0,2.0,INSERT_VALUES));
41: PetscCall(VecSetValue(v,1,-1.0,INSERT_VALUES));
42: PetscCall(VecSetValue(v,2,3.0,INSERT_VALUES));
43: PetscCall(VecSetValue(v,3,3.5,INSERT_VALUES));
44: }
45: if ((!rank && size==1) || (rank && size==2)) {
46: PetscCall(VecSetValue(v,4,1.2,INSERT_VALUES));
47: PetscCall(VecSetValue(v,5,1.8,INSERT_VALUES));
48: PetscCall(VecSetValue(v,6,-2.2,INSERT_VALUES));
49: PetscCall(VecSetValue(v,7,2.0,INSERT_VALUES));
50: }
51: PetscCall(VecAssemblyBegin(v));
52: PetscCall(VecAssemblyEnd(v));
53: PetscCall(VecDuplicate(v,&w));
54: PetscCall(VecSet(w,1.0));
55: PetscCall(VecDuplicate(v,&x));
56: PetscCall(VecDuplicate(v,&y));
57: if (!rank) PetscCall(VecSetValue(y,0,1.0,INSERT_VALUES));
58: PetscCall(VecAssemblyBegin(y));
59: PetscCall(VecAssemblyEnd(y));
61: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
62: Create veccomp vectors
63: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
65: PetscCall(VecCreate(PETSC_COMM_WORLD,&vparent));
66: PetscCall(VecSetSizes(vparent,4/size,4));
67: PetscCall(VecSetFromOptions(vparent));
69: /* create a veccomp vector with two subvectors */
70: PetscCall(VecDuplicate(vparent,&vchild[0]));
71: PetscCall(VecDuplicate(vparent,&vchild[1]));
72: if (!rank) {
73: PetscCall(VecSetValue(vchild[0],0,2.0,INSERT_VALUES));
74: PetscCall(VecSetValue(vchild[0],1,-1.0,INSERT_VALUES));
75: PetscCall(VecSetValue(vchild[1],0,1.2,INSERT_VALUES));
76: PetscCall(VecSetValue(vchild[1],1,1.8,INSERT_VALUES));
77: }
78: if ((!rank && size==1) || (rank && size==2)) {
79: PetscCall(VecSetValue(vchild[0],2,3.0,INSERT_VALUES));
80: PetscCall(VecSetValue(vchild[0],3,3.5,INSERT_VALUES));
81: PetscCall(VecSetValue(vchild[1],2,-2.2,INSERT_VALUES));
82: PetscCall(VecSetValue(vchild[1],3,2.0,INSERT_VALUES));
83: }
84: PetscCall(VecAssemblyBegin(vchild[0]));
85: PetscCall(VecAssemblyBegin(vchild[1]));
86: PetscCall(VecAssemblyEnd(vchild[0]));
87: PetscCall(VecAssemblyEnd(vchild[1]));
88: PetscCall(VecCreateCompWithVecs(vchild,2,vparent,&vc));
89: PetscCall(VecDestroy(&vchild[0]));
90: PetscCall(VecDestroy(&vchild[1]));
91: PetscCall(VecView(vc,NULL));
93: PetscCall(VecGetSize(vc,&k));
94: PetscCheck(k==8,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Vector global length should be 8");
96: /* create an empty veccomp vector with two subvectors */
97: Nx[0] = 4;
98: Nx[1] = 4;
99: PetscCall(VecCreateComp(PETSC_COMM_WORLD,Nx,2,VECSTANDARD,vparent,&wc));
100: PetscCall(VecCompGetSubVecs(wc,&n,&varray));
101: PetscCheck(n==2,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"n should be 2");
102: for (i=0;i<2;i++) PetscCall(VecSet(varray[i],1.0));
104: PetscCall(VecGetSize(wc,&k));
105: PetscCheck(k==8,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Vector global length should be 8");
107: /* duplicate a veccomp */
108: PetscCall(VecDuplicate(vc,&xc));
110: /* create a veccomp via VecSetType */
111: PetscCall(VecCreate(PETSC_COMM_WORLD,&yc));
112: PetscCall(VecSetType(yc,VECCOMP));
113: PetscCall(VecSetSizes(yc,8/size,8));
114: PetscCall(VecCompSetSubVecs(yc,2,NULL));
116: PetscCall(VecCompGetSubVecs(yc,&n,&varray));
117: PetscCheck(n==2,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"n should be 2");
118: if (!rank) PetscCall(VecSetValue(varray[0],0,1.0,INSERT_VALUES));
119: PetscCall(VecAssemblyBegin(varray[0]));
120: PetscCall(VecAssemblyEnd(varray[0]));
122: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
123: Operate with vectors
124: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
126: PetscCall(VecCopy(w,x));
127: PetscCall(VecAXPBY(x,1.0,-2.0,v));
128: PetscCall(VecNorm(x,NORM_2,&norm));
129: PetscCall(VecCopy(wc,xc));
130: PetscCall(VecAXPBY(xc,1.0,-2.0,vc));
131: PetscCall(VecNorm(xc,NORM_2,&normc));
132: PetscCheck(PetscAbsReal(norm-normc)<10*PETSC_MACHINE_EPSILON,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Norms are different");
134: PetscCall(VecCopy(w,x));
135: PetscCall(VecWAXPY(x,-2.0,w,v));
136: PetscCall(VecNorm(x,NORM_2,&norm));
137: PetscCall(VecCopy(wc,xc));
138: PetscCall(VecWAXPY(xc,-2.0,wc,vc));
139: PetscCall(VecNorm(xc,NORM_2,&normc));
140: PetscCheck(PetscAbsReal(norm-normc)<10*PETSC_MACHINE_EPSILON,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Norms are different");
142: PetscCall(VecAXPBYPCZ(y,3.0,-1.0,1.0,w,v));
143: PetscCall(VecNorm(y,NORM_2,&norm));
144: PetscCall(VecAXPBYPCZ(yc,3.0,-1.0,1.0,wc,vc));
145: PetscCall(VecNorm(yc,NORM_2,&normc));
146: PetscCheck(PetscAbsReal(norm-normc)<10*PETSC_MACHINE_EPSILON,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Norms are different");
148: PetscCall(VecMax(xc,NULL,&vmax));
149: PetscCall(VecMin(xc,NULL,&vmin));
150: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"xc has max value %g min value %g\n",(double)vmax,(double)vmin));
152: PetscCall(VecMaxPointwiseDivide(wc,xc,&vmax));
153: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"wc/xc has max value %g\n",(double)vmax));
155: PetscCall(VecDot(x,y,&dot[0]));
156: PetscCall(VecDot(xc,yc,&dotc[0]));
157: PetscCheck(PetscAbsScalar(dot[0]-dotc[0])<10*PETSC_MACHINE_EPSILON,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Dots are different");
158: PetscCall(VecTDot(x,y,&dot[0]));
159: PetscCall(VecTDot(xc,yc,&dotc[0]));
160: PetscCheck(PetscAbsScalar(dot[0]-dotc[0])<10*PETSC_MACHINE_EPSILON,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Dots are different");
162: vecs[0] = w; vecs[1] = y;
163: PetscCall(VecMDot(x,2,vecs,dot));
164: vecs[0] = wc; vecs[1] = yc;
165: PetscCall(VecMDot(xc,2,vecs,dotc));
166: PetscCheck(PetscAbsScalar(dot[0]-dotc[0])<10*PETSC_MACHINE_EPSILON || PetscAbsScalar(dot[1]-dotc[1])>10*PETSC_MACHINE_EPSILON,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Dots are different");
167: vecs[0] = w; vecs[1] = y;
168: PetscCall(VecMTDot(x,2,vecs,dot));
169: vecs[0] = wc; vecs[1] = yc;
170: PetscCall(VecMTDot(xc,2,vecs,dotc));
171: PetscCheck(PetscAbsScalar(dot[0]-dotc[0])<10*PETSC_MACHINE_EPSILON || PetscAbsScalar(dot[1]-dotc[1])>10*PETSC_MACHINE_EPSILON,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Dots are different");
173: PetscCall(VecDotNorm2(x,y,&dot[0],&norm));
174: PetscCall(VecDotNorm2(xc,yc,&dotc[0],&normc));
175: PetscCheck(PetscAbsScalar(dot[0]-dotc[0])<10*PETSC_MACHINE_EPSILON,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Dots are different");
176: PetscCheck(PetscAbsReal(norm-normc)<100*PETSC_MACHINE_EPSILON,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Norms are different");
178: PetscCall(VecAbs(w));
179: PetscCall(VecAbs(wc));
180: PetscCall(VecConjugate(x));
181: PetscCall(VecConjugate(xc));
182: PetscCall(VecShift(y,0.5));
183: PetscCall(VecShift(yc,0.5));
184: PetscCall(VecReciprocal(y));
185: PetscCall(VecReciprocal(yc));
186: PetscCall(VecExp(y));
187: PetscCall(VecExp(yc));
188: PetscCall(VecLog(y));
189: PetscCall(VecLog(yc));
190: PetscCall(VecNorm(y,NORM_1,&norm));
191: PetscCall(VecNorm(yc,NORM_1,&normc));
192: PetscCheck(PetscAbsReal(norm-normc)<10*PETSC_MACHINE_EPSILON,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Norms are different");
194: PetscCall(VecPointwiseMult(w,x,y));
195: PetscCall(VecPointwiseMult(wc,xc,yc));
196: PetscCall(VecNorm(w,NORM_INFINITY,&norm));
197: PetscCall(VecNorm(wc,NORM_INFINITY,&normc));
198: PetscCheck(PetscAbsReal(norm-normc)<10*PETSC_MACHINE_EPSILON,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Norms are different");
200: PetscCall(VecPointwiseMax(w,x,y));
201: PetscCall(VecPointwiseMax(wc,xc,yc));
202: PetscCall(VecNorm(w,NORM_INFINITY,&norm));
203: PetscCall(VecNorm(wc,NORM_INFINITY,&normc));
204: PetscCheck(PetscAbsReal(norm-normc)<10*PETSC_MACHINE_EPSILON,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Norms are different");
206: PetscCall(VecSwap(x,y));
207: PetscCall(VecSwap(xc,yc));
208: PetscCall(VecPointwiseDivide(w,x,y));
209: PetscCall(VecPointwiseDivide(wc,xc,yc));
210: PetscCall(VecScale(w,0.3));
211: PetscCall(VecScale(wc,0.3));
212: PetscCall(VecSqrtAbs(w));
213: PetscCall(VecSqrtAbs(wc));
214: PetscCall(VecNorm(w,NORM_1_AND_2,norm12));
215: PetscCall(VecNorm(wc,NORM_1_AND_2,norm12c));
216: PetscCheck(PetscAbsReal(norm12[0]-norm12c[0])<10*PETSC_MACHINE_EPSILON || PetscAbsReal(norm12[1]-norm12c[1])>10*PETSC_MACHINE_EPSILON,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Norms are different");
218: PetscCall(VecPointwiseMin(w,x,y));
219: PetscCall(VecPointwiseMin(wc,xc,yc));
220: PetscCall(VecPointwiseMaxAbs(x,y,w));
221: PetscCall(VecPointwiseMaxAbs(xc,yc,wc));
222: PetscCall(VecNorm(x,NORM_INFINITY,&norm));
223: PetscCall(VecNorm(xc,NORM_INFINITY,&normc));
224: PetscCheck(PetscAbsReal(norm-normc)<10*PETSC_MACHINE_EPSILON,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Norms are different");
226: PetscCall(VecSetRandom(wc,NULL));
228: PetscCall(VecDestroy(&v));
229: PetscCall(VecDestroy(&w));
230: PetscCall(VecDestroy(&x));
231: PetscCall(VecDestroy(&y));
232: PetscCall(VecDestroy(&vparent));
233: PetscCall(VecDestroy(&vc));
234: PetscCall(VecDestroy(&wc));
235: PetscCall(VecDestroy(&xc));
236: PetscCall(VecDestroy(&yc));
237: PetscCall(SlepcFinalize());
238: return 0;
239: }
241: /*TEST
243: test:
244: suffix: 1
246: test:
247: suffix: 2
248: nsize: 2
250: TEST*/