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