1 /*
2 ARPACK++ v1.2 2/20/2000
3 c++ interface to ARPACK code.
4
5 MODULE ARRSEig.h.
6 Arpack++ base class ARrcStdEig definition.
7 This class implements c++ version of the reverse
8 communication interface for standard problems.
9
10 ARPACK Authors
11 Richard Lehoucq
12 Danny Sorensen
13 Chao Yang
14 Dept. of Computational & Applied Mathematics
15 Rice University
16 Houston, Texas
17 */
18
19 #ifndef ARRSEIG_H
20 #define ARRSEIG_H
21
22 #include <new>
23 #include <cstddef>
24 #include <string>
25 #include "arch.h"
26 #include "arerror.h"
27 #include "debug.h"
28 #include "blas1c.h"
29
30
31 // "New" handler.
32
MemoryOverflow()33 void MemoryOverflow() { throw ArpackError(ArpackError::MEMORY_OVERFLOW); }
34
35 // ARrcStdEig class definition.
36
37 template<class ARFLOAT, class ARTYPE>
38 class ARrcStdEig {
39
40 protected:
41
42 // a) Protected variables:
43
44 // a.1) User defined parameters.
45
46 int n; // Dimension of the eigenproblem.
47 int nev; // Number of eigenvalues to be computed. 0 < nev < n-1.
48 int ncv; // Number of Arnoldi vectors generated at each iteration.
49 int maxit; // Maximum number of Arnoldi update iterations allowed.
50 std::string which; // Specify which of the Ritz values of OP to compute.
51 ARFLOAT tol; // Stopping criterion (relative accuracy of Ritz values).
52 ARFLOAT sigmaI; // Imaginary part of shift (for nonsymmetric problems).
53 ARTYPE sigmaR; // Shift (real part only if problem is nonsymmetric).
54 ARTYPE *resid; // Initial residual vector.
55
56
57 // a.2) Internal variables.
58
59 bool rvec; // Indicates if eigenvectors/Schur vectors were
60 // requested (or only eigenvalues will be determined).
61 bool newRes; // Indicates if a new "resid" vector was created.
62 bool newVal; // Indicates if a new "EigValR" vector was created.
63 bool newVec; // Indicates if a new "EigVec" vector was created.
64 bool PrepareOK; // Indicates if internal variables were correctly set.
65 bool BasisOK; // Indicates if an Arnoldi basis was found.
66 bool ValuesOK; // Indicates if eigenvalues were calculated.
67 bool VectorsOK; // Indicates if eigenvectors were determined.
68 bool SchurOK; // Indicates if Schur vectors were determined.
69 bool AutoShift; // Indicates if implicit shifts will be generated
70 // internally (or will be supplied by the user).
71 char bmat; // Indicates if the problem is a standard ('I') or
72 // generalized ('G") eigenproblem.
73 char HowMny; // Indicates if eigenvectors ('A') or Schur vectors ('P')
74 // were requested (not referenced if rvec = false).
75 int ido; // Original ARPACK reverse communication flag.
76 int info; // Original ARPACK error flag.
77 int mode; // Indicates the type of the eigenproblem (regular,
78 // shift and invert, etc).
79 int lworkl; // Dimension of array workl.
80 int lworkv; // Dimension of array workv.
81 int lrwork; // Dimension of array rwork.
82 int iparam[12]; // Vector that handles original ARPACK parameters.
83 int ipntr[15]; // Vector that handles original ARPACK pointers.
84 ARFLOAT *rwork; // Original ARPACK internal vector.
85 ARTYPE *workl; // Original ARPACK internal vector.
86 ARTYPE *workd; // Original ARPACK internal vector.
87 ARTYPE *workv; // Original ARPACK internal vector.
88 ARTYPE *V; // Arnoldi basis / Schur vectors.
89
90
91 // a.3) Pure output variables.
92
93 int nconv; // Number of "converged" Ritz values.
94 ARFLOAT *EigValI; // Imaginary part of eigenvalues (nonsymmetric problems).
95 ARTYPE *EigValR; // Eigenvalues (real part only if problem is nonsymmetric).
96 ARTYPE *EigVec; // Eigenvectors.
97
98
99 // b) Protected functions:
100
101 // b.1) Memory control functions.
102
OverV()103 bool OverV() { return (EigVec == &V[1]); }
104 // Indicates whether EigVec overrides V or no.
105
ValSize()106 virtual int ValSize() { return nev; }
107 // Provides the size of array EigVal.
108 // Redefined in ARrcNonSymStdEig.
109
110 void ClearFirst();
111 // Clears some boolean variables in order to define a entire new problem.
112
113 void ClearBasis();
114 // Clear some boolean variables in order to recalculate the arnoldi basis.
115
116 void ClearMem();
117 // Clears workspace.
118
119 virtual void ValAllocate();
120 // Creates arrays EigValR and EigValI.
121 // Redefined in ARrcNonSymStdEig.
122
123 virtual void VecAllocate(bool newV = true);
124 // Creates array EigVec.
125
126 virtual void WorkspaceAllocate();
127 // Function that must be defined by a derived class.
128 // Redefined in ARrc[Sym|NonSym|Complex]StdEig.
129
130
131 // b.2) Functions that call the original ARPACK FORTRAN code.
132
Aupp()133 virtual void Aupp() {
134 throw ArpackError(ArpackError::NOT_IMPLEMENTED, "Aupp");
135 }
136 // Interface to FORTRAN subroutines __AUPD.
137 // Must be defined by a derived class.
138 // Redefined in ARrc[Sym|NonSym|Complex]StdEig.
139
140 void AuppError();
141 // Handles errors occurred in function Aupp.
142
Eupp()143 virtual void Eupp() {
144 throw ArpackError(ArpackError::NOT_IMPLEMENTED, "Eupp");
145 }
146 // Interface to FORTRAN subroutines __EUPD.
147 // Must be defined by a derived class.
148 // Redefined in ARrc[Sym|NonSym|Complex]Eig.
149
150 void EuppError();
151 // Handles errors occurred in function Eupp.
152
153
154 // b.3) Functions that check user defined parameters.
155
156 int CheckN(int np);
157 // Does range checking on ncv.
158
159 int CheckNcv(int ncvp);
160 // Forces ncv to conform to its ranges.
161
162 virtual int CheckNev(int nevp);
163 // Does range checking on nev.
164 // Redefined in ARrcNonSymStdEig.
165
166 int CheckMaxit(int maxitp);
167 // Forces maxit to be greater than zero.
168
169 virtual std::string CheckWhich(const std::string& whichp);
170 // Determines if the value of variable "which" is valid.
171 // Redefined in ARrcSymStdEig.
172
173
174 // b.4) Functions that set internal variables.
175
176 void Restart();
177
178 virtual void Prepare();
179 // Defines internal variables and allocates working arrays.
180
181 virtual void Copy(const ARrcStdEig& other);
182 // Makes a deep copy of "other" over "this" object.
183 // Old values are not deleted (this function is to be used
184 // by the copy constructor and the assignment operator only).
185
186 public:
187
188 // c) Public functions:
189
190 // c.1) Function that stores user defined parameters.
191
192 virtual void DefineParameters(int np, int nevp, const std::string& whichp="LM",
193 int ncvp=0, ARFLOAT tolp=0.0, int maxitp=0,
194 ARTYPE* residp=NULL, bool ishiftp=true);
195 // Set values of problem parameters (also called by constructors).
196 // Redefined in ARrcStdEig and ARrcGenEig.
197
198
199 // c.2) Functions that detect if output data is ready.
200
ParametersDefined()201 bool ParametersDefined() { return PrepareOK; }
202 // Indicates if all internal variables and arrays were defined.
203
ArnoldiBasisFound()204 bool ArnoldiBasisFound() { return BasisOK; }
205 // Indicates if an Arnoldi basis is available.
206
EigenvaluesFound()207 bool EigenvaluesFound() { return ValuesOK; }
208 // Indicates if the requested eigenvalues are available.
209
EigenvectorsFound()210 bool EigenvectorsFound() { return VectorsOK; }
211 // Indicates if the requested eigenvectors are available.
212
SchurVectorsFound()213 bool SchurVectorsFound() { return SchurOK; }
214 // Indicates if the Schur vectors are available.
215
216
217 // c.3) Functions that provides access to internal variables' values.
218
ConvergedEigenvalues()219 int ConvergedEigenvalues() { return nconv; }
220 // Provides the number of "converged" eigenvalues found so far.
221
GetMaxit()222 int GetMaxit() { return maxit; }
223 // Returns the maximum number of Arnoldi update iterations allowed.
224
225 int GetIter();
226 // Returns the actual number of Arnoldi iterations taken.
227
GetTol()228 ARFLOAT GetTol() { return tol; }
229 // Returns the stopping criterion.
230
GetN()231 int GetN() { return n; }
232 // Returns the dimension of the problem.
233
GetNev()234 int GetNev() { return nev; }
235 // Returns the number of eigenvalues to be computed.
236
GetNcv()237 int GetNcv() { return ncv; }
238 // Returns the number of Arnoldi vectors generated at each iteration..
239
GetWhich()240 const std::string& GetWhich() { return which; }
241 // Returns "which".
242
GetShift()243 ARTYPE GetShift() { return sigmaR; }
244 // Returns the shift (when in shift and invert mode).
245
GetAutoShift()246 bool GetAutoShift() { return AutoShift; }
247 // Returns "AutoShift".
248
GetMode()249 int GetMode() { return mode; }
250 // Returns the type of problem (regular, shift and invert, etc).
251
252
253 // c.4) Functions that allow changes in problem parameters.
254
255 void ChangeMaxit(int maxitp);
256 // Changes the maximum number of Arnoldi update iterations allowed.
257
258 void ChangeTol(ARFLOAT tolp);
259 // Changes the stopping criterion.
260
261 virtual void ChangeNev(int nevp);
262 // Changes the number of eigenvalues to be computed.
263
264 virtual void ChangeNcv(int ncvp);
265 // Changes the number of Arnoldi vectors generated at each iteration..
266
267 virtual void ChangeWhich(const std::string& whichp);
268 // Changes "which".
269
270 virtual void ChangeShift(ARTYPE sigmaRp);
271 // Turns the problem to shift-and-invert mode with shift defined by sigmaRp.
272 // Redefined in many derived classes.
273
274 virtual void NoShift();
275 // Turns the problem to regular mode.
276 // Redefined in ARrcGenEig.
277
278 void InvertAutoShift();
279 // Inverts "AutoShift".
280
SetRegularMode()281 virtual void SetRegularMode() { NoShift(); }
282 // Turns problem to regular mode.
283
SetShiftInvertMode(ARTYPE sigmaRp)284 virtual void SetShiftInvertMode(ARTYPE sigmaRp) { ChangeShift(sigmaRp); }
285 // Turns problem to regular mode.
286
287 // c.5) Trace functions.
288
Trace()289 virtual void Trace() {
290 throw ArpackError(ArpackError::NOT_IMPLEMENTED, "Trace");
291 }
292 // Turns on trace mode.
293 // Must be defined by a derived class.
294 // Redefined in ARrc[Sym|NonSym|Complex]StdEig.
295
NoTrace()296 void NoTrace() { TraceOff(); }
297 // Turns off trace mode.
298
299
300 // c.6) Functions that permit step by step execution of ARPACK.
301
GetNp()302 int GetNp() { return iparam[8]; }
303 // Returns the number of shifts that must be supplied after a call to
304 // TakeStep() when shifts are provided by the user (AutoShift = false).
305
GetIdo()306 int GetIdo() { return ido; }
307 // Indicates the type of operation the user must perform between two
308 // successive calls to function TakeStep().
309 // ido = -1: compute y = OP*x, where GetVector() gives a pointer to the
310 // vector x, and PutVector() indicates where to store y.
311 // ido = 1: compute y = OP*x, where GetVector() gives a pointer to the
312 // vector x, and PutVector() indicates where to store y.
313 // When solving generalized problems, a pointer to the product
314 // B*x is also available by using GetProd().
315 // ido = 2: compute y = B*x, where GetVector() gives a pointer to the
316 // vector x, and PutVector() indicates where to store y.
317 // ido = 3: compute shifts, where PutVector() indicates where to store them.
318
319 virtual ARTYPE* GetVector();
320 // When ido = -1, 1 or 2 and the user must perform a product in the form
321 // y <- M*x, this function indicates where x is stored. When ido = 3, this
322 // function indicates where the eigenvalues of the current Hessenberg
323 // matrix are located.
324
325 ARTYPE* GetProd();
326 // When ido = 1 and the user must perform a product in the form
327 // y <- M*B*x, this function indicates where B*x is stored.
328
329 virtual ARTYPE* PutVector();
330 // When ido = -1, 1 or 2 and the user must perform a product in the form
331 // y <- M*x, this function indicates where to store y. When ido = 3, this
332 // function indicates where to store the shifts.
333 // Redefined in ARrcSymStdEig.
334
335 virtual int TakeStep();
336 // Calls Aupp once if there is no Arnoldi basis available.
337
338
339 // c.7) Functions that perform final calculations.
340
341 virtual int FindArnoldiBasis();
342 // Determines the Arnoldi basis related to the given problem.
343 // This function has no meaning for this class. It is maintained
344 // here for compatibility with all derived classes.
345 // Redefined in ARStdEig, ARGenEig and ARSymGenEig.
346
347 virtual int FindEigenvalues();
348 // Determines nev approximated eigenvalues of the given eigen-problem.
349 // Redefined in ARNonSymGenEig.
350
351 virtual int FindEigenvectors(bool schurp = false);
352 // Determines nev approximated eigenvectors of the given eigen-problem
353 // Optionally also determines nev Schur vectors that span the desired
354 // invariant subspace.
355 // Redefined in ARNonSymGenEig.
356
357 virtual int FindSchurVectors();
358 // Determines nev Schur vectors that span the desired invariant subspace.
359 // Redefined in ARrcSymStdEig and ARNonSymGenEig.
360
361
362 // c.8) Function that perform calculations using user supplied data structure.
363
364 virtual int Eigenvectors(ARTYPE* &EigVecp, bool ischur = false);
365 // Overrides array EigVecp sequentially with the eigenvectors of the
366 // given eigen-problem. Also calculates Schur vectors if requested.
367
368
369 // c.9) Functions that return elements of vectors and matrices.
370
371 ARTYPE ArnoldiBasisVector(int i, int j);
372 // Furnishes element j of the i-eth Arnoldi basis vector.
373
374 ARTYPE SchurVector(int i, int j);
375 // Furnishes element j of the i-eth Schur vector.
376
377 ARTYPE ResidualVector(int i);
378 // Furnishes element j of the residual vector.
379
380
381 // c.10) Functions that provide raw access to internal vectors and matrices.
382
383 ARTYPE* RawArnoldiBasisVectors();
384 // Provides raw access to Arnoldi basis vectors elements.
385
386 ARTYPE* RawArnoldiBasisVector(int i);
387 // Provides raw access to Arnoldi basis vector i.
388
389 ARTYPE* RawEigenvalues();
390 // Provides raw access to eigenvalues.
391
392 ARTYPE* RawEigenvectors();
393 // Provides raw access to eigenvectors elements.
394
395 ARTYPE* RawEigenvector(int i);
396 // Provides raw access to eigenvector i.
397
398 ARTYPE* RawSchurVectors();
399 // Provides raw access to Schur vectors elements.
400
401 ARTYPE* RawSchurVector(int i);
402 // Provides raw access to Schur vector i.
403
404 ARTYPE* RawResidualVector();
405 // Provides raw access to residual vector elements.
406
407
408 // c.11) Functions that use STL vector class.
409
410 #ifdef STL_VECTOR_H
411
412 vector<ARTYPE>* StlArnoldiBasisVectors();
413 // Returns a copy of the Arnoldi basis in a single STL vector.
414
415 vector<ARTYPE>* StlArnoldiBasisVector(int i);
416 // Returns the i-th Arnoldi basis vector in a STL vector.
417
418 vector<ARTYPE>* StlEigenvectors(bool ischur = false);
419 // Calculates the eigenvectors and stores them sequentially in a
420 // single STL vector. Also calculates Schur vectors if requested.
421
422 vector<ARTYPE>* StlSchurVectors();
423 // Returns a copy of the Schur vectors in a single STL vector.
424
425 vector<ARTYPE>* StlSchurVector(int i);
426 // Returns the i-th Schur vector in a STL vector.
427
428 vector<ARTYPE>* StlResidualVector();
429 // Returns the residual vector in a STL vector.
430
431 #endif // #ifdef STL_VECTOR_H.
432
433
434 // c.12) Constructors and destructor.
435
436 ARrcStdEig();
437 // Short constructor that does almost nothing.
438
ARrcStdEig(const ARrcStdEig & other)439 ARrcStdEig(const ARrcStdEig& other) { Copy(other); }
440 // Copy constructor.
441
~ARrcStdEig()442 virtual ~ARrcStdEig() { ClearMem(); }
443 // Very simple destructor.
444
445 // d) Operators:
446
447 ARrcStdEig& operator=(const ARrcStdEig& other);
448 // Assignment operator.
449
450 }; // class ARrcStdEig.
451
452
453 // ------------------------------------------------------------------------ //
454 // ARrcStdEig member functions definition. //
455 // ------------------------------------------------------------------------ //
456
457
458 template<class ARFLOAT, class ARTYPE>
ClearFirst()459 inline void ARrcStdEig<ARFLOAT, ARTYPE>::ClearFirst()
460 {
461
462 PrepareOK = newVal = newVec = newRes = false;
463
464 } // ClearFirst.
465
466
467 template<class ARFLOAT, class ARTYPE>
ClearBasis()468 inline void ARrcStdEig<ARFLOAT, ARTYPE>::ClearBasis()
469 {
470
471 BasisOK = ValuesOK = VectorsOK = SchurOK = false;
472
473 } // ClearBasis.
474
475
476 template<class ARFLOAT, class ARTYPE>
ClearMem()477 void ARrcStdEig<ARFLOAT, ARTYPE>::ClearMem()
478 {
479
480 // Deleting working arrays.
481
482 if (workl) delete[] workl;
483 if (workd) delete[] workd;
484 if (workv) delete[] workv;
485 if (rwork) delete[] rwork;
486 if (V) delete[] V;
487
488 workl = NULL;
489 workd = NULL;
490 workv = NULL;
491 rwork = NULL;
492 V = NULL;
493
494 // Deleting input and output arrays.
495
496 if (newRes) {
497 delete[] resid;
498 newRes = false;
499 resid = NULL; // Salwen. Mar 3, 2000.
500 }
501
502 if (newVal) {
503 delete[] EigValR;
504 delete[] EigValI;
505 newVal = false;
506 }
507 EigValR=NULL;
508 EigValI=NULL;
509
510 if (newVec) {
511 delete[] EigVec;
512 newVec = false;
513 }
514 EigVec=NULL;
515
516 // Adjusting boolean variables.
517
518 ClearFirst();
519
520 } // ClearMem.
521
522
523 template<class ARFLOAT, class ARTYPE>
ValAllocate()524 inline void ARrcStdEig<ARFLOAT, ARTYPE>::ValAllocate()
525 {
526
527 if (EigValR == NULL) { // Creating a new array EigValR.
528 EigValR = new ARTYPE[ValSize()];
529 newVal = true;
530 }
531
532 } // ValAllocate.
533
534 template<class ARFLOAT, class ARTYPE>
VecAllocate(bool newV)535 inline void ARrcStdEig<ARFLOAT, ARTYPE>::VecAllocate(bool newV)
536 {
537
538 if (EigVec == NULL) {
539 if (newV) { // Creating a new array EigVec.
540 EigVec = new ARTYPE[ValSize()*n];
541 newVec = true;
542 }
543 else { // Using V to store EigVec.
544 EigVec = &V[1];
545 }
546 }
547
548 } // VecAllocate.
549
550
551 template<class ARFLOAT, class ARTYPE>
WorkspaceAllocate()552 void ARrcStdEig<ARFLOAT, ARTYPE>::WorkspaceAllocate()
553 {
554
555 lworkl = 0;
556 lworkv = 0;
557 lrwork = 0;
558
559 } // WorkspaceAllocate.
560
561
562 template<class ARFLOAT, class ARTYPE>
AuppError()563 void ARrcStdEig<ARFLOAT, ARTYPE>::AuppError()
564 {
565
566 switch (info) {
567 case 0:
568 return;
569 case -8:
570 throw ArpackError(ArpackError::LAPACK_ERROR, "Aupp");
571 case -9:
572 throw ArpackError(ArpackError::START_RESID_ZERO, "Aupp");
573 case -9999:
574 throw ArpackError(ArpackError::ARNOLDI_NOT_BUILD, "Aupp");
575 case 1:
576 ArpackError(ArpackError::MAX_ITERATIONS, "Aupp");
577 return;
578 case 3:
579 ArpackError(ArpackError::NO_SHIFTS_APPLIED, "Aupp");
580 return;
581 default :
582 throw ArpackError(ArpackError::AUPP_ERROR, "Aupp");
583 }
584
585 } // AuppError.
586
587
588 template<class ARFLOAT, class ARTYPE>
EuppError()589 void ARrcStdEig<ARFLOAT, ARTYPE>::EuppError()
590 {
591
592 switch (info) {
593 case 0:
594 return;
595 case -8:
596 case -9:
597 throw ArpackError(ArpackError::LAPACK_ERROR, "Eupp");
598 case -14:
599 throw ArpackError(ArpackError::NOT_ACCURATE_EIG, "Eupp");
600 case 1:
601 throw ArpackError(ArpackError::REORDERING_ERROR, "Eupp");
602 default :
603 throw ArpackError(ArpackError::EUPP_ERROR, "Eupp");
604 }
605
606 } // EuppError.
607
608
609 template<class ARFLOAT, class ARTYPE>
CheckN(int np)610 inline int ARrcStdEig<ARFLOAT, ARTYPE>::CheckN(int np)
611 {
612
613 if (np < 2) {
614 throw ArpackError(ArpackError::N_SMALLER_THAN_2);
615 }
616 return np;
617
618 } // CheckN.
619
620
621 template<class ARFLOAT, class ARTYPE>
CheckNcv(int ncvp)622 inline int ARrcStdEig<ARFLOAT, ARTYPE>::CheckNcv(int ncvp)
623 {
624
625 // Adjusting ncv if ncv <= nev or ncv > n.
626
627 if (ncvp < nev+1) {
628 if (ncvp) ArpackError::Set(ArpackError::NCV_OUT_OF_BOUNDS);
629 return ((2*nev+1)>n)?n:(2*nev+1);
630 }
631 else if (ncvp > n) {
632 ArpackError::Set(ArpackError::NCV_OUT_OF_BOUNDS);
633 return n;
634 }
635 else {
636 return ncvp;
637 }
638
639 } // CheckNcv.
640
641
642 template<class ARFLOAT, class ARTYPE>
CheckNev(int nevp)643 inline int ARrcStdEig<ARFLOAT, ARTYPE>::CheckNev(int nevp)
644 {
645
646 if ((nevp<1)||(nevp>=n)) {
647 throw ArpackError(ArpackError::NEV_OUT_OF_BOUNDS);
648 }
649 return nevp;
650
651 } // CheckNev.
652
653
654 template<class ARFLOAT, class ARTYPE>
CheckMaxit(int maxitp)655 inline int ARrcStdEig<ARFLOAT, ARTYPE>::CheckMaxit(int maxitp)
656 {
657
658 if (maxitp >= 1) return maxitp;
659 if (maxitp < 0) ArpackError::Set(ArpackError::MAXIT_NON_POSITIVE);
660 return 100*nev;
661
662 } // CheckMaxit.
663
664 template<class ARFLOAT, class ARTYPE>
CheckWhich(const std::string & whichp)665 std::string ARrcStdEig<ARFLOAT, ARTYPE>::CheckWhich(const std::string& whichp)
666 {
667
668 switch (whichp[0]) { // The first ought to be S or L.
669 case 'S':
670 case 'L':
671 switch (whichp[1]) { // The second must be M, R or I.
672 case 'M':
673 case 'R':
674 case 'I':
675 return whichp;
676 }
677 default :
678 throw ArpackError(ArpackError::WHICH_UNDEFINED);
679 }
680
681 } // CheckWhich.
682
683
684 template<class ARFLOAT, class ARTYPE>
Restart()685 void ARrcStdEig<ARFLOAT, ARTYPE>::Restart()
686 {
687
688 nconv=0; // No eigenvalues found yet.
689 ido =0; // First call to AUPP.
690 iparam[1]=(int)AutoShift; // Shift strategy used.
691 iparam[3]=maxit; // Maximum number of Arnoldi iterations allowed.
692 iparam[4]=1; // Blocksize must be 1.
693 info =(int)(!newRes); // Starting vector used.
694 ClearBasis();
695
696 } // Restart.
697
698
699 template<class ARFLOAT, class ARTYPE>
Prepare()700 void ARrcStdEig<ARFLOAT, ARTYPE>::Prepare()
701 {
702
703 // Deleting old stuff.
704
705 ClearMem();
706
707 // Defining internal variables.
708
709 try {
710
711 if (resid == NULL) { // Using a random starting vector.
712 resid = new ARTYPE[n];
713 newRes = true;
714 }
715
716 // Setting dimensions of working arrays.
717
718 workd = new ARTYPE[3*n+1];
719 V = new ARTYPE[n*ncv+1];
720 WorkspaceAllocate();
721
722 }
723 catch (ArpackError) { // Returning from here if an error has occurred.
724 ArpackError(ArpackError::CANNOT_PREPARE, "Prepare");
725 return;
726 }
727
728 Restart();
729
730 // That's all.
731
732 PrepareOK = true;
733
734 } // Prepare.
735
736
737 template<class ARFLOAT, class ARTYPE>
Copy(const ARrcStdEig<ARFLOAT,ARTYPE> & other)738 void ARrcStdEig<ARFLOAT,ARTYPE>::Copy(const ARrcStdEig<ARFLOAT,ARTYPE>& other)
739 {
740
741 // Defining local variables.
742
743 int i;
744
745 // Copying variables that belong to fundamental types.
746
747 n = other.n;
748 nev = other.nev;
749 ncv = other.ncv;
750 maxit = other.maxit;
751 which = other.which;
752 tol = other.tol;
753 sigmaI = other.sigmaI;
754 sigmaR = other.sigmaR;
755 rvec = other.rvec;
756 newRes = other.newRes;
757 newVal = other.newVal;
758 newVec = other.newVec;
759 PrepareOK = other.PrepareOK;
760 BasisOK = other.BasisOK;
761 ValuesOK = other.ValuesOK;
762 VectorsOK = other.VectorsOK;
763 SchurOK = other.SchurOK;
764 AutoShift = other.AutoShift;
765 bmat = other.bmat;
766 HowMny = other.HowMny;
767 ido = other.ido;
768 info = other.info;
769 mode = other.mode;
770 nconv = other.nconv;
771
772 // Copying arrays with static dimension.
773
774 for (i=0; i<12; i++) iparam[i] = other.iparam[i];
775 for (i=0; i<15; i++) ipntr[i] = other.ipntr[i];
776
777 // Returning from here if "other" was not initialized.
778
779 if (!PrepareOK) return;
780
781 // Copying dynamic variables.
782
783 workd = new ARTYPE[3*n+1]; // workd.
784 copy(3*n+1,other.workd,1,workd,1);
785
786 V = new ARTYPE[n*ncv+1]; // V.
787 copy(n*ncv+1,other.V,1,V,1);
788
789 if (newRes) { // resid.
790 resid = new ARTYPE[n];
791 copy(n,other.resid,1,resid,1);
792 }
793 else {
794 resid = other.resid;
795 }
796
797 if (newVec) { // EigVec.
798 EigVec = new ARTYPE[ValSize()*n];
799 copy(ValSize()*n,other.EigVec,1,EigVec,1);
800 }
801 else if (other.EigVec == (&other.V[1])) {
802 EigVec = &V[1];
803 }
804 else {
805 EigVec = other.EigVec;
806 }
807
808 if (newVal) { // EigValR and EigValI.
809 EigValR = new ARTYPE[ValSize()];
810 copy(ValSize(),other.EigValR,1,EigValR,1);
811 if (other.EigValI != NULL) {
812 EigValI = new ARFLOAT[ValSize()];
813 copy(ValSize(),other.EigValI,1,EigValI,1);
814 }
815 else {
816 EigValI = NULL;
817 }
818 }
819 else {
820 EigValR = other.EigValR;
821 EigValI = other.EigValI;
822 }
823
824 WorkspaceAllocate(); // lworkl, workl, workv and rwork.
825 if (lworkl) copy(lworkl+1,other.workl,1,workl,1);
826 if (lworkv) copy(lworkv+1,other.workv,1,workv,1);
827 if (lrwork) copy(lrwork+1,other.rwork,1,rwork,1);
828
829 } // Copy.
830
831
832 template<class ARFLOAT, class ARTYPE>
833 void ARrcStdEig<ARFLOAT, ARTYPE>::
DefineParameters(int np,int nevp,const std::string & whichp,int ncvp,ARFLOAT tolp,int maxitp,ARTYPE * residp,bool ishiftp)834 DefineParameters(int np, int nevp, const std::string& whichp, int ncvp,
835 ARFLOAT tolp, int maxitp, ARTYPE* residp, bool ishiftp)
836
837 {
838
839 // Providing a "new" handler.
840
841 std::set_new_handler ( MemoryOverflow );
842
843 // Setting user defined parameters.
844
845 try {
846 n = CheckN(np);
847 nev = CheckNev(nevp);
848 ncv = CheckNcv(ncvp);
849 which = CheckWhich(whichp);
850 maxit = CheckMaxit(maxitp);
851 tol = tolp;
852 resid = residp;
853 AutoShift = ishiftp;
854 }
855
856 // Returning from here if an error has occurred.
857
858 catch (ArpackError) {
859 ArpackError(ArpackError::PARAMETER_ERROR, "DefineParameter");
860 return;
861 }
862
863 // Setting internal variables.
864
865 ClearFirst();
866 Prepare();
867
868 } // DefineParameters.
869
870
871 template<class ARFLOAT, class ARTYPE>
GetIter()872 int ARrcStdEig<ARFLOAT, ARTYPE>::GetIter()
873 {
874
875 if (BasisOK || ValuesOK || VectorsOK || SchurOK) {
876 return iparam[3];
877 }
878 else {
879 return 0;
880 }
881
882 } // GetIter.
883
884
885 template<class ARFLOAT, class ARTYPE>
ChangeMaxit(int maxitp)886 inline void ARrcStdEig<ARFLOAT, ARTYPE>::ChangeMaxit(int maxitp)
887 {
888
889 maxit = CheckMaxit(maxitp);
890 Restart();
891
892 } // ChangeMaxit.
893
894
895 template<class ARFLOAT, class ARTYPE>
ChangeTol(ARFLOAT tolp)896 inline void ARrcStdEig<ARFLOAT, ARTYPE>::ChangeTol(ARFLOAT tolp)
897 {
898
899 if (tolp < tol) Restart();
900 tol = tolp;
901
902 } // ChangeTol.
903
904
905 template<class ARFLOAT, class ARTYPE>
ChangeNev(int nevp)906 void ARrcStdEig<ARFLOAT, ARTYPE>::ChangeNev(int nevp)
907 {
908
909 try {
910 nev = CheckNev(nevp);
911 ncv = CheckNcv(ncv);
912 }
913 catch (ArpackError) { return; }
914 if (PrepareOK) Prepare();
915
916 } // ChangeNev.
917
918
919 template<class ARFLOAT, class ARTYPE>
ChangeNcv(int ncvp)920 inline void ARrcStdEig<ARFLOAT, ARTYPE>::ChangeNcv(int ncvp)
921 {
922
923 ncv = CheckNcv(ncvp);
924 if (PrepareOK) Prepare();
925
926 } // ChangeNcv.
927
928
929 template<class ARFLOAT, class ARTYPE>
ChangeWhich(const std::string & whichp)930 void ARrcStdEig<ARFLOAT, ARTYPE>::ChangeWhich(const std::string& whichp)
931 {
932
933 try { which = CheckWhich(whichp); }
934 catch (ArpackError) { return; }
935 Restart();
936
937 } // ChangeWhich.
938
939
940 template<class ARFLOAT, class ARTYPE>
ChangeShift(ARTYPE sigmaRp)941 inline void ARrcStdEig<ARFLOAT, ARTYPE>::ChangeShift(ARTYPE sigmaRp)
942 {
943
944 sigmaR = sigmaRp;
945 sigmaI = 0.0;
946 mode = 3;
947 iparam[7] = mode;
948 Restart();
949
950 } // ChangeShift.
951
952
953 template<class ARFLOAT, class ARTYPE>
NoShift()954 inline void ARrcStdEig<ARFLOAT, ARTYPE>::NoShift()
955 {
956
957 sigmaR = (ARTYPE)0;
958 sigmaI = 0.0;
959 mode = 1;
960 iparam[7] = mode;
961 Restart();
962
963 } // NoShift.
964
965
966 template<class ARFLOAT, class ARTYPE>
InvertAutoShift()967 void ARrcStdEig<ARFLOAT, ARTYPE>::InvertAutoShift()
968 {
969
970 AutoShift = !AutoShift;
971 Restart();
972
973 } // InvertAutoShift.
974
975
976 template<class ARFLOAT, class ARTYPE>
GetVector()977 ARTYPE* ARrcStdEig<ARFLOAT, ARTYPE>::GetVector()
978 {
979
980 switch (ido) {
981 case -1:
982 case 1:
983 case 2:
984 return &workd[ipntr[1]];
985 case 3:
986 return &workl[ipntr[6]];
987 default:
988 throw ArpackError(ArpackError::CANNOT_GET_VECTOR, "GetVector");
989 }
990
991 } // GetVector.
992
993
994 template<class ARFLOAT, class ARTYPE>
GetProd()995 ARTYPE* ARrcStdEig<ARFLOAT, ARTYPE>::GetProd()
996 {
997
998 if (ido != 1) {
999 throw ArpackError(ArpackError::CANNOT_GET_PROD, "GetProd");
1000 }
1001 return &workd[ipntr[3]];
1002
1003 } // GetProd.
1004
1005
1006 template<class ARFLOAT, class ARTYPE>
PutVector()1007 ARTYPE* ARrcStdEig<ARFLOAT, ARTYPE>::PutVector()
1008 {
1009
1010 switch (ido) {
1011 case -1:
1012 case 1:
1013 case 2:
1014 return &workd[ipntr[2]];
1015 case 3:
1016 return &workl[ipntr[14]];
1017 default:
1018 throw ArpackError(ArpackError::CANNOT_PUT_VECTOR, "PutVector");
1019 }
1020
1021 } // PutVector.
1022
1023
1024 template<class ARFLOAT, class ARTYPE>
TakeStep()1025 int ARrcStdEig<ARFLOAT, ARTYPE>::TakeStep()
1026 {
1027
1028 // Requiring the definition of all internal variables.
1029
1030 if (!PrepareOK) {
1031
1032 throw ArpackError(ArpackError::PREPARE_NOT_OK, "TakeStep");
1033
1034 }
1035 else if (!BasisOK) {
1036
1037 // Taking a step if the Arnoldi basis is not available.
1038
1039 Aupp();
1040
1041 // Checking if convergence was obtained.
1042
1043 if (ido==99) {
1044 nconv = iparam[5];
1045 AuppError();
1046 if (info >= 0) BasisOK = true;
1047 }
1048 }
1049
1050 return ido;
1051
1052 } // TakeStep.
1053
1054
1055 template<class ARFLOAT, class ARTYPE>
FindArnoldiBasis()1056 inline int ARrcStdEig<ARFLOAT, ARTYPE>::FindArnoldiBasis()
1057 {
1058
1059 if (!BasisOK) {
1060 throw ArpackError(ArpackError::CANNOT_FIND_BASIS, "FindArnoldiBasis");
1061 }
1062 return nconv;
1063
1064 } // FindArnoldiBasis.
1065
1066
1067 template<class ARFLOAT, class ARTYPE>
FindEigenvalues()1068 int ARrcStdEig<ARFLOAT, ARTYPE>::FindEigenvalues()
1069 {
1070
1071 // Determining eigenvalues if they are not available.
1072
1073 if (!ValuesOK) {
1074 try {
1075 ValAllocate();
1076 nconv = FindArnoldiBasis();
1077 rvec = false;
1078 if (nconv>0) {
1079 Eupp();
1080 EuppError();
1081 }
1082 }
1083 catch (ArpackError) {
1084 ArpackError(ArpackError::CANNOT_FIND_VALUES, "FindEigenvalues");
1085 return 0;
1086 }
1087 if (newVal) ValuesOK = true;
1088 }
1089 return nconv;
1090
1091 } // FindEigenvalues.
1092
1093
1094 template<class ARFLOAT, class ARTYPE>
FindEigenvectors(bool schurp)1095 int ARrcStdEig<ARFLOAT, ARTYPE>::FindEigenvectors(bool schurp)
1096 {
1097
1098 // Determining eigenvectors if they are not available.
1099
1100 if (!VectorsOK) {
1101 try {
1102 ValAllocate();
1103 VecAllocate(schurp);
1104 nconv = FindArnoldiBasis();
1105 rvec = true;
1106 HowMny = 'A';
1107 if (nconv>0) {
1108 Eupp();
1109 EuppError();
1110 }
1111 }
1112 catch (ArpackError) {
1113 ArpackError(ArpackError::CANNOT_FIND_VECTORS, "FindEigenvectors");
1114 return 0;
1115 }
1116 BasisOK = false;
1117 if (newVal) ValuesOK = true;
1118 if (newVec || OverV()) VectorsOK = true;
1119 if (!OverV()) SchurOK = true;
1120 }
1121 return nconv;
1122
1123 } // FindEigenvectors.
1124
1125
1126 template<class ARFLOAT, class ARTYPE>
FindSchurVectors()1127 int ARrcStdEig<ARFLOAT, ARTYPE>::FindSchurVectors()
1128 {
1129
1130 // Determining Schur vectors if they are not available.
1131
1132 if (!SchurOK) {
1133 try {
1134 ValAllocate();
1135 nconv = FindArnoldiBasis();
1136 rvec = true;
1137 HowMny = 'P';
1138 if (nconv>0) {
1139 Eupp();
1140 EuppError();
1141 }
1142 }
1143 catch (ArpackError) {
1144 ArpackError(ArpackError::CANNOT_FIND_SCHUR, "FindSchurVectors");
1145 return 0;
1146 }
1147 BasisOK = false;
1148 if (newVal) ValuesOK = true;
1149 SchurOK =true;
1150 }
1151 return nconv;
1152
1153 } // FindSchurVectors.
1154
1155
1156 template<class ARFLOAT, class ARTYPE>
1157 int ARrcStdEig<ARFLOAT, ARTYPE>::
Eigenvectors(ARTYPE * & EigVecp,bool ischur)1158 Eigenvectors(ARTYPE* &EigVecp, bool ischur)
1159 {
1160
1161 // Overriding EigVecp with the converged eigenvectors.
1162
1163 if (VectorsOK) { // Eigenvectors are available.
1164 if ((EigVecp == NULL) && (newVec)) { // Moving eigenvectors.
1165 EigVecp = EigVec;
1166 EigVec = NULL;
1167 newVec = false;
1168 VectorsOK = false;
1169 }
1170 else { // Copying eigenvectors.
1171 if (EigVecp == NULL) {
1172 try { EigVecp = new ARTYPE[ValSize()*n]; }
1173 catch (ArpackError) { return 0; }
1174 }
1175 copy(ValSize()*n,EigVec,1,EigVecp,1);
1176 }
1177 }
1178 else { // Eigenvectors are not available.
1179 if (newVec) {
1180 delete[] EigVec;
1181 newVec = false;
1182 }
1183 if (EigVecp == NULL) {
1184 try { EigVecp = new ARTYPE[ValSize()*n]; }
1185 catch (ArpackError) { return 0; }
1186 }
1187 EigVec = EigVecp;
1188 nconv = FindEigenvectors(ischur);
1189 EigVec = NULL;
1190 }
1191 return nconv;
1192
1193 } // Eigenvectors(EigVecp, ischur).
1194
1195
1196 template<class ARFLOAT, class ARTYPE>
ArnoldiBasisVector(int i,int j)1197 inline ARTYPE ARrcStdEig<ARFLOAT, ARTYPE>::ArnoldiBasisVector(int i, int j)
1198 {
1199
1200 // Returning element j of Arnoldi basis vector i.
1201
1202 if (!BasisOK) {
1203 throw ArpackError(ArpackError::BASIS_NOT_OK, "ArnoldiBasisVector(i,j)");
1204 }
1205 else if ((i>=ncv)||(i<0)||(j>=n)||(j<0)) {
1206 throw ArpackError(ArpackError::RANGE_ERROR,"ArnoldiBasisVector(i,j)");
1207 }
1208 return V[i*n+j+1];
1209
1210 } // ArnoldiBasisVector(i,j).
1211
1212
1213 template<class ARFLOAT, class ARTYPE>
SchurVector(int i,int j)1214 inline ARTYPE ARrcStdEig<ARFLOAT, ARTYPE>::SchurVector(int i, int j)
1215 {
1216
1217 // Returning element j of Schur vector i.
1218
1219 if (!SchurOK) {
1220 throw ArpackError(ArpackError::SCHUR_NOT_OK, "SchurVector(i,j)");
1221 }
1222 else if ((i>=nconv)||(i<0)||(j>=n)||(j<0)) {
1223 throw ArpackError(ArpackError::RANGE_ERROR, "SchurVector(i,j)");
1224 }
1225 return V[i*n+j+1];
1226
1227 } // SchurVector(i,j).
1228
1229
1230 template<class ARFLOAT, class ARTYPE>
ResidualVector(int i)1231 inline ARTYPE ARrcStdEig<ARFLOAT, ARTYPE>::ResidualVector(int i)
1232 {
1233
1234 // Returning element i of the residual vector.
1235
1236 if ((!newRes)||(!(BasisOK||ValuesOK||VectorsOK||SchurOK))) {
1237 throw ArpackError(ArpackError::RESID_NOT_OK, "ResidualVector(i)");
1238 }
1239 else if ((i>=n)||(i<0)) {
1240 throw ArpackError(ArpackError::RANGE_ERROR, "ResidualVector(i)");
1241 }
1242 return resid[i];
1243
1244 } // ResidualVector(i).
1245
1246
1247 template<class ARFLOAT, class ARTYPE>
RawArnoldiBasisVectors()1248 inline ARTYPE* ARrcStdEig<ARFLOAT, ARTYPE>::RawArnoldiBasisVectors()
1249 {
1250
1251 // Returning a constant pointer to Arnoldi basis.
1252
1253 if (!BasisOK) {
1254 throw ArpackError(ArpackError::BASIS_NOT_OK, "RawArnoldiBasisVectors");
1255 }
1256 return &V[1];
1257
1258 } // RawArnoldiBasisVectors.
1259
1260
1261 template<class ARFLOAT, class ARTYPE>
RawArnoldiBasisVector(int i)1262 inline ARTYPE* ARrcStdEig<ARFLOAT, ARTYPE>::RawArnoldiBasisVector(int i)
1263 {
1264
1265 // Returning a constant pointer to Arnoldi basis vector i.
1266
1267 if (!BasisOK) {
1268 throw ArpackError(ArpackError::BASIS_NOT_OK, "RawArnoldiBasisVector(i)");
1269 }
1270 else if ((i>=ncv)||(i<0)) {
1271 throw ArpackError(ArpackError::RANGE_ERROR,"RawArnoldiBasisVector(i)");
1272 }
1273 return &V[i*n+1];
1274
1275 } // RawArnoldiBasisVector(i).
1276
1277
1278 template<class ARFLOAT, class ARTYPE>
RawEigenvalues()1279 inline ARTYPE* ARrcStdEig<ARFLOAT, ARTYPE>::RawEigenvalues()
1280 {
1281
1282 if (!ValuesOK) {
1283 throw ArpackError(ArpackError::VALUES_NOT_OK, "RawEigenvalues");
1284 }
1285 return EigValR;
1286
1287 } // RawEigenvalues.
1288
1289
1290 template<class ARFLOAT, class ARTYPE>
RawEigenvectors()1291 inline ARTYPE* ARrcStdEig<ARFLOAT, ARTYPE>::RawEigenvectors()
1292 {
1293
1294 if (!VectorsOK) {
1295 throw ArpackError(ArpackError::VECTORS_NOT_OK, "RawEigenvectors");
1296 }
1297 return EigVec;
1298
1299 } // RawEigenvectors.
1300
1301
1302 template<class ARFLOAT, class ARTYPE>
RawEigenvector(int i)1303 inline ARTYPE* ARrcStdEig<ARFLOAT, ARTYPE>::RawEigenvector(int i)
1304 {
1305
1306 // Returning a constant pointer to eigenvector i.
1307
1308 if (!VectorsOK) {
1309 throw ArpackError(ArpackError::VECTORS_NOT_OK, "RawEigenvector(i)");
1310 }
1311 else if ((i>=ValSize())||(i<0)) {
1312 throw ArpackError(ArpackError::RANGE_ERROR, "RawEigenvector(i)");
1313 }
1314 return &EigVec[i*n];
1315
1316 } // RawEigenvector(i).
1317
1318
1319 template<class ARFLOAT, class ARTYPE>
RawSchurVectors()1320 inline ARTYPE* ARrcStdEig<ARFLOAT, ARTYPE>::RawSchurVectors()
1321 {
1322
1323 if (!SchurOK) {
1324 throw ArpackError(ArpackError::SCHUR_NOT_OK, "RawSchurVectors");
1325 }
1326 return &V[1];
1327
1328 } // RawSchurVectors.
1329
1330
1331 template<class ARFLOAT, class ARTYPE>
RawSchurVector(int i)1332 inline ARTYPE* ARrcStdEig<ARFLOAT, ARTYPE>::RawSchurVector(int i)
1333 {
1334
1335 // Returning a constant pointer to Schur vector i.
1336
1337 if (!SchurOK) {
1338 throw ArpackError(ArpackError::SCHUR_NOT_OK, "RawSchurVector(i)");
1339 }
1340 else if ((i>=nev)||(i<0)) {
1341 throw ArpackError(ArpackError::RANGE_ERROR, "RawSchurVector(i)");
1342 }
1343 return &V[i*n+1];
1344
1345 } // RawSchurVector(i).
1346
1347
1348 template<class ARFLOAT, class ARTYPE>
RawResidualVector()1349 inline ARTYPE* ARrcStdEig<ARFLOAT, ARTYPE>::RawResidualVector()
1350 {
1351
1352 if (!newRes) {
1353 throw ArpackError(ArpackError::RESID_NOT_OK, "RawResidualVector");
1354 }
1355 return resid;
1356
1357 } // RawResidualVector.
1358
1359
1360 #ifdef STL_VECTOR_H // Defining some functions that use STL vector class.
1361
1362 template<class ARFLOAT, class ARTYPE>
StlArnoldiBasisVectors()1363 inline vector<ARTYPE>* ARrcStdEig<ARFLOAT, ARTYPE>::StlArnoldiBasisVectors()
1364 {
1365
1366 // Returning the Arnoldi basis in a single STL vector.
1367
1368 vector<ARTYPE>* StlBasis;
1369
1370 if (!BasisOK) {
1371 nconv = FindArnoldiBasis();
1372 }
1373 try {
1374 StlBasis = new vector<ARTYPE>(&V[1], &V[n*ncv+1]);
1375 }
1376 catch (ArpackError) { return NULL; }
1377 return StlBasis;
1378
1379 } // StlArnoldiBasisVectors.
1380
1381
1382 template<class ARFLOAT, class ARTYPE>
StlArnoldiBasisVector(int i)1383 inline vector<ARTYPE>* ARrcStdEig<ARFLOAT,ARTYPE>::StlArnoldiBasisVector(int i)
1384 {
1385
1386 // Returning the i-th Arnoldi basis vector in a STL vector.
1387
1388 vector<ARTYPE>* StlBasis;
1389
1390 if (!BasisOK) {
1391 throw ArpackError(ArpackError::BASIS_NOT_OK, "StlArnoldiBasisVector(i)");
1392 }
1393 else if ((i>=ncv)||(i<0)) {
1394 throw ArpackError(ArpackError::RANGE_ERROR,"StlArnoldiBasisVector(i)");
1395 }
1396 try {
1397 StlBasis = new vector<ARTYPE>(&V[i*n+1], &V[(i+1)*n+1]);
1398 }
1399 catch (ArpackError) { return NULL; }
1400 return StlBasis;
1401
1402 } // StlArnoldiBasisVector(i).
1403
1404
1405 template<class ARFLOAT, class ARTYPE>
StlEigenvectors(bool ischur)1406 inline vector<ARTYPE>* ARrcStdEig<ARFLOAT,ARTYPE>::StlEigenvectors(bool ischur)
1407 {
1408
1409 // Returning the eigenvector in a single STL vector.
1410
1411 vector<ARTYPE>* StlEigVec;
1412 ARTYPE* VecPtr;
1413
1414 try { StlEigVec = new vector<ARTYPE>(ValSize()*n); }
1415 catch (ArpackError) { return NULL; }
1416 VecPtr = StlEigVec->begin();
1417 nconv = Eigenvectors(VecPtr, ischur);
1418 return StlEigVec;
1419
1420 } // StlEigenvectors.
1421
1422
1423 template<class ARFLOAT, class ARTYPE>
StlSchurVectors()1424 inline vector<ARTYPE>* ARrcStdEig<ARFLOAT, ARTYPE>::StlSchurVectors()
1425 {
1426
1427 vector<ARTYPE>* StlSchurVec;
1428
1429 if (!SchurOK) {
1430 nconv = FindSchurVectors();
1431 }
1432 try {
1433 StlSchurVec = new vector<ARTYPE>(&V[1], &V[nev*n+1]);
1434 }
1435 catch (ArpackError) { return NULL; }
1436 return StlSchurVec;
1437
1438 } // StlSchurVectors.
1439
1440
1441 template<class ARFLOAT, class ARTYPE>
StlSchurVector(int i)1442 inline vector<ARTYPE>* ARrcStdEig<ARFLOAT, ARTYPE>::StlSchurVector(int i)
1443 {
1444
1445 // Returning the i-th Schur vector in a STL vector.
1446
1447 vector<ARTYPE>* StlSchurVec;
1448
1449 if (!SchurOK) {
1450 throw ArpackError(ArpackError::SCHUR_NOT_OK, "StlSchurVector(i)");
1451 }
1452 else if ((i>=nev)||(i<0)) {
1453 throw ArpackError(ArpackError::RANGE_ERROR, "StlSchurVector(i)");
1454 }
1455 try {
1456 StlSchurVec = new vector<ARTYPE>(&V[i*n+1], &V[(i+1)*n+1]);
1457 }
1458 catch (ArpackError) { return NULL; }
1459 return StlSchurVec;
1460
1461 } // StlSchurVector(i).
1462
1463
1464 template<class ARFLOAT, class ARTYPE>
StlResidualVector()1465 inline vector<ARTYPE>* ARrcStdEig<ARFLOAT, ARTYPE>::StlResidualVector()
1466 {
1467
1468 // Returning the residual vector in a STL vector.
1469
1470 vector<ARTYPE>* StlResid;
1471
1472 if (!newRes) {
1473 throw ArpackError(ArpackError::RESID_NOT_OK, "StlResidualVector");
1474 }
1475 try {
1476 StlResid = new vector<ARTYPE>(resid, &resid[n]);
1477 }
1478 catch (ArpackError) { return NULL; }
1479 return StlResid;
1480
1481 } // StlResidualVector.
1482
1483 #endif // #ifdef STL_VECTOR_H.
1484
1485
1486 template<class ARFLOAT, class ARTYPE>
ARrcStdEig()1487 inline ARrcStdEig<ARFLOAT, ARTYPE>::ARrcStdEig()
1488 {
1489
1490 resid = NULL;
1491 rwork = NULL;
1492 workl = NULL;
1493 workd = NULL;
1494 workv = NULL;
1495 V = NULL;
1496 EigValR = NULL;
1497 EigValI = NULL;
1498 EigVec = NULL;
1499 bmat = 'I'; // This is a standard problem.
1500 ClearFirst();
1501 NoShift();
1502 NoTrace();
1503
1504 } // Short constructor.
1505
1506
1507 template<class ARFLOAT, class ARTYPE>
1508 ARrcStdEig<ARFLOAT, ARTYPE>& ARrcStdEig<ARFLOAT, ARTYPE>::
1509 operator=(const ARrcStdEig<ARFLOAT, ARTYPE>& other)
1510 {
1511
1512 if (this != &other) { // Stroustrup suggestion.
1513 ClearMem();
1514 Copy(other);
1515 }
1516 return *this;
1517
1518 } // operator=.
1519
1520
1521 #endif // ARRSEIG_H
1522
1523