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 16 : int main(int argc,char **argv)
21 : {
22 16 : Mat A,B; /* matrices */
23 16 : SVD svd; /* singular value problem solver context */
24 16 : PetscInt i,m,n,p,Istart,Iend,col[3];
25 16 : PetscScalar vals[3];
26 16 : char filename[PETSC_MAX_PATH_LEN];
27 16 : PetscViewer viewer;
28 16 : PetscBool flg,terse;
29 :
30 16 : PetscFunctionBeginUser;
31 16 : PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
32 :
33 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
34 : Load matrices that define the generalized singular value problem
35 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
36 :
37 16 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nGeneralized singular value problem stored in file.\n\n"));
38 16 : PetscCall(PetscOptionsGetString(NULL,NULL,"-f1",filename,sizeof(filename),&flg));
39 16 : 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 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Reading COMPLEX matrix from a binary file...\n"));
43 : #else
44 16 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Reading REAL matrix from a binary file...\n"));
45 : #endif
46 16 : PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer));
47 16 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
48 16 : PetscCall(MatSetFromOptions(A));
49 16 : PetscCall(MatLoad(A,viewer));
50 16 : PetscCall(PetscViewerDestroy(&viewer));
51 :
52 16 : PetscCall(MatGetSize(A,&m,&n));
53 :
54 16 : PetscCall(PetscOptionsGetString(NULL,NULL,"-f2",filename,sizeof(filename),&flg));
55 16 : 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 16 : PetscCall(PetscStrcmp(filename,"identity",&flg));
57 16 : if (flg) {
58 8 : p = n;
59 8 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-p",&p,&flg));
60 8 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Using B=I with p=%" PetscInt_FMT "\n\n",p));
61 8 : PetscCall(MatCreate(PETSC_COMM_WORLD,&B));
62 8 : PetscCall(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,p,n));
63 8 : PetscCall(MatSetFromOptions(B));
64 8 : PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY));
65 8 : PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY));
66 8 : PetscCall(MatShift(B,1.0));
67 : } else {
68 8 : PetscCall(PetscStrcmp(filename,"bidiagonal",&flg));
69 8 : if (flg) {
70 0 : p = n+1;
71 0 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-p",&p,&flg));
72 0 : vals[0]=-1; vals[1]=1;
73 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Using B=bidiag(1,-1) with p=%" PetscInt_FMT "\n\n",p));
74 0 : PetscCall(MatCreate(PETSC_COMM_WORLD,&B));
75 0 : PetscCall(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,p,n));
76 0 : PetscCall(MatSetFromOptions(B));
77 0 : PetscCall(MatGetOwnershipRange(B,&Istart,&Iend));
78 0 : for (i=Istart;i<Iend;i++) {
79 0 : col[0]=i-1; col[1]=i;
80 0 : if (i==0) PetscCall(MatSetValue(B,i,col[1],vals[1],INSERT_VALUES));
81 0 : else if (i<n) PetscCall(MatSetValues(B,1,&i,2,col,vals,INSERT_VALUES));
82 0 : else if (i==n) PetscCall(MatSetValue(B,i,col[0],vals[0],INSERT_VALUES));
83 : }
84 0 : PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY));
85 0 : PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY));
86 : } else {
87 8 : PetscCall(PetscStrcmp(filename,"tridiagonal",&flg));
88 8 : 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 8 : PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer));
105 8 : PetscCall(MatCreate(PETSC_COMM_WORLD,&B));
106 8 : PetscCall(MatSetFromOptions(B));
107 8 : PetscCall(MatLoad(B,viewer));
108 8 : 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 16 : PetscCall(SVDCreate(PETSC_COMM_WORLD,&svd));
121 :
122 : /*
123 : Set operators of GSVD problem
124 : */
125 16 : PetscCall(SVDSetOperators(svd,A,B));
126 16 : PetscCall(SVDSetProblemType(svd,SVD_GENERALIZED));
127 :
128 : /*
129 : Set solver parameters at runtime
130 : */
131 16 : PetscCall(SVDSetFromOptions(svd));
132 :
133 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
134 : Solve the problem and print solution
135 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
136 :
137 16 : PetscCall(SVDSolve(svd));
138 :
139 : /* show detailed info unless -terse option is given by user */
140 16 : PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
141 16 : 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 16 : PetscCall(SVDDestroy(&svd));
149 16 : PetscCall(MatDestroy(&A));
150 16 : PetscCall(MatDestroy(&B));
151 16 : 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 : TEST*/
|
