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[] = "Computes a partial generalized singular value decomposition (GSVD).\n"
12 : "The command line options are:\n"
13 : " -m <m>, where <m> = number of rows of A.\n"
14 : " -n <n>, where <n> = number of columns of A.\n"
15 : " -p <p>, where <p> = number of rows of B.\n\n";
16 :
17 : #include <slepcsvd.h>
18 :
19 42 : int main(int argc,char **argv)
20 : {
21 42 : Mat A,B; /* operator matrices */
22 42 : Vec u,v,x; /* singular vectors */
23 42 : SVD svd; /* singular value problem solver context */
24 42 : SVDType type;
25 42 : Vec uv,aux[2],*U,*V;
26 42 : PetscReal error,tol,sigma,lev1=0.0,lev2=0.0;
27 42 : PetscInt m=100,n,p=14,i,j,d,Istart,Iend,nsv,maxit,its,nconv;
28 42 : PetscBool flg,skiporth=PETSC_FALSE;
29 :
30 42 : PetscFunctionBeginUser;
31 42 : PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
32 :
33 42 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL));
34 42 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,&flg));
35 42 : if (!flg) n = m;
36 42 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-p",&p,&flg));
37 42 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nGeneralized singular value decomposition, (%" PetscInt_FMT "+%" PetscInt_FMT ")x%" PetscInt_FMT "\n\n",m,p,n));
38 42 : PetscCall(PetscOptionsGetBool(NULL,NULL,"-skiporth",&skiporth,NULL));
39 :
40 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
41 : Build the matrices
42 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
43 :
44 42 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
45 42 : PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m,n));
46 42 : PetscCall(MatSetFromOptions(A));
47 :
48 42 : PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
49 852 : for (i=Istart;i<Iend;i++) {
50 810 : if (i>0 && i-1<n) PetscCall(MatSetValue(A,i,i-1,-1.0,INSERT_VALUES));
51 810 : if (i+1<n) PetscCall(MatSetValue(A,i,i+1,-1.0,INSERT_VALUES));
52 810 : if (i<n) PetscCall(MatSetValue(A,i,i,2.0,INSERT_VALUES));
53 810 : if (i>n) PetscCall(MatSetValue(A,i,n-1,1.0,INSERT_VALUES));
54 : }
55 42 : PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
56 42 : PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
57 :
58 42 : PetscCall(MatCreate(PETSC_COMM_WORLD,&B));
59 42 : PetscCall(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,p,n));
60 42 : PetscCall(MatSetFromOptions(B));
61 :
62 42 : PetscCall(MatGetOwnershipRange(B,&Istart,&Iend));
63 42 : d = PetscMax(0,n-p);
64 898 : for (i=Istart;i<Iend;i++) {
65 9912 : for (j=0;j<=PetscMin(i,n-1);j++) PetscCall(MatSetValue(B,i,j+d,1.0,INSERT_VALUES));
66 : }
67 42 : PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY));
68 42 : PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY));
69 :
70 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
71 : Create the singular value solver and set various options
72 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
73 :
74 : /*
75 : Create singular value solver context
76 : */
77 42 : PetscCall(SVDCreate(PETSC_COMM_WORLD,&svd));
78 :
79 : /*
80 : Set operators and problem type
81 : */
82 42 : PetscCall(SVDSetOperators(svd,A,B));
83 42 : PetscCall(SVDSetProblemType(svd,SVD_GENERALIZED));
84 :
85 : /*
86 : Set solver parameters at runtime
87 : */
88 42 : PetscCall(SVDSetFromOptions(svd));
89 :
90 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
91 : Solve the singular value system
92 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
93 :
94 42 : PetscCall(SVDSolve(svd));
95 42 : PetscCall(SVDGetIterationNumber(svd,&its));
96 42 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %" PetscInt_FMT "\n",its));
97 :
98 : /*
99 : Optional: Get some information from the solver and display it
100 : */
101 42 : PetscCall(SVDGetType(svd,&type));
102 42 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type));
103 42 : PetscCall(SVDGetDimensions(svd,&nsv,NULL,NULL));
104 42 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of requested generalized singular values: %" PetscInt_FMT "\n",nsv));
105 42 : PetscCall(SVDGetTolerances(svd,&tol,&maxit));
106 42 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%" PetscInt_FMT "\n",(double)tol,maxit));
107 :
108 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
109 : Display solution and clean up
110 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
111 :
112 : /*
113 : Get number of converged singular triplets
114 : */
115 42 : PetscCall(SVDGetConverged(svd,&nconv));
116 42 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of converged approximate singular triplets: %" PetscInt_FMT "\n\n",nconv));
117 :
118 42 : if (nconv>0) {
119 : /*
120 : Create vectors. The interface returns u and v as stacked on top of each other
121 : [u;v] so need to create a special vector (VecNest) to extract them
122 : */
123 42 : PetscCall(MatCreateVecs(A,&x,&u));
124 42 : PetscCall(MatCreateVecs(B,NULL,&v));
125 42 : aux[0] = u;
126 42 : aux[1] = v;
127 42 : PetscCall(VecCreateNest(PETSC_COMM_WORLD,2,NULL,aux,&uv));
128 :
129 42 : PetscCall(VecDuplicateVecs(u,nconv,&U));
130 42 : PetscCall(VecDuplicateVecs(v,nconv,&V));
131 :
132 : /*
133 : Display singular values and errors relative to the norms
134 : */
135 42 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,
136 : " sigma ||r||/||[A;B]||\n"
137 : " --------------------- ------------------\n"));
138 243 : for (i=0;i<nconv;i++) {
139 : /*
140 : Get converged singular triplets: i-th singular value is stored in sigma
141 : */
142 201 : PetscCall(SVDGetSingularTriplet(svd,i,&sigma,uv,x));
143 :
144 : /* at this point, u and v can be used normally as individual vectors */
145 201 : PetscCall(VecCopy(u,U[i]));
146 201 : PetscCall(VecCopy(v,V[i]));
147 :
148 : /*
149 : Compute the error associated to each singular triplet
150 : */
151 201 : PetscCall(SVDComputeError(svd,i,SVD_ERROR_NORM,&error));
152 :
153 201 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," % 6f ",(double)sigma));
154 201 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," % 12g\n",(double)error));
155 : }
156 42 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n"));
157 :
158 42 : if (!skiporth) {
159 42 : PetscCall(VecCheckOrthonormality(U,nconv,NULL,nconv,NULL,NULL,&lev1));
160 42 : PetscCall(VecCheckOrthonormality(V,nconv,NULL,nconv,NULL,NULL,&lev2));
161 : }
162 42 : if (lev1+lev2<20*tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality below the tolerance\n"));
163 0 : else PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality: %g (U) %g (V)\n",(double)lev1,(double)lev2));
164 42 : PetscCall(VecDestroyVecs(nconv,&U));
165 42 : PetscCall(VecDestroyVecs(nconv,&V));
166 42 : PetscCall(VecDestroy(&x));
167 42 : PetscCall(VecDestroy(&u));
168 42 : PetscCall(VecDestroy(&v));
169 42 : PetscCall(VecDestroy(&uv));
170 : }
171 :
172 : /*
173 : Free work space
174 : */
175 42 : PetscCall(SVDDestroy(&svd));
176 42 : PetscCall(MatDestroy(&A));
177 42 : PetscCall(MatDestroy(&B));
178 42 : PetscCall(SlepcFinalize());
179 : return 0;
180 : }
181 :
182 : /*TEST
183 :
184 : testset:
185 : filter: grep -v "Solution method" | grep -v "Number of iterations" | sed -e "s/, maxit=1[0]*$//" | sed -e "s/[0-9]\.[0-9]*e[+-]\([0-9]*\)/removed/g"
186 : requires: double
187 : test:
188 : args: -svd_type lapack -m 20 -n 10 -p 6
189 : suffix: 1
190 : test:
191 : args: -svd_type lapack -m 15 -n 20 -p 10 -svd_smallest
192 : suffix: 2
193 : test:
194 : args: -svd_type lapack -m 15 -n 20 -p 21
195 : suffix: 3
196 : test:
197 : args: -svd_type lapack -m 20 -n 15 -p 21
198 : suffix: 4
199 :
200 : testset:
201 : args: -m 25 -n 20 -p 21 -svd_smallest -svd_nsv 2
202 : filter: grep -v "Solution method" | grep -v "Number of iterations" | sed -e "s/, maxit=1[0]*$//" | sed -e "s/[0-9]\.[0-9]*e[+-]\([0-9]*\)/removed/g"
203 : requires: double
204 : output_file: output/ex45_5.out
205 : test:
206 : args: -svd_type trlanczos -svd_ncv 8 -svd_trlanczos_gbidiag {{upper lower}} -svd_trlanczos_oneside {{0 1}}
207 : suffix: 5
208 : test:
209 : args: -svd_type cross -svd_ncv 10 -svd_cross_explicitmatrix
210 : suffix: 5_cross
211 : test:
212 : args: -svd_type cross -svd_ncv 10 -svd_cross_eps_type krylovschur -svd_cross_st_type sinvert -svd_cross_eps_target 0 -svd_cross_st_ksp_type gmres -svd_cross_st_pc_type jacobi
213 : suffix: 5_cross_implicit
214 : test:
215 : args: -svd_type cyclic -svd_ncv 12 -svd_cyclic_explicitmatrix {{0 1}}
216 : suffix: 5_cyclic
217 : requires: !complex
218 :
219 : testset:
220 : args: -m 15 -n 20 -p 21 -svd_nsv 4 -svd_ncv 9
221 : filter: grep -v "Solution method" | grep -v "Number of iterations" | sed -e "s/7.884967/7.884968/" | sed -e "s/, maxit=1[0]*$//" | sed -e "s/[0-9]\.[0-9]*e[+-]\([0-9]*\)/removed/g"
222 : requires: double
223 : output_file: output/ex45_6.out
224 : test:
225 : args: -svd_type trlanczos -svd_trlanczos_gbidiag {{single upper lower}} -svd_trlanczos_locking {{0 1}} -svd_trlanczos_oneside {{0 1}}
226 : suffix: 6
227 : test:
228 : args: -svd_type cross -svd_cross_explicitmatrix {{0 1}}
229 : suffix: 6_cross
230 :
231 : test:
232 : args: -m 15 -n 20 -p 21 -svd_nsv 4 -svd_ncv 9 -svd_type cyclic -svd_cyclic_explicitmatrix {{0 1}}
233 : filter: grep -v "Number of iterations" | sed -e "s/7.884967/7.884968/" | sed -e "s/[0-9]\.[0-9]*e[+-]\([0-9]*\)/removed/g"
234 : requires: double
235 : suffix: 6_cyclic
236 : output_file: output/ex45_6_cyclic.out
237 :
238 : testset:
239 : args: -m 20 -n 15 -p 21 -svd_nsv 4 -svd_ncv 9
240 : filter: grep -v "Solution method" | grep -v "Number of iterations" | sed -e "s/, maxit=1[0]*$//" | sed -e "s/[0-9]\.[0-9]*e[+-]\([0-9]*\)/removed/g"
241 : requires: double
242 : output_file: output/ex45_7.out
243 : test:
244 : args: -svd_type trlanczos -svd_trlanczos_gbidiag {{single upper lower}} -svd_trlanczos_restart 0.4 -svd_trlanczos_oneside {{0 1}}
245 : suffix: 7
246 : test:
247 : args: -svd_type cross -svd_cross_explicitmatrix {{0 1}}
248 : suffix: 7_cross
249 :
250 : test:
251 : args: -m 20 -n 15 -p 21 -svd_nsv 4 -svd_ncv 16 -svd_type cyclic -svd_cyclic_explicitmatrix {{0 1}}
252 : filter: grep -v "Number of iterations" | sed -e "s/[0-9]\.[0-9]*e[+-]\([0-9]*\)/removed/g"
253 : requires: double
254 : suffix: 7_cyclic
255 : output_file: output/ex45_7_cyclic.out
256 :
257 : test:
258 : args: -m 25 -n 20 -p 21 -svd_smallest -svd_nsv 2 -svd_ncv 5 -svd_type trlanczos -svd_trlanczos_gbidiag {{upper lower}} -svd_trlanczos_scale {{0.1 -20}}
259 : filter: grep -v "Solution method" | grep -v "Number of iterations" | grep -v "Stopping condition" | sed -e "s/, maxit=1[0]*$//" | sed -e "s/[0-9]\.[0-9]*e[+-]\([0-9]*\)/removed/g"
260 : suffix: 8
261 :
262 : TEST*/
|
