Actual source code: test1.c

slepc-3.17.0 2022-03-31
Report Typos and Errors
  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:   SlepcInitialize(&argc,&argv,(char*)0,help);
 25:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 26:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 28:   PetscPrintf(PETSC_COMM_WORLD,"VecComp test\n");

 30:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 31:      Create standard vectors
 32:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 34:   VecCreate(PETSC_COMM_WORLD,&v);
 35:   VecSetSizes(v,8/size,8);
 36:   VecSetFromOptions(v);

 38:   if (!rank) {
 39:     VecSetValue(v,0,2.0,INSERT_VALUES);
 40:     VecSetValue(v,1,-1.0,INSERT_VALUES);
 41:     VecSetValue(v,2,3.0,INSERT_VALUES);
 42:     VecSetValue(v,3,3.5,INSERT_VALUES);
 43:   }
 44:   if ((!rank && size==1) || (rank && size==2)) {
 45:     VecSetValue(v,4,1.2,INSERT_VALUES);
 46:     VecSetValue(v,5,1.8,INSERT_VALUES);
 47:     VecSetValue(v,6,-2.2,INSERT_VALUES);
 48:     VecSetValue(v,7,2.0,INSERT_VALUES);
 49:   }
 50:   VecAssemblyBegin(v);
 51:   VecAssemblyEnd(v);
 52:   VecDuplicate(v,&w);
 53:   VecSet(w,1.0);
 54:   VecDuplicate(v,&x);
 55:   VecDuplicate(v,&y);
 56:   if (!rank) VecSetValue(y,0,1.0,INSERT_VALUES);
 57:   VecAssemblyBegin(y);
 58:   VecAssemblyEnd(y);

 60:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 61:      Create veccomp vectors
 62:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 64:   VecCreate(PETSC_COMM_WORLD,&vparent);
 65:   VecSetSizes(vparent,4/size,4);
 66:   VecSetFromOptions(vparent);

 68:   /* create a veccomp vector with two subvectors */
 69:   VecDuplicate(vparent,&vchild[0]);
 70:   VecDuplicate(vparent,&vchild[1]);
 71:   if (!rank) {
 72:     VecSetValue(vchild[0],0,2.0,INSERT_VALUES);
 73:     VecSetValue(vchild[0],1,-1.0,INSERT_VALUES);
 74:     VecSetValue(vchild[1],0,1.2,INSERT_VALUES);
 75:     VecSetValue(vchild[1],1,1.8,INSERT_VALUES);
 76:   }
 77:   if ((!rank && size==1) || (rank && size==2)) {
 78:     VecSetValue(vchild[0],2,3.0,INSERT_VALUES);
 79:     VecSetValue(vchild[0],3,3.5,INSERT_VALUES);
 80:     VecSetValue(vchild[1],2,-2.2,INSERT_VALUES);
 81:     VecSetValue(vchild[1],3,2.0,INSERT_VALUES);
 82:   }
 83:   VecAssemblyBegin(vchild[0]);
 84:   VecAssemblyBegin(vchild[1]);
 85:   VecAssemblyEnd(vchild[0]);
 86:   VecAssemblyEnd(vchild[1]);
 87:   VecCreateCompWithVecs(vchild,2,vparent,&vc);
 88:   VecDestroy(&vchild[0]);
 89:   VecDestroy(&vchild[1]);
 90:   VecView(vc,NULL);

 92:   VecGetSize(vc,&k);

 95:   /* create an empty veccomp vector with two subvectors */
 96:   Nx[0] = 4;
 97:   Nx[1] = 4;
 98:   VecCreateComp(PETSC_COMM_WORLD,Nx,2,VECSTANDARD,vparent,&wc);
 99:   VecCompGetSubVecs(wc,&n,&varray);
101:   for (i=0;i<2;i++) VecSet(varray[i],1.0);

103:   VecGetSize(wc,&k);

106:   /* duplicate a veccomp */
107:   VecDuplicate(vc,&xc);

109:   /* create a veccomp via VecSetType */
110:   VecCreate(PETSC_COMM_WORLD,&yc);
111:   VecSetType(yc,VECCOMP);
112:   VecSetSizes(yc,8/size,8);
113:   VecCompSetSubVecs(yc,2,NULL);

115:   VecCompGetSubVecs(yc,&n,&varray);
117:   if (!rank) VecSetValue(varray[0],0,1.0,INSERT_VALUES);
118:   VecAssemblyBegin(varray[0]);
119:   VecAssemblyEnd(varray[0]);

121:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
122:      Operate with vectors
123:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

