1 /*
2 ARPACK++ v1.2 2/20/2000
3 c++ interface to ARPACK code.
4
5 MODULE ARBSPen.h.
6 Arpack++ class ARbdSymPencil 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 ARBSPEN_H
18 #define ARBSPEN_H
19
20 #include "arch.h"
21 #include "arerror.h"
22 #include "blas1c.h"
23 #include "lapackc.h"
24 #include "arbsmat.h"
25
26
27 template<class ARTYPE>
28 class ARbdSymPencil
29 {
30
31 protected:
32
33 ARbdSymMatrix<ARTYPE>* A;
34 ARbdSymMatrix<ARTYPE>* B;
35 ARbdSymMatrix<ARTYPE> AsB;
36
max(int a,int b)37 int max(int a, int b) { return (a>b)?a:b; }
38
min(int a,int b)39 int min(int a, int b) { return (a<b)?a:b; }
40
41 virtual void Copy(const ARbdSymPencil& other);
42
43 void SubtractAsB(ARTYPE sigma);
44
45 public:
46
IsFactored()47 bool IsFactored() { return AsB.IsFactored(); }
48
49 void FactorAsB(ARTYPE sigma);
50
MultAv(ARTYPE * v,ARTYPE * w)51 void MultAv(ARTYPE* v, ARTYPE* w) { A->MultMv(v,w); }
52
MultBv(ARTYPE * v,ARTYPE * w)53 void MultBv(ARTYPE* v, ARTYPE* w) { B->MultMv(v,w); }
54
55 void MultInvBAv(ARTYPE* v, ARTYPE* w);
56
MultInvAsBv(ARTYPE * v,ARTYPE * w)57 void MultInvAsBv(ARTYPE* v, ARTYPE* w) { AsB.MultInvv(v,w); }
58
59 void DefineMatrices(ARbdSymMatrix<ARTYPE>& Ap, ARbdSymMatrix<ARTYPE>& Bp);
60
ARbdSymPencil()61 ARbdSymPencil() { AsB.factored = false; }
62 // Short constructor that does nothing.
63
64 ARbdSymPencil(ARbdSymMatrix<ARTYPE>& Ap, ARbdSymMatrix<ARTYPE>& Bp);
65 // Long constructor.
66
ARbdSymPencil(const ARbdSymPencil & other)67 ARbdSymPencil(const ARbdSymPencil& other) { Copy(other); }
68 // Copy constructor.
69
~ARbdSymPencil()70 virtual ~ARbdSymPencil() { }
71 // Destructor.
72
73 ARbdSymPencil& operator=(const ARbdSymPencil& other);
74 // Assignment operator.
75
76 };
77
78 // ------------------------------------------------------------------------ //
79 // ARbdSymPencil member functions definition. //
80 // ------------------------------------------------------------------------ //
81
82
83 template<class ARTYPE>
Copy(const ARbdSymPencil<ARTYPE> & other)84 inline void ARbdSymPencil<ARTYPE>::Copy(const ARbdSymPencil<ARTYPE>& other)
85 {
86
87 A = other.A;
88 B = other.B;
89 AsB = other.AsB;
90
91 } // Copy.
92
93
94 template<class ARTYPE>
SubtractAsB(ARTYPE sigma)95 void ARbdSymPencil<ARTYPE>::SubtractAsB(ARTYPE sigma)
96 {
97
98 int i, n, minD, ndA, ndB, ndAsB, lda;
99 ARTYPE negsig;
100
101 negsig = -sigma;
102 n = AsB.n;
103 ndA = A->nsdiag;
104 ndB = B->nsdiag;
105 ndAsB = AsB.nsdiag;
106 lda = AsB.lda;
107
108 // Expanding A.
109
110 if (A->uplo == 'U') {
111
112 // Copying the main diagonal of A.
113
114 copy(n, &A->A[ndA], ndA+1, &AsB.Ainv[2*ndAsB], lda);
115
116 // Copying the superdiagonals of A.
117
118 for (i = 0; i < ndA; i++) {
119 copy(n, &A->A[i], ndA+1, &AsB.Ainv[2*ndAsB-ndA+i], lda);
120 copy(n-ndA+i, &A->A[i+(ndA-i)*(ndA+1)], ndA+1,
121 &AsB.Ainv[2*ndAsB+ndA-i], lda);
122 }
123
124 }
125 else {
126
127 // Copying the main diagonal of A to Ainv.
128
129 copy(n, &A->A[0], ndA+1, &AsB.Ainv[2*ndAsB], lda);
130
131 // Copying the subdiagonals of A to Ainv.
132
133 for (i = 1; i <= ndA; i++) {
134 copy(n, &A->A[i], ndA+1, &AsB.Ainv[2*ndAsB+i], lda);
135 copy(n-i, &A->A[i], ndA+1, &AsB.Ainv[2*ndAsB-i+i*lda], lda);
136 }
137
138 }
139
140 // Transferring part of B (*(-sigma)) if AsB.nsdiag > A->nsdiag.
141
142 if (A->nsdiag < AsB.nsdiag) {
143
144 if (B->uplo == 'U') {
145
146 for (i = 0; i < ndAsB-ndA; i++) {
147 copy(n, &B->A[i], ndB+1, &AsB.Ainv[ndAsB+i], lda);
148 scal(n, negsig, &AsB.Ainv[ndAsB+i], lda);
149 copy(n-ndAsB+i, &AsB.Ainv[ndAsB+i+(ndAsB-i)*lda], lda,
150 &AsB.Ainv[lda-i-1], lda);
151 }
152
153 }
154 else {
155
156 for (i = ndA+1; i <= ndAsB; i++) {
157 copy(n, &B->A[i], ndB+1, &AsB.Ainv[2*ndAsB+i], lda);
158 scal(n, negsig, &AsB.Ainv[2*ndAsB+i], lda);
159 copy(n-i, &AsB.Ainv[2*ndAsB+i], lda,
160 &AsB.Ainv[2*ndAsB-i+i*lda], lda);
161 }
162
163 }
164
165 }
166
167 // Subtracting sigma*B from A.
168
169 minD = min(ndA, ndB);
170
171 if (B->uplo == 'U') {
172
173 // Subtracting the main diagonal of B.
174
175 axpy(n, negsig, &B->A[ndB], ndB+1, &AsB.Ainv[2*ndAsB], lda);
176
177 // Subtracting the superdiagonals.
178
179 for (i = 0; i < minD; i++) {
180 axpy(n, negsig, &B->A[ndB-minD+i], ndB+1,
181 &AsB.Ainv[2*ndAsB-minD+i], lda);
182 copy(n-minD+i, &AsB.Ainv[2*ndAsB-minD+i+(minD-i)*lda], lda,
183 &AsB.Ainv[2*ndAsB+minD-i], lda);
184 }
185
186 }
187 else {
188
189 // Subtracting the main diagonal of B.
190
191 axpy(n, negsig, &B->A[0], ndB+1, &AsB.Ainv[2*ndAsB], lda);
192
193 // Subtracting the subdiagonals.
194
195 for (i = 1; i <= minD; i++) {
196 axpy(n, negsig, &B->A[i], ndB+1, &AsB.Ainv[2*ndAsB+i], lda);
197 copy(n-i, &AsB.Ainv[2*ndAsB+i], lda,
198 &AsB.Ainv[2*ndAsB-i+i*lda], lda);
199 }
200
201 }
202
203 } // SubtractAsB (ARTYPE shift).
204
205
206 template<class ARTYPE>
FactorAsB(ARTYPE sigma)207 void ARbdSymPencil<ARTYPE>::FactorAsB(ARTYPE sigma)
208 {
209
210 // Quitting the function if A and B were not defined.
211
212 if (!(A->IsDefined()&&B->IsDefined())) {
213 throw ArpackError(ArpackError::DATA_UNDEFINED,
214 "ARbdSymPencil::FactorAsB");
215 }
216
217 // Copying A to AsB if sigma = 0.
218
219 if (sigma == (ARTYPE)0) {
220
221 AsB = *A;
222 if (!AsB.IsFactored()) AsB.FactorA();
223 return;
224
225 }
226
227 // Defining matrix AsB.
228
229 if (!AsB.IsDefined()) {
230 AsB.DefineMatrix(A->ncols(), max(A->nsdiag, B->nsdiag), A->A);
231 }
232
233 // Reserving memory for some vectors used in matrix decomposition.
234
235 AsB.CreateStructure();
236
237 // Subtracting sigma*B from A and storing the result on AsB.
238
239 SubtractAsB(sigma);
240
241 // Decomposing AsB.
242
243 gbtrf(AsB.n, AsB.n, AsB.nsdiag, AsB.nsdiag,
244 AsB.Ainv, AsB.lda, AsB.ipiv, AsB.info);
245
246 // Handling errors.
247
248 AsB.ThrowError();
249
250 AsB.factored = true;
251
252 } // FactorAsB (ARTYPE shift).
253
254
255 template<class ARTYPE>
MultInvBAv(ARTYPE * v,ARTYPE * w)256 void ARbdSymPencil<ARTYPE>::MultInvBAv(ARTYPE* v, ARTYPE* w)
257 {
258
259 if (!B->IsFactored()) B->FactorA();
260
261 A->MultMv(v, w);
262 copy(A->ncols(), w, 1, v, 1);
263 B->MultInvv(w, w);
264
265 } // MultInvBAv.
266
267
268 template<class ARTYPE>
269 inline void ARbdSymPencil<ARTYPE>::
DefineMatrices(ARbdSymMatrix<ARTYPE> & Ap,ARbdSymMatrix<ARTYPE> & Bp)270 DefineMatrices(ARbdSymMatrix<ARTYPE>& Ap, ARbdSymMatrix<ARTYPE>& Bp)
271 {
272
273 A = &Ap;
274 B = &Bp;
275
276 } // DefineMatrices.
277
278
279 template<class ARTYPE>
280 inline ARbdSymPencil<ARTYPE>::
ARbdSymPencil(ARbdSymMatrix<ARTYPE> & Ap,ARbdSymMatrix<ARTYPE> & Bp)281 ARbdSymPencil(ARbdSymMatrix<ARTYPE>& Ap, ARbdSymMatrix<ARTYPE>& Bp)
282 {
283
284 AsB.factored = false;
285 DefineMatrices(Ap, Bp);
286
287 } // Long constructor.
288
289
290 template<class ARTYPE>
291 ARbdSymPencil<ARTYPE>& ARbdSymPencil<ARTYPE>::
292 operator=(const ARbdSymPencil<ARTYPE>& other)
293 {
294
295 if (this != &other) { // Stroustrup suggestion.
296 Copy(other);
297 }
298 return *this;
299
300 } // operator=.
301
302
303 #endif // ARBSPEN_H
304