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 DSNHEP.\n\n";
12 :
13 : #include <slepcds.h>
14 :
15 2 : int main(int argc,char **argv)
16 : {
17 2 : DS ds;
18 2 : SlepcSC sc;
19 2 : DSType type;
20 2 : DSStateType state;
21 2 : PetscScalar *A,*X,*Q,*wr,*wi,d;
22 2 : PetscReal re,im,rnorm,aux;
23 2 : PetscInt i,j,n=10,ld,method;
24 2 : PetscViewer viewer;
25 2 : PetscBool verbose,extrarow;
26 :
27 2 : PetscFunctionBeginUser;
28 2 : PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
29 2 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
30 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Solve a Dense System of type NHEP - dimension %" PetscInt_FMT ".\n",n));
31 2 : PetscCall(PetscOptionsHasName(NULL,NULL,"-verbose",&verbose));
32 2 : PetscCall(PetscOptionsHasName(NULL,NULL,"-extrarow",&extrarow));
33 :
34 : /* Create DS object */
35 2 : PetscCall(DSCreate(PETSC_COMM_WORLD,&ds));
36 2 : PetscCall(DSSetType(ds,DSNHEP));
37 2 : PetscCall(DSSetFromOptions(ds));
38 2 : ld = n+2; /* test leading dimension larger than n */
39 2 : PetscCall(DSAllocate(ds,ld));
40 2 : PetscCall(DSSetDimensions(ds,n,0,0));
41 2 : PetscCall(DSSetExtraRow(ds,extrarow));
42 :
43 : /* Set up viewer */
44 2 : PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&viewer));
45 2 : PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO_DETAIL));
46 2 : PetscCall(DSView(ds,viewer));
47 2 : PetscCall(PetscViewerPopFormat(viewer));
48 2 : if (verbose) PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB));
49 :
50 : /* Fill with Grcar matrix */
51 2 : PetscCall(DSGetArray(ds,DS_MAT_A,&A));
52 20 : for (i=1;i<n;i++) A[i+(i-1)*ld]=-1.0;
53 10 : for (j=0;j<4;j++) {
54 76 : for (i=0;i<n-j;i++) A[i+(i+j)*ld]=1.0;
55 : }
56 2 : if (extrarow) A[n+(n-1)*ld]=-1.0;
57 2 : PetscCall(DSRestoreArray(ds,DS_MAT_A,&A));
58 2 : PetscCall(DSSetState(ds,DS_STATE_INTERMEDIATE));
59 2 : if (verbose) {
60 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Initial - - - - - - - - -\n"));
61 0 : PetscCall(DSView(ds,viewer));
62 : }
63 :
64 : /* Solve */
65 2 : PetscCall(PetscMalloc2(n,&wr,n,&wi));
66 2 : PetscCall(DSGetSlepcSC(ds,&sc));
67 2 : sc->comparison = SlepcCompareLargestMagnitude;
68 2 : sc->comparisonctx = NULL;
69 2 : sc->map = NULL;
70 2 : sc->mapobj = NULL;
71 2 : PetscCall(DSSolve(ds,wr,wi));
72 2 : PetscCall(DSSort(ds,wr,wi,NULL,NULL,NULL));
73 2 : if (extrarow) PetscCall(DSUpdateExtraRow(ds));
74 :
75 2 : PetscCall(DSGetType(ds,&type));
76 2 : PetscCall(DSGetMethod(ds,&method));
77 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"DS of type %s, method used=%" PetscInt_FMT "\n",type,method));
78 2 : PetscCall(DSGetState(ds,&state));
79 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"State after solve: %s\n",DSStateTypes[state]));
80 :
81 2 : if (verbose) {
82 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"After solve - - - - - - - - -\n"));
83 0 : PetscCall(DSView(ds,viewer));
84 : }
85 :
86 : /* Print eigenvalues */
87 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Computed eigenvalues =\n"));
88 22 : for (i=0;i<n;i++) {
89 : #if defined(PETSC_USE_COMPLEX)
90 : re = PetscRealPart(wr[i]);
91 : im = PetscImaginaryPart(wr[i]);
92 : #else
93 20 : re = wr[i];
94 20 : im = wi[i];
95 : #endif
96 20 : if (PetscAbs(im)<1e-10) PetscCall(PetscViewerASCIIPrintf(viewer," %.5f\n",(double)re));
97 20 : else PetscCall(PetscViewerASCIIPrintf(viewer," %.5f%+.5fi\n",(double)re,(double)im));
98 : }
99 :
100 2 : if (extrarow) {
101 : /* Check that extra row is correct */
102 1 : PetscCall(DSGetArray(ds,DS_MAT_A,&A));
103 1 : PetscCall(DSGetArray(ds,DS_MAT_Q,&Q));
104 : d = 0.0;
105 11 : for (i=0;i<n;i++) d += A[n+i*ld]+Q[n-1+i*ld];
106 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)));
107 1 : PetscCall(DSRestoreArray(ds,DS_MAT_A,&A));
108 1 : PetscCall(DSRestoreArray(ds,DS_MAT_Q,&Q));
109 : }
110 :
111 : /* Eigenvectors */
112 2 : j = 2;
113 2 : PetscCall(DSVectors(ds,DS_MAT_X,&j,&rnorm)); /* third eigenvector */
114 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Value of rnorm for 3rd vector = %.3f\n",(double)rnorm));
115 2 : PetscCall(DSVectors(ds,DS_MAT_X,NULL,NULL)); /* all eigenvectors */
116 2 : j = 0;
117 2 : rnorm = 0.0;
118 2 : PetscCall(DSGetArray(ds,DS_MAT_X,&X));
119 22 : for (i=0;i<n;i++) {
120 : #if defined(PETSC_USE_COMPLEX)
121 : aux = PetscAbsScalar(X[i+j*ld]);
122 : #else
123 20 : if (PetscAbs(wi[j])==0.0) aux = PetscAbsScalar(X[i+j*ld]);
124 20 : else aux = SlepcAbsEigenvalue(X[i+j*ld],X[i+(j+1)*ld]);
125 : #endif
126 20 : rnorm += aux*aux;
127 : }
128 2 : PetscCall(DSRestoreArray(ds,DS_MAT_X,&X));
129 2 : rnorm = PetscSqrtReal(rnorm);
130 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Norm of 1st vector = %.3f\n",(double)rnorm));
131 2 : if (verbose) {
132 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"After vectors - - - - - - - - -\n"));
133 0 : PetscCall(DSView(ds,viewer));
134 : }
135 :
136 2 : PetscCall(PetscFree2(wr,wi));
137 2 : PetscCall(DSDestroy(&ds));
138 2 : PetscCall(SlepcFinalize());
139 : return 0;
140 : }
141 :
142 : /*TEST
143 :
144 : testset:
145 : filter: sed -e "s/[+-]\([0-9]\.[0-9]*i\)/+-\\1/" | sed -e "s/extrarow//"
146 : output_file: output/test1_1.out
147 : requires: !single
148 : test:
149 : suffix: 1
150 : test:
151 : suffix: 2
152 : args: -extrarow
153 :
154 : TEST*/
|
