Actual source code: test6.c
slepc-3.17.2 2022-08-09
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 BV orthogonalization functions with constraints.\n\n";
13: #include <slepcbv.h>
15: int main(int argc,char **argv)
16: {
17: BV X;
18: Mat M;
19: Vec v,t,*C;
20: PetscInt i,j,n=20,k=8,nc=2,nloc;
21: PetscViewer view;
22: PetscBool verbose;
23: PetscReal norm;
24: PetscScalar alpha;
26: SlepcInitialize(&argc,&argv,(char*)0,help);
27: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
28: PetscOptionsGetInt(NULL,NULL,"-k",&k,NULL);
29: PetscOptionsGetInt(NULL,NULL,"-nc",&nc,NULL);
30: PetscOptionsHasName(NULL,NULL,"-verbose",&verbose);
31: PetscPrintf(PETSC_COMM_WORLD,"Test BV orthogonalization with %" PetscInt_FMT " columns + %" PetscInt_FMT " constraints, of length %" PetscInt_FMT ".\n",k,nc,n);
33: /* Create template vector */
34: VecCreate(PETSC_COMM_WORLD,&t);
35: VecSetSizes(t,PETSC_DECIDE,n);
36: VecSetFromOptions(t);
37: VecGetLocalSize(t,&nloc);
39: /* Create BV object X */
40: BVCreate(PETSC_COMM_WORLD,&X);
41: PetscObjectSetName((PetscObject)X,"X");
42: BVSetSizes(X,nloc,n,k);
43: BVSetFromOptions(X);
45: /* Generate constraints and attach them to X */
46: if (nc>0) {
47: VecDuplicateVecs(t,nc,&C);
48: for (j=0;j<nc;j++) {
49: for (i=0;i<=j;i++) VecSetValue(C[j],i,1.0,INSERT_VALUES);
50: VecAssemblyBegin(C[j]);
51: VecAssemblyEnd(C[j]);
52: }
53: BVInsertConstraints(X,&nc,C);
54: VecDestroyVecs(nc,&C);
55: }
57: /* Set up viewer */
58: PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&view);
59: if (verbose) PetscViewerPushFormat(view,PETSC_VIEWER_ASCII_MATLAB);
61: /* Fill X entries */
62: for (j=0;j<k;j++) {
63: BVGetColumn(X,j,&v);
64: VecSet(v,0.0);
65: for (i=0;i<=n/2;i++) {
66: if (i+j<n) {
67: alpha = (3.0*i+j-2)/(2*(i+j+1));
68: VecSetValue(v,i+j,alpha,INSERT_VALUES);
69: }
70: }
71: VecAssemblyBegin(v);
72: VecAssemblyEnd(v);
73: BVRestoreColumn(X,j,&v);
74: }
75: if (verbose) BVView(X,view);
77: /* Test BVOrthogonalizeColumn */
78: for (j=0;j<k;j++) {
79: BVOrthogonalizeColumn(X,j,NULL,&norm,NULL);
80: alpha = 1.0/norm;
81: BVScaleColumn(X,j,alpha);
82: }
83: if (verbose) BVView(X,view);
85: /* Check orthogonality */
86: MatCreateSeqDense(PETSC_COMM_SELF,k,k,NULL,&M);
87: BVDot(X,X,M);
88: MatShift(M,-1.0);
89: MatNorm(M,NORM_1,&norm);
90: if (norm<100*PETSC_MACHINE_EPSILON) PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality < 100*eps\n");
91: else PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality: %g\n",(double)norm);
93: MatDestroy(&M);
94: BVDestroy(&X);
95: VecDestroy(&t);
96: SlepcFinalize();
97: return 0;
98: }
100: /*TEST
102: testset:
103: output_file: output/test6_1.out
104: test:
105: suffix: 1
106: args: -bv_type {{vecs contiguous svec mat}shared output}
107: test:
108: suffix: 1_cuda
109: args: -bv_type svec -vec_type cuda
110: requires: cuda
111: test:
112: suffix: 2
113: nsize: 2
114: args: -bv_type {{vecs contiguous svec mat}shared output}
115: test:
116: suffix: 3
117: args: -bv_type {{vecs contiguous svec mat}shared output} -bv_orthog_type mgs
118: test:
119: suffix: 3_cuda
120: args: -bv_type svec -vec_type cuda -bv_orthog_type mgs
121: requires: cuda
123: TEST*/