Actual source code: test5.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 BV operations with indefinite inner product.\n\n";
13: #include <slepcbv.h>
15: int main(int argc,char **argv)
16: {
18: Vec t,v,w,omega;
19: Mat B,M;
20: BV X,Y;
21: PetscInt i,j,n=10,k=5,l,Istart,Iend;
22: PetscScalar alpha;
23: PetscReal nrm;
24: PetscViewer view;
25: PetscBool verbose;
27: SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
28: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
29: PetscOptionsGetInt(NULL,NULL,"-k",&k,NULL);
30: PetscOptionsHasName(NULL,NULL,"-verbose",&verbose);
31: PetscPrintf(PETSC_COMM_WORLD,"Test BV with indefinite inner product (n=%D, k=%D).\n",n,k);
33: /* Create inner product matrix (standard involutionary permutation) */
34: MatCreate(PETSC_COMM_WORLD,&B);
35: MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,n,n);
36: MatSetFromOptions(B);
37: MatSetUp(B);
38: PetscObjectSetName((PetscObject)B,"B");
40: MatGetOwnershipRange(B,&Istart,&Iend);
41: for (i=Istart;i<Iend;i++) {
42: MatSetValue(B,i,n-i-1,1.0,INSERT_VALUES);
43: }
44: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
45: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
46: MatCreateVecs(B,&t,NULL);
48: /* Create BV object X */
49: BVCreate(PETSC_COMM_WORLD,&X);
50: PetscObjectSetName((PetscObject)X,"X");
51: BVSetSizesFromVec(X,t,k);
52: BVSetFromOptions(X);
53: BVSetMatrix(X,B,PETSC_TRUE);
55: /* Set up viewer */
56: PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&view);
57: if (verbose) {
58: PetscViewerPushFormat(view,PETSC_VIEWER_ASCII_MATLAB);
59: }
61: /* Fill X entries */
62: l = -3;
63: for (j=0;j<k;j++) {
64: BVGetColumn(X,j,&v);
65: VecSet(v,-1.0);
66: for (i=0;i<n/2;i++) {
67: if (i+j<n) {
68: l = (l + 3*i+j-2) % n;
69: VecSetValue(v,i+j,(PetscScalar)l,INSERT_VALUES);
70: }
71: }
72: VecAssemblyBegin(v);
73: VecAssemblyEnd(v);
74: BVRestoreColumn(X,j,&v);
75: }
76: if (verbose) {
77: MatView(B,view);
78: BVView(X,view);
79: }
81: /* Test BVNormColumn */
82: BVNormColumn(X,0,NORM_2,&nrm);
83: PetscPrintf(PETSC_COMM_WORLD,"B-Norm of X[0] = %g\n",(double)nrm);
85: /* Test BVOrthogonalizeColumn */
86: for (j=0;j<k;j++) {
87: BVOrthogonalizeColumn(X,j,NULL,&nrm,NULL);
88: alpha = 1.0/nrm;
89: BVScaleColumn(X,j,alpha);
90: }
91: if (verbose) {
92: BVView(X,view);
93: }
95: /* Create a copy on Y */
96: BVDuplicate(X,&Y);
97: PetscObjectSetName((PetscObject)Y,"Y");
98: BVCopy(X,Y);
100: /* Check orthogonality */
101: MatCreateSeqDense(PETSC_COMM_SELF,k,k,NULL,&M);
102: BVDot(Y,Y,M);
103: VecCreateSeq(PETSC_COMM_SELF,k,&omega);
104: BVGetSignature(Y,omega);
105: VecScale(omega,-1.0);
106: MatDiagonalSet(M,omega,ADD_VALUES);
107: MatNorm(M,NORM_1,&nrm);
108: if (nrm<100*PETSC_MACHINE_EPSILON) {
109: PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality < 100*eps\n");
110: } else {
111: PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality: %g\n",(double)nrm);
112: }
114: /* Test BVSetSignature */
115: VecScale(omega,-1.0);
116: BVSetSignature(Y,omega);
117: VecDestroy(&omega);
119: /* Test BVApplyMatrix */
120: VecDuplicate(t,&w);
121: BVGetColumn(X,0,&v);
122: BVApplyMatrix(X,v,w);
123: BVApplyMatrix(X,w,t);
124: VecAXPY(t,-1.0,v);
125: BVRestoreColumn(X,0,&v);
126: VecNorm(t,NORM_2,&nrm);
127: if (PetscAbsReal(nrm)>10*PETSC_MACHINE_EPSILON) SETERRQ1(PETSC_COMM_WORLD,1,"Wrong value, nrm = %g\n",(double)nrm);
129: BVApplyMatrixBV(X,Y);
130: BVGetColumn(Y,0,&v);
131: VecAXPY(w,-1.0,v);
132: BVRestoreColumn(Y,0,&v);
133: VecNorm(w,NORM_2,&nrm);
134: if (PetscAbsReal(nrm)>10*PETSC_MACHINE_EPSILON) SETERRQ1(PETSC_COMM_WORLD,1,"Wrong value, nrm = %g\n",(double)nrm);
136: BVDestroy(&X);
137: BVDestroy(&Y);
138: MatDestroy(&M);
139: MatDestroy(&B);
140: VecDestroy(&w);
141: VecDestroy(&t);
142: SlepcFinalize();
143: return ierr;
144: }
146: /*TEST
148: test:
149: suffix: 1
150: nsize: 1
151: args: -bv_orthog_refine always -bv_type {{vecs contiguous svec mat}shared output}
152: output_file: output/test5_1.out
154: test:
155: suffix: 1_cuda
156: nsize: 1
157: args: -bv_orthog_refine always -bv_type svec -mat_type aijcusparse
158: requires: cuda
159: output_file: output/test5_1.out
161: test:
162: suffix: 2
163: nsize: 1
164: args: -bv_orthog_refine always -bv_type {{vecs contiguous svec mat}shared output} -bv_orthog_type mgs
165: output_file: output/test5_1.out
167: test:
168: suffix: 2_cuda
169: nsize: 1
170: args: -bv_orthog_refine always -bv_type svec -mat_type aijcusparse -bv_orthog_type mgs
171: requires: cuda
172: output_file: output/test5_1.out
175: TEST*/