1 /*
2 ARPACK++ v1.2 2/20/2000
3 c++ interface to ARPACK code.
4
5 MODULE ARBNSPen.h.
6 Arpack++ class ARbdNonSymPencil 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 ARBNSPEN_H
18 #define ARBNSPEN_H
19
20 #include "arch.h"
21 #include "arerror.h"
22 #include "blas1c.h"
23 #include "lapackc.h"
24 #include "arbnsmat.h"
25
26
27 template<class ARTYPE, class ARFLOAT>
28 class ARbdNonSymPencil
29 {
30
31 protected:
32
33 char part;
34 ARbdNonSymMatrix<ARTYPE, ARFLOAT>* A;
35 ARbdNonSymMatrix<ARTYPE, ARFLOAT>* B;
36 ARbdNonSymMatrix<ARTYPE, ARFLOAT> AsB;
37 #ifdef ARCOMP_H
38 ARbdNonSymMatrix<arcomplex<ARFLOAT>, ARFLOAT> AsBc;
39 #endif
40
max(int a,int b)41 int max(int a, int b) { return (a>b)?a:b; }
42
min(int a,int b)43 int min(int a, int b) { return (a<b)?a:b; }
44
45 void ComplexCopy(int n, ARFLOAT dx[], int incx,
46 arcomplex<ARFLOAT> dy[], int incy);
47
48 void ComplexAxpy(int n, arcomplex<ARFLOAT> da, ARTYPE dx[],
49 int incx, arcomplex<ARFLOAT> dy[], int incy);
50
51 virtual void Copy(const ARbdNonSymPencil& other);
52
53 void SubtractAsB(ARTYPE sigma);
54
55 #ifdef ARCOMP_H
56 void SubtractAsB(ARFLOAT sigmaR, ARFLOAT sigmaI);
57 #endif
58
59 public:
60
61 #ifdef ARCOMP_H
IsFactored()62 bool IsFactored() { return (AsB.IsFactored()||AsBc.IsFactored()); }
63 #else
64 bool IsFactored() { return AsB.IsFactored(); }
65 #endif
66
67 void FactorAsB(ARTYPE sigma);
68
69 #ifdef ARCOMP_H
70 void FactorAsB(ARFLOAT sigmaR, ARFLOAT sigmaI, char partp = 'R');
71 #endif
72
MultAv(ARTYPE * v,ARTYPE * w)73 void MultAv(ARTYPE* v, ARTYPE* w) { A->MultMv(v,w); }
74
MultBv(ARTYPE * v,ARTYPE * w)75 void MultBv(ARTYPE* v, ARTYPE* w) { B->MultMv(v,w); }
76
77 void MultInvBAv(ARTYPE* v, ARTYPE* w);
78
79 #ifdef ARCOMP_H
80 void MultInvAsBv(arcomplex<ARFLOAT>* v, arcomplex<ARFLOAT>* w);
81 #endif
82
83 void MultInvAsBv(ARFLOAT* v, ARFLOAT* w);
84
85 void DefineMatrices(ARbdNonSymMatrix<ARTYPE, ARFLOAT>& Ap,
86 ARbdNonSymMatrix<ARTYPE, ARFLOAT>& Bp);
87
ARbdNonSymPencil()88 ARbdNonSymPencil() { part = 'N'; }
89 // Short constructor that does nothing.
90
91 ARbdNonSymPencil(ARbdNonSymMatrix<ARTYPE, ARFLOAT>& Ap,
92 ARbdNonSymMatrix<ARTYPE, ARFLOAT>& Bp);
93 // Long constructor.
94
ARbdNonSymPencil(const ARbdNonSymPencil & other)95 ARbdNonSymPencil(const ARbdNonSymPencil& other) { Copy(other); }
96 // Copy constructor.
97
~ARbdNonSymPencil()98 virtual ~ARbdNonSymPencil() { }
99 // Destructor.
100
101 ARbdNonSymPencil& operator=(const ARbdNonSymPencil& other);
102 // Assignment operator.
103
104 };
105
106 // ------------------------------------------------------------------------ //
107 // ARbdNonSymPencil member functions definition. //
108 // ------------------------------------------------------------------------ //
109
110
111 template<class ARTYPE, class ARFLOAT>
112 inline void ARbdNonSymPencil<ARTYPE, ARFLOAT>::
Copy(const ARbdNonSymPencil<ARTYPE,ARFLOAT> & other)113 Copy(const ARbdNonSymPencil<ARTYPE, ARFLOAT>& other)
114 {
115
116 part = other.part;
117 A = other.A;
118 B = other.B;
119 AsB = other.AsB;
120 #ifdef ARCOMP_H
121 AsBc = other.AsBc;
122 #endif
123
124 } // Copy.
125
126
127 template<class ARTYPE, class ARFLOAT>
128 void ARbdNonSymPencil<ARTYPE, ARFLOAT>::
ComplexCopy(int n,ARFLOAT dx[],int incx,arcomplex<ARFLOAT> dy[],int incy)129 ComplexCopy(int n, ARFLOAT dx[], int incx, arcomplex<ARFLOAT> dy[], int incy)
130 {
131
132 for (int ix=0, iy=0; ix<(n*incx); ix+=incx, iy+=incy) {
133 dy[iy] = arcomplex<ARFLOAT>(dx[ix], 0.0);
134 }
135
136 } // ComplexCopy.
137
138
139 template<class ARTYPE, class ARFLOAT>
140 void ARbdNonSymPencil<ARTYPE, ARFLOAT>::
ComplexAxpy(int n,arcomplex<ARFLOAT> da,ARTYPE dx[],int incx,arcomplex<ARFLOAT> dy[],int incy)141 ComplexAxpy(int n, arcomplex<ARFLOAT> da, ARTYPE dx[], int incx,
142 arcomplex<ARFLOAT> dy[], int incy)
143 {
144
145 for (int ix=0, iy=0; ix<(n*incx); ix+=incx, iy+=incy) {
146 dy[iy] += da*dx[ix];
147 }
148
149 } // ComplexAxpy.
150
151
152 template<class ARTYPE, class ARFLOAT>
SubtractAsB(ARTYPE sigma)153 void ARbdNonSymPencil<ARTYPE, ARFLOAT>::SubtractAsB(ARTYPE sigma)
154 {
155
156 int i, inca, incb, minL, minU, begB, begAsB;
157 ARTYPE negsig;
158
159 inca = A->ndiagL+A->ndiagU+1;
160 incb = B->ndiagL+B->ndiagU+1;
161 negsig = -sigma;
162
163 // Expanding A.
164
165 begAsB = AsB.ndiagL+AsB.ndiagU-A->ndiagU;
166 for (i = 0; i < inca; i++) {
167 copy(AsB.n, &A->A[i], inca, &AsB.Ainv[begAsB+i], AsB.lda);
168 }
169
170 // Transferring part of B (*(-sigma)) if AsB.ndiagU > A->ndiagU.
171
172 if (A->ndiagU < AsB.ndiagU) {
173 for (i = 0; i < AsB.ndiagU-A->ndiagU; i++) {
174 copy(AsB.n, &B->A[i], incb, &AsB.Ainv[AsB.ndiagL+i], AsB.lda);
175 scal(AsB.n, negsig, &AsB.Ainv[AsB.ndiagL+i], AsB.lda);
176 }
177 }
178
179 // Subtracting sigma*B from A.
180
181 minL = min(A->ndiagL, B->ndiagL);
182 minU = min(A->ndiagU, B->ndiagU);
183 begB = B->ndiagU-minU;
184 begAsB = AsB.ndiagL+AsB.ndiagU-minU;
185
186 for (i = 0; i < minL+minU+1; i++) {
187 axpy(AsB.n, -sigma, &B->A[begB+i], incb, &AsB.Ainv[begAsB+i], AsB.lda);
188 }
189
190 // Transferring part of B (*(-sigma)) if AsB.ndiagL > A->ndiagL.
191
192 if (A->ndiagL < AsB.ndiagL) {
193 begB = B->ndiagU+1+minL;
194 begAsB = AsB.ndiagL+AsB.ndiagU+1+minL;
195 for (i = 0; i < AsB.ndiagL-A->ndiagL; i++) {
196 copy(AsB.n, &B->A[begB+i], incb, &AsB.Ainv[begAsB+i], AsB.lda);
197 scal(AsB.n, negsig, &AsB.Ainv[begAsB+i], AsB.lda);
198 }
199 }
200
201 } // SubtractAsB (ARTYPE shift).
202
203
204 #ifdef ARCOMP_H
205 template<class ARTYPE, class ARFLOAT>
206 void ARbdNonSymPencil<ARTYPE, ARFLOAT>::
SubtractAsB(ARFLOAT sigmaR,ARFLOAT sigmaI)207 SubtractAsB(ARFLOAT sigmaR, ARFLOAT sigmaI)
208 {
209
210 int i, inca, incb, minL, minU, begB, begAsB;
211 arcomplex<ARFLOAT> sigma;
212
213 inca = A->ndiagL+A->ndiagU+1;
214 incb = B->ndiagL+B->ndiagU+1;
215 sigma = arcomplex<ARFLOAT>(sigmaR, sigmaI);
216
217 // Expanding A.
218
219 begAsB = AsBc.ndiagL+AsBc.ndiagU-A->ndiagU;
220 for (i = 0; i < inca; i++) {
221 ComplexCopy(AsBc.n,(ARFLOAT*)(&A->A[i]),inca,&AsBc.Ainv[begAsB+i],AsBc.lda);
222 }
223
224 // Transferring part of B (*(-sigma)) if AsBc.ndiagU > A->ndiagU.
225
226 if (A->ndiagU < AsBc.ndiagU) {
227 for (i = 0; i < AsBc.ndiagU-A->ndiagU; i++) {
228 ComplexCopy(AsBc.n, (ARFLOAT*)(&B->A[i]), incb,
229 &AsBc.Ainv[AsBc.ndiagL+i], AsBc.lda);
230 scal(AsBc.n, -sigma, &AsBc.Ainv[AsBc.ndiagL+i], AsBc.lda);
231 }
232 }
233
234 // Subtracting sigma*B from A.
235
236 minL = min(A->ndiagL, B->ndiagL);
237 minU = min(A->ndiagU, B->ndiagU);
238 begB = B->ndiagU-minU;
239 begAsB = AsBc.ndiagL+AsBc.ndiagU-minU;
240
241 for (i = 0; i < minL+minU+1; i++) {
242 ComplexAxpy(AsBc.n, -sigma, &B->A[begB+i], incb,
243 &AsBc.Ainv[begAsB+i], AsBc.lda);
244 }
245
246 // Transferring part of B (*(-sigma)) if AsBc.ndiagL > A->ndiagL.
247
248 if (A->ndiagL < AsBc.ndiagL) {
249 begB = B->ndiagU+1+minL;
250 begAsB = AsBc.ndiagL+AsBc.ndiagU+1+minL;
251 for (i = 0; i < AsBc.ndiagL-A->ndiagL; i++) {
252 ComplexCopy(AsBc.n, (ARFLOAT*)(&B->A[begB+i]), incb,
253 &AsBc.Ainv[begAsB+i], AsBc.lda);
254 scal(AsBc.n, -sigma, &AsBc.Ainv[begAsB+i], AsBc.lda);
255 }
256 }
257
258 } // SubtractAsB (arcomplex<ARFLOAT> shift).
259 #endif // ARCOMP_H
260
261
262 template<class ARTYPE, class ARFLOAT>
FactorAsB(ARTYPE sigma)263 void ARbdNonSymPencil<ARTYPE, ARFLOAT>::FactorAsB(ARTYPE sigma)
264 {
265
266 // Quitting the function if A and B were not defined.
267
268 if (!(A->IsDefined()&&B->IsDefined())) {
269 throw ArpackError(ArpackError::DATA_UNDEFINED,
270 "ARbdNonSymPencil::FactorAsB");
271 }
272
273 // Quitting the function if A and B are not square.
274
275 if ((A->nrows() != A->ncols()) || (B->nrows() != B->ncols())) {
276 throw ArpackError(ArpackError::NOT_SQUARE_MATRIX,
277 "ARbdNonSymPencil::FactorAsB");
278 }
279
280 // Copying A to AsB if sigma = 0.
281
282 if (sigma == (ARTYPE)0) {
283
284 AsB = *A;
285 if (!AsB.IsFactored()) AsB.FactorA();
286 return;
287
288 }
289
290 // Defining matrix AsB.
291
292 if (!AsB.IsDefined()) {
293 AsB.DefineMatrix(A->ncols(), max(A->ndiagL, B->ndiagL),
294 max(A->ndiagU, B->ndiagU), A->A);
295 }
296
297 // Reserving memory for some vectors used in matrix decomposition.
298
299 AsB.CreateStructure();
300
301 // Subtracting sigma*B from A and storing the result on AsB.
302
303 SubtractAsB(sigma);
304
305 // Decomposing AsB.
306
307 gbtrf(AsB.n, AsB.n, AsB.ndiagL, AsB.ndiagU,
308 AsB.Ainv, AsB.lda, AsB.ipiv, AsB.info);
309
310 // Handling errors.
311
312 AsB.ThrowError();
313
314 AsB.factored = true;
315
316 } // FactorAsB (ARTYPE shift).
317
318
319 #ifdef ARCOMP_H
320 template<class ARTYPE, class ARFLOAT>
321 void ARbdNonSymPencil<ARTYPE, ARFLOAT>::
FactorAsB(ARFLOAT sigmaR,ARFLOAT sigmaI,char partp)322 FactorAsB(ARFLOAT sigmaR, ARFLOAT sigmaI, char partp)
323 {
324
325 // Quitting the function if A and B were not defined.
326
327 if (!(A->IsDefined()&&B->IsDefined())) {
328 throw ArpackError(ArpackError::DATA_UNDEFINED,
329 "ARbdNonSymPencil::FactorAsB");
330 }
331
332 // Quitting the function if A and B are not square.
333
334 if ((A->nrows() != A->ncols()) || (B->nrows() != B->ncols())) {
335 throw ArpackError(ArpackError::NOT_SQUARE_MATRIX,
336 "ARbdNonSymPencil::FactorAsB");
337 }
338
339 // Defining matrix AsBc.
340
341 if (!AsBc.IsDefined()) {
342 part = partp;
343 AsBc.DefineMatrix(A->ncols(), max(A->ndiagL,B->ndiagL),
344 max(A->ndiagU,B->ndiagU), 0);
345 }
346
347 // Reserving memory for some vectors used in matrix decomposition.
348
349 AsBc.CreateStructure();
350
351 // Subtracting sigma*B from A and storing the result on AsBc.
352
353 SubtractAsB(sigmaR, sigmaI);
354
355 // Decomposing AsBc.
356
357 gbtrf(AsBc.n, AsBc.n, AsBc.ndiagL, AsBc.ndiagU,
358 AsBc.Ainv, AsBc.lda, AsBc.ipiv, AsBc.info);
359
360 // Handling errors.
361
362 AsBc.ThrowError();
363
364 AsBc.factored = true;
365
366 } // FactorAsB (arcomplex<ARFLOAT> shift).
367 #endif // ARCOMP_H.
368
369
370 template<class ARTYPE, class ARFLOAT>
MultInvBAv(ARTYPE * v,ARTYPE * w)371 void ARbdNonSymPencil<ARTYPE, ARFLOAT>::MultInvBAv(ARTYPE* v, ARTYPE* w)
372 {
373
374 if (!B->IsFactored()) B->FactorA();
375
376 A->MultMv(v, w);
377 B->MultInvv(w, w);
378
379 } // MultInvBAv.
380
381
382 #ifdef ARCOMP_H
383 template<class ARTYPE, class ARFLOAT>
384 void ARbdNonSymPencil<ARTYPE, ARFLOAT>::
MultInvAsBv(arcomplex<ARFLOAT> * v,arcomplex<ARFLOAT> * w)385 MultInvAsBv(arcomplex<ARFLOAT>* v, arcomplex<ARFLOAT>* w)
386 {
387
388 AsB.MultInvv((ARTYPE*)v, (ARTYPE*)w);
389
390 } // MultInvAsBv (arcomplex<ARFLOAT>).
391 #endif // ARCOMP_H.
392
393
394 template<class ARTYPE, class ARFLOAT>
MultInvAsBv(ARFLOAT * v,ARFLOAT * w)395 void ARbdNonSymPencil<ARTYPE, ARFLOAT>::MultInvAsBv(ARFLOAT* v, ARFLOAT* w)
396 {
397
398 if (part == 'N') { // shift is real.
399
400 AsB.MultInvv((ARTYPE*)v, (ARTYPE*)w);
401
402 }
403 else { // shift is complex.
404
405 #ifdef ARCOMP_H
406
407 int i;
408 arcomplex<ARFLOAT> *tv, *tw;
409
410 tv = new arcomplex<ARFLOAT>[AsBc.ncols()];
411 tw = new arcomplex<ARFLOAT>[AsBc.ncols()];
412
413 for (i=0; i!=AsBc.ncols(); i++) tv[i] = arcomplex<ARFLOAT>(v[i], 0.0);
414
415 AsBc.MultInvv(tv, tw);
416
417 if (part=='I') {
418 for (i=0; i!=AsBc.ncols(); i++) w[i] = imag(tw[i]);
419 }
420 else {
421 for (i=0; i!=AsBc.ncols(); i++) w[i] = real(tw[i]);
422 }
423
424 delete[] tv;
425 delete[] tw;
426
427 #endif // ARCOMP_H.
428
429 }
430
431 } // MultInvAsBv (ARFLOAT).
432
433
434 template<class ARTYPE, class ARFLOAT>
435 inline void ARbdNonSymPencil<ARTYPE, ARFLOAT>::
DefineMatrices(ARbdNonSymMatrix<ARTYPE,ARFLOAT> & Ap,ARbdNonSymMatrix<ARTYPE,ARFLOAT> & Bp)436 DefineMatrices(ARbdNonSymMatrix<ARTYPE, ARFLOAT>& Ap,
437 ARbdNonSymMatrix<ARTYPE, ARFLOAT>& Bp)
438 {
439
440 A = &Ap;
441 B = &Bp;
442
443 } // DefineMatrices.
444
445
446 template<class ARTYPE, class ARFLOAT>
447 inline ARbdNonSymPencil<ARTYPE, ARFLOAT>::
ARbdNonSymPencil(ARbdNonSymMatrix<ARTYPE,ARFLOAT> & Ap,ARbdNonSymMatrix<ARTYPE,ARFLOAT> & Bp)448 ARbdNonSymPencil(ARbdNonSymMatrix<ARTYPE, ARFLOAT>& Ap,
449 ARbdNonSymMatrix<ARTYPE, ARFLOAT>& Bp)
450 {
451
452 DefineMatrices(Ap, Bp);
453
454 } // Long constructor.
455
456
457 template<class ARTYPE, class ARFLOAT>
458 ARbdNonSymPencil<ARTYPE, ARFLOAT>& ARbdNonSymPencil<ARTYPE, ARFLOAT>::
459 operator=(const ARbdNonSymPencil<ARTYPE, ARFLOAT>& other)
460 {
461
462 if (this != &other) { // Stroustrup suggestion.
463 Copy(other);
464 }
465 return *this;
466
467 } // operator=.
468
469
470 #endif // ARBNSPEN_H
471