Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
6 : This file is part of SLEPc.
7 : SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9 : */
10 :
11 : static char help[] = "Test DSSVD with compact storage.\n\n";
12 :
13 : #include <slepcds.h>
14 :
15 3 : int main(int argc,char **argv)
16 : {
17 3 : DS ds;
18 3 : SlepcSC sc;
19 3 : PetscReal *T,sigma;
20 3 : PetscScalar *U,*w,d;
21 3 : PetscInt i,n=10,m,l=2,k=5,ld;
22 3 : PetscViewer viewer;
23 3 : PetscBool verbose,extrarow;
24 :
25 3 : PetscFunctionBeginUser;
26 3 : PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
27 3 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
28 3 : m = n;
29 3 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Solve a Dense System of type SVD with compact storage - dimension %" PetscInt_FMT "x%" PetscInt_FMT ".\n",n,m));
30 3 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-l",&l,NULL));
31 3 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-k",&k,NULL));
32 3 : PetscCheck(l<=n && k<=n && l<=k,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Wrong value of dimensions");
33 3 : PetscCall(PetscOptionsHasName(NULL,NULL,"-verbose",&verbose));
34 3 : PetscCall(PetscOptionsHasName(NULL,NULL,"-extrarow",&extrarow));
35 :
36 : /* Create DS object */
37 3 : PetscCall(DSCreate(PETSC_COMM_WORLD,&ds));
38 3 : PetscCall(DSSetType(ds,DSSVD));
39 3 : PetscCall(DSSetFromOptions(ds));
40 3 : ld = n+2; /* test leading dimension larger than n */
41 3 : PetscCall(DSAllocate(ds,ld));
42 3 : PetscCall(DSSetDimensions(ds,n,l,k));
43 3 : PetscCall(DSSVDSetDimensions(ds,m));
44 3 : PetscCall(DSSetCompact(ds,PETSC_TRUE));
45 3 : PetscCall(DSSetExtraRow(ds,extrarow));
46 :
47 : /* Set up viewer */
48 3 : PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&viewer));
49 3 : PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO_DETAIL));
50 3 : PetscCall(DSView(ds,viewer));
51 3 : PetscCall(PetscViewerPopFormat(viewer));
52 3 : PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB));
53 :
54 : /* Fill upper arrow-bidiagonal matrix */
55 3 : PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
56 33 : for (i=0;i<n;i++) T[i] = (PetscReal)(i+1);
57 26 : for (i=l;i<n-1;i++) T[i+ld] = 1.0;
58 3 : if (extrarow) T[n-1+ld] = 1.0;
59 3 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
60 3 : if (l==0 && k==0) PetscCall(DSSetState(ds,DS_STATE_INTERMEDIATE));
61 2 : else PetscCall(DSSetState(ds,DS_STATE_RAW));
62 3 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Initial - - - - - - - - -\n"));
63 3 : PetscCall(DSView(ds,viewer));
64 :
65 : /* Solve */
66 3 : PetscCall(PetscMalloc1(n,&w));
67 3 : PetscCall(DSGetSlepcSC(ds,&sc));
68 3 : sc->comparison = SlepcCompareLargestReal;
69 3 : sc->comparisonctx = NULL;
70 3 : sc->map = NULL;
71 3 : sc->mapobj = NULL;
72 3 : PetscCall(DSSolve(ds,w,NULL));
73 3 : PetscCall(DSSort(ds,w,NULL,NULL,NULL,NULL));
74 3 : if (extrarow) PetscCall(DSUpdateExtraRow(ds));
75 3 : if (verbose) {
76 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"After solve - - - - - - - - -\n"));
77 0 : PetscCall(DSView(ds,viewer));
78 : }
79 :
80 : /* Print singular values */
81 3 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Computed singular values =\n"));
82 33 : for (i=0;i<n;i++) {
83 30 : sigma = PetscRealPart(w[i]);
84 30 : PetscCall(PetscViewerASCIIPrintf(viewer," %.5f\n",(double)sigma));
85 : }
86 :
87 3 : if (extrarow) {
88 : /* Check that extra row is correct */
89 1 : PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
90 1 : PetscCall(DSGetArray(ds,DS_MAT_U,&U));
91 1 : d = 0.0;
92 9 : for (i=l;i<n;i++) d += T[i+ld]-U[n-1+i*ld];
93 1 : if (PetscAbsScalar(d)>10*PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Warning: there is a mismatch in the extra row of %g\n",(double)PetscAbsScalar(d)));
94 1 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
95 1 : PetscCall(DSRestoreArray(ds,DS_MAT_U,&U));
96 : }
97 3 : PetscCall(PetscFree(w));
98 3 : PetscCall(DSDestroy(&ds));
99 3 : PetscCall(SlepcFinalize());
100 : return 0;
101 : }
102 :
103 : /*TEST
104 :
105 : test:
106 : suffix: 1
107 : requires: !single
108 :
109 : test:
110 : args: -l 0 -k 0
111 : suffix: 2
112 : requires: !single
113 :
114 : test:
115 : args: -extrarow
116 : suffix: 3
117 : requires: !single
118 :
119 : TEST*/
|