Actual source code: test1.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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];
25: SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
26: MPI_Comm_size(PETSC_COMM_WORLD,&size);
27: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
28: if (size > 2) SETERRQ(PETSC_COMM_WORLD,1,"This test needs one or two processes");
29: PetscPrintf(PETSC_COMM_WORLD,"VecComp test\n");
31: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
32: Create standard vectors
33: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
35: VecCreate(PETSC_COMM_WORLD,&v);
36: VecSetSizes(v,8/size,8);
37: VecSetFromOptions(v);
39: if (!rank) {
40: VecSetValue(v,0,2.0,INSERT_VALUES);
41: VecSetValue(v,1,-1.0,INSERT_VALUES);
42: VecSetValue(v,2,3.0,INSERT_VALUES);
43: VecSetValue(v,3,3.5,INSERT_VALUES);
44: }
45: if ((!rank && size==1) || (rank && size==2)) {
46: VecSetValue(v,4,1.2,INSERT_VALUES);
47: VecSetValue(v,5,1.8,INSERT_VALUES);
48: VecSetValue(v,6,-2.2,INSERT_VALUES);
49: VecSetValue(v,7,2.0,INSERT_VALUES);
50: }
51: VecAssemblyBegin(v);
52: VecAssemblyEnd(v);
53: VecDuplicate(v,&w);
54: VecSet(w,1.0);
55: VecDuplicate(v,&x);
56: VecDuplicate(v,&y);
57: if (!rank) {
58: VecSetValue(y,0,1.0,INSERT_VALUES);
59: }
60: VecAssemblyBegin(y);
61: VecAssemblyEnd(y);
63: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
64: Create veccomp vectors
65: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
67: VecCreate(PETSC_COMM_WORLD,&vparent);
68: VecSetSizes(vparent,4/size,4);
69: VecSetFromOptions(vparent);
71: /* create a veccomp vector with two subvectors */
72: VecDuplicate(vparent,&vchild[0]);
73: VecDuplicate(vparent,&vchild[1]);
74: if (!rank) {
75: VecSetValue(vchild[0],0,2.0,INSERT_VALUES);
76: VecSetValue(vchild[0],1,-1.0,INSERT_VALUES);
77: VecSetValue(vchild[1],0,1.2,INSERT_VALUES);
78: VecSetValue(vchild[1],1,1.8,INSERT_VALUES);
79: }
80: if ((!rank && size==1) || (rank && size==2)) {
81: VecSetValue(vchild[0],2,3.0,INSERT_VALUES);
82: VecSetValue(vchild[0],3,3.5,INSERT_VALUES);
83: VecSetValue(vchild[1],2,-2.2,INSERT_VALUES);
84: VecSetValue(vchild[1],3,2.0,INSERT_VALUES);
85: }
86: VecAssemblyBegin(vchild[0]);
87: VecAssemblyBegin(vchild[1]);
88: VecAssemblyEnd(vchild[0]);
89: VecAssemblyEnd(vchild[1]);
90: VecCreateCompWithVecs(vchild,2,vparent,&vc);
91: VecDestroy(&vchild[0]);
92: VecDestroy(&vchild[1]);
94: VecGetSize(vc,&k);
95: if (k!=8) SETERRQ(PETSC_COMM_WORLD,1,"Vector global length should be 8");
97: /* create an empty veccomp vector with two subvectors */
98: Nx[0] = 4;
99: Nx[1] = 4;
100: VecCreateComp(PETSC_COMM_WORLD,Nx,2,VECSTANDARD,vparent,&wc);
101: VecCompGetSubVecs(wc,&n,&varray);
102: if (n!=2) SETERRQ(PETSC_COMM_WORLD,1,"n should be 2");
103: for (i=0;i<2;i++) {
104: VecSet(varray[i],1.0);
105: }
107: VecGetSize(wc,&k);
108: if (k!=8) SETERRQ(PETSC_COMM_WORLD,1,"Vector global length should be 8");
110: /* duplicate a veccomp */
111: VecDuplicate(vc,&xc);
113: /* create a veccomp via VecSetType */
114: VecCreate(PETSC_COMM_WORLD,&yc);
115: VecSetType(yc,VECCOMP);
116: VecSetSizes(yc,8/size,8);
117: VecCompSetSubVecs(yc,2,NULL);
119: VecCompGetSubVecs(yc,&n,&varray);
120: if (n!=2) SETERRQ(PETSC_COMM_WORLD,1,"n should be 2");
121: if (!rank) {
122: VecSetValue(varray[0],0,1.0,INSERT_VALUES);
123: }
124: VecAssemblyBegin(varray[0]);
125: VecAssemblyEnd(varray[0]);
127: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128: Operate with vectors
129: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
131: VecCopy(w,x);
132: VecAXPBY(x,1.0,-2.0,v);
133: VecNorm(x,NORM_2,&norm);
134: VecCopy(wc,xc);
135: VecAXPBY(xc,1.0,-2.0,vc);
136: VecNorm(xc,NORM_2,&normc);
137: if (PetscAbsReal(norm-normc)>10*PETSC_MACHINE_EPSILON) SETERRQ(PETSC_COMM_WORLD,1,"Norms are different");
139: VecCopy(w,x);
140: VecWAXPY(x,-2.0,w,v);
141: VecNorm(x,NORM_2,&norm);
142: VecCopy(wc,xc);
143: VecWAXPY(xc,-2.0,wc,vc);
144: VecNorm(xc,NORM_2,&normc);
145: if (PetscAbsReal(norm-normc)>10*PETSC_MACHINE_EPSILON) SETERRQ(PETSC_COMM_WORLD,1,"Norms are different");
147: VecAXPBYPCZ(y,3.0,-1.0,1.0,w,v);
148: VecNorm(y,NORM_2,&norm);
149: VecAXPBYPCZ(yc,3.0,-1.0,1.0,wc,vc);
150: VecNorm(yc,NORM_2,&normc);
151: if (PetscAbsReal(norm-normc)>10*PETSC_MACHINE_EPSILON) SETERRQ(PETSC_COMM_WORLD,1,"Norms are different");
153: VecMax(xc,NULL,&vmax);
154: VecMin(xc,NULL,&vmin);
155: PetscPrintf(PETSC_COMM_WORLD,"xc has max value %g min value %g\n",(double)vmax,(double)vmin);
157: VecMaxPointwiseDivide(wc,xc,&vmax);
158: PetscPrintf(PETSC_COMM_WORLD,"wc/xc has max value %g\n",(double)vmax);
160: VecDot(x,y,&dot[0]);
161: VecDot(xc,yc,&dotc[0]);
162: if (PetscAbsScalar(dot[0]-dotc[0])>10*PETSC_MACHINE_EPSILON) SETERRQ(PETSC_COMM_WORLD,1,"Dots are different");
163: VecTDot(x,y,&dot[0]);
164: VecTDot(xc,yc,&dotc[0]);
165: if (PetscAbsScalar(dot[0]-dotc[0])>10*PETSC_MACHINE_EPSILON) SETERRQ(PETSC_COMM_WORLD,1,"Dots are different");
167: vecs[0] = w; vecs[1] = y;
168: VecMDot(x,2,vecs,dot);
169: vecs[0] = wc; vecs[1] = yc;
170: VecMDot(xc,2,vecs,dotc);
171: if (PetscAbsScalar(dot[0]-dotc[0])>10*PETSC_MACHINE_EPSILON || PetscAbsScalar(dot[1]-dotc[1])>10*PETSC_MACHINE_EPSILON) SETERRQ(PETSC_COMM_WORLD,1,"Dots are different");
172: vecs[0] = w; vecs[1] = y;
173: VecMTDot(x,2,vecs,dot);
174: vecs[0] = wc; vecs[1] = yc;
175: VecMTDot(xc,2,vecs,dotc);
176: if (PetscAbsScalar(dot[0]-dotc[0])>10*PETSC_MACHINE_EPSILON || PetscAbsScalar(dot[1]-dotc[1])>10*PETSC_MACHINE_EPSILON) SETERRQ(PETSC_COMM_WORLD,1,"Dots are different");
178: VecDotNorm2(x,y,&dot[0],&norm);
179: VecDotNorm2(xc,yc,&dotc[0],&normc);
180: if (PetscAbsScalar(dot[0]-dotc[0])>10*PETSC_MACHINE_EPSILON) SETERRQ(PETSC_COMM_WORLD,1,"Dots are different");
181: if (PetscAbsReal(norm-normc)>100*PETSC_MACHINE_EPSILON) SETERRQ(PETSC_COMM_WORLD,1,"Norms are different");
183: VecAbs(w);
184: VecAbs(wc);
185: VecConjugate(x);
186: VecConjugate(xc);
187: VecShift(y,0.5);
188: VecShift(yc,0.5);
189: VecReciprocal(y);
190: VecReciprocal(yc);
191: VecNorm(y,NORM_1,&norm);
192: VecNorm(yc,NORM_1,&normc);
193: if (PetscAbsReal(norm-normc)>10*PETSC_MACHINE_EPSILON) SETERRQ(PETSC_COMM_WORLD,1,"Norms are different");
195: VecPointwiseMult(w,x,y);
196: VecPointwiseMult(wc,xc,yc);
197: VecNorm(w,NORM_INFINITY,&norm);
198: VecNorm(wc,NORM_INFINITY,&normc);
199: if (PetscAbsReal(norm-normc)>10*PETSC_MACHINE_EPSILON) SETERRQ(PETSC_COMM_WORLD,1,"Norms are different");
201: VecSwap(x,y);
202: VecSwap(xc,yc);
203: VecPointwiseDivide(w,x,y);
204: VecPointwiseDivide(wc,xc,yc);
205: VecScale(w,0.3);
206: VecScale(wc,0.3);
207: VecNorm(w,NORM_1_AND_2,norm12);
208: VecNorm(wc,NORM_1_AND_2,norm12c);
209: if (PetscAbsReal(norm12[0]-norm12c[0])>10*PETSC_MACHINE_EPSILON || PetscAbsReal(norm12[1]-norm12c[1])>10*PETSC_MACHINE_EPSILON) SETERRQ(PETSC_COMM_WORLD,1,"Norms are different");
211: VecDestroy(&v);
212: VecDestroy(&w);
213: VecDestroy(&x);
214: VecDestroy(&y);
215: VecDestroy(&vparent);
216: VecDestroy(&vc);
217: VecDestroy(&wc);
218: VecDestroy(&xc);
219: VecDestroy(&yc);
220: SlepcFinalize();
221: return ierr;
222: }
224: /*TEST
226: test:
227: suffix: 1
228: nsize: {{1 2}}
230: TEST*/