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[] = "Solves a GSVD problem with matrices loaded from a file.\n"
12 : "The command line options are:\n"
13 : " -f1 <filename>, PETSc binary file containing matrix A.\n"
14 : " -f2 <filename>, PETSc binary file containing matrix B (optional). Instead of"
15 : " a file it is possible to specify one of 'identity', 'bidiagonal' or 'tridiagonal'"
16 : " -p <p>, in case B is not taken from a file.\n\n";
17 :
18 : #include <slepcsvd.h>
19 :
20 9 : int main(int argc,char **argv)
21 : {
22 9 : Mat A,B; /* matrices */
23 9 : SVD svd; /* singular value problem solver context */
24 9 : PetscInt i,m,n,p,Istart,Iend,col[3];
25 9 : PetscScalar vals[3];
26 9 : char filename[PETSC_MAX_PATH_LEN];
27 9 : PetscViewer viewer;
28 9 : PetscBool flg,terse;
29 :
30 9 : PetscFunctionBeginUser;
31 9 : PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
32 :
33 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
34 : Load matrices that define the generalized singular value problem
35 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
36 :
37 9 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nGeneralized singular value problem stored in file.\n\n"));
38 9 : PetscCall(PetscOptionsGetString(NULL,NULL,"-f1",filename,sizeof(filename),&flg));
39 9 : PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Must indicate a file name for matrix A with the -f1 option");
40 :
41 : #if defined(PETSC_USE_COMPLEX)
42 9 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Reading COMPLEX matrix from a binary file...\n"));
43 : #else
44 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Reading REAL matrix from a binary file...\n"));
45 : #endif
46 9 : PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer));
47 9 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
48 9 : PetscCall(MatSetFromOptions(A));
49 9 : PetscCall(MatLoad(A,viewer));
50 9 : PetscCall(PetscViewerDestroy(&viewer));
51 :
52 9 : PetscCall(MatGetSize(A,&m,&n));
53 :
54 9 : PetscCall(PetscOptionsGetString(NULL,NULL,"-f2",filename,sizeof(filename),&flg));
55 9 : PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Must indicate a file name for matrix B with the -f2 option, or alternatively the strings 'identity', 'bidiagonal', or 'tridiagonal'");
56 9 : PetscCall(PetscStrcmp(filename,"identity",&flg));
57 9 : if (flg) {
58 0 : p = n;
59 0 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-p",&p,&flg));
60 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Using B=I with p=%" PetscInt_FMT "\n\n",p));
61 0 : PetscCall(MatCreate(PETSC_COMM_WORLD,&B));
62 0 : PetscCall(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,p,n));
63 0 : PetscCall(MatSetFromOptions(B));
64 0 : PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY));
65 0 : PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY));
66 0 : PetscCall(MatShift(B,1.0));
67 : } else {
68 9 : PetscCall(PetscStrcmp(filename,"bidiagonal",&flg));
69 9 : if (flg) {
70 9 : p = n+1;
71 9 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-p",&p,&flg));
72 9 : vals[0]=-1; vals[1]=1;
73 9 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Using B=bidiag(1,-1) with p=%" PetscInt_FMT "\n\n",p));
74 9 : PetscCall(MatCreate(PETSC_COMM_WORLD,&B));
75 9 : PetscCall(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,p,n));
76 9 : PetscCall(MatSetFromOptions(B));
77 9 : PetscCall(MatGetOwnershipRange(B,&Istart,&Iend));
78 2929 : for (i=Istart;i<Iend;i++) {
79 2920 : col[0]=i-1; col[1]=i;
80 2920 : if (i==0) PetscCall(MatSetValue(B,i,col[1],vals[1],INSERT_VALUES));
81 2911 : else if (i<n) PetscCall(MatSetValues(B,1,&i,2,col,vals,INSERT_VALUES));
82 2920 : else if (i==n) PetscCall(MatSetValue(B,i,col[0],vals[0],INSERT_VALUES));
83 : }
84 9 : PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY));
85 9 : PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY));
86 : } else {
87 0 : PetscCall(PetscStrcmp(filename,"tridiagonal",&flg));
88 0 : if (flg) {
89 0 : p = n-2;
90 0 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-p",&p,&flg));
91 0 : vals[0]=-1; vals[1]=2; vals[2]=-1;
92 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Using B=tridiag(-1,2,-1) with p=%" PetscInt_FMT "\n\n",p));
93 0 : PetscCall(MatCreate(PETSC_COMM_WORLD,&B));
94 0 : PetscCall(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,p,n));
95 0 : PetscCall(MatSetFromOptions(B));
96 0 : PetscCall(MatGetOwnershipRange(B,&Istart,&Iend));
97 0 : for (i=Istart;i<Iend;i++) {
98 0 : col[0]=i; col[1]=i+1; col[2]=i+2;
99 0 : PetscCall(MatSetValues(B,1,&i,3,col,vals,INSERT_VALUES));
100 : }
101 0 : PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY));
102 0 : PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY));
103 : } else { /* load file */
104 0 : PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer));
105 0 : PetscCall(MatCreate(PETSC_COMM_WORLD,&B));
106 0 : PetscCall(MatSetFromOptions(B));
107 0 : PetscCall(MatLoad(B,viewer));
108 0 : PetscCall(PetscViewerDestroy(&viewer));
109 : }
110 : }
111 : }
112 :
113 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114 : Create the singular value solver and set various options
115 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
116 :
117 : /*
118 : Create singular value solver context
119 : */
120 9 : PetscCall(SVDCreate(PETSC_COMM_WORLD,&svd));
121 :
122 : /*
123 : Set operators of GSVD problem
124 : */
125 9 : PetscCall(SVDSetOperators(svd,A,B));
126 9 : PetscCall(SVDSetProblemType(svd,SVD_GENERALIZED));
127 :
128 : /*
129 : Set solver parameters at runtime
130 : */
131 9 : PetscCall(SVDSetFromOptions(svd));
132 :
133 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
134 : Solve the problem and print solution
135 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
136 :
137 9 : PetscCall(SVDSolve(svd));
138 :
139 : /* show detailed info unless -terse option is given by user */
140 9 : PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
141 9 : if (terse) PetscCall(SVDErrorView(svd,SVD_ERROR_NORM,NULL));
142 : else {
143 0 : PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
144 0 : PetscCall(SVDConvergedReasonView(svd,PETSC_VIEWER_STDOUT_WORLD));
145 0 : PetscCall(SVDErrorView(svd,SVD_ERROR_NORM,PETSC_VIEWER_STDOUT_WORLD));
146 0 : PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
147 : }
148 9 : PetscCall(SVDDestroy(&svd));
149 9 : PetscCall(MatDestroy(&A));
150 9 : PetscCall(MatDestroy(&B));
151 9 : PetscCall(SlepcFinalize());
152 : return 0;
153 : }
154 : /*TEST
155 :
156 : testset:
157 : requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
158 : args: -f1 ${SLEPC_DIR}/share/slepc/datafiles/matrices/bfw62a.petsc -f2 ${SLEPC_DIR}/share/slepc/datafiles/matrices/bfw62b.petsc -svd_nsv 3 -terse
159 : output_file: output/ex48_1.out
160 : test:
161 : suffix: 1
162 : args: -svd_type trlanczos -svd_trlanczos_explicitmatrix {{0 1}} -svd_trlanczos_scale 1e5 -svd_trlanczos_ksp_rtol 1e-15
163 : test:
164 : suffix: 1_spqr
165 : args: -svd_type trlanczos -svd_trlanczos_explicitmatrix -svd_trlanczos_pc_type qr -svd_trlanczos_scale 1e5 -svd_trlanczos_oneside {{0 1}}
166 : requires: suitesparse
167 : test:
168 : suffix: 1_autoscale
169 : args: -svd_type trlanczos -svd_trlanczos_gbidiag {{lower upper}} -svd_trlanczos_scale -5 -svd_trlanczos_ksp_rtol 1e-16 -svd_trlanczos_oneside {{0 1}}
170 : test:
171 : suffix: 1_cross
172 : args: -svd_type cross -svd_cross_explicitmatrix
173 : test:
174 : suffix: 1_cyclic
175 : args: -svd_type cyclic -svd_cyclic_explicitmatrix
176 :
177 : testset:
178 : requires: double complex datafilespath !defined(PETSC_USE_64BIT_INDICES)
179 : args: -f1 ${DATAFILESPATH}/matrices/complex/qc324.petsc -f2 bidiagonal -svd_nsv 3 -terse
180 : output_file: output/ex48_2.out
181 : filter: sed -e "s/30749/30748/"
182 : timeoutfactor: 2
183 : test:
184 : suffix: 2
185 : args: -svd_type trlanczos -svd_trlanczos_explicitmatrix {{0 1}} -svd_trlanczos_ksp_rtol 1e-10 -svd_trlanczos_scale 100
186 : requires: !defined(PETSCTEST_VALGRIND)
187 : test:
188 : suffix: 2_spqr
189 : args: -svd_type trlanczos -svd_trlanczos_explicitmatrix -svd_trlanczos_pc_type qr -svd_trlanczos_ksp_rtol 1e-10
190 : requires: suitesparse
191 : test:
192 : suffix: 2_cross
193 : args: -svd_type cross -svd_cross_explicitmatrix
194 : test:
195 : suffix: 2_cyclic
196 : args: -svd_type cyclic -svd_cyclic_explicitmatrix
197 :
198 : test:
199 : requires: double complex datafilespath !defined(PETSC_USE_64BIT_INDICES) !defined(PETSCTEST_VALGRIND)
200 : args: -f1 ${DATAFILESPATH}/matrices/complex/qc324.petsc -f2 bidiagonal -p 320 -svd_nsv 3 -svd_type trlanczos -svd_trlanczos_ksp_rtol 1e-14 -svd_trlanczos_scale 100 -terse
201 : timeoutfactor: 2
202 : suffix: 3
203 :
204 : testset:
205 : requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
206 : args: -f1 ${SLEPC_DIR}/share/slepc/datafiles/matrices/rdb200.petsc -f2 identity -svd_nsv 3 -svd_ncv 24 -svd_smallest -terse
207 : output_file: output/ex48_4.out
208 : test:
209 : suffix: 4
210 : args: -svd_type trlanczos
211 : test:
212 : suffix: 4_spqr
213 : args: -svd_type trlanczos -svd_trlanczos_explicitmatrix -svd_trlanczos_pc_type qr
214 : requires: suitesparse
215 : test:
216 : suffix: 4_cross
217 : args: -svd_type cross -svd_cross_explicitmatrix
218 : test:
219 : suffix: 4_cross_implicit
220 : args: -svd_type cross -svd_cross_eps_type lobpcg -svd_cross_st_ksp_type cg -svd_cross_st_pc_type jacobi -svd_max_it 1000
221 : test:
222 : suffix: 4_cyclic
223 : args: -svd_type cyclic -svd_cyclic_explicitmatrix
224 : test:
225 : suffix: 4_hpddm
226 : nsize: 4
227 : args: -svd_type trlanczos -svd_trlanczos_explicitmatrix -svd_trlanczos_pc_type hpddm
228 : args: -prefix_push svd_trlanczos_pc_hpddm_ -levels_1_st_share_sub_ksp -levels_1_eps_nev 10 -levels_1_eps_threshold 0.005 -levels_1_pc_asm_type basic -define_subdomains -levels_1_pc_asm_sub_mat_type sbaij -levels_1_sub_pc_type cholesky -prefix_pop
229 : requires: hpddm
230 :
231 : testset:
232 : args: -f1 ${SLEPC_DIR}/share/slepc/datafiles/matrices/rdb200.petsc -f2 identity -svd_threshold_relative 0.87 -terse
233 : output_file: output/ex48_5.out
234 : requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
235 : test:
236 : suffix: 5
237 : args: -svd_type trlanczos -svd_trlanczos_explicitmatrix -svd_trlanczos_scale 100 -svd_trlanczos_gbidiag {{lower upper single}}
238 : test:
239 : suffix: 5_cross
240 : args: -svd_type cross -svd_cross_explicitmatrix
241 :
242 : testset:
243 : args: -f1 ${DATAFILESPATH}/matrices/complex/qc324.petsc -f2 bidiagonal -svd_threshold_relative 0.3 -svd_ncv 8 -svd_tol 1e-7 -terse
244 : output_file: output/ex48_5_complex.out
245 : requires: double complex datafilespath !defined(PETSC_USE_64BIT_INDICES)
246 : test:
247 : suffix: 5_complex
248 : args: -svd_type trlanczos -svd_trlanczos_explicitmatrix -svd_trlanczos_scale 100 -svd_trlanczos_ksp_rtol 1e-6 -svd_trlanczos_gbidiag {{lower upper single}}
249 : test:
250 : suffix: 5_complex_cross
251 : args: -svd_type cross -svd_cross_explicitmatrix
252 :
253 : TEST*/
|
