Actual source code: test1.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  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*/