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: }