1 /*
2 ARPACK++ v1.2 2/20/2000
3 c++ interface to ARPACK code.
4
5 MODULE ARLNSPen.h.
6 Arpack++ class ARluNonSymPencil 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 ARLNSPEN_H
18 #define ARLNSPEN_H
19
20 #include <cstddef>
21
22 #include "arch.h"
23 #include "arerror.h"
24 #include "blas1c.h"
25 #include "superluc.h"
26 #include "arlspdef.h"
27 #include "arlutil.h"
28 #include "arlnsmat.h"
29
30
31 template<class ARTYPE, class ARFLOAT>
32 class ARluNonSymPencil
33 {
34
35 protected:
36
37 bool factored;
38 int* permc;
39 int* permr;
40 char part;
41 ARluNonSymMatrix<ARTYPE, ARFLOAT>* A;
42 ARluNonSymMatrix<ARTYPE, ARFLOAT>* B;
43 SuperMatrix L;
44 SuperMatrix U;
45 SuperLUStat_t stat;
46
47 virtual void Copy(const ARluNonSymPencil& other);
48
49 void ClearMem();
50
51 void SparseSaxpy(ARTYPE a, ARTYPE x[], int xind[], int nx, ARTYPE y[],
52 int yind[], int ny, ARTYPE z[], int zind[], int& nz);
53
54 #ifdef ARCOMP_H
55 void SparseSaxpy(arcomplex<ARFLOAT> a, ARFLOAT x[], int xind[], int nx,
56 ARFLOAT y[], int yind[], int ny, arcomplex<ARFLOAT> z[],
57 int zind[], int& nz);
58 #endif
59
60 void SubtractAsB(int n, ARTYPE sigma, NCformat& A,
61 NCformat& B, NCformat& AsB);
62
63 #ifdef ARCOMP_H
64 void SubtractAsB(int n, ARFLOAT sigmaR, ARFLOAT sigmaI,
65 NCformat& A, NCformat& B, NCformat& AsB);
66 #endif
67
68 public:
69
IsFactored()70 bool IsFactored() { return factored; }
71
72 void FactorAsB(ARTYPE sigma);
73
74 #ifdef ARCOMP_H
75 void FactorAsB(ARFLOAT sigmaR, ARFLOAT sigmaI, char partp = 'R');
76 #endif
77
MultAv(ARTYPE * v,ARTYPE * w)78 void MultAv(ARTYPE* v, ARTYPE* w) { A->MultMv(v,w); }
79
MultBv(ARTYPE * v,ARTYPE * w)80 void MultBv(ARTYPE* v, ARTYPE* w) { B->MultMv(v,w); }
81
82 void MultInvBAv(ARTYPE* v, ARTYPE* w);
83
84 #ifdef ARCOMP_H
85 void MultInvAsBv(arcomplex<ARFLOAT>* v, arcomplex<ARFLOAT>* w);
86 #endif
87
88 void MultInvAsBv(ARFLOAT* v, ARFLOAT* w);
89
90 void DefineMatrices(ARluNonSymMatrix<ARTYPE, ARFLOAT>& Ap,
91 ARluNonSymMatrix<ARTYPE, ARFLOAT>& Bp);
92
93 ARluNonSymPencil();
94 // Short constructor that does nothing.
95
96 ARluNonSymPencil(ARluNonSymMatrix<ARTYPE, ARFLOAT>& Ap,
97 ARluNonSymMatrix<ARTYPE, ARFLOAT>& Bp);
98 // Long constructor.
99
ARluNonSymPencil(const ARluNonSymPencil & other)100 ARluNonSymPencil(const ARluNonSymPencil& other) { Copy(other); }
101 // Copy constructor.
102
~ARluNonSymPencil()103 virtual ~ARluNonSymPencil() { ClearMem(); }
104 // Destructor.
105
106 ARluNonSymPencil& operator=(const ARluNonSymPencil& other);
107 // Assignment operator.
108
109 };
110
111 // ------------------------------------------------------------------------ //
112 // ARluNonSymPencil member functions definition. //
113 // ------------------------------------------------------------------------ //
114
115
116 template<class ARTYPE, class ARFLOAT>
117 inline void ARluNonSymPencil<ARTYPE, ARFLOAT>::
Copy(const ARluNonSymPencil<ARTYPE,ARFLOAT> & other)118 Copy(const ARluNonSymPencil<ARTYPE, ARFLOAT>& other)
119 {
120
121 factored = other.factored;
122 part = other.part;
123 A = other.A;
124 B = other.B;
125
126 // Throwing the original factorization away (this procedure
127 // is really awkward, but it is necessary because there
128 // is no copy function for matrices L and U in the SuperLU
129 // library and it is not a good idea to do this kind of deep
130 // copy here).
131
132 if (factored) {
133 ArpackError(ArpackError::DISCARDING_FACTORS, "ARluNonSymPencil");
134 factored = false;
135 }
136
137 } // Copy.
138
139
140 template<class ARTYPE, class ARFLOAT>
ClearMem()141 void ARluNonSymPencil<ARTYPE, ARFLOAT>::ClearMem()
142 {
143
144 if (factored) {
145 Destroy_SuperNode_Matrix(&L);
146 Destroy_CompCol_Matrix(&U);
147 StatFree(&stat);
148 delete[] permc;
149 delete[] permr;
150 permc = NULL;
151 permr = NULL;
152 }
153
154 } // ClearMem.
155
156
157 template<class ARTYPE, class ARFLOAT>
158 void ARluNonSymPencil<ARTYPE, ARFLOAT>::
SparseSaxpy(ARTYPE a,ARTYPE x[],int xind[],int nx,ARTYPE y[],int yind[],int ny,ARTYPE z[],int zind[],int & nz)159 SparseSaxpy(ARTYPE a, ARTYPE x[], int xind[], int nx, ARTYPE y[],
160 int yind[], int ny, ARTYPE z[], int zind[], int& nz)
161 // A strongly sequential (and inefficient) sparse saxpy algorithm.
162 {
163
164 int ix, iy;
165
166 nz = 0;
167 if ((nx == 0) || (a == (ARTYPE)0)) {
168 copy(ny,y,1,z,1);
169 for (iy=0; iy!=ny; iy++) zind[iy] = yind[iy];
170 nz = ny;
171 return;
172 }
173 if (ny == 0) {
174 copy(nx,x,1,z,1);
175 scal(nx,a,z,1);
176 for (ix=0; ix!=nx; ix++) zind[ix] = xind[ix];
177 nz = nx;
178 return;
179 }
180 ix = 0;
181 iy = 0;
182 while (true) {
183 if (xind[ix] == yind[iy]) {
184 zind[nz] = xind[ix];
185 z[nz++] = a*x[ix++]+y[iy++];
186 if ((ix == nx)||(iy == ny)) break;
187 }
188 else if (xind[ix] < yind[iy]) {
189 zind[nz] = xind[ix];
190 z[nz++] = a*x[ix++];
191 if (ix == nx) break;
192 }
193 else {
194 zind[nz] = yind[iy];
195 z[nz++] = y[iy++];
196 if (iy == ny) break;
197 }
198 }
199 while (iy < ny) {
200 zind[nz] = yind[iy];
201 z[nz++] = y[iy++];
202 }
203 while (ix < nx) {
204 zind[nz] = xind[ix];
205 z[nz++] = x[ix++];
206 }
207
208 } // SparseSaxpy (ARTYPE).
209
210
211 #ifdef ARCOMP_H
212 template<class ARTYPE, class ARFLOAT>
213 void ARluNonSymPencil<ARTYPE, ARFLOAT>::
SparseSaxpy(arcomplex<ARFLOAT> a,ARFLOAT x[],int xind[],int nx,ARFLOAT y[],int yind[],int ny,arcomplex<ARFLOAT> z[],int zind[],int & nz)214 SparseSaxpy(arcomplex<ARFLOAT> a, ARFLOAT x[], int xind[], int nx, ARFLOAT y[],
215 int yind[], int ny, arcomplex<ARFLOAT> z[], int zind[], int& nz)
216 // A strongly sequential (and inefficient) sparse saxpy algorithm.
217 {
218
219 int ix, iy;
220
221 nz = 0;
222 if ((nx == 0) || (a == arcomplex<ARFLOAT>(0.0,0.0))) {
223 for (iy=0; iy!=ny; iy++) {
224 z[iy] = arcomplex<ARFLOAT>(y[iy],0.0);
225 zind[iy] = yind[iy];
226 }
227 nz = ny;
228 return;
229 }
230 if (ny == 0) {
231 for (ix=0; ix!=ny; ix++) {
232 z[ix] = a*arcomplex<ARFLOAT>(x[ix],0.0);
233 zind[ix] = xind[ix];
234 }
235 nz = nx;
236 return;
237 }
238 ix = 0;
239 iy = 0;
240 while (true) {
241 if (xind[ix] == yind[iy]) {
242 zind[nz] = xind[ix];
243 z[nz++] = a*x[ix++]+y[iy++];
244 if ((ix == nx)||(iy == ny)) break;
245 }
246 else if (xind[ix] < yind[iy]) {
247 zind[nz] = xind[ix];
248 z[nz++] = a*x[ix++];
249 if (ix == nx) break;
250 }
251 else {
252 zind[nz] = yind[iy];
253 z[nz++] = arcomplex<ARFLOAT>(y[iy++],0.0);
254 if (iy == ny) break;
255 }
256 }
257 while (iy < ny) {
258 zind[nz] = yind[iy];
259 z[nz++] = arcomplex<ARFLOAT>(y[iy++],0.0);
260 }
261 while (ix < nx) {
262 zind[nz] = xind[ix];
263 z[nz++] = arcomplex<ARFLOAT>(x[ix++],0.0);
264 }
265
266 } // SparseSaxpy (arcomplex<ARFLOAT>).
267 #endif // ARCOMP_H.
268
269
270 template<class ARTYPE, class ARFLOAT>
271 void ARluNonSymPencil<ARTYPE, ARFLOAT>::
SubtractAsB(int n,ARTYPE sigma,NCformat & A,NCformat & B,NCformat & AsB)272 SubtractAsB(int n, ARTYPE sigma, NCformat& A, NCformat& B, NCformat& AsB)
273 {
274
275 int i, acol, bcol, asbcol, scol;
276 ARTYPE* anzval;
277 ARTYPE* bnzval;
278 ARTYPE* asbnzval;
279
280 // Telling the compiler that nzval must ve viewed as a vector of ARTYPE.
281
282 anzval = (ARTYPE*)A.nzval;
283 bnzval = (ARTYPE*)B.nzval;
284 asbnzval = (ARTYPE*)AsB.nzval;
285
286 // Subtracting sigma*B from A.
287
288 AsB.colptr[0] = 0;
289 asbcol = 0;
290
291 for (i=0; i!=n; i++) {
292 bcol = B.colptr[i];
293 acol = A.colptr[i];
294 SparseSaxpy(-sigma, &bnzval[bcol], &B.rowind[bcol], B.colptr[i+1]-bcol,
295 &anzval[acol], &A.rowind[acol], A.colptr[i+1]-acol,
296 &asbnzval[asbcol], &AsB.rowind[asbcol], scol);
297 asbcol += scol;
298 AsB.colptr[i+1] = asbcol;
299 }
300
301 AsB.nnz = AsB.colptr[n];
302
303 } // SubtractAsB (ARTYPE shift).
304
305
306 #ifdef ARCOMP_H
307 template<class ARTYPE, class ARFLOAT>
308 void ARluNonSymPencil<ARTYPE, ARFLOAT>::
SubtractAsB(int n,ARFLOAT sigmaR,ARFLOAT sigmaI,NCformat & A,NCformat & B,NCformat & AsB)309 SubtractAsB(int n, ARFLOAT sigmaR, ARFLOAT sigmaI,
310 NCformat& A, NCformat& B, NCformat& AsB)
311 {
312
313 int i, acol, bcol, asbcol, scol;
314 ARTYPE* anzval;
315 ARTYPE* bnzval;
316 arcomplex<ARFLOAT>* asbnzval;
317 arcomplex<ARFLOAT> sigma;
318
319 // Telling the compiler that nzval must ve viewed as a vector of ARTYPE.
320
321 anzval = (ARTYPE*)A.nzval;
322 bnzval = (ARTYPE*)B.nzval;
323 asbnzval = (arcomplex<ARFLOAT>*)AsB.nzval;
324
325 // Subtracting sigma*B from A.
326
327 sigma = arcomplex<ARFLOAT>(sigmaR, sigmaI);
328 AsB.colptr[0] = 0;
329 asbcol = 0;
330
331 for (i=0; i!=n; i++) {
332 bcol = B.colptr[i];
333 acol = A.colptr[i];
334 SparseSaxpy(-sigma, &bnzval[bcol], &B.rowind[bcol], B.colptr[i+1]-bcol,
335 &anzval[acol], &A.rowind[acol], A.colptr[i+1]-acol,
336 &asbnzval[asbcol], &AsB.rowind[asbcol], scol);
337 asbcol += scol;
338 AsB.colptr[i+1] = asbcol;
339 }
340
341 AsB.nnz = AsB.colptr[n];
342
343 } // SubtractAsB (arcomplex<ARFLOAT> shift).
344 #endif // ARCOMP_H.
345
346
347 template<class ARTYPE, class ARFLOAT>
FactorAsB(ARTYPE sigma)348 void ARluNonSymPencil<ARTYPE, ARFLOAT>::FactorAsB(ARTYPE sigma)
349 {
350
351 // Quitting the function if A and B were not defined.
352
353 if (!(A->IsDefined()&&B->IsDefined())) {
354 throw ArpackError(ArpackError::DATA_UNDEFINED,
355 "ARluNonSymPencil::FactorAsB");
356 }
357
358 // Quitting the function if A and B are not square.
359
360 if ((A->nrows() != A->ncols()) || (B->nrows() != B->ncols())) {
361 throw ArpackError(ArpackError::NOT_SQUARE_MATRIX,
362 "ARluNonSymPencil::FactorAsB");
363 }
364
365 // Defining local variables.
366
367 int nnzi, info;
368 int* etree;
369 int* irowi;
370 int* pcoli;
371 ARTYPE* asb;
372 SuperMatrix AsB;
373 SuperMatrix AC;
374 NCformat* Astore;
375 NCformat* Bstore;
376 NCformat* AsBstore;
377
378 // Deleting old versions of L, U, perm_r and perm_c.
379
380 ClearMem();
381
382 // Setting default values for gstrf parameters.
383
384 int panel_size = sp_ienv(1);
385 int relax = sp_ienv(2);
386 superlu_options_t options;
387 /* Set the default input options:
388 options.Fact = DOFACT;
389 options.Equil = YES;
390 options.ColPerm = COLAMD;
391 options.DiagPivotThresh = 1.0;
392 options.Trans = NOTRANS;
393 options.IterRefine = NOREFINE;
394 options.SymmetricMode = NO;
395 options.PivotGrowth = NO;
396 options.ConditionNumber = NO;
397 options.PrintStat = YES;
398 */
399 set_default_options(&options);
400 options.DiagPivotThresh = A->threshold;
401
402 // Defining A and B format.
403
404 Astore = (NCformat*)A->A.Store;
405 Bstore = (NCformat*)B->A.Store;
406
407 // Creating a temporary matrix AsB.
408
409 nnzi = Astore->nnz+Bstore->nnz;
410 irowi = new int[nnzi];
411 pcoli = new int[A->ncols()+1];
412 asb = new ARTYPE[nnzi];
413 Create_CompCol_Matrix(&AsB, A->nrows(), A->ncols(), nnzi, asb,
414 irowi, pcoli, SLU_NC, SLU_GE);
415
416 // Subtracting sigma*B from A and storing the result on AsB.
417
418 AsBstore = (NCformat*)AsB.Store;
419 SubtractAsB(A->ncols(), sigma, *Astore, *Bstore, *AsBstore);
420
421 // Reserving memory for some vectors used in matrix decomposition.
422
423 etree = new int[A->ncols()];
424 if (permc == NULL) permc = new int[A->ncols()];
425 if (permr == NULL) permr = new int[A->ncols()];
426
427 // Defining LUStat.
428
429 // StatInit(panel_size, relax);
430 SuperLUStat_t stat;
431 StatInit(&stat);
432
433 // Defining the column permutation of matrix AsB
434 // (using minimum degree ordering on AsB'*AsB).
435
436 get_perm_c(A->order, &AsB, permc);
437
438 // Permuting columns of AsB and
439 // creating the elimination tree of AsB'*AsB.
440
441 // sp_preorder("N", &AsB, permc, etree, &AC);
442 sp_preorder(&options, &AsB, permc, etree, &AC);
443
444 // Decomposing AsB.
445
446 // gstrf("N",&AC, A->threshold, drop_tol, relax, panel_size, etree,
447 // NULL, 0, permr, permc, &L, &U, &info);
448 gstrf(&options, &AC, relax, panel_size, etree,
449 NULL, 0, permc, permr, &L, &U, &stat, &info);
450
451 // Deleting AC, AsB and etree.
452
453 Destroy_CompCol_Permuted(&AC);
454 Destroy_CompCol_Matrix(&AsB);
455 delete[] etree;
456
457 factored = (info == 0);
458
459 // Handling errors.
460
461 if (info < 0) { // Illegal argument.
462 throw ArpackError(ArpackError::PARAMETER_ERROR,
463 "ARluNonSymPencil::FactorAsB");
464 }
465 else if (info > A->ncols()) { // Memory is not sufficient.
466 throw ArpackError(ArpackError::MEMORY_OVERFLOW,
467 "ARluNonSymPencil::FactorAsB");
468 }
469 else if (info > 0) { // Matrix is singular.
470 throw ArpackError(ArpackError::MATRIX_IS_SINGULAR,
471 "ARluNonSymPencil::FactorAsB");
472 }
473
474 } // FactorAsB (ARTYPE shift).
475
476
477 #ifdef ARCOMP_H
478 template<class ARTYPE, class ARFLOAT>
479 void ARluNonSymPencil<ARTYPE, ARFLOAT>::
FactorAsB(ARFLOAT sigmaR,ARFLOAT sigmaI,char partp)480 FactorAsB(ARFLOAT sigmaR, ARFLOAT sigmaI, char partp)
481 {
482
483 // Quitting the function if A and B were not defined.
484
485 if (!(A->IsDefined()&&B->IsDefined())) {
486 throw ArpackError(ArpackError::DATA_UNDEFINED,
487 "ARluNonSymPencil::FactorAsB");
488 }
489
490 // Quitting the function if A and B are not square.
491
492 if ((A->nrows() != A->ncols()) || (B->nrows() != B->ncols())) {
493 throw ArpackError(ArpackError::NOT_SQUARE_MATRIX,
494 "ARluNonSymPencil::FactorAsB");
495 }
496
497 // Defining local variables.
498
499 int nnzi, info;
500 int* etree;
501 int* irowi;
502 int* pcoli;
503 arcomplex<ARFLOAT>* asb;
504 SuperMatrix AsB;
505 SuperMatrix AC;
506 NCformat* Astore;
507 NCformat* Bstore;
508 NCformat* AsBstore;
509
510 // Deleting old versions of L, U, perm_r and perm_c.
511
512 ClearMem();
513
514 // Setting default values for gstrf parameters.
515
516 int panel_size = sp_ienv(1);
517 int relax = sp_ienv(2);
518 superlu_options_t options;
519 /* Set the default input options:
520 options.Fact = DOFACT;
521 options.Equil = YES;
522 options.ColPerm = COLAMD;
523 options.DiagPivotThresh = 1.0;
524 options.Trans = NOTRANS;
525 options.IterRefine = NOREFINE;
526 options.SymmetricMode = NO;
527 options.PivotGrowth = NO;
528 options.ConditionNumber = NO;
529 options.PrintStat = YES;
530 */
531 set_default_options(&options);
532 options.DiagPivotThresh = A->threshold;
533
534
535 // Defining A and B format.
536
537 Astore = (NCformat*)A->A.Store;
538 Bstore = (NCformat*)B->A.Store;
539
540 // Creating a temporary matrix AsB.
541
542 part = partp;
543 nnzi = Astore->nnz+Bstore->nnz;
544 irowi = new int[nnzi];
545 pcoli = new int[A->ncols()+1];
546 asb = new arcomplex<ARFLOAT>[nnzi];
547 Create_CompCol_Matrix(&AsB, A->nrows(), A->ncols(), nnzi, asb,
548 irowi, pcoli, SLU_NC, SLU_GE);
549
550 // Subtracting sigma*B from A and storing the result on AsB.
551
552 AsBstore = (NCformat*)AsB.Store;
553 SubtractAsB(A->ncols(), sigmaR, sigmaI, *Astore, *Bstore, *AsBstore);
554
555 // Reserving memory for some vectors used in matrix decomposition.
556
557 etree = new int[A->ncols()];
558 if (permc == NULL) permc = new int[A->ncols()];
559 if (permr == NULL) permr = new int[A->ncols()];
560
561 // Defining LUStat.
562
563 // StatInit(panel_size, relax);
564 SuperLUStat_t stat;
565 StatInit(&stat);
566
567 // Defining the column permutation of matrix AsB
568 // (using minimum degree ordering on AsB'*AsB).
569
570 get_perm_c(A->order, &AsB, permc);
571
572 // Permuting columns of AsB and
573 // creating the elimination tree of AsB'*AsB.
574
575 //sp_preorder("N", &AsB, permc, etree, &AC);
576 sp_preorder(&options, &AsB, permc, etree, &AC);
577
578 // Decomposing AsB.
579
580 // gstrf("N",&AC, A->threshold, drop_tol, relax, panel_size, etree, NULL,
581 // 0, permr, permc, &L, &U, &info);
582 gstrf(&options, &AC, relax, panel_size, etree,
583 NULL, 0, permc, permr, &L, &U, &stat, &info);
584
585 // Deleting AC, AsB and etree.
586
587 Destroy_CompCol_Permuted(&AC);
588 Destroy_CompCol_Matrix(&AsB);
589 delete[] etree;
590
591 factored = (info == 0);
592
593 // Handling errors.
594
595 if (info < 0) { // Illegal argument.
596 throw ArpackError(ArpackError::PARAMETER_ERROR,
597 "ARluNonSymPencil::FactorAsB");
598 }
599 else if (info > A->ncols()) { // Memory is not sufficient.
600 throw ArpackError(ArpackError::MEMORY_OVERFLOW,
601 "ARluNonSymPencil::FactorAsB");
602 }
603 else if (info > 0) { // Matrix is singular.
604 throw ArpackError(ArpackError::MATRIX_IS_SINGULAR,
605 "ARluNonSymPencil::FactorAsB");
606 }
607
608 } // FactorAsB (arcomplex<ARFLOAT> shift).
609 #endif // ARCOMP_H.
610
611
612 template<class ARTYPE, class ARFLOAT>
MultInvBAv(ARTYPE * v,ARTYPE * w)613 void ARluNonSymPencil<ARTYPE, ARFLOAT>::MultInvBAv(ARTYPE* v, ARTYPE* w)
614 {
615
616 if (!B->IsFactored()) B->FactorA();
617
618 A->MultMv(v, w);
619 B->MultInvv(w, w);
620
621 } // MultInvBAv.
622
623
624 #ifdef ARCOMP_H
625
626 template<class ARTYPE, class ARFLOAT>
627 void ARluNonSymPencil<ARTYPE, ARFLOAT>::
MultInvAsBv(arcomplex<ARFLOAT> * v,arcomplex<ARFLOAT> * w)628 MultInvAsBv(arcomplex<ARFLOAT>* v, arcomplex<ARFLOAT>* w)
629 {
630
631 // Quitting the function if AsB was not factored.
632
633 if (!IsFactored()) {
634 throw ArpackError(ArpackError::NOT_FACTORED_MATRIX,
635 "ARluNonSymPencil::MultInvAsBv");
636 }
637
638 // Solving AsB.w = v.
639
640 int info;
641 SuperMatrix RHS;
642
643 copy(A->nrows(), v, 1, w, 1);
644 Create_Dense_Matrix(&RHS, A->nrows(), 1, w, A->nrows(), SLU_DN, SLU_GE);
645 // gstrs("N", &L, &U, permr, permc, &RHS, &info);
646 trans_t trans = NOTRANS;
647 StatInit(&stat);
648
649 gstrs(trans, &L, &U, permc, permr, &RHS, &stat, &info);
650
651 Destroy_SuperMatrix_Store(&RHS); // delete RHS.Store;
652
653 } // MultInvAsBv (arcomplex<ARFLOAT>).
654
655 #endif
656
657
658 template<class ARTYPE, class ARFLOAT>
MultInvAsBv(ARFLOAT * v,ARFLOAT * w)659 void ARluNonSymPencil<ARTYPE, ARFLOAT>::MultInvAsBv(ARFLOAT* v, ARFLOAT* w)
660 {
661
662 // Quitting the function if AsB was not factored.
663
664 if (!IsFactored()) {
665 throw ArpackError(ArpackError::NOT_FACTORED_MATRIX,
666 "ARluNonSymPencil::MultInvAsBv");
667 }
668
669 // Solving AsB.w = v.
670
671 int info;
672 SuperMatrix RHS;
673
674 if (part == 'N') { // shift is real.
675
676 copy(A->nrows(), v, 1, w, 1);
677 Create_Dense_Matrix(&RHS, A->nrows(), 1, w, A->nrows(), SLU_DN, SLU_GE);
678 //gstrs("N", &L, &U, permr, permc, &RHS, &info);
679 trans_t trans = NOTRANS;
680 StatInit(&stat);
681 gstrs(trans, &L, &U, permc, permr, &RHS, &stat, &info);
682
683 }
684 else { // shift is complex.
685
686 #ifdef ARCOMP_H
687
688 int i;
689 arcomplex<ARFLOAT> *tv = new arcomplex<ARFLOAT>[A->ncols()];
690
691 for (i=0; i!=A->ncols(); i++) tv[i] = arcomplex<ARFLOAT>(v[i],0.0);
692 Create_Dense_Matrix(&RHS, A->ncols(), 1, tv, A->ncols(), SLU_DN, SLU_GE);
693 //gstrs("N", &L, &U, permr, permc, &RHS, &info);
694 trans_t trans = NOTRANS;
695 StatInit(&stat);
696 gstrs(trans, &L, &U, permc, permr, &RHS, &stat, &info);
697
698
699 if (part=='I') {
700 for (i=0; i!=A->ncols(); i++) w[i] = imag(tv[i]);
701 }
702 else {
703 for (i=0; i!=A->ncols(); i++) w[i] = real(tv[i]);
704 }
705
706 delete[] tv;
707
708 #endif
709
710 }
711
712 Destroy_SuperMatrix_Store(&RHS); // delete RHS.Store;
713
714 } // MultInvAsBv (ARFLOAT).
715
716
717 template<class ARTYPE, class ARFLOAT>
718 inline void ARluNonSymPencil<ARTYPE, ARFLOAT>::
DefineMatrices(ARluNonSymMatrix<ARTYPE,ARFLOAT> & Ap,ARluNonSymMatrix<ARTYPE,ARFLOAT> & Bp)719 DefineMatrices(ARluNonSymMatrix<ARTYPE, ARFLOAT>& Ap,
720 ARluNonSymMatrix<ARTYPE, ARFLOAT>& Bp)
721 {
722
723 A = &Ap;
724 B = &Bp;
725 permc = NULL;
726 permr = NULL;
727
728 if ((A->n != B->n)||(A->m != B->m)) {
729 throw ArpackError(ArpackError::INCOMPATIBLE_SIZES,
730 "ARluNonSymMatrix::DefineMatrices");
731 }
732
733 } // DefineMatrices.
734
735
736 template<class ARTYPE, class ARFLOAT>
ARluNonSymPencil()737 inline ARluNonSymPencil<ARTYPE, ARFLOAT>::ARluNonSymPencil()
738 {
739
740 factored = false;
741 part = 'N';
742 permr = NULL;
743 permc = NULL;
744
745 } // Short constructor.
746
747
748 template<class ARTYPE, class ARFLOAT>
749 inline ARluNonSymPencil<ARTYPE, ARFLOAT>::
ARluNonSymPencil(ARluNonSymMatrix<ARTYPE,ARFLOAT> & Ap,ARluNonSymMatrix<ARTYPE,ARFLOAT> & Bp)750 ARluNonSymPencil(ARluNonSymMatrix<ARTYPE, ARFLOAT>& Ap,
751 ARluNonSymMatrix<ARTYPE, ARFLOAT>& Bp)
752 {
753
754 factored = false;
755 DefineMatrices(Ap, Bp);
756
757 } // Long constructor.
758
759
760 template<class ARTYPE, class ARFLOAT>
761 ARluNonSymPencil<ARTYPE, ARFLOAT>& ARluNonSymPencil<ARTYPE, ARFLOAT>::
762 operator=(const ARluNonSymPencil<ARTYPE, ARFLOAT>& other)
763 {
764
765 if (this != &other) { // Stroustrup suggestion.
766 ClearMem();
767 Copy(other);
768 }
769 return *this;
770
771 } // operator=.
772
773
774 #endif // ARLNSPEN_H
775