| 1294 |
slepc |
1 |
|
|
|
2 |
static char help[] = "Solves a singular value problem with the matrix loaded from a file.\n"
|
|
|
3 |
"This example works for both real and complex numbers.\n\n"
|
|
|
4 |
"The command line options are:\n"
|
|
|
5 |
" -file <filename>, where <filename> = matrix file in PETSc binary form.\n\n";
|
|
|
6 |
|
|
|
7 |
#include "slepcsvd.h"
|
|
|
8 |
|
|
|
9 |
#undef __FUNCT__
|
|
|
10 |
#define __FUNCT__ "main"
|
|
|
11 |
int main( int argc, char **argv )
|
|
|
12 |
{
|
|
|
13 |
Mat A; /* operator matrix */
|
|
|
14 |
SVD svd; /* singular value problem solver context */
|
|
|
15 |
SVDType type;
|
| 1317 |
slepc |
16 |
PetscReal error, tol, sigma;
|
| 1294 |
slepc |
17 |
PetscErrorCode ierr;
|
|
|
18 |
int nsv, maxit, i, its, nconv;
|
|
|
19 |
char filename[256];
|
|
|
20 |
PetscViewer viewer;
|
|
|
21 |
PetscTruth flg;
|
|
|
22 |
|
|
|
23 |
|
|
|
24 |
SlepcInitialize(&argc,&argv,(char*)0,help);
|
|
|
25 |
|
|
|
26 |
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
|
27 |
Load the operator matrix that defines the singular value problem
|
|
|
28 |
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
|
|
|
29 |
|
|
|
30 |
ierr = PetscPrintf(PETSC_COMM_WORLD,"\nSingular value problem stored in file.\n\n");CHKERRQ(ierr);
|
|
|
31 |
ierr = PetscOptionsGetString(PETSC_NULL,"-file",filename,256,&flg);CHKERRQ(ierr);
|
|
|
32 |
if (!flg) {
|
|
|
33 |
SETERRQ(1,"Must indicate a file name with the -file option.");
|
|
|
34 |
}
|
|
|
35 |
|
|
|
36 |
#if defined(PETSC_USE_COMPLEX)
|
|
|
37 |
ierr = PetscPrintf(PETSC_COMM_WORLD," Reading COMPLEX matrix from a binary file...\n");CHKERRQ(ierr);
|
|
|
38 |
#else
|
|
|
39 |
ierr = PetscPrintf(PETSC_COMM_WORLD," Reading REAL matrix from a binary file...\n");CHKERRQ(ierr);
|
|
|
40 |
#endif
|
|
|
41 |
ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);CHKERRQ(ierr);
|
|
|
42 |
ierr = MatLoad(viewer,MATAIJ,&A);CHKERRQ(ierr);
|
|
|
43 |
ierr = PetscViewerDestroy(viewer);CHKERRQ(ierr);
|
|
|
44 |
|
|
|
45 |
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
|
46 |
Create the singular value solver and set various options
|
|
|
47 |
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
|
|
|
48 |
|
|
|
49 |
/*
|
|
|
50 |
Create singular value solver context
|
|
|
51 |
*/
|
|
|
52 |
ierr = SVDCreate(PETSC_COMM_WORLD,&svd);CHKERRQ(ierr);
|
|
|
53 |
|
|
|
54 |
/*
|
|
|
55 |
Set operator
|
|
|
56 |
*/
|
|
|
57 |
ierr = SVDSetOperator(svd,A);CHKERRQ(ierr);
|
|
|
58 |
|
|
|
59 |
/*
|
|
|
60 |
Set solver parameters at runtime
|
|
|
61 |
*/
|
|
|
62 |
ierr = SVDSetFromOptions(svd);CHKERRQ(ierr);
|
|
|
63 |
|
|
|
64 |
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
|
65 |
Solve the singular value system
|
|
|
66 |
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
|
|
|
67 |
|
|
|
68 |
ierr = SVDSolve(svd);CHKERRQ(ierr);
|
|
|
69 |
ierr = SVDGetIterationNumber(svd, &its);CHKERRQ(ierr);
|
|
|
70 |
ierr = PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %d\n",its);CHKERRQ(ierr);
|
|
|
71 |
|
|
|
72 |
/*
|
|
|
73 |
Optional: Get some information from the solver and display it
|
|
|
74 |
*/
|
|
|
75 |
ierr = SVDGetType(svd,&type);CHKERRQ(ierr);
|
|
|
76 |
ierr = PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);CHKERRQ(ierr);
|
|
|
77 |
ierr = SVDGetDimensions(svd,&nsv,PETSC_NULL);CHKERRQ(ierr);
|
|
|
78 |
ierr = PetscPrintf(PETSC_COMM_WORLD," Number of requested singular values: %d\n",nsv);CHKERRQ(ierr);
|
|
|
79 |
ierr = SVDGetTolerances(svd,&tol,&maxit);CHKERRQ(ierr);
|
|
|
80 |
ierr = PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%d\n",tol,maxit);CHKERRQ(ierr);
|
|
|
81 |
|
|
|
82 |
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
|
83 |
Display solution and clean up
|
|
|
84 |
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
|
|
|
85 |
|
|
|
86 |
/*
|
|
|
87 |
Get number of converged singular triplets
|
|
|
88 |
*/
|
|
|
89 |
ierr = SVDGetConverged(svd,&nconv);CHKERRQ(ierr);
|
|
|
90 |
ierr = PetscPrintf(PETSC_COMM_WORLD," Number of converged approximate singular triplets: %d\n\n",nconv);CHKERRQ(ierr);
|
|
|
91 |
|
|
|
92 |
if (nconv>0) {
|
|
|
93 |
/*
|
|
|
94 |
Display singular values and relative errors
|
|
|
95 |
*/
|
|
|
96 |
ierr = PetscPrintf(PETSC_COMM_WORLD,
|
|
|
97 |
" sigma residual norm\n"
|
|
|
98 |
" --------------------- ------------------\n" );CHKERRQ(ierr);
|
|
|
99 |
for( i=0; i<nconv; i++ ) {
|
|
|
100 |
/*
|
|
|
101 |
Get converged singular triplets: i-th singular value is stored in sigma
|
|
|
102 |
*/
|
|
|
103 |
ierr = SVDGetSingularTriplet(svd,i,&sigma,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
|
|
|
104 |
|
|
|
105 |
/*
|
|
|
106 |
Compute the error associated to each singular triplet
|
|
|
107 |
*/
|
| 1317 |
slepc |
108 |
ierr = SVDComputeRelativeError(svd,i,&error);CHKERRQ(ierr);
|
| 1294 |
slepc |
109 |
ierr = PetscPrintf(PETSC_COMM_WORLD," % 6f ",sigma); CHKERRQ(ierr);
|
|
|
110 |
ierr = PetscPrintf(PETSC_COMM_WORLD," % 12g\n",error);CHKERRQ(ierr);
|
|
|
111 |
}
|
|
|
112 |
ierr = PetscPrintf(PETSC_COMM_WORLD,"\n" );CHKERRQ(ierr);
|
|
|
113 |
}
|
|
|
114 |
|
|
|
115 |
/*
|
|
|
116 |
Free work space
|
|
|
117 |
*/
|
|
|
118 |
ierr = SVDDestroy(svd);CHKERRQ(ierr);
|
|
|
119 |
ierr = MatDestroy(A);CHKERRQ(ierr);
|
|
|
120 |
ierr = SlepcFinalize();CHKERRQ(ierr);
|
|
|
121 |
return 0;
|
|
|
122 |
}
|