Actual source code: ex8.c
2: 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. "
3: "The matrix is a Grcar matrix.\n\n"
4: "The command line options are:\n"
5: " -n <n>, where <n> = matrix dimension.\n\n";
7: #include slepceps.h
9: /*
10: This example computes the singular values of A by computing the eigenvalues
11: of A^T*A, where A^T denotes the transpose of A.
13: An nxn Grcar matrix is a nonsymmetric Toeplitz matrix:
15: | 1 1 1 1 |
16: | -1 1 1 1 1 |
17: | -1 1 1 1 1 |
18: | . . . . . |
19: A = | . . . . . |
20: | -1 1 1 1 1 |
21: | -1 1 1 1 |
22: | -1 1 1 |
23: | -1 1 |
25: Note that working with A^T*A can lead to poor accuracy of the computed
26: singular values when A is ill-conditioned (which is not the case here).
27: Another alternative would be to compute the eigenvalues of
29: | 0 A |
30: | A^T 0 |
32: but this significantly increases the cost of the solution process.
34: */
37: /*
38: Matrix multiply routine
39: */
42: PetscErrorCode MatSVD_Mult(Mat H,Vec x,Vec y)
43: {
44: Mat A;
45: Vec w;
49: MatShellGetContext(H,(void**)&A);
50: MatGetVecs(A,PETSC_NULL,&w);
51: MatMult(A,x,w);
52: MatMultTranspose(A,w,y);
53: VecDestroy(w);
54: return(0);
55: }
59: int main( int argc, char **argv )
60: {
62: Mat A; /* Grcar matrix */
63: Mat H; /* eigenvalue problem matrix, H=A^T*A */
64: EPS eps; /* eigenproblem solver context */
65: PetscInt N=30, n, Istart, Iend, i, col[5];
66: int nconv1, nconv2;
67: PetscScalar kl, ks, value[] = { -1, 1, 1, 1, 1 };
68: PetscReal sigma_1, sigma_n;
70: SlepcInitialize(&argc,&argv,(char*)0,help);
72: PetscOptionsGetInt(PETSC_NULL,"-n",&N,PETSC_NULL);
73: PetscPrintf(PETSC_COMM_WORLD,"\nEstimate the condition number of a Grcar matrix, n=%d\n\n",N);
75: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
76: Generate the matrix
77: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
79: MatCreate(PETSC_COMM_WORLD,&A);
80: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N);
81: MatSetFromOptions(A);
83: MatGetOwnershipRange(A,&Istart,&Iend);
84: for( i=Istart; i<Iend; i++ ) {
85: col[0]=i-1; col[1]=i; col[2]=i+1; col[3]=i+2; col[4]=i+3;
86: if (i==0) {
87: MatSetValues(A,1,&i,4,col+1,value+1,INSERT_VALUES);
88: }
89: else {
90: MatSetValues(A,1,&i,PetscMin(5,N-i+1),col,value,INSERT_VALUES);
91: }
92: }
94: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
95: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
97: /*
98: Now create a symmetric shell matrix H=A^T*A
99: */
100: MatGetLocalSize(A,PETSC_NULL,&n);
101: MatCreateShell(PETSC_COMM_WORLD,n,n,PETSC_DETERMINE,PETSC_DETERMINE,(void*)A,&H);
102: MatShellSetOperation(H,MATOP_MULT,(void(*)())MatSVD_Mult);
103: MatShellSetOperation(H,MATOP_MULT_TRANSPOSE,(void(*)())MatSVD_Mult);
105: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106: Create the eigensolver and set the solution method
107: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
109: /*
110: Create eigensolver context
111: */
112: EPSCreate(PETSC_COMM_WORLD,&eps);
114: /*
115: Set operators. In this case, it is a standard symmetric eigenvalue problem
116: */
117: EPSSetOperators(eps,H,PETSC_NULL);
118: EPSSetProblemType(eps,EPS_HEP);
120: /*
121: Set solver parameters at runtime
122: */
123: EPSSetFromOptions(eps);
124: EPSSetDimensions(eps,1,PETSC_DEFAULT);
125: EPSSetTolerances(eps,PETSC_DEFAULT,1000);
127: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128: Solve the eigensystem
129: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
131: /*
132: First request an eigenvalue from one end of the spectrum
133: */
134: EPSSetWhichEigenpairs(eps,EPS_LARGEST_REAL);
135: EPSSolve(eps);
136: /*
137: Get number of converged eigenpairs
138: */
139: EPSGetConverged(eps,&nconv1);
140: /*
141: Get converged eigenpairs: largest eigenvalue is stored in kl. In this
142: example, we are not interested in the eigenvectors
143: */
144: if (nconv1 > 0) {
145: EPSGetEigenpair(eps,0,&kl,PETSC_NULL,PETSC_NULL,PETSC_NULL);
146: }
148: /*
149: Request an eigenvalue from the other end of the spectrum
150: */
151: EPSSetWhichEigenpairs(eps,EPS_SMALLEST_REAL);
152: EPSSolve(eps);
153: /*
154: Get number of converged eigenpairs
155: */
156: EPSGetConverged(eps,&nconv2);
157: /*
158: Get converged eigenpairs: smallest eigenvalue is stored in ks. In this
159: example, we are not interested in the eigenvectors
160: */
161: if (nconv2 > 0) {
162: EPSGetEigenpair(eps,0,&ks,PETSC_NULL,PETSC_NULL,PETSC_NULL);
163: }
165: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
166: Display solution and clean up
167: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
168: if (nconv1 > 0 && nconv2 > 0) {
169: sigma_1 = sqrt(PetscRealPart(kl));
170: sigma_n = sqrt(PetscRealPart(ks));
172: PetscPrintf(PETSC_COMM_WORLD," Computed singular values: sigma_1=%6f, sigma_n=%6f\n",sigma_1,sigma_n);
173: PetscPrintf(PETSC_COMM_WORLD," Estimated condition number: sigma_1/sigma_n=%6f\n\n",sigma_1/sigma_n);
174: } else {
175: PetscPrintf(PETSC_COMM_WORLD," Process did not converge! Try running with a larger value for -eps_ncv\n\n");
176: }
177:
178: /*
179: Free work space
180: */
181: EPSDestroy(eps);
182: MatDestroy(A);
183: MatDestroy(H);
184: SlepcFinalize();
185: return 0;
186: }