Actual source code: ex8.c
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2010, Universidad Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7:
8: SLEPc is free software: you can redistribute it and/or modify it under the
9: terms of version 3 of the GNU Lesser General Public License as published by
10: the Free Software Foundation.
12: SLEPc is distributed in the hope that it will be useful, but WITHOUT ANY
13: WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14: FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
15: more details.
17: You should have received a copy of the GNU Lesser General Public License
18: along with SLEPc. If not, see <http://www.gnu.org/licenses/>.
19: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
20: */
22: static char help[] = "Estimates the 2-norm condition number of a matrix A, that is, the ratio of the largest to the smallest singular values of A. "
23: "The matrix is a Grcar matrix.\n\n"
24: "The command line options are:\n"
25: " -n <n>, where <n> = matrix dimension.\n\n";
27: #include slepcsvd.h
29: /*
30: This example computes the singular values of an nxn Grcar matrix,
31: which is a nonsymmetric Toeplitz matrix:
33: | 1 1 1 1 |
34: | -1 1 1 1 1 |
35: | -1 1 1 1 1 |
36: | . . . . . |
37: A = | . . . . . |
38: | -1 1 1 1 1 |
39: | -1 1 1 1 |
40: | -1 1 1 |
41: | -1 1 |
43: */
47: int main( int argc, char **argv )
48: {
50: Mat A; /* Grcar matrix */
51: SVD svd; /* singular value solver context */
52: PetscInt N=30, Istart, Iend, i, col[5], nconv1, nconv2;
53: PetscScalar value[] = { -1, 1, 1, 1, 1 };
54: PetscReal sigma_1, sigma_n;
56: SlepcInitialize(&argc,&argv,(char*)0,help);
58: PetscOptionsGetInt(PETSC_NULL,"-n",&N,PETSC_NULL);
59: PetscPrintf(PETSC_COMM_WORLD,"\nEstimate the condition number of a Grcar matrix, n=%d\n\n",N);
61: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
62: Generate the matrix
63: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
65: MatCreate(PETSC_COMM_WORLD,&A);
66: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N);
67: MatSetFromOptions(A);
69: MatGetOwnershipRange(A,&Istart,&Iend);
70: for( i=Istart; i<Iend; i++ ) {
71: col[0]=i-1; col[1]=i; col[2]=i+1; col[3]=i+2; col[4]=i+3;
72: if (i==0) {
73: MatSetValues(A,1,&i,4,col+1,value+1,INSERT_VALUES);
74: }
75: else {
76: MatSetValues(A,1,&i,PetscMin(5,N-i+1),col,value,INSERT_VALUES);
77: }
78: }
80: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
81: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
83: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
84: Create the singular value solver and set the solution method
85: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
87: /*
88: Create singular value context
89: */
90: SVDCreate(PETSC_COMM_WORLD,&svd);
92: /*
93: Set operator
94: */
95: SVDSetOperator(svd,A);
97: /*
98: Set solver parameters at runtime
99: */
100: SVDSetFromOptions(svd);
101: SVDSetDimensions(svd,1,PETSC_IGNORE,PETSC_IGNORE);
103: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
104: Solve the eigensystem
105: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
107: /*
108: First request an eigenvalue from one end of the spectrum
109: */
110: SVDSetWhichSingularTriplets(svd,SVD_LARGEST);
111: SVDSolve(svd);
112: /*
113: Get number of converged singular values
114: */
115: SVDGetConverged(svd,&nconv1);
116: /*
117: Get converged singular values: largest singular value is stored in sigma_1.
118: In this example, we are not interested in the singular vectors
119: */
120: if (nconv1 > 0) {
121: SVDGetSingularTriplet(svd,0,&sigma_1,PETSC_NULL,PETSC_NULL);
122: } else {
123: PetscPrintf(PETSC_COMM_WORLD," Unable to compute large singular value!\n\n");
124: }
126: /*
127: Request an eigenvalue from the other end of the spectrum
128: */
129: SVDSetWhichSingularTriplets(svd,SVD_SMALLEST);
130: SVDSolve(svd);
131: /*
132: Get number of converged eigenpairs
133: */
134: SVDGetConverged(svd,&nconv2);
135: /*
136: Get converged singular values: smallest singular value is stored in sigma_n.
137: As before, we are not interested in the singular vectors
138: */
139: if (nconv2 > 0) {
140: SVDGetSingularTriplet(svd,0,&sigma_n,PETSC_NULL,PETSC_NULL);
141: } else {
142: PetscPrintf(PETSC_COMM_WORLD," Unable to compute small singular value!\n\n");
143: }
145: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146: Display solution and clean up
147: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
148: if (nconv1 > 0 && nconv2 > 0) {
149: PetscPrintf(PETSC_COMM_WORLD," Computed singular values: sigma_1=%6f, sigma_n=%6f\n",sigma_1,sigma_n);
150: PetscPrintf(PETSC_COMM_WORLD," Estimated condition number: sigma_1/sigma_n=%6f\n\n",sigma_1/sigma_n);
151: }
152:
153: /*
154: Free work space
155: */
156: SVDDestroy(svd);
157: MatDestroy(A);
158: SlepcFinalize();
159: return 0;
160: }