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 DSHSVD with compact storage.\n\n";
12 :
13 : #include <slepcds.h>
14 :
15 4 : int main(int argc,char **argv)
16 : {
17 4 : DS ds;
18 4 : SlepcSC sc;
19 4 : PetscReal *T,*D,sigma;
20 4 : PetscScalar *U,*w,d;
21 4 : PetscInt i,n=10,m,l=2,k=5,p=1,ld;
22 4 : PetscViewer viewer;
23 4 : PetscBool verbose,extrarow,reorthog,flg;
24 :
25 4 : PetscFunctionBeginUser;
26 4 : PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
27 4 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
28 4 : m = n;
29 4 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Solve a Dense System of type HSVD with compact storage - dimension %" PetscInt_FMT "x%" PetscInt_FMT ".\n",n,m));
30 4 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-l",&l,NULL));
31 4 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-k",&k,NULL));
32 4 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-p",&p,NULL));
33 4 : PetscCheck(l<=n && k<=n && l<=k && p>=0 && p<=n-l,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Wrong value of dimensions");
34 4 : PetscCall(PetscOptionsHasName(NULL,NULL,"-verbose",&verbose));
35 4 : PetscCall(PetscOptionsHasName(NULL,NULL,"-extrarow",&extrarow));
36 4 : PetscCall(PetscOptionsHasName(NULL,NULL,"-reorthog",&reorthog));
37 :
38 : /* Create DS object */
39 4 : PetscCall(DSCreate(PETSC_COMM_WORLD,&ds));
40 4 : PetscCall(DSSetType(ds,DSHSVD));
41 4 : PetscCall(DSSetFromOptions(ds));
42 4 : ld = n+2; /* test leading dimension larger than n */
43 4 : PetscCall(DSAllocate(ds,ld));
44 4 : PetscCall(DSSetDimensions(ds,n,l,k));
45 4 : PetscCall(DSHSVDSetDimensions(ds,m));
46 4 : PetscCall(DSSetCompact(ds,PETSC_TRUE));
47 4 : PetscCall(DSSetExtraRow(ds,extrarow));
48 4 : PetscCall(DSHSVDSetReorthogonalize(ds,reorthog));
49 4 : PetscCall(DSHSVDGetReorthogonalize(ds,&flg));
50 4 : if (flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"reorthogonalizing\n"));
51 :
52 : /* Set up viewer */
53 4 : PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&viewer));
54 4 : PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO_DETAIL));
55 4 : PetscCall(DSView(ds,viewer));
56 4 : PetscCall(PetscViewerPopFormat(viewer));
57 4 : PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB));
58 :
59 : /* Fill upper arrow-bidiagonal matrix and signature matrix */
60 4 : PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
61 44 : for (i=0;i<n;i++) T[i] = (PetscReal)(i+1);
62 34 : for (i=l;i<n-1;i++) T[i+ld] = 1.0;
63 4 : if (extrarow) T[n-1+ld] = 1.0;
64 4 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
65 4 : PetscCall(DSGetArrayReal(ds,DS_MAT_D,&D));
66 44 : for (i=0;i<n;i++) D[i] = (i>=l && i<l+p)? -1.0: 1.0;
67 4 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&D));
68 4 : if (l==0 && k==0) PetscCall(DSSetState(ds,DS_STATE_INTERMEDIATE));
69 3 : else PetscCall(DSSetState(ds,DS_STATE_RAW));
70 4 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Initial - - - - - - - - -\n"));
71 4 : PetscCall(DSView(ds,viewer));
72 :
73 : /* Solve */
74 4 : PetscCall(PetscMalloc1(n,&w));
75 4 : PetscCall(DSGetSlepcSC(ds,&sc));
76 4 : sc->comparison = SlepcCompareLargestReal;
77 4 : sc->comparisonctx = NULL;
78 4 : sc->map = NULL;
79 4 : sc->mapobj = NULL;
80 4 : PetscCall(DSSolve(ds,w,NULL));
81 4 : PetscCall(DSSort(ds,w,NULL,NULL,NULL,NULL));
82 4 : if (extrarow) PetscCall(DSUpdateExtraRow(ds));
83 4 : if (verbose) {
84 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"After solve - - - - - - - - -\n"));
85 0 : PetscCall(DSView(ds,viewer));
86 : }
87 :
88 : /* Print singular values */
89 4 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Computed singular values =\n"));
90 44 : for (i=0;i<n;i++) {
91 40 : sigma = PetscRealPart(w[i]);
92 40 : PetscCall(PetscViewerASCIIPrintf(viewer," %.5f\n",(double)sigma));
93 : }
94 :
95 4 : if (extrarow) {
96 : /* Check that extra row is correct */
97 2 : PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
98 2 : PetscCall(DSGetArray(ds,DS_MAT_U,&U));
99 2 : PetscCall(DSGetArrayReal(ds,DS_MAT_D,&D));
100 2 : d = 0.0;
101 18 : for (i=l;i<n;i++) d += T[i+ld]-D[i]*U[n-1+i*ld];
102 2 : 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)));
103 2 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
104 2 : PetscCall(DSRestoreArray(ds,DS_MAT_U,&U));
105 2 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&D));
106 : }
107 4 : PetscCall(PetscFree(w));
108 4 : PetscCall(DSDestroy(&ds));
109 4 : PetscCall(SlepcFinalize());
110 : return 0;
111 : }
112 :
113 : /*TEST
114 :
115 : test:
116 : suffix: 1
117 : requires: !single
118 :
119 : test:
120 : args: -l 0 -k 0
121 : suffix: 2
122 : requires: !single
123 :
124 : test:
125 : args: -extrarow -reorthog {{0 1}}
126 : suffix: 3
127 : requires: !single
128 : filter: grep -v reorthogonalizing
129 :
130 : TEST*/
|
