Actual source code: ex52.c
slepc-3.21.1 2024-04-26
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
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";
17: #include <slepcsvd.h>
19: int main(int argc,char **argv)
20: {
21: Mat A,Omega; /* operator matrix, signature matrix */
22: SVD svd; /* singular value problem solver context */
23: Mat At;
24: Vec u,v,vomega,*U,*V;
25: MatType Atype;
26: PetscReal tol,lev1=0.0,lev2=0.0;
27: PetscInt M,N,p=0,i,Istart,Iend,nconv,nsv;
28: char filename[PETSC_MAX_PATH_LEN];
29: PetscViewer viewer;
30: PetscBool flg,terse,skiporth=PETSC_FALSE,transpose=PETSC_FALSE;
32: PetscFunctionBeginUser;
33: PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
35: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
36: Load matrix that defines the hyperbolic singular value problem
37: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
39: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nHyperbolic singular value problem stored in file.\n\n"));
40: PetscCall(PetscOptionsGetString(NULL,NULL,"-file",filename,sizeof(filename),&flg));
41: PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Must indicate a file name for matrix A with the -file option");
43: #if defined(PETSC_USE_COMPLEX)
44: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Reading COMPLEX matrix from a binary file...\n"));
45: #else
46: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Reading REAL matrix from a binary file...\n"));
47: #endif
48: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer));
49: PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
50: PetscCall(MatSetFromOptions(A));
51: PetscCall(MatLoad(A,viewer));
52: PetscCall(PetscViewerDestroy(&viewer));
54: /* transpose the matrix if requested */
55: PetscCall(PetscOptionsGetBool(NULL,NULL,"-transpose",&transpose,NULL));
56: if (transpose) {
57: PetscCall(MatHermitianTranspose(A,MAT_INITIAL_MATRIX,&At));
58: PetscCall(MatDestroy(&A));
59: A = At;
60: }
62: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
63: Create the signature
64: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
66: PetscCall(MatGetSize(A,&M,&N));
67: PetscCall(PetscOptionsGetInt(NULL,NULL,"-p",&p,&flg));
68: PetscCheck(p>=0,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Parameter p cannot be negative");
69: PetscCheck(p<M,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Parameter p cannot be larger than the number of rows of A");
70: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n Matrix dimensions: %" PetscInt_FMT "x%" PetscInt_FMT,M,N));
71: PetscCall(PetscPrintf(PETSC_COMM_WORLD,", using signature Omega=blkdiag(-I_%" PetscInt_FMT ",I_%" PetscInt_FMT ")\n\n",p,M-p));
73: PetscCall(MatCreateVecs(A,NULL,&vomega));
74: PetscCall(VecSet(vomega,1.0));
75: PetscCall(VecGetOwnershipRange(vomega,&Istart,&Iend));
76: for (i=Istart;i<Iend;i++) {
77: if (i<p) PetscCall(VecSetValue(vomega,i,-1.0,INSERT_VALUES));
78: }
79: PetscCall(VecAssemblyBegin(vomega));
80: PetscCall(VecAssemblyEnd(vomega));
82: PetscCall(MatGetType(A,&Atype));
83: PetscCall(MatCreate(PETSC_COMM_WORLD,&Omega));
84: PetscCall(MatSetSizes(Omega,PETSC_DECIDE,PETSC_DECIDE,M,M));
85: PetscCall(MatSetType(Omega,Atype));
86: PetscCall(MatDiagonalSet(Omega,vomega,INSERT_VALUES));
88: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
89: Create the singular value solver and set various options
90: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
92: PetscCall(SVDCreate(PETSC_COMM_WORLD,&svd));
94: PetscCall(SVDSetOperators(svd,A,NULL));
95: PetscCall(SVDSetSignature(svd,vomega));
96: PetscCall(SVDSetProblemType(svd,SVD_HYPERBOLIC));
98: PetscCall(SVDSetFromOptions(svd));
100: PetscCall(SVDIsHyperbolic(svd,&flg));
101: PetscCheck(flg,PetscObjectComm((PetscObject)(svd)),PETSC_ERR_COR,"Problem should be hyperbolic");
103: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
104: Solve the problem, display solution
105: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
107: PetscCall(MatCreateVecs(A,&v,&u));
108: PetscCall(VecSet(u,1.0));
109: PetscCall(VecSet(v,1.0));
110: PetscCall(SVDSolve(svd));
112: /* show detailed info unless -terse option is given by user */
113: PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
114: if (terse) PetscCall(SVDErrorView(svd,SVD_ERROR_NORM,NULL));
115: else {
116: PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
117: PetscCall(SVDConvergedReasonView(svd,PETSC_VIEWER_STDOUT_WORLD));
118: PetscCall(SVDErrorView(svd,SVD_ERROR_NORM,PETSC_VIEWER_STDOUT_WORLD));
119: PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
120: }
122: /* check orthogonality */
123: PetscCall(PetscOptionsGetBool(NULL,NULL,"-skiporth",&skiporth,NULL));
124: PetscCall(SVDGetConverged(svd,&nconv));
125: PetscCall(SVDGetDimensions(svd,&nsv,NULL,NULL));
126: nconv = PetscMin(nconv,nsv);
127: if (nconv>0 && !skiporth) {
128: PetscCall(SVDGetTolerances(svd,&tol,NULL));
129: PetscCall(VecDuplicateVecs(u,nconv,&U));
130: PetscCall(VecDuplicateVecs(v,nconv,&V));
131: for (i=0;i<nconv;i++) PetscCall(SVDGetSingularTriplet(svd,i,NULL,U[i],V[i]));
132: PetscCall(VecCheckOrthonormality(U,nconv,NULL,nconv,Omega,NULL,&lev1));
133: PetscCall(VecCheckOrthonormality(V,nconv,NULL,nconv,NULL,NULL,&lev2));
134: if (lev1+lev2<20*tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality below the tolerance\n"));
135: else PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Level of orthogonality: %g (U) %g (V)\n",(double)lev1,(double)lev2));
136: PetscCall(VecDestroyVecs(nconv,&U));
137: PetscCall(VecDestroyVecs(nconv,&V));
138: }
139: PetscCall(VecDestroy(&u));
140: PetscCall(VecDestroy(&v));
142: /* free work space */
143: PetscCall(SVDDestroy(&svd));
144: PetscCall(MatDestroy(&A));
145: PetscCall(MatDestroy(&Omega));
146: PetscCall(VecDestroy(&vomega));
147: PetscCall(SlepcFinalize());
148: return 0;
149: }
151: /*TEST
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
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
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
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
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 30 -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
237: TEST*/