Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
6 : This file is part of SLEPc.
7 : SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9 : */
10 :
11 : static char help[] = "Test block orthogonalization on a rank-deficient BV.\n\n";
12 :
13 : #include <slepcbv.h>
14 :
15 4 : int main(int argc,char **argv)
16 : {
17 4 : BV X,Z;
18 4 : Mat M,R;
19 4 : Vec v,w,t;
20 4 : PetscInt i,j,n=20,k=8;
21 4 : PetscViewer view;
22 4 : PetscBool verbose;
23 4 : PetscReal norm;
24 4 : PetscScalar alpha;
25 :
26 4 : PetscFunctionBeginUser;
27 4 : PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
28 4 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
29 4 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-k",&k,NULL));
30 4 : PetscCall(PetscOptionsHasName(NULL,NULL,"-verbose",&verbose));
31 4 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test BV block orthogonalization (length %" PetscInt_FMT ", k=%" PetscInt_FMT ").\n",n,k));
32 4 : PetscCheck(k>5,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"k must be at least 6");
33 :
34 : /* Create template vector */
35 4 : PetscCall(VecCreate(PETSC_COMM_WORLD,&t));
36 4 : PetscCall(VecSetSizes(t,PETSC_DECIDE,n));
37 4 : PetscCall(VecSetFromOptions(t));
38 :
39 : /* Create BV object X */
40 4 : PetscCall(BVCreate(PETSC_COMM_WORLD,&X));
41 4 : PetscCall(PetscObjectSetName((PetscObject)X,"X"));
42 4 : PetscCall(BVSetSizesFromVec(X,t,k));
43 4 : PetscCall(BVSetFromOptions(X));
44 :
45 : /* Set up viewer */
46 4 : PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&view));
47 4 : if (verbose) PetscCall(PetscViewerPushFormat(view,PETSC_VIEWER_ASCII_MATLAB));
48 :
49 : /* Fill X entries (first half) */
50 20 : for (j=0;j<k/2;j++) {
51 16 : PetscCall(BVGetColumn(X,j,&v));
52 16 : PetscCall(VecSet(v,0.0));
53 192 : for (i=0;i<=n/2;i++) {
54 176 : if (i+j<n) {
55 176 : alpha = (3.0*i+j-2)/(2*(i+j+1));
56 176 : PetscCall(VecSetValue(v,i+j,alpha,INSERT_VALUES));
57 : }
58 : }
59 16 : PetscCall(VecAssemblyBegin(v));
60 16 : PetscCall(VecAssemblyEnd(v));
61 16 : PetscCall(BVRestoreColumn(X,j,&v));
62 : }
63 :
64 : /* make middle column linearly dependent wrt columns 0 and 1 */
65 4 : PetscCall(BVCopyColumn(X,0,j));
66 4 : PetscCall(BVGetColumn(X,j,&v));
67 4 : PetscCall(BVGetColumn(X,1,&w));
68 4 : PetscCall(VecAXPY(v,0.5,w));
69 4 : PetscCall(BVRestoreColumn(X,1,&w));
70 4 : PetscCall(BVRestoreColumn(X,j,&v));
71 4 : j++;
72 :
73 : /* Fill X entries (second half) */
74 12 : for (;j<k-1;j++) {
75 8 : PetscCall(BVGetColumn(X,j,&v));
76 8 : PetscCall(VecSet(v,0.0));
77 96 : for (i=0;i<=n/2;i++) {
78 88 : if (i+j<n) {
79 88 : alpha = (3.0*i+j-2)/(2*(i+j+1));
80 88 : PetscCall(VecSetValue(v,i+j,alpha,INSERT_VALUES));
81 : }
82 : }
83 8 : PetscCall(VecAssemblyBegin(v));
84 8 : PetscCall(VecAssemblyEnd(v));
85 8 : PetscCall(BVRestoreColumn(X,j,&v));
86 : }
87 :
88 : /* make middle column linearly dependent wrt columns 1 and k/2+1 */
89 4 : PetscCall(BVCopyColumn(X,1,j));
90 4 : PetscCall(BVGetColumn(X,j,&v));
91 4 : PetscCall(BVGetColumn(X,k/2+1,&w));
92 4 : PetscCall(VecAXPY(v,-1.2,w));
93 4 : PetscCall(BVRestoreColumn(X,k/2+1,&w));
94 4 : PetscCall(BVRestoreColumn(X,j,&v));
95 :
96 4 : if (verbose) PetscCall(BVView(X,view));
97 :
98 : /* Create a copy on Z */
99 4 : PetscCall(BVDuplicate(X,&Z));
100 4 : PetscCall(PetscObjectSetName((PetscObject)Z,"Z"));
101 4 : PetscCall(BVCopy(X,Z));
102 :
103 : /* Test BVOrthogonalize */
104 4 : PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,k,k,NULL,&R));
105 4 : PetscCall(PetscObjectSetName((PetscObject)R,"R"));
106 4 : PetscCall(BVOrthogonalize(X,R));
107 4 : if (verbose) {
108 0 : PetscCall(BVView(X,view));
109 0 : PetscCall(MatView(R,view));
110 : }
111 :
112 : /* Check orthogonality */
113 4 : PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,k,k,NULL,&M));
114 4 : PetscCall(MatShift(M,1.0)); /* set leading part to identity */
115 4 : PetscCall(BVDot(X,X,M));
116 4 : PetscCall(MatShift(M,-1.0));
117 4 : PetscCall(MatNorm(M,NORM_1,&norm));
118 4 : if (norm<100*PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality < 100*eps\n"));
119 0 : else PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality: %g\n",(double)norm));
120 :
121 : /* Check residual */
122 4 : PetscCall(BVMult(Z,-1.0,1.0,X,R));
123 4 : PetscCall(BVNorm(Z,NORM_FROBENIUS,&norm));
124 4 : if (norm<100*PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Residual ||X-QR|| < 100*eps\n"));
125 0 : else PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Residual ||X-QR||: %g\n",(double)norm));
126 :
127 4 : PetscCall(MatDestroy(&R));
128 4 : PetscCall(MatDestroy(&M));
129 4 : PetscCall(BVDestroy(&X));
130 4 : PetscCall(BVDestroy(&Z));
131 4 : PetscCall(VecDestroy(&t));
132 4 : PetscCall(SlepcFinalize());
133 : return 0;
134 : }
135 :
136 : /*TEST
137 :
138 : test:
139 : suffix: 1
140 : nsize: 1
141 : args: -bv_orthog_block gs -bv_type {{vecs contiguous svec mat}shared output}
142 : output_file: output/test12_1.out
143 :
144 : TEST*/
|
