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[] = "Partial hyperbolic singular value decomposition (HSVD) from a file.\n"
12 : "The command line options are:\n"
13 : " -file <filename>, PETSc binary file containing matrix A.\n"
14 : " -p <p>, where <p> = number of -1's in signature.\n"
15 : " -transpose, to transpose the matrix before doing the computation.\n\n";
16 :
17 : #include <slepcsvd.h>
18 :
19 25 : int main(int argc,char **argv)
20 : {
21 25 : Mat A,Omega; /* operator matrix, signature matrix */
22 25 : SVD svd; /* singular value problem solver context */
23 25 : Mat At;
24 25 : Vec u,v,vomega,*U,*V;
25 25 : MatType Atype;
26 25 : PetscReal tol,lev1=0.0,lev2=0.0;
27 25 : PetscInt M,N,p=0,i,Istart,Iend,nconv,nsv;
28 25 : char filename[PETSC_MAX_PATH_LEN];
29 25 : PetscViewer viewer;
30 25 : PetscBool flg,terse,skiporth=PETSC_FALSE,transpose=PETSC_FALSE;
31 :
32 25 : PetscFunctionBeginUser;
33 25 : PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
34 :
35 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
36 : Load matrix that defines the hyperbolic singular value problem
37 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
38 :
39 25 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nHyperbolic singular value problem stored in file.\n\n"));
40 25 : PetscCall(PetscOptionsGetString(NULL,NULL,"-file",filename,sizeof(filename),&flg));
41 25 : PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Must indicate a file name for matrix A with the -file option");
42 :
43 : #if defined(PETSC_USE_COMPLEX)
44 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Reading COMPLEX matrix from a binary file...\n"));
45 : #else
46 25 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Reading REAL matrix from a binary file...\n"));
47 : #endif
48 25 : PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer));
49 25 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
50 25 : PetscCall(MatSetFromOptions(A));
51 25 : PetscCall(MatLoad(A,viewer));
52 25 : PetscCall(PetscViewerDestroy(&viewer));
53 :
54 : /* transpose the matrix if requested */
55 25 : PetscCall(PetscOptionsGetBool(NULL,NULL,"-transpose",&transpose,NULL));
56 25 : if (transpose) {
57 10 : PetscCall(MatHermitianTranspose(A,MAT_INITIAL_MATRIX,&At));
58 10 : PetscCall(MatDestroy(&A));
59 10 : A = At;
60 : }
61 :
62 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
63 : Create the signature
64 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
65 :
66 25 : PetscCall(MatGetSize(A,&M,&N));
67 25 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-p",&p,&flg));
68 25 : PetscCheck(p>=0,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Parameter p cannot be negative");
69 25 : PetscCheck(p<M,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Parameter p cannot be larger than the number of rows of A");
70 25 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n Matrix dimensions: %" PetscInt_FMT "x%" PetscInt_FMT,M,N));
71 25 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,", using signature Omega=blkdiag(-I_%" PetscInt_FMT ",I_%" PetscInt_FMT ")\n\n",p,M-p));
72 :
73 25 : PetscCall(MatCreateVecs(A,NULL,&vomega));
74 25 : PetscCall(VecSet(vomega,1.0));
75 25 : PetscCall(VecGetOwnershipRange(vomega,&Istart,&Iend));
76 13722 : for (i=Istart;i<Iend;i++) {
77 13697 : if (i<p) PetscCall(VecSetValue(vomega,i,-1.0,INSERT_VALUES));
78 : }
79 25 : PetscCall(VecAssemblyBegin(vomega));
80 25 : PetscCall(VecAssemblyEnd(vomega));
81 :
82 25 : PetscCall(MatGetType(A,&Atype));
83 25 : PetscCall(MatCreate(PETSC_COMM_WORLD,&Omega));
84 25 : PetscCall(MatSetSizes(Omega,PETSC_DECIDE,PETSC_DECIDE,M,M));
85 25 : PetscCall(MatSetType(Omega,Atype));
86 25 : PetscCall(MatDiagonalSet(Omega,vomega,INSERT_VALUES));
87 :
88 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
89 : Create the singular value solver and set various options
90 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
91 :
92 25 : PetscCall(SVDCreate(PETSC_COMM_WORLD,&svd));
93 :
94 25 : PetscCall(SVDSetOperators(svd,A,NULL));
95 25 : PetscCall(SVDSetSignature(svd,vomega));
96 25 : PetscCall(SVDSetProblemType(svd,SVD_HYPERBOLIC));
97 :
98 25 : PetscCall(SVDSetFromOptions(svd));
99 :
100 25 : PetscCall(SVDIsHyperbolic(svd,&flg));
101 25 : PetscCheck(flg,PetscObjectComm((PetscObject)(svd)),PETSC_ERR_COR,"Problem should be hyperbolic");
102 :
103 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
104 : Solve the problem, display solution
105 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
106 :
107 25 : PetscCall(MatCreateVecs(A,&v,&u));
108 25 : PetscCall(VecSet(u,1.0));
109 25 : PetscCall(VecSet(v,1.0));
110 25 : PetscCall(SVDSolve(svd));
111 :
112 : /* show detailed info unless -terse option is given by user */
113 25 : PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
114 25 : if (terse) PetscCall(SVDErrorView(svd,SVD_ERROR_NORM,NULL));
115 : else {
116 0 : PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
117 0 : PetscCall(SVDConvergedReasonView(svd,PETSC_VIEWER_STDOUT_WORLD));
118 0 : PetscCall(SVDErrorView(svd,SVD_ERROR_NORM,PETSC_VIEWER_STDOUT_WORLD));
119 0 : PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
120 : }
121 :
122 : /* check orthogonality */
123 25 : PetscCall(PetscOptionsGetBool(NULL,NULL,"-skiporth",&skiporth,NULL));
124 25 : PetscCall(SVDGetConverged(svd,&nconv));
125 25 : PetscCall(SVDGetDimensions(svd,&nsv,NULL,NULL));
126 25 : if (nsv) nconv = PetscMin(nconv,nsv);
127 25 : if (nconv>0 && !skiporth) {
128 25 : PetscCall(SVDGetTolerances(svd,&tol,NULL));
129 25 : PetscCall(VecDuplicateVecs(u,nconv,&U));
130 25 : PetscCall(VecDuplicateVecs(v,nconv,&V));
131 657 : for (i=0;i<nconv;i++) PetscCall(SVDGetSingularTriplet(svd,i,NULL,U[i],V[i]));
132 25 : PetscCall(VecCheckOrthonormality(U,nconv,NULL,nconv,Omega,NULL,&lev1));
133 25 : PetscCall(VecCheckOrthonormality(V,nconv,NULL,nconv,NULL,NULL,&lev2));
134 25 : if (lev1+lev2<20*tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality below the tolerance\n"));
135 0 : else PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality: %g (U) %g (V)\n",(double)lev1,(double)lev2));
136 25 : PetscCall(VecDestroyVecs(nconv,&U));
137 25 : PetscCall(VecDestroyVecs(nconv,&V));
138 : }
139 25 : PetscCall(VecDestroy(&u));
140 25 : PetscCall(VecDestroy(&v));
141 :
142 : /* free work space */
143 25 : PetscCall(SVDDestroy(&svd));
144 25 : PetscCall(MatDestroy(&A));
145 25 : PetscCall(MatDestroy(&Omega));
146 25 : PetscCall(VecDestroy(&vomega));
147 25 : PetscCall(SlepcFinalize());
148 : return 0;
149 : }
150 :
151 : /*TEST
152 :
153 : testset:
154 : args: -file ${DATAFILESPATH}/matrices/real/illc1033.petsc -svd_nsv 62 -p 40 -terse
155 : requires: double !complex datafilespath !defined(PETSC_USE_64BIT_INDICES)
156 : filter: grep -v Reading
157 : output_file: output/ex52_1.out
158 : test:
159 : args: -svd_type cross -svd_cross_explicitmatrix {{0 1}} -svd_implicittranspose {{0 1}}
160 : suffix: 1_cross
161 : test:
162 : args: -svd_type cyclic -svd_cyclic_explicitmatrix {{0 1}} -svd_ncv 300
163 : suffix: 1_cyclic
164 : test:
165 : args: -svd_type trlanczos -svd_trlanczos_explicitmatrix {{0 1}}
166 : suffix: 1_trlanczos
167 : test:
168 : args: -svd_type lapack
169 : suffix: 1_lapack
170 :
171 : testset:
172 : args: -file ${DATAFILESPATH}/matrices/real/illc1033.petsc -transpose -svd_nsv 6 -p 130 -terse
173 : requires: double !complex datafilespath !defined(PETSC_USE_64BIT_INDICES)
174 : filter: grep -v Reading
175 : output_file: output/ex52_2.out
176 : test:
177 : args: -svd_type cross -svd_cross_explicitmatrix {{0 1}} -svd_implicittranspose {{0 1}} -svd_ncv 100
178 : suffix: 2_cross
179 : test:
180 : args: -svd_type cyclic -svd_cyclic_explicitmatrix {{0 1}} -svd_ncv 25
181 : suffix: 2_cyclic
182 : test:
183 : args: -svd_type trlanczos -svd_trlanczos_explicitmatrix {{0 1}} -svd_implicittranspose {{0 1}}
184 : suffix: 2_trlanczos
185 :
186 : testset:
187 : args: -file ${DATAFILESPATH}/matrices/complex/illc1033.petsc -svd_nsv 62 -p 40 -terse
188 : requires: double complex datafilespath !defined(PETSC_USE_64BIT_INDICES)
189 : filter: grep -v Reading
190 : output_file: output/ex52_1.out
191 : test:
192 : args: -svd_type cross -svd_cross_explicitmatrix {{0 1}}
193 : suffix: 3_cross
194 : test:
195 : args: -svd_type cyclic -svd_cyclic_explicitmatrix {{0 1}} -svd_cyclic_bv_definite_tol 1e-13 -svd_cyclic_st_ksp_type gcr -svd_cyclic_st_pc_type jacobi -svd_ncv 250
196 : suffix: 3_cyclic
197 : test:
198 : args: -svd_type trlanczos -svd_trlanczos_explicitmatrix {{0 1}} -bv_definite_tol 1e-13
199 : suffix: 3_trlanczos
200 : test:
201 : args: -svd_type lapack
202 : suffix: 3_lapack
203 :
204 : testset:
205 : args: -file ${DATAFILESPATH}/matrices/complex/illc1033.petsc -transpose -svd_nsv 6 -p 130 -terse
206 : requires: double complex datafilespath !defined(PETSC_USE_64BIT_INDICES)
207 : filter: grep -v Reading
208 : output_file: output/ex52_2.out
209 : test:
210 : args: -svd_type cross -svd_cross_explicitmatrix {{0 1}} -svd_ncv 100 -svd_cross_bv_definite_tol 1e-14
211 : suffix: 4_cross
212 : test:
213 : args: -svd_type cyclic -svd_cyclic_explicitmatrix {{0 1}} -svd_ncv 26
214 : suffix: 4_cyclic
215 : test:
216 : args: -svd_type trlanczos -svd_trlanczos_explicitmatrix {{0 1}}
217 : suffix: 4_trlanczos
218 :
219 : testset:
220 : args: -file ${SLEPC_DIR}/share/slepc/datafiles/matrices/rdb200.petsc -svd_smallest -svd_nsv 3 -p 1 -terse
221 : requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
222 : filter: grep -v Reading
223 : output_file: output/ex52_5.out
224 : test:
225 : args: -svd_type cross -svd_max_it 1000 -svd_cross_bv_definite_tol 1e-14
226 : suffix: 5_cross
227 : test:
228 : args: -svd_type cyclic -svd_max_it 4000 -svd_ncv 32 -svd_cyclic_st_ksp_type preonly -svd_cyclic_st_pc_type jacobi -svd_cyclic_bv_definite_tol 1e-14
229 : suffix: 5_cyclic
230 : test:
231 : args: -svd_type trlanczos -svd_max_it 4000 -svd_ncv 28 -bv_definite_tol 1e-14
232 : suffix: 5_trlanczos
233 : test:
234 : args: -svd_type lapack
235 : suffix: 5_lapack
236 :
237 : testset:
238 : args: -svd_type {{trlanczos cross}} -terse -p 100
239 : test:
240 : suffix: 6
241 : filter: sed -e "s/27.29445, 27.29445/27.29445/"
242 : args: -file ${SLEPC_DIR}/share/slepc/datafiles/matrices/rdb200.petsc -svd_threshold_relative 0.8
243 : requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
244 : test:
245 : suffix: 6_complex
246 : args: -file ${DATAFILESPATH}/matrices/complex/qc324.petsc -svd_threshold_relative 0.6
247 : requires: double complex datafilespath !defined(PETSC_USE_64BIT_INDICES)
248 :
249 : TEST*/
|
