1 /*
2 ARPACK++ v1.2 2/20/2000
3 c++ interface to ARPACK code.
4
5 MODULE ARDNSPen.h.
6 Arpack++ class ARdsNonSymPencil definition.
7
8 ARPACK Authors
9 Richard Lehoucq
10 Danny Sorensen
11 Chao Yang
12 Dept. of Computational & Applied Mathematics
13 Rice University
14 Houston, Texas
15 */
16
17 #ifndef ARDNSPEN_H
18 #define ARDNSPEN_H
19
20 #include "arch.h"
21 #include "arerror.h"
22 #include "blas1c.h"
23 #include "lapackc.h"
24 #include "ardnsmat.h"
25
26
27 template<class ARTYPE, class ARFLOAT>
28 class ARdsNonSymPencil
29 {
30
31 protected:
32
33 char part;
34 ARdsNonSymMatrix<ARTYPE, ARFLOAT>* A;
35 ARdsNonSymMatrix<ARTYPE, ARFLOAT>* B;
36 ARdsNonSymMatrix<ARTYPE, ARFLOAT> AsB;
37 #ifdef ARCOMP_H
38 ARdsNonSymMatrix<arcomplex<ARFLOAT>, ARFLOAT> AsBc;
39 #endif
40
41 virtual void Copy(const ARdsNonSymPencil& other);
42
43 public:
44
45 #ifdef ARCOMP_H
IsFactored()46 bool IsFactored() { return (AsB.IsFactored()||AsBc.IsFactored()); }
47 #else
48 bool IsFactored() { return AsB.IsFactored(); }
49 #endif
50
51 void FactorAsB(ARTYPE sigma);
52
53 #ifdef ARCOMP_H
54 void FactorAsB(ARFLOAT sigmaR, ARFLOAT sigmaI, char partp = 'R');
55 #endif
56
MultAv(ARTYPE * v,ARTYPE * w)57 void MultAv(ARTYPE* v, ARTYPE* w) { A->MultMv(v,w); }
58
MultBv(ARTYPE * v,ARTYPE * w)59 void MultBv(ARTYPE* v, ARTYPE* w) { B->MultMv(v,w); }
60
61 void MultInvBAv(ARTYPE* v, ARTYPE* w);
62
63 #ifdef ARCOMP_H
64 void MultInvAsBv(arcomplex<ARFLOAT>* v, arcomplex<ARFLOAT>* w);
65 #endif
66
67 void MultInvAsBv(ARFLOAT* v, ARFLOAT* w);
68
69 void DefineMatrices(ARdsNonSymMatrix<ARTYPE, ARFLOAT>& Ap,
70 ARdsNonSymMatrix<ARTYPE, ARFLOAT>& Bp);
71
ARdsNonSymPencil()72 ARdsNonSymPencil() { part = 'N'; }
73 // Short constructor that does nothing.
74
75 ARdsNonSymPencil(ARdsNonSymMatrix<ARTYPE, ARFLOAT>& Ap,
76 ARdsNonSymMatrix<ARTYPE, ARFLOAT>& Bp);
77 // Long constructor.
78
ARdsNonSymPencil(const ARdsNonSymPencil & other)79 ARdsNonSymPencil(const ARdsNonSymPencil& other) { Copy(other); }
80 // Copy constructor.
81
~ARdsNonSymPencil()82 virtual ~ARdsNonSymPencil() { }
83 // Destructor.
84
85 ARdsNonSymPencil& operator=(const ARdsNonSymPencil& other);
86 // Assignment operator.
87
88 };
89
90 // ------------------------------------------------------------------------ //
91 // ARdsNonSymPencil member functions definition. //
92 // ------------------------------------------------------------------------ //
93
94
95 template<class ARTYPE, class ARFLOAT>
96 inline void ARdsNonSymPencil<ARTYPE, ARFLOAT>::
Copy(const ARdsNonSymPencil<ARTYPE,ARFLOAT> & other)97 Copy(const ARdsNonSymPencil<ARTYPE, ARFLOAT>& other)
98 {
99
100 part = other.part;
101 A = other.A;
102 B = other.B;
103 AsB = other.AsB;
104 #ifdef ARCOMP_H
105 AsBc = other.AsBc;
106 #endif
107
108 } // Copy.
109
110
111 template<class ARTYPE, class ARFLOAT>
FactorAsB(ARTYPE sigma)112 void ARdsNonSymPencil<ARTYPE, ARFLOAT>::FactorAsB(ARTYPE sigma)
113 {
114
115 // Quitting the function if A and B were not defined.
116
117 if (!(A->IsDefined()&&B->IsDefined())) {
118 throw ArpackError(ArpackError::DATA_UNDEFINED,
119 "ARdsNonSymPencil::FactorAsB");
120 }
121
122 // Quitting the function if A and B are not square.
123
124 if ((A->nrows() != A->ncols()) || (B->nrows() != B->ncols())) {
125 throw ArpackError(ArpackError::NOT_SQUARE_MATRIX,
126 "ARdsNonSymPencil::FactorAsB");
127 }
128
129 // Copying A to AsB if sigma = 0.
130
131 if (sigma == (ARTYPE)0) {
132
133 AsB = *A;
134 if (!AsB.IsFactored()) AsB.FactorA();
135 return;
136
137 }
138
139 // Defining matrix AsB.
140
141 if (!AsB.IsDefined()) {
142 AsB.DefineMatrix(A->ncols(), A->A);
143 }
144
145 // Reserving memory for some vectors used in matrix decomposition.
146
147 AsB.CreateStructure();
148
149 // Subtracting sigma*B from A and storing the result on AsB.
150
151 ::copy(A->m*A->n, A->A, 1, AsB.Ainv, 1);
152 axpy(A->m*A->n, -sigma, B->A, 1, AsB.Ainv, 1);
153
154 // Decomposing AsB.
155
156 getrf(AsB.m, AsB.n, AsB.Ainv, AsB.m, AsB.ipiv, AsB.info);
157
158 // Handling errors.
159
160 AsB.ThrowError();
161
162 AsB.factored = true;
163
164 } // FactorAsB (ARTYPE shift).
165
166
167 #ifdef ARCOMP_H
168 template<class ARTYPE, class ARFLOAT>
169 void ARdsNonSymPencil<ARTYPE, ARFLOAT>::
FactorAsB(ARFLOAT sigmaR,ARFLOAT sigmaI,char partp)170 FactorAsB(ARFLOAT sigmaR, ARFLOAT sigmaI, char partp)
171 {
172
173 // Quitting the function if A and B were not defined.
174
175 if (!(A->IsDefined()&&B->IsDefined())) {
176 throw ArpackError(ArpackError::DATA_UNDEFINED,
177 "ARdsNonSymPencil::FactorAsB");
178 }
179
180 // Quitting the function if A and B are not square.
181
182 if ((A->nrows() != A->ncols()) || (B->nrows() != B->ncols())) {
183 throw ArpackError(ArpackError::NOT_SQUARE_MATRIX,
184 "ARdsNonSymPencil::FactorAsB");
185 }
186
187 // Defining matrix AsB.
188
189 if (!AsBc.IsDefined()) {
190 part = partp;
191 AsBc.DefineMatrix(A->ncols(), 0);
192 }
193
194 // Reserving memory for some vectors used in matrix decomposition.
195
196 AsBc.CreateStructure();
197
198 // Subtracting sigma*B from A and storing the result on AsBc.
199
200 arcomplex<ARFLOAT> sigma(sigmaR, sigmaI);
201 for (int i=0; i<(A->m*A->n); i++) AsBc.Ainv[i] = A->A[i]-sigma*B->A[i];
202
203 // Decomposing AsBc.
204
205 getrf(AsBc.m, AsBc.n, AsBc.Ainv, AsBc.m, AsBc.ipiv, AsBc.info);
206
207 // Handling errors.
208
209 AsBc.ThrowError();
210
211 AsBc.factored = true;
212
213 } // FactorAsB (arcomplex<ARFLOAT> shift).
214 #endif // ARCOMP_H.
215
216
217 template<class ARTYPE, class ARFLOAT>
MultInvBAv(ARTYPE * v,ARTYPE * w)218 void ARdsNonSymPencil<ARTYPE, ARFLOAT>::MultInvBAv(ARTYPE* v, ARTYPE* w)
219 {
220
221 if (!B->IsFactored()) B->FactorA();
222
223 A->MultMv(v, w);
224 B->MultInvv(w, w);
225
226 } // MultInvBAv.
227
228
229 #ifdef ARCOMP_H
230
231 template<class ARTYPE, class ARFLOAT>
232 void ARdsNonSymPencil<ARTYPE, ARFLOAT>::
MultInvAsBv(arcomplex<ARFLOAT> * v,arcomplex<ARFLOAT> * w)233 MultInvAsBv(arcomplex<ARFLOAT>* v, arcomplex<ARFLOAT>* w)
234 {
235
236 AsB.MultInvv((ARTYPE*)v,(ARTYPE*)w);
237
238 } // MultInvAsBv (arcomplex<ARFLOAT>).
239
240 #endif // ARCOMP_H.
241
242
243 template<class ARTYPE, class ARFLOAT>
MultInvAsBv(ARFLOAT * v,ARFLOAT * w)244 void ARdsNonSymPencil<ARTYPE, ARFLOAT>::MultInvAsBv(ARFLOAT* v, ARFLOAT* w)
245 {
246
247 if (part == 'N') { // shift is real.
248
249 AsB.MultInvv((ARTYPE*)v,(ARTYPE*)w);
250
251 }
252 else { // shift is complex.
253
254 #ifdef ARCOMP_H
255
256 int i;
257 arcomplex<ARFLOAT> *tv, *tw;
258
259 tv = new arcomplex<ARFLOAT>[AsBc.ncols()];
260 tw = new arcomplex<ARFLOAT>[AsBc.ncols()];
261
262 for (i=0; i!=AsBc.ncols(); i++) tv[i] = arcomplex<ARFLOAT>(v[i], 0.0);
263
264 AsBc.MultInvv(tv, tw);
265
266 if (part=='I') {
267 for (i=0; i!=AsBc.ncols(); i++) w[i] = imag(tw[i]);
268 }
269 else {
270 for (i=0; i!=AsBc.ncols(); i++) w[i] = real(tw[i]);
271 }
272
273 delete[] tv;
274 delete[] tw;
275
276 #endif // ARCOMP_H.
277
278 }
279
280 } // MultInvAsBv (ARFLOAT).
281
282
283 template<class ARTYPE, class ARFLOAT>
284 inline void ARdsNonSymPencil<ARTYPE, ARFLOAT>::
DefineMatrices(ARdsNonSymMatrix<ARTYPE,ARFLOAT> & Ap,ARdsNonSymMatrix<ARTYPE,ARFLOAT> & Bp)285 DefineMatrices(ARdsNonSymMatrix<ARTYPE, ARFLOAT>& Ap,
286 ARdsNonSymMatrix<ARTYPE, ARFLOAT>& Bp)
287 {
288
289 A = &Ap;
290 B = &Bp;
291
292 if ((A->n != B->n)||(A->m != B->m)) {
293 throw ArpackError(ArpackError::INCOMPATIBLE_SIZES,
294 "ARdsNonSymMatrix::DefineMatrices");
295 }
296
297 } // DefineMatrices.
298
299
300 template<class ARTYPE, class ARFLOAT>
301 inline ARdsNonSymPencil<ARTYPE, ARFLOAT>::
ARdsNonSymPencil(ARdsNonSymMatrix<ARTYPE,ARFLOAT> & Ap,ARdsNonSymMatrix<ARTYPE,ARFLOAT> & Bp)302 ARdsNonSymPencil(ARdsNonSymMatrix<ARTYPE, ARFLOAT>& Ap,
303 ARdsNonSymMatrix<ARTYPE, ARFLOAT>& Bp)
304 {
305
306 DefineMatrices(Ap, Bp);
307
308 } // Long constructor.
309
310
311 template<class ARTYPE, class ARFLOAT>
312 ARdsNonSymPencil<ARTYPE, ARFLOAT>& ARdsNonSymPencil<ARTYPE, ARFLOAT>::
313 operator=(const ARdsNonSymPencil<ARTYPE, ARFLOAT>& other)
314 {
315
316 if (this != &other) { // Stroustrup suggestion.
317 Copy(other);
318 }
319 return *this;
320
321 } // operator=.
322
323
324 #endif // ARDNSPEN_H
325