Actual source code: test2.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.\n\n";
13: #include <slepcbv.h>
15: int main(int argc,char **argv)
16: {
17: BV X,Y,Z;
18: Mat M,R;
19: Vec v,t,e;
20: PetscInt i,j,n=20,k=8;
21: PetscViewer view;
22: PetscBool verbose;
23: PetscReal norm,condn=1.0;
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: PetscOptionsGetReal(NULL,NULL,"-condn",&condn,NULL);
31: PetscOptionsHasName(NULL,NULL,"-verbose",&verbose);
32: PetscPrintf(PETSC_COMM_WORLD,"Test BV orthogonalization with %" PetscInt_FMT " columns of length %" PetscInt_FMT ".\n",k,n);
33: if (condn>1.0) PetscPrintf(PETSC_COMM_WORLD," - Using a random BV with condition number = %g\n",(double)condn);
35: /* Create template vector */
36: VecCreate(PETSC_COMM_WORLD,&t);
37: VecSetSizes(t,PETSC_DECIDE,n);
38: VecSetFromOptions(t);
40: /* Create BV object X */
41: BVCreate(PETSC_COMM_WORLD,&X);
42: PetscObjectSetName((PetscObject)X,"X");
43: BVSetSizesFromVec(X,t,k);
44: BVSetFromOptions(X);
46: /* Set up viewer */
47: PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&view);
48: if (verbose) PetscViewerPushFormat(view,PETSC_VIEWER_ASCII_MATLAB);
50: /* Fill X entries */
51: if (condn==1.0) {
52: for (j=0;j<k;j++) {
53: BVGetColumn(X,j,&v);
54: VecSet(v,0.0);
55: for (i=0;i<=n/2;i++) {
56: if (i+j<n) {
57: alpha = (3.0*i+j-2)/(2*(i+j+1));
58: VecSetValue(v,i+j,alpha,INSERT_VALUES);
59: }
60: }
61: VecAssemblyBegin(v);
62: VecAssemblyEnd(v);
63: BVRestoreColumn(X,j,&v);
64: }
65: } else BVSetRandomCond(X,condn);
66: if (verbose) BVView(X,view);
68: /* Create copies on Y and Z */
69: BVDuplicate(X,&Y);
70: PetscObjectSetName((PetscObject)Y,"Y");
71: BVCopy(X,Y);
72: BVDuplicate(X,&Z);
73: PetscObjectSetName((PetscObject)Z,"Z");
74: BVCopy(X,Z);
76: /* Test BVOrthogonalizeColumn */
77: for (j=0;j<k;j++) {
78: BVOrthogonalizeColumn(X,j,NULL,&norm,NULL);
79: alpha = 1.0/norm;
80: BVScaleColumn(X,j,alpha);
81: }
82: if (verbose) BVView(X,view);
84: /* Check orthogonality */
85: MatCreateSeqDense(PETSC_COMM_SELF,k,k,NULL,&M);
86: BVDot(X,X,M);
87: MatShift(M,-1.0);
88: MatNorm(M,NORM_1,&norm);
89: if (norm<100*PETSC_MACHINE_EPSILON) PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality < 100*eps\n");
90: else PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality: %g\n",(double)norm);
92: /* Test BVOrthogonalize */
93: MatCreateSeqDense(PETSC_COMM_SELF,k,k,NULL,&R);
94: PetscObjectSetName((PetscObject)R,"R");
95: BVOrthogonalize(Y,R);
96: if (verbose) {
97: BVView(Y,view);
98: MatView(R,view);
99: }
101: /* Check orthogonality */
102: BVDot(Y,Y,M);
103: MatShift(M,-1.0);
104: MatNorm(M,NORM_1,&norm);
105: if (norm<100*PETSC_MACHINE_EPSILON) PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality < 100*eps\n");
106: else PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality: %g\n",(double)norm);
108: /* Check residual */
109: BVMult(Z,-1.0,1.0,Y,R);
110: BVNorm(Z,NORM_FROBENIUS,&norm);
111: if (norm<100*PETSC_MACHINE_EPSILON) PetscPrintf(PETSC_COMM_WORLD,"Residual ||X-QR|| < 100*eps\n");
112: else PetscPrintf(PETSC_COMM_WORLD,"Residual ||X-QR||: %g\n",(double)norm);
114: /* Test BVOrthogonalizeVec */
115: VecDuplicate(t,&e);
116: VecSet(e,1.0);
117: BVOrthogonalizeVec(X,e,NULL,&norm,NULL);
118: PetscPrintf(PETSC_COMM_WORLD,"Norm of ones(n,1) after orthogonalizing against X: %g\n",(double)norm);
120: MatDestroy(&M);
121: MatDestroy(&R);
122: BVDestroy(&X);
123: BVDestroy(&Y);
124: BVDestroy(&Z);
125: VecDestroy(&e);
126: VecDestroy(&t);
127: SlepcFinalize();
128: return 0;
129: }
131: /*TEST
133: testset:
134: output_file: output/test2_1.out
135: test:
136: suffix: 1
137: args: -bv_type {{vecs contiguous svec mat}shared output} -bv_orthog_type cgs
138: test:
139: suffix: 1_cuda
140: args: -bv_type svec -vec_type cuda -bv_orthog_type cgs
141: requires: cuda
142: test:
143: suffix: 2
144: args: -bv_type {{vecs contiguous svec mat}shared output} -bv_orthog_type mgs
145: test:
146: suffix: 2_cuda
147: args: -bv_type svec -vec_type cuda -bv_orthog_type mgs
148: requires: cuda
150: test:
151: suffix: 3
152: nsize: 1
153: args: -bv_type {{vecs contiguous svec mat}shared output} -condn 1e8
154: requires: !single
155: filter: grep -v "against"
156: output_file: output/test2_3.out
158: TEST*/