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[] = "Test the solution of a SVD without calling SVDSetFromOptions (based on ex8.c).\n\n"
12 : "The command line options are:\n"
13 : " -n <n>, where <n> = matrix dimension.\n"
14 : " -type <svd_type> = svd type to test.\n\n";
15 :
16 : #include <slepcsvd.h>
17 :
18 : /*
19 : This example computes the singular values of an nxn Grcar matrix,
20 : which is a nonsymmetric Toeplitz matrix:
21 :
22 : | 1 1 1 1 |
23 : | -1 1 1 1 1 |
24 : | -1 1 1 1 1 |
25 : | . . . . . |
26 : A = | . . . . . |
27 : | -1 1 1 1 1 |
28 : | -1 1 1 1 |
29 : | -1 1 1 |
30 : | -1 1 |
31 :
32 : */
33 :
34 8 : int main(int argc,char **argv)
35 : {
36 8 : Mat A; /* Grcar matrix */
37 8 : SVD svd; /* singular value solver context */
38 8 : PetscInt N=30,Istart,Iend,i,col[5],nconv1,nconv2;
39 8 : PetscScalar value[] = { -1, 1, 1, 1, 1 };
40 8 : PetscReal sigma_1,sigma_n;
41 8 : char svdtype[30] = "cross",epstype[30] = "";
42 8 : PetscBool flg;
43 8 : EPS eps;
44 :
45 8 : PetscFunctionBeginUser;
46 8 : PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
47 :
48 8 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&N,NULL));
49 8 : PetscCall(PetscOptionsGetString(NULL,NULL,"-type",svdtype,sizeof(svdtype),NULL));
50 8 : PetscCall(PetscOptionsGetString(NULL,NULL,"-epstype",epstype,sizeof(epstype),&flg));
51 8 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nEstimate the condition number of a Grcar matrix, n=%" PetscInt_FMT,N));
52 8 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n\n"));
53 :
54 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
55 : Generate the matrix
56 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
57 :
58 8 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
59 8 : PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N));
60 8 : PetscCall(MatSetFromOptions(A));
61 :
62 8 : PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
63 248 : for (i=Istart;i<Iend;i++) {
64 240 : col[0]=i-1; col[1]=i; col[2]=i+1; col[3]=i+2; col[4]=i+3;
65 240 : if (i==0) PetscCall(MatSetValues(A,1,&i,4,col+1,value+1,INSERT_VALUES));
66 240 : else PetscCall(MatSetValues(A,1,&i,PetscMin(5,N-i+1),col,value,INSERT_VALUES));
67 : }
68 :
69 8 : PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
70 8 : PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
71 :
72 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
73 : Create the singular value solver and set the solution method
74 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
75 :
76 : /*
77 : Create singular value context
78 : */
79 8 : PetscCall(SVDCreate(PETSC_COMM_WORLD,&svd));
80 :
81 : /*
82 : Set operator
83 : */
84 8 : PetscCall(SVDSetOperators(svd,A,NULL));
85 :
86 : /*
87 : Set solver parameters at runtime
88 : */
89 8 : PetscCall(SVDSetType(svd,svdtype));
90 8 : if (flg) {
91 2 : PetscCall(PetscObjectTypeCompare((PetscObject)svd,SVDCROSS,&flg));
92 2 : if (flg) {
93 1 : PetscCall(SVDCrossGetEPS(svd,&eps));
94 1 : PetscCall(EPSSetType(eps,epstype));
95 : }
96 2 : PetscCall(PetscObjectTypeCompare((PetscObject)svd,SVDCYCLIC,&flg));
97 2 : if (flg) {
98 1 : PetscCall(SVDCyclicGetEPS(svd,&eps));
99 1 : PetscCall(EPSSetType(eps,epstype));
100 : }
101 : }
102 8 : PetscCall(SVDSetDimensions(svd,1,PETSC_DETERMINE,PETSC_DETERMINE));
103 8 : PetscCall(SVDSetTolerances(svd,1e-6,1000));
104 :
105 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106 : Compute the singular values
107 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
108 :
109 : /*
110 : First request the largest singular value
111 : */
112 8 : PetscCall(SVDSetWhichSingularTriplets(svd,SVD_LARGEST));
113 8 : PetscCall(SVDSolve(svd));
114 : /*
115 : Get number of converged singular values
116 : */
117 8 : PetscCall(SVDGetConverged(svd,&nconv1));
118 : /*
119 : Get converged singular values: largest singular value is stored in sigma_1.
120 : In this example, we are not interested in the singular vectors
121 : */
122 8 : if (nconv1 > 0) PetscCall(SVDGetSingularTriplet(svd,0,&sigma_1,NULL,NULL));
123 0 : else PetscCall(PetscPrintf(PETSC_COMM_WORLD," Unable to compute large singular value!\n\n"));
124 :
125 : /*
126 : Request the smallest singular value
127 : */
128 8 : PetscCall(SVDSetWhichSingularTriplets(svd,SVD_SMALLEST));
129 8 : PetscCall(SVDSolve(svd));
130 : /*
131 : Get number of converged triplets
132 : */
133 8 : PetscCall(SVDGetConverged(svd,&nconv2));
134 : /*
135 : Get converged singular values: smallest singular value is stored in sigma_n.
136 : As before, we are not interested in the singular vectors
137 : */
138 8 : if (nconv2 > 0) PetscCall(SVDGetSingularTriplet(svd,0,&sigma_n,NULL,NULL));
139 0 : else PetscCall(PetscPrintf(PETSC_COMM_WORLD," Unable to compute small singular value!\n\n"));
140 :
141 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
142 : Display solution and clean up
143 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
144 8 : if (nconv1 > 0 && nconv2 > 0) {
145 8 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Computed singular values: sigma_1=%.4f, sigma_n=%.4f\n",(double)sigma_1,(double)sigma_n));
146 8 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Estimated condition number: sigma_1/sigma_n=%.4f\n\n",(double)(sigma_1/sigma_n)));
147 : }
148 :
149 : /*
150 : Free work space
151 : */
152 8 : PetscCall(SVDDestroy(&svd));
153 8 : PetscCall(MatDestroy(&A));
154 8 : PetscCall(SlepcFinalize());
155 : return 0;
156 : }
157 :
158 : /*TEST
159 :
160 : test:
161 : suffix: 1
162 : args: -type {{lanczos trlanczos cross cyclic lapack}}
163 :
164 : test:
165 : suffix: 1_cross_gd
166 : args: -type cross -epstype gd
167 : output_file: output/test1_1.out
168 :
169 : test:
170 : suffix: 1_cyclic_gd
171 : args: -type cyclic -epstype gd
172 : output_file: output/test1_1.out
173 : requires: !single
174 :
175 : test:
176 : suffix: 1_primme
177 : args: -type primme
178 : requires: primme
179 : output_file: output/test1_1.out
180 :
181 : TEST*/
|
