1 /*
2 ARPACK++ v1.2 2/18/2000
3 c++ interface to ARPACK code.
4
5 MODULE RNSymGSC.cc.
6 Example program that illustrates how to solve a real
7 nonsymmetric generalized eigenvalue problem in complex shift
8 and invert mode using the ARrcNonSymGenEig class.
9
10 1) Problem description:
11
12 In this example we try to solve A*x = B*x*lambda in shift and
13 invert mode, where A is the tridiagonal matrix with 2 on the
14 diagonal, -2 on the subdiagonal and 3 on the superdiagonal, and
15 B is the tridiagonal matrix with 4 on the diagonal and 1 on the
16 off-diagonals.
17 The shift sigma is a complex number.
18
19 2) Data structure used to represent matrix A:
20
21 To obtain the eigenvalues of the above problem, the user is
22 required to provide a way to perform the matrix-vector products
23 w = OP*Bv = real{inv(A-sigma*B)}*B*v, w = A*v and w = B*v. In
24 this example, a class called NonSymGenProblemC was created with
25 this purpose. NonSymGenProblemC contains a member function,
26 MultOPv(v,w), that takes a vector v and returns the product OPv
27 in w. It also contains two objects, A and B, that
28 store matrices A and B, respectively. The product Bv is
29 performed by MultMv, a member function of B, and Av is obtained
30 by calling A.MultMv.
31
32 3) The reverse communication interface:
33
34 This example uses the reverse communication interface, which
35 means that the desired eigenvalues cannot be obtained directly
36 from an ARPACK++ class.
37 Here, the overall process of finding eigenvalues by using the
38 Arnoldi method is splitted into two parts. In the first, a
39 sequence of calls to a function called TakeStep is combined
40 with matrix-vector products in order to find an Arnoldi basis.
41 In the second part, an ARPACK++ function like FindEigenvectors
42 (or EigenValVectors) is used to extract eigenvalues and
43 eigenvectors.
44
45 4) Included header files:
46
47 File Contents
48 ----------- -------------------------------------------
49 ngenprbc.h The NonSymGenProblemC class definition.
50 arrgnsym.h The ARrcNonSymGenEig class definition.
51 rnsymgsl.h The Solution function.
52
53 5) ARPACK Authors:
54
55 Richard Lehoucq
56 Kristyn Maschhoff
57 Danny Sorensen
58 Chao Yang
59 Dept. of Computational & Applied Mathematics
60 Rice University
61 Houston, Texas
62 */
63
64 #include "ngenprbc.h"
65 #include "rnsymgsl.h"
66 #include "arrgnsym.h"
67
68
69 template<class ARFLOAT>
RecoverEigenvalues(long nconv,NonSymGenProblemC<ARFLOAT> & P,ARFLOAT EigVec[],ARFLOAT EigValR[],ARFLOAT EigValI[])70 void RecoverEigenvalues(long nconv, NonSymGenProblemC<ARFLOAT>& P,
71 ARFLOAT EigVec[], ARFLOAT EigValR[], ARFLOAT EigValI[])
72 /*
73 This function is used to recover the eigenvalues of the original
74 problem A.x = B.x.lambda, after calling ARPACK++.
75 When a complex shift is used to solve a nonsymmetric problem
76 defined by the ARRCNonSymGenEig class, the Rayleigh quotient
77 lambda = x'Ax/x'Bx must be formed by the user to obtain the
78 desired eigenvalues.
79 The Rayleigh quotient cannot be calculated automatically by
80 ARPACK++ because ARRCNonSymGenEig do not handle matrix information.
81 Other classes such as ARNonSymGenEig and ARLUNonSymGenEig do
82 not require the user to define this eigenvalue transformation.
83 */
84
85 {
86
87 int j, n, ColJ, ColJp1;
88 ARFLOAT numr, numi, denr, deni;
89 ARFLOAT* Ax;
90
91 n = P.A.ncols();
92 Ax = new ARFLOAT[n];
93
94 for (j=0; j<nconv; j++) {
95
96 ColJ = j*n;
97 ColJp1 = ColJ+n;
98
99 if (EigValI[j] == (ARFLOAT)0.0) {
100
101 // Eigenvalue is real. Computing EigVal = x'(Ax)/x'(Mx).
102
103 P.A.MultMv(&EigVec[ColJ], Ax);
104 numr = dot(n, &EigVec[ColJ], 1, Ax, 1);
105 P.B.MultMv(&EigVec[ColJ], Ax);
106 denr = dot(n, &EigVec[ColJ], 1, Ax, 1);
107 EigValR[j] = numr / denr;
108
109 }
110 else {
111
112 // Eigenvalue is complex.
113
114 // Computing x'(Ax).
115
116 P.A.MultMv(&EigVec[ColJ], Ax);
117 numr = dot(n, &EigVec[ColJ], 1, Ax, 1);
118 numi = dot(n, &EigVec[ColJp1], 1, Ax, 1);
119 P.A.MultMv(&EigVec[ColJp1], Ax);
120 numr = numr + dot(n, &EigVec[ColJp1], 1, Ax, 1);
121 numi = -numi + dot(n, &EigVec[ColJ], 1, Ax, 1);
122
123 // Computing x'(Mx).
124
125 P.B.MultMv(&EigVec[ColJ], Ax);
126 denr = dot(n, &EigVec[ColJ], 1, Ax, 1);
127 deni = dot(n, &EigVec[ColJp1], 1, Ax, 1);
128 P.B.MultMv(&EigVec[ColJp1], Ax);
129 denr = denr + dot(n, &EigVec[ColJp1], 1, Ax, 1);
130 deni = -deni + dot(n, &EigVec[ColJ], 1, Ax, 1);
131
132 // Computing the first eigenvalue of the conjugate pair.
133
134 EigValR[j] = (numr*denr+numi*deni) / lapy2(denr, deni);
135 EigValI[j] = (numi*denr-numr*deni) / lapy2(denr, deni);
136
137 // Getting the second eigenvalue of the conjugate pair by taking
138 // the conjugate of the first.
139
140 EigValR[j+1] = EigValR[j];
141 EigValI[j+1] = -EigValI[j];
142 j++;
143
144 }
145
146 }
147
148 delete[] Ax;
149
150 } // RecoverEigenvalues.
151
152
153 template<class T>
Test(T type)154 void Test(T type)
155 {
156
157 long nconv;
158
159 // Defining a temporary vector.
160
161 T temp[101];
162
163 // Creating a pencil.
164
165 NonSymGenProblemC<T> P(100, 0.4, 0.6); // n = 100, sigma = (0.4, 0.6).
166
167 // Creating a nonsymmetric eigenvalue problem. 'R' indicates that
168 // we will use the real part of OPv. P.A.ncols() furnishes
169 // the dimension of the problem. 4 is the number of eigenvalues
170 // sought and 0.4 + 0.6I is the shift.
171
172 ARrcNonSymGenEig<T> prob(P.A.ncols(), 4L, 'R', 0.4, 0.6);
173
174 // Finding an Arnoldi basis.
175
176 while (!prob.ArnoldiBasisFound()) {
177
178 // Calling ARPACK FORTRAN code. Almost all work needed to
179 // find an Arnoldi basis is performed by TakeStep.
180
181 prob.TakeStep();
182
183 switch (prob.GetIdo()) {
184 case -1:
185
186 // Performing w <- Re{OP*B*v} for the first time.
187 // This product must be performed only if GetIdo is equal to
188 // -1. GetVector supplies a pointer to the input vector, v,
189 // and PutVector a pointer to the output vector, w.
190
191 P.B.MultMv(prob.GetVector(), temp);
192 P.MultOPvRe(temp, prob.PutVector());
193 break;
194
195 case 1:
196
197 // Performing w <- Real{OP*B*v} when Bv is available.
198 // This product must be performed whenever GetIdo is equal to
199 // 1. GetProd supplies a pointer to the previously calculated
200 // product Bv and PutVector a pointer to the output vector w.
201
202 P.MultOPvRe(prob.GetProd(), prob.PutVector());
203 break;
204
205 case 2:
206
207 // Performing w <- B*v.
208
209 P.B.MultMv(prob.GetVector(), prob.PutVector());
210
211 }
212 }
213
214 // Finding eigenvalues and eigenvectors.
215
216 nconv = prob.FindEigenvectors();
217
218 // Recovering eigenvalues of the original problem
219 // using the Rayleigh quotient.
220
221 RecoverEigenvalues(nconv, P, prob.RawEigenvectors(),
222 prob.RawEigenvalues(), prob.RawEigenvaluesImag());
223
224 // Printing solution.
225
226 Solution(prob);
227
228 } // Test.
229
230
main()231 int main()
232 {
233
234 // Solving a single precision problem with n = 100.
235
236 #ifndef __SUNPRO_CC
237
238 Test((float)0.0);
239
240 #endif
241
242 // Solving a double precision problem with n = 100.
243
244 Test((double)0.0);
245
246 } // main
247
248