Actual source code: test1.c
slepc-3.18.1 2022-11-02
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
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";
16: #include <slepcsvd.h>
18: /*
19: This example computes the singular values of an nxn Grcar matrix,
20: which is a nonsymmetric Toeplitz matrix:
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 |
32: */
34: int main(int argc,char **argv)
35: {
36: Mat A; /* Grcar matrix */
37: SVD svd; /* singular value solver context */
38: PetscInt N=30,Istart,Iend,i,col[5],nconv1,nconv2;
39: PetscScalar value[] = { -1, 1, 1, 1, 1 };
40: PetscReal sigma_1,sigma_n;
41: char svdtype[30] = "cross",epstype[30] = "";
42: PetscBool flg;
43: EPS eps;
46: SlepcInitialize(&argc,&argv,(char*)0,help);
48: PetscOptionsGetInt(NULL,NULL,"-n",&N,NULL);
49: PetscOptionsGetString(NULL,NULL,"-type",svdtype,sizeof(svdtype),NULL);
50: PetscOptionsGetString(NULL,NULL,"-epstype",epstype,sizeof(epstype),&flg);
51: PetscPrintf(PETSC_COMM_WORLD,"\nEstimate the condition number of a Grcar matrix, n=%" PetscInt_FMT,N);
52: PetscPrintf(PETSC_COMM_WORLD,"\n\n");
54: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
55: Generate the matrix
56: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
58: MatCreate(PETSC_COMM_WORLD,&A);
59: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N);
60: MatSetFromOptions(A);
61: MatSetUp(A);
63: MatGetOwnershipRange(A,&Istart,&Iend);
64: for (i=Istart;i<Iend;i++) {
65: col[0]=i-1; col[1]=i; col[2]=i+1; col[3]=i+2; col[4]=i+3;
66: if (i==0) MatSetValues(A,1,&i,4,col+1,value+1,INSERT_VALUES);
67: else MatSetValues(A,1,&i,PetscMin(5,N-i+1),col,value,INSERT_VALUES);
68: }
70: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
71: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
73: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
74: Create the singular value solver and set the solution method
75: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
77: /*
78: Create singular value context
79: */
80: SVDCreate(PETSC_COMM_WORLD,&svd);
82: /*
83: Set operator
84: */
85: SVDSetOperators(svd,A,NULL);
87: /*
88: Set solver parameters at runtime
89: */
90: SVDSetType(svd,svdtype);
91: if (flg) {
92: PetscObjectTypeCompare((PetscObject)svd,SVDCROSS,&flg);
93: if (flg) {
94: SVDCrossGetEPS(svd,&eps);
95: EPSSetType(eps,epstype);
96: }
97: PetscObjectTypeCompare((PetscObject)svd,SVDCYCLIC,&flg);
98: if (flg) {
99: SVDCyclicGetEPS(svd,&eps);
100: EPSSetType(eps,epstype);
101: }
102: }
103: SVDSetDimensions(svd,1,PETSC_DEFAULT,PETSC_DEFAULT);
104: SVDSetTolerances(svd,1e-6,1000);
106: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107: Compute the singular values
108: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
110: /*
111: First request the largest singular value
112: */
113: SVDSetWhichSingularTriplets(svd,SVD_LARGEST);
114: SVDSolve(svd);
115: /*
116: Get number of converged singular values
117: */
118: SVDGetConverged(svd,&nconv1);
119: /*
120: Get converged singular values: largest singular value is stored in sigma_1.
121: In this example, we are not interested in the singular vectors
122: */
123: if (nconv1 > 0) SVDGetSingularTriplet(svd,0,&sigma_1,NULL,NULL);
124: else PetscPrintf(PETSC_COMM_WORLD," Unable to compute large singular value!\n\n");
126: /*
127: Request the smallest singular value
128: */
129: SVDSetWhichSingularTriplets(svd,SVD_SMALLEST);
130: SVDSolve(svd);
131: /*
132: Get number of converged triplets
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) SVDGetSingularTriplet(svd,0,&sigma_n,NULL,NULL);
140: else PetscPrintf(PETSC_COMM_WORLD," Unable to compute small singular value!\n\n");
142: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
143: Display solution and clean up
144: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
145: if (nconv1 > 0 && nconv2 > 0) {
146: PetscPrintf(PETSC_COMM_WORLD," Computed singular values: sigma_1=%.4f, sigma_n=%.4f\n",(double)sigma_1,(double)sigma_n);
147: PetscPrintf(PETSC_COMM_WORLD," Estimated condition number: sigma_1/sigma_n=%.4f\n\n",(double)(sigma_1/sigma_n));
148: }
150: /*
151: Free work space
152: */
153: SVDDestroy(&svd);
154: MatDestroy(&A);
155: SlepcFinalize();
156: return 0;
157: }
159: /*TEST
161: test:
162: suffix: 1
163: args: -type {{lanczos trlanczos cross cyclic lapack}}
165: test:
166: suffix: 1_cross_gd
167: args: -type cross -epstype gd
168: output_file: output/test1_1.out
170: test:
171: suffix: 1_cyclic_gd
172: args: -type cyclic -epstype gd
173: output_file: output/test1_1.out
174: requires: !single
176: test:
177: suffix: 1_primme
178: args: -type primme
179: requires: primme
180: output_file: output/test1_1.out
182: TEST*/