125:   VecCopy(w,x);
126:   VecAXPBY(x,1.0,-2.0,v);
127:   VecNorm(x,NORM_2,&norm);
128:   VecCopy(wc,xc);
129:   VecAXPBY(xc,1.0,-2.0,vc);
130:   VecNorm(xc,NORM_2,&normc);

133:   VecCopy(w,x);
134:   VecWAXPY(x,-2.0,w,v);
135:   VecNorm(x,NORM_2,&norm);
136:   VecCopy(wc,xc);
137:   VecWAXPY(xc,-2.0,wc,vc);
138:   VecNorm(xc,NORM_2,&normc);

141:   VecAXPBYPCZ(y,3.0,-1.0,1.0,w,v);
142:   VecNorm(y,NORM_2,&norm);
143:   VecAXPBYPCZ(yc,3.0,-1.0,1.0,wc,vc);
144:   VecNorm(yc,NORM_2,&normc);

147:   VecMax(xc,NULL,&vmax);
148:   VecMin(xc,NULL,&vmin);
149:   PetscPrintf(PETSC_COMM_WORLD,"xc has max value %g min value %g\n",(double)vmax,(double)vmin);

151:   VecMaxPointwiseDivide(wc,xc,&vmax);
152:   PetscPrintf(PETSC_COMM_WORLD,"wc/xc has max value %g\n",(double)vmax);

154:   VecDot(x,y,&dot[0]);
155:   VecDot(xc,yc,&dotc[0]);
157:   VecTDot(x,y,&dot[0]);
158:   VecTDot(xc,yc,&dotc[0]);

161:   vecs[0] = w; vecs[1] = y;
162:   VecMDot(x,2,vecs,dot);
163:   vecs[0] = wc; vecs[1] = yc;
164:   VecMDot(xc,2,vecs,dotc);
166:   vecs[0] = w; vecs[1] = y;
167:   VecMTDot(x,2,vecs,dot);
168:   vecs[0] = wc; vecs[1] = yc;
169:   VecMTDot(xc,2,vecs,dotc);

172:   VecDotNorm2(x,y,&dot[0],&norm);
173:   VecDotNorm2(xc,yc,&dotc[0],&normc);

177:   VecAbs(w);
178:   VecAbs(wc);
179:   VecConjugate(x);
180:   VecConjugate(xc);
181:   VecShift(y,0.5);
182:   VecShift(yc,0.5);
183:   VecReciprocal(y);
184:   VecReciprocal(yc);
185:   VecExp(y);
186:   VecExp(yc);
187:   VecLog(y);
188:   VecLog(yc);
189:   VecNorm(y,NORM_1,&norm);
190:   VecNorm(yc,NORM_1,&normc);

193:   VecPointwiseMult(w,x,y);
194:   VecPointwiseMult(wc,xc,yc);
195:   VecNorm(w,NORM_INFINITY,&norm);
196:   VecNorm(wc,NORM_INFINITY,&normc);

199:   VecPointwiseMax(w,x,y);
200:   VecPointwiseMax(wc,xc,yc);
201:   VecNorm(w,NORM_INFINITY,&norm);
202:   VecNorm(wc,NORM_INFINITY,&normc);

205:   VecSwap(x,y);
206:   VecSwap(xc,yc);
207:   VecPointwiseDivide(w,x,y);
208:   VecPointwiseDivide(wc,xc,yc);
209:   VecScale(w,0.3);
210:   VecScale(wc,0.3);
211:   VecSqrtAbs(w);
212:   VecSqrtAbs(wc);
213:   VecNorm(w,NORM_1_AND_2,norm12);
214:   VecNorm(wc,NORM_1_AND_2,norm12c);

217:   VecPointwiseMin(w,x,y);
218:   VecPointwiseMin(wc,xc,yc);
219:   VecPointwiseMaxAbs(x,y,w);
220:   VecPointwiseMaxAbs(xc,yc,wc);
221:   VecNorm(x,NORM_INFINITY,&norm);
222:   VecNorm(xc,NORM_INFINITY,&normc);

225:   VecSetRandom(wc,NULL);

227:   VecDestroy(&v);
228:   VecDestroy(&w);
229:   VecDestroy(&x);
230:   VecDestroy(&y);
231:   VecDestroy(&vparent);
232:   VecDestroy(&vc);
233:   VecDestroy(&wc);
234:   VecDestroy(&xc);
235:   VecDestroy(&yc);
236:   SlepcFinalize();
237:   return 0;
238: }

240: /*TEST

242:    test:
243:       suffix: 1

245:    test:
246:       suffix: 2
247:       nsize: 2

249: TEST*/