Actual source code: test21.c
slepc-3.18.0 2022-10-01
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 DSGSVD.\n\n";
13: #include <slepcds.h>
15: int main(int argc,char **argv)
16: {
17: DS ds;
18: SlepcSC sc;
19: Mat X;
20: Vec x0;
21: PetscReal sigma,rnorm,cond;
22: PetscScalar *A,*B,*w;
23: PetscInt i,j,k,n=15,m=10,p=10,m1,p1,ld;
24: PetscViewer viewer;
25: PetscBool verbose;
28: SlepcInitialize(&argc,&argv,(char*)0,help);
29: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
30: PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
31: PetscOptionsGetInt(NULL,NULL,"-p",&p,NULL);
32: PetscPrintf(PETSC_COMM_WORLD,"Solve a Dense System of type GSVD - dimension (%" PetscInt_FMT "+%" PetscInt_FMT ")x%" PetscInt_FMT ".\n",n,p,m);
33: PetscOptionsHasName(NULL,NULL,"-verbose",&verbose);
35: /* Create DS object */
36: DSCreate(PETSC_COMM_WORLD,&ds);
37: DSSetType(ds,DSGSVD);
38: DSSetFromOptions(ds);
39: ld = PetscMax(PetscMax(p,m),n)+2; /* test leading dimension larger than n */
40: DSAllocate(ds,ld);
41: DSSetDimensions(ds,n,0,0);
42: DSGSVDSetDimensions(ds,m,p);
43: DSGSVDGetDimensions(ds,&m1,&p1);
46: /* Set up viewer */
47: PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&viewer);
48: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO_DETAIL);
49: DSView(ds,viewer);
50: PetscViewerPopFormat(viewer);
52: k = PetscMin(n,m);
53: /* Fill A with a rectangular Toeplitz matrix */
54: DSGetArray(ds,DS_MAT_A,&A);
55: for (i=0;i<k;i++) A[i+i*ld]=1.0;
56: for (j=1;j<3;j++) {
57: for (i=0;i<n-j;i++) { if ((i+j)<m) A[i+(i+j)*ld]=(PetscScalar)(j+1); }
58: }
59: for (j=1;j<n/2;j++) {
60: for (i=0;i<n-j;i++) { if ((i+j)<n && i<m) A[(i+j)+i*ld]=-1.0; }
61: }
62: DSRestoreArray(ds,DS_MAT_A,&A);
64: k = PetscMin(p,m);
65: /* Fill B with a shifted bidiagonal matrix */
66: DSGetArray(ds,DS_MAT_B,&B);
67: for (i=m-k;i<m;i++) {
68: B[i-m+k+i*ld]=2.0-1.0/(PetscScalar)(i+1);
69: if (i) B[i-1-m+k+i*ld]=1.0;
70: }
71: DSRestoreArray(ds,DS_MAT_B,&B);
73: DSSetState(ds,DS_STATE_RAW);
74: if (verbose) {
75: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
76: PetscPrintf(PETSC_COMM_WORLD,"Initial - - - - - - - - -\n");
77: }
78: DSView(ds,viewer);
80: /* Condition number */
81: DSCond(ds,&cond);
82: PetscPrintf(PETSC_COMM_WORLD,"Condition number = %.3f\n",(double)cond);
84: /* Solve */
85: PetscMalloc1(m,&w);
86: DSGetSlepcSC(ds,&sc);
87: sc->comparison = SlepcCompareLargestReal;
88: sc->comparisonctx = NULL;
89: sc->map = NULL;
90: sc->mapobj = NULL;
91: DSSolve(ds,w,NULL);
92: DSSort(ds,w,NULL,NULL,NULL,NULL);
93: DSSynchronize(ds,w,NULL);
94: if (verbose) {
95: PetscPrintf(PETSC_COMM_WORLD,"After solve - - - - - - - - -\n");
96: DSView(ds,viewer);
97: }
98: /* Print singular values */
99: PetscPrintf(PETSC_COMM_WORLD,"Computed singular values =\n");
100: DSGetDimensions(ds,NULL,NULL,NULL,&k);
101: for (i=0;i<k;i++) {
102: sigma = PetscRealPart(w[i]);
103: PetscViewerASCIIPrintf(viewer," %g\n",(double)sigma);
104: }
106: /* Singular vectors */
107: DSVectors(ds,DS_MAT_X,NULL,NULL); /* all singular vectors */
108: DSGetMat(ds,DS_MAT_X,&X);
109: MatCreateVecs(X,NULL,&x0);
110: MatGetColumnVector(X,x0,0);
111: VecNorm(x0,NORM_2,&rnorm);
112: DSRestoreMat(ds,DS_MAT_X,&X);
113: VecDestroy(&x0);
114: PetscPrintf(PETSC_COMM_WORLD,"Norm of 1st X vector = %.3f\n",(double)rnorm);
116: DSGetMat(ds,DS_MAT_U,&X);
117: MatCreateVecs(X,NULL,&x0);
118: MatGetColumnVector(X,x0,0);
119: VecNorm(x0,NORM_2,&rnorm);
120: DSRestoreMat(ds,DS_MAT_U,&X);
121: VecDestroy(&x0);
122: if (PetscAbs(rnorm-1.0)>10*PETSC_MACHINE_EPSILON) PetscPrintf(PETSC_COMM_WORLD,"Warning: the 1st U vector has norm %g\n",(double)rnorm);
124: DSGetMat(ds,DS_MAT_V,&X);
125: MatCreateVecs(X,NULL,&x0);
126: MatGetColumnVector(X,x0,0);
127: VecNorm(x0,NORM_2,&rnorm);
128: DSRestoreMat(ds,DS_MAT_V,&X);
129: VecDestroy(&x0);
130: if (PetscAbs(rnorm-1.0)>10*PETSC_MACHINE_EPSILON) PetscPrintf(PETSC_COMM_WORLD,"Warning: the 1st V vector has norm %g\n",(double)rnorm);
132: PetscFree(w);
133: DSDestroy(&ds);
134: SlepcFinalize();
135: return 0;
136: }
138: /*TEST
140: testset:
141: output_file: output/test21_1.out
142: requires: !single
143: nsize: {{1 2 3}}
144: filter: grep -v "parallel operation mode" | grep -v " MPI process"
145: test:
146: suffix: 1
147: args: -ds_parallel redundant
148: test:
149: suffix: 2
150: args: -ds_parallel synchronized
152: TEST*/