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[] = "Singular value decomposition of the Lauchli matrix.\n"
12 : "The command line options are:\n"
13 : " -n <n>, where <n> = matrix dimension.\n"
14 : " -mu <mu>, where <mu> = subdiagonal value.\n\n";
15 :
16 : #include <slepcsvd.h>
17 :
18 1 : int main(int argc,char **argv)
19 : {
20 1 : Mat A; /* operator matrix */
21 1 : Vec u,v; /* left and right singular vectors */
22 1 : SVD svd; /* singular value problem solver context */
23 1 : SVDType type;
24 1 : PetscReal error,tol,sigma,mu=PETSC_SQRT_MACHINE_EPSILON;
25 1 : PetscInt n=100,i,j,Istart,Iend,nsv,maxit,its,nconv;
26 :
27 1 : PetscFunctionBeginUser;
28 1 : PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
29 :
30 1 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
31 1 : PetscCall(PetscOptionsGetReal(NULL,NULL,"-mu",&mu,NULL));
32 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nLauchli singular value decomposition, (%" PetscInt_FMT " x %" PetscInt_FMT ") mu=%g\n\n",n+1,n,(double)mu));
33 :
34 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
35 : Build the Lauchli matrix
36 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
37 :
38 1 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
39 1 : PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n+1,n));
40 1 : PetscCall(MatSetFromOptions(A));
41 :
42 1 : PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
43 102 : for (i=Istart;i<Iend;i++) {
44 101 : if (i == 0) {
45 101 : for (j=0;j<n;j++) PetscCall(MatSetValue(A,0,j,1.0,INSERT_VALUES));
46 101 : } else PetscCall(MatSetValue(A,i,i-1,mu,INSERT_VALUES));
47 : }
48 :
49 1 : PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
50 1 : PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
51 1 : PetscCall(MatCreateVecs(A,&v,&u));
52 :
53 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
54 : Create the singular value solver and set various options
55 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
56 :
57 : /*
58 : Create singular value solver context
59 : */
60 1 : PetscCall(SVDCreate(PETSC_COMM_WORLD,&svd));
61 :
62 : /*
63 : Set operators and problem type
64 : */
65 1 : PetscCall(SVDSetOperators(svd,A,NULL));
66 1 : PetscCall(SVDSetProblemType(svd,SVD_STANDARD));
67 :
68 : /*
69 : Use thick-restart Lanczos as default solver
70 : */
71 1 : PetscCall(SVDSetType(svd,SVDTRLANCZOS));
72 :
73 : /*
74 : Set solver parameters at runtime
75 : */
76 1 : PetscCall(SVDSetFromOptions(svd));
77 :
78 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
79 : Solve the singular value system
80 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
81 :
82 1 : PetscCall(SVDSolve(svd));
83 1 : PetscCall(SVDGetIterationNumber(svd,&its));
84 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %" PetscInt_FMT "\n",its));
85 :
86 : /*
87 : Optional: Get some information from the solver and display it
88 : */
89 1 : PetscCall(SVDGetType(svd,&type));
90 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type));
91 1 : PetscCall(SVDGetDimensions(svd,&nsv,NULL,NULL));
92 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of requested singular values: %" PetscInt_FMT "\n",nsv));
93 1 : PetscCall(SVDGetTolerances(svd,&tol,&maxit));
94 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%" PetscInt_FMT "\n",(double)tol,maxit));
95 :
96 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
97 : Display solution and clean up
98 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
99 :
100 : /*
101 : Get number of converged singular triplets
102 : */
103 1 : PetscCall(SVDGetConverged(svd,&nconv));
104 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of converged approximate singular triplets: %" PetscInt_FMT "\n\n",nconv));
105 :
106 1 : if (nconv>0) {
107 : /*
108 : Display singular values and relative errors
109 : */
110 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,
111 : " sigma relative error\n"
112 : " --------------------- ------------------\n"));
113 11 : for (i=0;i<nconv;i++) {
114 : /*
115 : Get converged singular triplets: i-th singular value is stored in sigma
116 : */
117 10 : PetscCall(SVDGetSingularTriplet(svd,i,&sigma,u,v));
118 :
119 : /*
120 : Compute the error associated to each singular triplet
121 : */
122 10 : PetscCall(SVDComputeError(svd,i,SVD_ERROR_RELATIVE,&error));
123 :
124 10 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," % 6f ",(double)sigma));
125 10 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," % 12g\n",(double)error));
126 : }
127 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n"));
128 : }
129 :
130 : /*
131 : Free work space
132 : */
133 1 : PetscCall(SVDDestroy(&svd));
134 1 : PetscCall(MatDestroy(&A));
135 1 : PetscCall(VecDestroy(&u));
136 1 : PetscCall(VecDestroy(&v));
137 1 : PetscCall(SlepcFinalize());
138 : return 0;
139 : }
140 :
141 : /*TEST
142 :
143 : testset:
144 : filter: sed -e "s/[0-9]\.[0-9]*e[+-]\([0-9]*\)/removed/g"
145 : requires: double
146 : test:
147 : suffix: 1
148 : test:
149 : suffix: 1_scalapack
150 : nsize: {{1 2}}
151 : args: -svd_type scalapack
152 : requires: scalapack
153 :
154 : TEST*/
|
