GNU Octave 11.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 
Loading...
Searching...
No Matches
eigs-base.cc
Go to the documentation of this file.
1////////////////////////////////////////////////////////////////////////
2//
3// Copyright (C) 2005-2026 The Octave Project Developers
4//
5// See the file COPYRIGHT.md in the top-level directory of this
6// distribution or <https://octave.org/copyright/>.
7//
8// This file is part of Octave.
9//
10// Octave is free software: you can redistribute it and/or modify it
11// under the terms of the GNU General Public License as published by
12// the Free Software Foundation, either version 3 of the License, or
13// (at your option) any later version.
14//
15// Octave is distributed in the hope that it will be useful, but
16// WITHOUT ANY WARRANTY; without even the implied warranty of
17// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18// GNU General Public License for more details.
19//
20// You should have received a copy of the GNU General Public License
21// along with Octave; see the file COPYING. If not, see
22// <https://www.gnu.org/licenses/>.
23//
24////////////////////////////////////////////////////////////////////////
25
26#if defined (HAVE_CONFIG_H)
27# include "config.h"
28#endif
29
30#include <cmath>
31
32#include <ostream>
33
34#include "Array-oct.h"
35#include "CSparse.h"
36#include "MatrixType.h"
37#include "PermMatrix.h"
38#include "arpack-proto.h"
39#include "blas-proto.h"
40#include "chol.h"
41#include "dSparse.h"
42#include "eigs-base.h"
43#include "lo-ieee.h"
44#include "lo-utils.h"
45#include "lu.h"
46#include "mx-ops.h"
47#include "oct-error.h"
48#include "oct-locbuf.h"
49#include "oct-rand.h"
50#include "sparse-chol.h"
51#include "sparse-lu.h"
52
53#if defined (HAVE_ARPACK)
54
55static void
56warn_convergence ()
57{
58 (*current_liboctave_warning_with_id_handler)
59 ("Octave:convergence",
60 "eigs: 'A - sigma*B' is singular, indicating sigma is exactly "
61 "an eigenvalue so convergence is not guaranteed");
62}
63
64// Conversion from error number to strings
65std::string
66arpack_errno2str (const octave_idx_type& errnum, const std::string& fcn_name)
67{
68 std::string msg;
69 std::string bug_msg = "\nThis should not happen. Please, see https://www.gnu.org/software/octave/bugs.html, and file a bug report";
70
71 switch (errnum)
72 {
73 case -1:
74 msg = "N must be positive";
75 break;
76
77 case -2:
78 msg = "NEV must be positive";
79 break;
80
81 case -3:
82 msg = "NCV-NEV >= 2 and less than or equal to N";
83 break;
84
85 case -4:
86 msg = "The maximum number of Arnoldi update iterations allowed must be greater than zero";
87 break;
88
89 case -5:
90 msg = "WHICH must be one of 'LM', 'SM', 'LR', 'SR', 'LI', 'SI'";
91 break;
92
93 case -6:
94 msg = "BMAT must be one of 'I' or 'G'";
95 break;
96
97 case -7:
98 msg = "Length of private work WORKL array is insufficient";
99 break;
100
101 case -8:
102 msg = "Error return from LAPACK eigenvalue calculation";
103 break;
104
105 case -9:
106 if (fcn_name.compare ("zneupd") == 0)
107 msg = "Error return from calculation of eigenvectors. Informational error from LAPACK routine ztrevc";
108 else if (fcn_name.compare ("dneupd") == 0)
109 msg = "Error return from calculation of eigenvectors. Informational error from LAPACK routine dtrevc";
110 else
111 msg = "Starting vector is zero";
112
113 break;
114
115 case -10:
116 if (fcn_name.compare ("dneupd") == 0
117 || fcn_name.compare ("dnaupd") == 0)
118 msg = "IPARAM(7) must be 1,2,3,4";
119 else if (fcn_name.compare ("zneupd") == 0
120 || fcn_name.compare ("znaupd") == 0)
121 msg = "IPARAM(7) must be 1,2,3";
122 else
123 msg = "IPARAM(7) must be 1,2,3,4,5";
124
125 break;
126
127 case -11:
128 msg = "IPARAM(7) = 1 and BMAT = 'G' are incompatible";
129 break;
130
131 case -12:
132 if (fcn_name.compare ("dnaupd") == 0
133 || fcn_name.compare ("znaupd") == 0
134 || fcn_name.compare ("dsaupd") == 0)
135 msg = std::string ("IPARAM(1) must be equal to 0 or 1");
136 else if (fcn_name.compare ("dneupd") == 0
137 || fcn_name.compare ("zneupd") == 0)
138 msg = "HOWMNY = 'S' not yet implemented";
139 else
140 msg = "NEV and WHICH = 'BE' are incompatible";
141
142 break;
143
144 case -13:
145 if (fcn_name.compare ("dneupd") == 0
146 || fcn_name.compare ("zneupd") == 0)
147 msg = "HOWMNY must be one of 'A' or 'P' if RVEC = .true.";
148 else if (fcn_name.compare ("dsaupd") == 0)
149 msg = "NEV and WHICH = 'BE' are incompatible";
150
151 break;
152
153 case -14:
154 if (fcn_name.compare ("dneupd") == 0)
155 msg = "DNAUPD did not find any eigenvalues to sufficient accuracy.";
156 else if (fcn_name.compare ("zneupd") == 0)
157 msg = "ZNAUPD did not find any eigenvalues to sufficient accuracy.";
158 else if (fcn_name.compare ("dseupd") == 0)
159 msg = "DSAUPD did not find any eigenvalues to sufficient accuracy.";
160
161 msg += " Consider changing tolerance (TOL), maximum iterations (MAXIT), number of Lanzcos basis vectors (P), or starting vector (V0) in OPTS structure.";
162
163 break;
164
165 case -15:
166 if (fcn_name.compare ("dseupd") == 0)
167 msg = "HOWMNY must be one of 'A' or 'S' if RVEC = .true.";
168
169 break;
170
171 case -16:
172 if (fcn_name.compare ("dseupd") == 0)
173 msg = "HOWMNY = 'S' not yet implemented";
174
175 break;
176
177 case -9999:
178 if (fcn_name.compare ("dnaupd") == 0)
179 msg = "Could not build an Arnoldi factorization. IPARAM(5) returns the size of the current Arnoldi factorization";
180
181 break;
182
183 case 1:
184 if (fcn_name.compare ("dneupd") == 0)
185 msg = "The Schur form computed by LAPACK routine dlahqr could not be reordered by LAPACK routine dtrsen. Re-enter subroutine DNEUPD with IPARAM(5)=NCV and increase the size of the arrays DR and DI to have dimension at least dimension NCV and allocate at least NCV columns for Z. NOTE: Not necessary if Z and V share the same space. Please notify the authors if this error occurs.";
186 else if (fcn_name.compare ("dnaupd") == 0
187 || fcn_name.compare ("znaupd") == 0
188 || fcn_name.compare ("dsaupd") == 0)
189 msg = "Maximum number of iterations taken. All possible eigenvalues of OP has been found. IPARAM(5) returns the number of wanted converged Ritz values";
190 else if (fcn_name.compare ("znaupd") == 0)
191 msg = "The Schur form computed by LAPACK routine csheqr could not be reordered by LAPACK routine ztrsen. Re-enter subroutine ZNEUPD with IPARAM(5)=NCV and increase the size of the array D to have dimension at least dimension NCV and allocate at least NCV columns for Z. NOTE: Not necessary if Z and V share the same space. Please notify the authors if this error occurs.";
192
193 break;
194
195 case 2:
196 if (fcn_name.compare ("dnaupd") == 0
197 || fcn_name.compare ("znaupd") == 0
198 || fcn_name.compare ("dsaupd") == 0)
199 msg = "No longer an informational error. Deprecated starting with release 2 of ARPACK.";
200
201 break;
202
203 case 3:
204 if (fcn_name.compare ("dnaupd") == 0
205 || fcn_name.compare ("znaupd") == 0
206 || fcn_name.compare ("dsaupd") == 0)
207 msg = "No shifts could be applied during a cycle of the implicitly restarted Arnoldi iteration. One possibility is to increase the size of NCV relative to NEV.";
208
209 break;
210
211 }
212
213 if ((errnum != -9) && (errnum != -14) && (errnum != -9999))
214 msg.append (bug_msg); // This is a bug in Octave interface to ARPACK
215
216 return msg;
217}
218
219template <typename M, typename SM>
220static octave_idx_type
221lusolve (const SM& L, const SM& U, M& m)
222{
223 octave_idx_type err = 0;
224 double rcond;
225 MatrixType utyp (MatrixType::Upper);
226
227 // Sparse L is lower triangular, Dense L is permuted lower triangular!!!
228 MatrixType ltyp (MatrixType::Lower);
229 m = L.solve (ltyp, m, err, rcond, nullptr);
230 if (err)
231 return err;
232
233 m = U.solve (utyp, m, err, rcond, nullptr);
234
235 return err;
236}
237
238template <typename SM, typename M>
239static M
240ltsolve (const SM& L, const ColumnVector& Q, const M& m)
241{
242 // Solve (Q_mat * L) * x = m, that is L * x = Q_mat' * m = m(Q)
243 octave_idx_type n = L.cols ();
244 octave_idx_type b_nc = m.cols ();
245 octave_idx_type err = 0;
246 double rcond;
247 MatrixType ltyp (MatrixType::Lower);
248 M retval (n, b_nc);
249 const double *qv = Q.data ();
250 for (octave_idx_type j = 0; j < b_nc; j++)
251 {
252 for (octave_idx_type i = 0; i < n; i++)
253 retval.elem (i, j) = m.elem (static_cast<octave_idx_type> (qv[i]), j);
254 }
255 return L.solve (ltyp, retval, err, rcond, nullptr);
256}
257
258template <typename SM, typename M>
259static M
260utsolve (const SM& U, const ColumnVector& Q, const M& m)
261{
262 // Solve (U * Q_mat') * x = m by U * tmp = m, x(Q) = tmp (Q_mat * tmp = x)
263 octave_idx_type n = U.cols ();
264 octave_idx_type b_nc = m.cols ();
265 octave_idx_type err = 0;
266 double rcond;
267 MatrixType utyp (MatrixType::Upper);
268 M tmp = U.solve (utyp, m, err, rcond, nullptr);
269 M retval;
270 const double *qv = Q.data ();
271
272 if (! err)
273 {
274 retval.resize (n, b_nc);
275 for (octave_idx_type j = 0; j < b_nc; j++)
276 {
277 for (octave_idx_type i = 0; i < n; i++)
278 retval.elem (static_cast<octave_idx_type> (qv[i]), j)
279 = tmp.elem (i, j);
280 }
281 }
282
283 return retval;
284}
285
286static bool
287vector_product (const SparseMatrix& m, const double *x, double *y)
288{
289 octave_idx_type nc = m.cols ();
290
291 for (octave_idx_type j = 0; j < nc; j++)
292 y[j] = 0.;
293
294 for (octave_idx_type j = 0; j < nc; j++)
295 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
296 y[m.ridx (i)] += m.data (i) * x[j];
297
298 return true;
299}
300
301static bool
302vector_product (const Matrix& m, const double *x, double *y)
303{
304 F77_INT nr = octave::to_f77_int (m.rows ());
305 F77_INT nc = octave::to_f77_int (m.cols ());
306
307 F77_XFCN (dgemv, DGEMV, (F77_CONST_CHAR_ARG2 ("N", 1),
308 nr, nc, 1.0, m.data (), nr,
309 x, 1, 0.0, y, 1
310 F77_CHAR_ARG_LEN (1)));
311
312 return true;
313}
314
315static bool
316vector_product (const SparseComplexMatrix& m, const Complex *x,
317 Complex *y)
318{
319 octave_idx_type nc = m.cols ();
320
321 for (octave_idx_type j = 0; j < nc; j++)
322 y[j] = 0.;
323
324 for (octave_idx_type j = 0; j < nc; j++)
325 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
326 y[m.ridx (i)] += m.data (i) * x[j];
327
328 return true;
329}
330
331static bool
332vector_product (const ComplexMatrix& m, const Complex *x, Complex *y)
333{
334 F77_INT nr = octave::to_f77_int (m.rows ());
335 F77_INT nc = octave::to_f77_int (m.cols ());
336
337 F77_XFCN (zgemv, ZGEMV, (F77_CONST_CHAR_ARG2 ("N", 1),
338 nr, nc, 1.0, F77_CONST_DBLE_CMPLX_ARG (m.data ()),
339 nr,
340 F77_CONST_DBLE_CMPLX_ARG (x), 1, 0.0,
341 F77_DBLE_CMPLX_ARG (y), 1
342 F77_CHAR_ARG_LEN (1)));
343
344 return true;
345}
346
347static bool
348make_cholb (Matrix& b, Matrix& bt, ColumnVector& permB)
349{
350 octave_idx_type info;
351 octave::math::chol<Matrix> fact (b, info);
352 octave_idx_type n = b.cols ();
353
354 if (info != 0)
355 return false;
356 else
357 {
358 bt = fact.chol_matrix (); // upper triangular
359 b = bt.transpose ();
360 permB = ColumnVector (n);
361 for (octave_idx_type i = 0; i < n; i++)
362 permB(i) = i;
363 return true;
364 }
365}
366
367static bool
368make_cholb (SparseMatrix& b, SparseMatrix& bt, ColumnVector& permB)
369{
370 octave_idx_type info;
371 octave::math::sparse_chol<SparseMatrix> fact (b, info, false);
372
373 if (info != 0)
374 return false;
375 else
376 {
377 b = fact.L (); // lower triangular
378 bt = b.transpose ();
379 permB = fact.perm () - 1.0;
380 return true;
381 }
382}
383
384static bool
385make_cholb (ComplexMatrix& b, ComplexMatrix& bt, ColumnVector& permB)
386{
387 octave_idx_type info;
388 octave::math::chol<ComplexMatrix> fact (b, info);
389 octave_idx_type n = b.cols ();
390
391 if (info != 0)
392 return false;
393 else
394 {
395 bt = fact.chol_matrix (); // upper triangular
396 b = bt.hermitian ();
397 permB = ColumnVector (n);
398 for (octave_idx_type i = 0; i < n; i++)
399 permB(i) = i;
400 return true;
401 }
402}
403
404static bool
405make_cholb (SparseComplexMatrix& b, SparseComplexMatrix& bt,
406 ColumnVector& permB)
407{
408 octave_idx_type info;
409 octave::math::sparse_chol<SparseComplexMatrix> fact (b, info, false);
410
411 if (info != 0)
412 return false;
413 else
414 {
415 b = fact.L (); // lower triangular
416 bt = b.hermitian ();
417 permB = fact.perm () - 1.0;
418 return true;
419 }
420}
421
422static bool
423LuAminusSigmaB (const SparseMatrix& m, const SparseMatrix& b,
424 bool cholB, const ColumnVector& permB, double sigma,
425 SparseMatrix& L, SparseMatrix& U, octave_idx_type *P,
426 octave_idx_type *Q, ColumnVector& r)
427{
428 bool have_b = ! b.isempty ();
429 octave_idx_type n = m.rows ();
430
431 // Calculate LU decomposition of 'M = A - sigma * B'
432 // P * (R \ M) * Q = L * U
433 SparseMatrix AminusSigmaB (m);
434
435 if (sigma != 0.0)
436 {
437 if (have_b)
438 {
439 if (cholB)
440 {
441 if (permB.numel ())
442 {
443 SparseMatrix tmp (n, n, n);
444 for (octave_idx_type i = 0; i < n; i++)
445 {
446 tmp.xcidx (i) = i;
447 tmp.xridx (i) = static_cast<octave_idx_type> (permB(i));
448 tmp.xdata (i) = 1;
449 }
450 tmp.xcidx (n) = n;
451
452 AminusSigmaB -= sigma * tmp *
453 b.transpose () * b * tmp.transpose ();
454 }
455 else
456 AminusSigmaB -= sigma * b.transpose () * b;
457 }
458 else
459 AminusSigmaB -= sigma * b;
460 }
461 else
462 {
463 SparseMatrix sigmat (n, n, n);
464
465 // Create sigma * speye (n,n)
466 sigmat.xcidx (0) = 0;
467 for (octave_idx_type i = 0; i < n; i++)
468 {
469 sigmat.xdata (i) = sigma;
470 sigmat.xridx (i) = i;
471 sigmat.xcidx (i+1) = i + 1;
472 }
473
474 AminusSigmaB -= sigmat;
475 }
476 }
477
478 octave::math::sparse_lu<SparseMatrix> fact (AminusSigmaB, Matrix (), true);
479
480 L = fact.L ();
481 U = fact.U ();
482 SparseMatrix R = fact.R ();
483 for (octave_idx_type i = 0; i < n; i++)
484 r(i) = R.xdata(i);
485
486 const octave_idx_type *P2 = fact.row_perm ();
487 const octave_idx_type *Q2 = fact.col_perm ();
488
489 for (octave_idx_type j = 0; j < n; j++)
490 {
491 P[j] = P2[j];
492 Q[j] = Q2[j];
493 }
494
495 // Test condition number of LU decomposition
496 double minU = octave::numeric_limits<double>::NaN ();
497 double maxU = octave::numeric_limits<double>::NaN ();
498 for (octave_idx_type j = 0; j < n; j++)
499 {
500 double d = 0.;
501 if (U.xcidx (j+1) > U.xcidx (j)
502 && U.xridx (U.xcidx (j+1)-1) == j)
503 d = std::abs (U.xdata (U.xcidx (j+1)-1));
504
505 if (octave::math::isnan (minU) || d < minU)
506 minU = d;
507
508 if (octave::math::isnan (maxU) || d > maxU)
509 maxU = d;
510 }
511
512 double rcond = (minU / maxU);
513 // Prevent use of extra precision.
514 double rcond_plus_one = rcond + 1.0;
515
516 if (rcond_plus_one == 1.0 || octave::math::isnan (rcond))
517 warn_convergence ();
518
519 return true;
520}
521
522static bool
523LuAminusSigmaB (const Matrix& m, const Matrix& b,
524 bool cholB, const ColumnVector& permB, double sigma,
525 Matrix& L, Matrix& U, octave_idx_type *P, octave_idx_type *Q,
526 ColumnVector& r)
527{
528 bool have_b = ! b.isempty ();
529 octave_idx_type n = m.cols ();
530
531 // Calculate LU decomposition of 'M = A - sigma * B'
532 // P * M = L * U
533 Matrix AminusSigmaB (m);
534
535 if (sigma != 0.0)
536 {
537 if (have_b)
538 {
539 if (cholB)
540 {
541 Matrix tmp = sigma * b.transpose () * b;
542 const double *pB = permB.data ();
543 double *p = AminusSigmaB.rwdata ();
544
545 if (permB.numel ())
546 {
547 for (octave_idx_type j = 0;
548 j < b.cols (); j++)
549 for (octave_idx_type i = 0;
550 i < b.rows (); i++)
551 *p++ -= tmp.xelem (static_cast<octave_idx_type> (pB[i]),
552 static_cast<octave_idx_type> (pB[j]));
553 }
554 else
555 AminusSigmaB -= tmp;
556 }
557 else
558 AminusSigmaB -= sigma * b;
559 }
560 else
561 {
562 double *p = AminusSigmaB.rwdata ();
563
564 for (octave_idx_type i = 0; i < n; i++)
565 p[i*(n+1)] -= sigma;
566 }
567 }
568
569 octave::math::lu<Matrix> fact (AminusSigmaB);
570
571 L = fact.L ();
572 U = fact.U ();
573 ColumnVector P2 = fact.P_vec();
574
575 for (octave_idx_type j = 0; j < n; j++)
576 {
577 Q[j] = j;
578 P[j] = P2(j) - 1;
579 r(j) = 1.;
580 }
581
582 // Test condition number of LU decomposition
583 double minU = octave::numeric_limits<double>::NaN ();
584 double maxU = octave::numeric_limits<double>::NaN ();
585 for (octave_idx_type j = 0; j < n; j++)
586 {
587 double d = std::abs (U.xelem (j, j));
588 if (octave::math::isnan (minU) || d < minU)
589 minU = d;
590
591 if (octave::math::isnan (maxU) || d > maxU)
592 maxU = d;
593 }
594
595 double rcond = (minU / maxU);
596 // Prevent use of extra precision.
597 double rcond_plus_one = rcond + 1.0;
598
599 if (rcond_plus_one == 1.0 || octave::math::isnan (rcond))
600 warn_convergence ();
601
602 return true;
603}
604
605static bool
606LuAminusSigmaB (const SparseComplexMatrix& m, const SparseComplexMatrix& b,
607 bool cholB, const ColumnVector& permB, Complex sigma,
608 SparseComplexMatrix& L, SparseComplexMatrix& U,
609 octave_idx_type *P, octave_idx_type *Q, ColumnVector& r)
610{
611 bool have_b = ! b.isempty ();
612 octave_idx_type n = m.rows ();
613
614 // Calculate LU decomposition of 'M = A - sigma * B'
615 // P * (R \ M) * Q = L * U
616 SparseComplexMatrix AminusSigmaB (m);
617
618 if (std::real (sigma) != 0.0 || std::imag (sigma) != 0.0)
619 {
620 if (have_b)
621 {
622 if (cholB)
623 {
624 if (permB.numel ())
625 {
626 SparseMatrix tmp (n, n, n);
627 for (octave_idx_type i = 0; i < n; i++)
628 {
629 tmp.xcidx (i) = i;
630 tmp.xridx (i) = static_cast<octave_idx_type> (permB(i));
631 tmp.xdata (i) = 1;
632 }
633 tmp.xcidx (n) = n;
634
635 AminusSigmaB -= tmp * b.hermitian () * b *
636 tmp.transpose () * sigma;
637 }
638 else
639 AminusSigmaB -= sigma * b.hermitian () * b;
640 }
641 else
642 AminusSigmaB -= sigma * b;
643 }
644 else
645 {
646 SparseComplexMatrix sigmat (n, n, n);
647
648 // Create sigma * speye (n,n)
649 sigmat.xcidx (0) = 0;
650 for (octave_idx_type i = 0; i < n; i++)
651 {
652 sigmat.xdata (i) = sigma;
653 sigmat.xridx (i) = i;
654 sigmat.xcidx (i+1) = i + 1;
655 }
656
657 AminusSigmaB -= sigmat;
658 }
659 }
660
661 octave::math::sparse_lu<SparseComplexMatrix> fact (AminusSigmaB, Matrix(),
662 true);
663
664 L = fact.L ();
665 U = fact.U ();
666 SparseMatrix R = fact.R ();
667 for (octave_idx_type i = 0; i < n; i++)
668 r(i) = R.xdata(i);
669
670 const octave_idx_type *P2 = fact.row_perm ();
671 const octave_idx_type *Q2 = fact.col_perm ();
672
673 for (octave_idx_type j = 0; j < n; j++)
674 {
675 P[j] = P2[j];
676 Q[j] = Q2[j];
677 }
678
679 // Test condition number of LU decomposition
680 double minU = octave::numeric_limits<double>::NaN ();
681 double maxU = octave::numeric_limits<double>::NaN ();
682 for (octave_idx_type j = 0; j < n; j++)
683 {
684 double d = 0.;
685 if (U.xcidx (j+1) > U.xcidx (j)
686 && U.xridx (U.xcidx (j+1)-1) == j)
687 d = std::abs (U.xdata (U.xcidx (j+1)-1));
688
689 if (octave::math::isnan (minU) || d < minU)
690 minU = d;
691
692 if (octave::math::isnan (maxU) || d > maxU)
693 maxU = d;
694 }
695
696 double rcond = (minU / maxU);
697 // Prevent use of extra precision.
698 double rcond_plus_one = rcond + 1.0;
699
700 if (rcond_plus_one == 1.0 || octave::math::isnan (rcond))
701 warn_convergence ();
702
703 return true;
704}
705
706static bool
707LuAminusSigmaB (const ComplexMatrix& m, const ComplexMatrix& b,
708 bool cholB, const ColumnVector& permB, Complex sigma,
709 ComplexMatrix& L, ComplexMatrix& U, octave_idx_type *P,
710 octave_idx_type *Q, ColumnVector& r)
711{
712 bool have_b = ! b.isempty ();
713 octave_idx_type n = m.cols ();
714
715 // Calculate LU decomposition of 'M = A - sigma * B'
716 // P * M = L * U
717 ComplexMatrix AminusSigmaB (m);
718
719 if (std::real (sigma) != 0.0 || std::imag (sigma) != 0.0)
720 {
721 if (have_b)
722 {
723 if (cholB)
724 {
725 ComplexMatrix tmp = sigma * b.hermitian () * b;
726 const double *pB = permB.data ();
727 Complex *p = AminusSigmaB.rwdata ();
728
729 if (permB.numel ())
730 {
731 for (octave_idx_type j = 0;
732 j < b.cols (); j++)
733 for (octave_idx_type i = 0;
734 i < b.rows (); i++)
735 *p++ -= tmp.xelem (static_cast<octave_idx_type> (pB[i]),
736 static_cast<octave_idx_type> (pB[j]));
737 }
738 else
739 AminusSigmaB -= tmp;
740 }
741 else
742 AminusSigmaB -= sigma * b;
743 }
744 else
745 {
746 Complex *p = AminusSigmaB.rwdata ();
747
748 for (octave_idx_type i = 0; i < n; i++)
749 p[i*(n+1)] -= sigma;
750 }
751 }
752
753 octave::math::lu<ComplexMatrix> fact (AminusSigmaB);
754
755 L = fact.L ();
756 U = fact.U ();
757 ColumnVector P2 = fact.P_vec ();
758
759 for (octave_idx_type j = 0; j < n; j++)
760 {
761 Q[j] = j;
762 P[j] = P2(j) - 1;
763 r(j) = 1.;
764 }
765
766 // Test condition number of LU decomposition
767 double minU = octave::numeric_limits<double>::NaN ();
768 double maxU = octave::numeric_limits<double>::NaN ();
769 for (octave_idx_type j = 0; j < n; j++)
770 {
771 double d = std::abs (U.xelem (j, j));
772 if (octave::math::isnan (minU) || d < minU)
773 minU = d;
774
775 if (octave::math::isnan (maxU) || d > maxU)
776 maxU = d;
777 }
778
779 double rcond = (minU / maxU);
780 // Prevent use of extra precision.
781 double rcond_plus_one = rcond + 1.0;
782
783 if (rcond_plus_one == 1.0 || octave::math::isnan (rcond))
784 warn_convergence ();
785
786 return true;
787}
788
789template <typename M>
790octave_idx_type
791EigsRealSymmetricMatrix (const M& m, const std::string typ,
792 octave_idx_type k_arg, octave_idx_type p_arg,
793 octave_idx_type& info, Matrix& eig_vec,
794 ColumnVector& eig_val, const M& _b,
795 ColumnVector& permB, ColumnVector& resid,
796 std::ostream& os, double tol, bool rvec,
797 bool cholB, int disp, int maxit)
798{
799 F77_INT k = octave::to_f77_int (k_arg);
800 F77_INT p = octave::to_f77_int (p_arg);
801 M b(_b);
802 F77_INT n = octave::to_f77_int (m.cols ());
803 F77_INT mode = 1;
804 bool have_b = ! b.isempty ();
805 bool note3 = false;
806 char bmat = 'I';
807 double sigma = 0.;
808 M bt;
809
810 if (m.rows () != m.cols ())
811 (*current_liboctave_error_handler) ("eigs: A must be square");
812 if (have_b && (m.rows () != b.rows () || m.rows () != b.cols ()))
813 (*current_liboctave_error_handler)
814 ("eigs: B must be square and the same size as A");
815
816 if (resid.isempty ())
817 {
818 std::string rand_dist = octave::rand::distribution ();
819 octave::rand::distribution ("uniform");
820 resid = ColumnVector (octave::rand::vector (n));
821 octave::rand::distribution (rand_dist);
822 }
823 else if (m.cols () != resid.numel ())
824 (*current_liboctave_error_handler) ("eigs: opts.v0 must be n-by-1");
825
826 if (n < 3)
827 (*current_liboctave_error_handler) ("eigs: n must be at least 3");
828
829 if (p < 0)
830 {
831 p = k * 2;
832
833 if (p < 20)
834 p = 20;
835
836 if (p > n)
837 p = n;
838 }
839
840 if (k < 1 || k > n - 2)
841 (*current_liboctave_error_handler)
842 ("eigs: Invalid number of eigenvalues to extract"
843 " (must be 0 < k < n-1-1).\n"
844 " Use 'eig (full (A))' instead");
845
846 if (p <= k || p > n)
847 (*current_liboctave_error_handler)
848 ("eigs: opts.p must be greater than k and less than or equal to n");
849
850 if (have_b && cholB && ! permB.isempty ())
851 {
852 // Check the we really have a permutation vector
853 if (permB.numel () != n)
854 (*current_liboctave_error_handler) ("eigs: permB vector invalid");
855
856 Array<bool> checked (dim_vector (n, 1), false);
857 for (F77_INT i = 0; i < n; i++)
858 {
859 octave_idx_type bidx = static_cast<octave_idx_type> (permB(i));
860
861 if (checked(bidx) || bidx < 0 || bidx >= n)
862 (*current_liboctave_error_handler) ("eigs: permB vector invalid");
863 }
864 }
865
866 if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA"
867 && typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI"
868 && typ != "SI")
869 (*current_liboctave_error_handler) ("eigs: unrecognized sigma value");
870
871 if (typ == "LI" || typ == "SI" || typ == "LR" || typ == "SR")
872 (*current_liboctave_error_handler)
873 ("eigs: invalid sigma value for real symmetric problem");
874
875 if (have_b)
876 {
877 // See Note 3 dsaupd
878 note3 = true;
879 if (cholB)
880 {
881 bt = b;
882 b = b.transpose ();
883 if (permB.isempty ())
884 {
885 permB = ColumnVector (n);
886 for (F77_INT i = 0; i < n; i++)
887 permB(i) = i;
888 }
889 }
890 else
891 {
892 if (! make_cholb (b, bt, permB))
893 (*current_liboctave_error_handler)
894 ("eigs: The matrix B is not positive definite");
895 }
896 }
897
898 Array<F77_INT> ip (dim_vector (11, 1));
899 F77_INT *iparam = ip.rwdata ();
900
901 ip(0) = 1; //ishift
902 ip(1) = 0; // ip(1) not referenced
903 ip(2) = maxit; // mxiter, maximum number of iterations
904 ip(3) = 1; // NB blocksize in recurrence
905 ip(4) = 0; // nconv, number of Ritz values that satisfy convergence
906 ip(5) = 0; //ip(5) not referenced
907 ip(6) = mode; // mode
908 ip(7) = 0;
909 ip(8) = 0;
910 ip(9) = 0;
911 ip(10) = 0;
912 // ip(7) to ip(10) return values
913
914 Array<F77_INT> iptr (dim_vector (14, 1));
915 F77_INT *ipntr = iptr.rwdata ();
916
917 F77_INT ido = 0;
918 int iter = 0;
919 F77_INT elems;
920 F77_INT lwork;
921
922 if (octave::math::int_multiply_overflow (n, p, &elems))
923 (*current_liboctave_error_handler)
924 ("eigs: array too large for Fortran integers");
925 if (octave::math::int_multiply_overflow (p, p + 8, &lwork))
926 (*current_liboctave_error_handler)
927 ("eigs: array too large for Fortran integers");
928
929 OCTAVE_LOCAL_BUFFER (double, v, elems);
930 OCTAVE_LOCAL_BUFFER (double, workl, lwork);
931 OCTAVE_LOCAL_BUFFER (double, workd, 3 * n);
932 double *presid = resid.rwdata ();
933
934 do
935 {
936 F77_INT tmp_info = octave::to_f77_int (info);
937
938 F77_FUNC (dsaupd, DSAUPD)
939 (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
940 F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2),
941 k, tol, presid, p, v, n, iparam,
942 ipntr, workd, workl, lwork, tmp_info
943 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
944
945 info = tmp_info;
946
947 if (disp > 0 && ! octave::math::isnan (workl[iptr (5)-1]))
948 {
949 if (iter++)
950 {
951 os << "Iteration " << iter - 1 <<
952 ": a few Ritz values of the " << p << "-by-" <<
953 p << " matrix\n";
954 if (ido == 99) // convergence
955 {
956 for (F77_INT i = 0; i < k; i++)
957 os << " " << workl[iptr(5)+i-1] << "\n";
958 }
959 else
960 {
961 // the wanted Ritz estimates are at the end
962 for (F77_INT i = p - k; i < p; i++)
963 os << " " << workl[iptr(5)+i-1] << "\n";
964 }
965 }
966
967 // This is a kludge, as ARPACK doesn't give its
968 // iteration pointer. But as workl[iptr(5)-1] is
969 // an output value updated at each iteration, setting
970 // a value in this array to NaN and testing for it
971 // is a way of obtaining the iteration counter.
972 if (ido != 99)
973 workl[iptr(5)-1] = octave::numeric_limits<double>::NaN ();
974 }
975
976 if (ido == -1 || ido == 1 || ido == 2)
977 {
978 if (have_b)
979 {
980 Matrix mtmp (n, 1);
981 for (F77_INT i = 0; i < n; i++)
982 mtmp(i, 0) = workd[i + iptr(0) - 1];
983
984 mtmp = ltsolve (b, permB, m * utsolve (bt, permB, mtmp));
985
986 for (F77_INT i = 0; i < n; i++)
987 workd[i+iptr(1)-1] = mtmp(i, 0);
988 }
989 else if (! vector_product (m, workd + iptr(0) - 1,
990 workd + iptr(1) - 1))
991 break;
992 }
993 else
994 {
995 if (info < 0)
996 (*current_liboctave_error_handler)
997 ("eigs: error in dsaupd: %s",
998 arpack_errno2str (info, "dsaupd").c_str ());
999
1000 break;
1001 }
1002 }
1003 while (1);
1004
1005 F77_INT info2;
1006
1007 // We have a problem in that the size of the C++ bool
1008 // type relative to the fortran logical type. It appears
1009 // that fortran uses 4- or 8-bytes per logical and C++ 1-byte
1010 // per bool, though this might be system dependent. As
1011 // long as the HOWMNY arg is not "S", the logical array
1012 // is just workspace for ARPACK, so use int type to
1013 // avoid problems.
1014 Array<F77_INT> s (dim_vector (p, 1));
1015 F77_INT *sel = s.rwdata ();
1016
1017 eig_vec.resize (n, k);
1018 double *z = eig_vec.rwdata ();
1019
1020 eig_val.resize (k);
1021 double *d = eig_val.rwdata ();
1022
1023 F77_FUNC (dseupd, DSEUPD)
1024 (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma,
1025 F77_CONST_CHAR_ARG2 (&bmat, 1), n,
1026 F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), k, tol, presid, p, v, n, iparam,
1027 ipntr, workd, workl, lwork, info2 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1)
1028 F77_CHAR_ARG_LEN(2));
1029
1030 if (info2 == 0)
1031 {
1032 for (F77_INT i = ip(4); i < k; i++)
1033 d[i] = octave::numeric_limits<double>::NaN ();
1034 F77_INT k2 = ip(4) / 2;
1035 if (typ != "SM" && typ != "BE" && ! (typ == "SA" && rvec))
1036 {
1037 for (F77_INT i = 0; i < k2; i++)
1038 {
1039 double dtmp = d[i];
1040 d[i] = d[ip(4) - i - 1];
1041 d[ip(4) - i - 1] = dtmp;
1042 }
1043 }
1044
1045 if (rvec)
1046 {
1047 for (F77_INT i = ip(4); i < k; i++)
1048 {
1049 F77_INT off1 = i * n;
1050 for (F77_INT j = 0; j < n; j++)
1051 z[off1 + j] = octave::numeric_limits<double>::NaN ();
1052 }
1053 if (typ != "SM" && typ != "BE" && typ != "SA")
1054 {
1055 OCTAVE_LOCAL_BUFFER (double, dtmp, n);
1056
1057 for (F77_INT i = 0; i < k2; i++)
1058 {
1059 F77_INT off1 = i * n;
1060 F77_INT off2 = (ip(4) - i - 1) * n;
1061
1062 if (off1 == off2)
1063 continue;
1064
1065 for (F77_INT j = 0; j < n; j++)
1066 dtmp[j] = z[off1 + j];
1067
1068 for (F77_INT j = 0; j < n; j++)
1069 z[off1 + j] = z[off2 + j];
1070
1071 for (F77_INT j = 0; j < n; j++)
1072 z[off2 + j] = dtmp[j];
1073 }
1074 }
1075
1076 if (note3)
1077 eig_vec = utsolve (bt, permB, eig_vec);
1078 }
1079 }
1080 else
1081 (*current_liboctave_error_handler)
1082 ("eigs: error in dseupd: %s",
1083 arpack_errno2str (info2, "dseupd").c_str ());
1084
1085 return ip(4);
1086}
1087
1088template <typename M>
1089octave_idx_type
1090EigsRealSymmetricMatrixShift (const M& m, double sigma,
1091 octave_idx_type k_arg, octave_idx_type p_arg,
1092 octave_idx_type& info, Matrix& eig_vec,
1093 ColumnVector& eig_val, const M& _b,
1094 ColumnVector& permB, ColumnVector& resid,
1095 std::ostream& os, double tol, bool rvec,
1096 bool cholB, int disp, int maxit)
1097{
1098 F77_INT k = octave::to_f77_int (k_arg);
1099 F77_INT p = octave::to_f77_int (p_arg);
1100 M b(_b);
1101 F77_INT n = octave::to_f77_int (m.cols ());
1102 F77_INT mode = 3;
1103 bool have_b = ! b.isempty ();
1104 std::string typ = "LM";
1105
1106 if (m.rows () != m.cols ())
1107 (*current_liboctave_error_handler) ("eigs: A must be square");
1108 if (have_b && (m.rows () != b.rows () || m.rows () != b.cols ()))
1109 (*current_liboctave_error_handler)
1110 ("eigs: B must be square and the same size as A");
1111
1112 // FIXME: The "SM" type for mode 1 seems unstable though faster!!
1113 //if (! std::abs (sigma))
1114 // return EigsRealSymmetricMatrix (m, "SM", k, p, info, eig_vec, eig_val,
1115 // _b, permB, resid, os, tol, rvec, cholB,
1116 // disp, maxit);
1117
1118 if (resid.isempty ())
1119 {
1120 std::string rand_dist = octave::rand::distribution ();
1121 octave::rand::distribution ("uniform");
1122 resid = ColumnVector (octave::rand::vector (n));
1123 octave::rand::distribution (rand_dist);
1124 }
1125 else if (m.cols () != resid.numel ())
1126 (*current_liboctave_error_handler) ("eigs: opts.v0 must be n-by-1");
1127
1128 if (n < 3)
1129 (*current_liboctave_error_handler) ("eigs: n must be at least 3");
1130
1131 if (k <= 0 || k >= n - 1)
1132 (*current_liboctave_error_handler)
1133 ("eigs: Invalid number of eigenvalues to extract"
1134 " (must be 0 < k < n-1-1).\n"
1135 " Use 'eig (full (A))' instead");
1136
1137 if (p < 0)
1138 {
1139 p = k * 2;
1140
1141 if (p < 20)
1142 p = 20;
1143
1144 if (p > n)
1145 p = n;
1146 }
1147
1148 if (p <= k || p > n)
1149 (*current_liboctave_error_handler)
1150 ("eigs: opts.p must be greater than k and less than or equal to n");
1151
1152 if (have_b && cholB && ! permB.isempty ())
1153 {
1154 // Check the we really have a permutation vector
1155 if (permB.numel () != n)
1156 (*current_liboctave_error_handler) ("eigs: permB vector invalid");
1157
1158 Array<bool> checked (dim_vector (n, 1), false);
1159 for (F77_INT i = 0; i < n; i++)
1160 {
1161 octave_idx_type bidx = static_cast<octave_idx_type> (permB(i));
1162
1163 if (checked(bidx) || bidx < 0 || bidx >= n)
1164 (*current_liboctave_error_handler) ("eigs: permB vector invalid");
1165 }
1166 }
1167
1168 char bmat = 'I';
1169 if (have_b)
1170 bmat = 'G';
1171
1172 Array<F77_INT> ip (dim_vector (11, 1));
1173 F77_INT *iparam = ip.rwdata ();
1174
1175 ip(0) = 1; //ishift
1176 ip(1) = 0; // ip(1) not referenced
1177 ip(2) = maxit; // mxiter, maximum number of iterations
1178 ip(3) = 1; // NB blocksize in recurrence
1179 ip(4) = 0; // nconv, number of Ritz values that satisfy convergence
1180 ip(5) = 0; //ip(5) not referenced
1181 ip(6) = mode; // mode
1182 ip(7) = 0;
1183 ip(8) = 0;
1184 ip(9) = 0;
1185 ip(10) = 0;
1186 // ip(7) to ip(10) return values
1187
1188 Array<F77_INT> iptr (dim_vector (14, 1));
1189 F77_INT *ipntr = iptr.rwdata ();
1190
1191 F77_INT ido = 0;
1192 int iter = 0;
1193 M L, U;
1194 ColumnVector r(n);
1195
1196 OCTAVE_LOCAL_BUFFER (octave_idx_type, P, (have_b ? b.rows () : m.rows ()));
1197 OCTAVE_LOCAL_BUFFER (octave_idx_type, Q, (have_b ? b.cols () : m.cols ()));
1198
1199 if (! LuAminusSigmaB (m, b, cholB, permB, sigma, L, U, P, Q, r))
1200 return -1;
1201
1202 F77_INT elems;
1203 F77_INT lwork;
1204
1205 if (octave::math::int_multiply_overflow (n, p, &elems))
1206 (*current_liboctave_error_handler)
1207 ("eigs: array too large for Fortran integers");
1208 if (octave::math::int_multiply_overflow (p, p + 8, &lwork))
1209 (*current_liboctave_error_handler)
1210 ("eigs: array too large for Fortran integers");
1211
1212 OCTAVE_LOCAL_BUFFER (double, v, elems);
1213 OCTAVE_LOCAL_BUFFER (double, workl, lwork);
1214 OCTAVE_LOCAL_BUFFER (double, workd, 3 * n);
1215 double *presid = resid.rwdata ();
1216
1217 do
1218 {
1219 F77_INT tmp_info = octave::to_f77_int (info);
1220
1221 F77_FUNC (dsaupd, DSAUPD)
1222 (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
1223 F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2),
1224 k, tol, presid, p, v, n, iparam,
1225 ipntr, workd, workl, lwork, tmp_info
1226 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
1227 info = tmp_info;
1228
1229 if (disp > 0 && ! octave::math::isnan (workl[iptr (5)-1]))
1230 {
1231 if (iter++)
1232 {
1233 os << "Iteration " << iter - 1 <<
1234 ": a few Ritz values of the " << p << "-by-" <<
1235 p << " matrix\n";
1236 if (ido == 99) // convergence
1237 {
1238 for (F77_INT i = 0; i < k; i++)
1239 os << " " << workl[iptr(5)+i-1] << "\n";
1240 }
1241 else
1242 {
1243 // the wanted Ritz estimates are at the end
1244 for (F77_INT i = p - k; i < p; i++)
1245 os << " " << workl[iptr(5)+i-1] << "\n";
1246 }
1247 }
1248
1249 // This is a kludge, as ARPACK doesn't give its
1250 // iteration pointer. But as workl[iptr(5)-1] is
1251 // an output value updated at each iteration, setting
1252 // a value in this array to NaN and testing for it
1253 // is a way of obtaining the iteration counter.
1254 if (ido != 99)
1255 workl[iptr(5)-1] = octave::numeric_limits<double>::NaN ();
1256 }
1257
1258 if (ido == -1 || ido == 1 || ido == 2)
1259 {
1260 if (have_b)
1261 {
1262 if (ido == -1)
1263 {
1264 OCTAVE_LOCAL_BUFFER (double, dtmp, n);
1265
1266 vector_product (b, workd+iptr(0)-1, dtmp);
1267
1268 Matrix tmp (n, 1);
1269
1270 for (F77_INT i = 0; i < n; i++)
1271 tmp(i, 0) = dtmp[P[i]] / r(P[i]);
1272
1273 lusolve (L, U, tmp);
1274
1275 double *ip2 = workd+iptr(1)-1;
1276 for (F77_INT i = 0; i < n; i++)
1277 ip2[Q[i]] = tmp(i, 0);
1278 }
1279 else if (ido == 2)
1280 vector_product (b, workd+iptr(0)-1, workd+iptr(1)-1);
1281 else
1282 {
1283 double *ip2 = workd+iptr(2)-1;
1284 Matrix tmp (n, 1);
1285
1286 for (F77_INT i = 0; i < n; i++)
1287 tmp(i, 0) = ip2[P[i]] / r(P[i]);
1288
1289 lusolve (L, U, tmp);
1290
1291 ip2 = workd+iptr(1)-1;
1292 for (F77_INT i = 0; i < n; i++)
1293 ip2[Q[i]] = tmp(i, 0);
1294 }
1295 }
1296 else
1297 {
1298 // ido cannot be 2 for non-generalized problems (see dsaupd2).
1299 double *ip2 = workd+iptr(0)-1;
1300 Matrix tmp (n, 1);
1301
1302 for (F77_INT i = 0; i < n; i++)
1303 tmp(i, 0) = ip2[P[i]] / r(P[i]);
1304
1305 lusolve (L, U, tmp);
1306
1307 ip2 = workd+iptr(1)-1;
1308 for (F77_INT i = 0; i < n; i++)
1309 ip2[Q[i]] = tmp(i, 0);
1310 }
1311 }
1312 else
1313 {
1314 if (info < 0)
1315 (*current_liboctave_error_handler)
1316 ("eigs: error in dsaupd: %s",
1317 arpack_errno2str (info, "dsaupd").c_str ());
1318
1319 break;
1320 }
1321 }
1322 while (1);
1323
1324 F77_INT info2;
1325
1326 // We have a problem in that the size of the C++ bool
1327 // type relative to the fortran logical type. It appears
1328 // that fortran uses 4- or 8-bytes per logical and C++ 1-byte
1329 // per bool, though this might be system dependent. As
1330 // long as the HOWMNY arg is not "S", the logical array
1331 // is just workspace for ARPACK, so use int type to
1332 // avoid problems.
1333 Array<F77_INT> s (dim_vector (p, 1));
1334 F77_INT *sel = s.rwdata ();
1335
1336 eig_vec.resize (n, k);
1337 double *z = eig_vec.rwdata ();
1338
1339 eig_val.resize (k);
1340 double *d = eig_val.rwdata ();
1341
1342 F77_FUNC (dseupd, DSEUPD)
1343 (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma,
1344 F77_CONST_CHAR_ARG2 (&bmat, 1), n,
1345 F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2),
1346 k, tol, presid, p, v, n, iparam, ipntr, workd, workl, lwork, info2
1347 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
1348
1349 if (info2 == 0)
1350 {
1351 for (F77_INT i = ip(4); i < k; i++)
1352 d[i] = octave::numeric_limits<double>::NaN ();
1353 F77_INT k2 = ip(4) / 2;
1354 for (F77_INT i = 0; i < k2; i++)
1355 {
1356 double dtmp = d[i];
1357 d[i] = d[ip(4) - i - 1];
1358 d[ip(4) - i - 1] = dtmp;
1359 }
1360
1361 if (rvec)
1362 {
1363 OCTAVE_LOCAL_BUFFER (double, dtmp, n);
1364
1365 for (F77_INT i = ip(4); i < k; i++)
1366 {
1367 F77_INT off1 = i * n;
1368 for (F77_INT j = 0; j < n; j++)
1369 z[off1 + j] = octave::numeric_limits<double>::NaN ();
1370 }
1371 for (F77_INT i = 0; i < k2; i++)
1372 {
1373 F77_INT off1 = i * n;
1374 F77_INT off2 = (ip(4) - i - 1) * n;
1375
1376 if (off1 == off2)
1377 continue;
1378
1379 for (F77_INT j = 0; j < n; j++)
1380 dtmp[j] = z[off1 + j];
1381
1382 for (F77_INT j = 0; j < n; j++)
1383 z[off1 + j] = z[off2 + j];
1384
1385 for (F77_INT j = 0; j < n; j++)
1386 z[off2 + j] = dtmp[j];
1387 }
1388 }
1389 }
1390 else
1391 (*current_liboctave_error_handler)
1392 ("eigs: error in dseupd: %s",
1393 arpack_errno2str (info2, "dseupd").c_str ());
1394
1395 return ip(4);
1396}
1397
1398template <typename M>
1399octave_idx_type
1400EigsRealSymmetricFunc (EigsFunc fcn, octave_idx_type n_arg,
1401 const std::string& _typ, double sigma,
1402 octave_idx_type k_arg, octave_idx_type p_arg,
1403 octave_idx_type& info, Matrix& eig_vec,
1404 ColumnVector& eig_val, const M& _b,
1405 ColumnVector& permB, ColumnVector& resid,
1406 std::ostream& os, double tol, bool rvec,
1407 bool cholB, int disp, int maxit)
1408{
1409 F77_INT n = octave::to_f77_int (n_arg);
1410 F77_INT k = octave::to_f77_int (k_arg);
1411 F77_INT p = octave::to_f77_int (p_arg);
1412 M b(_b);
1413 std::string typ (_typ);
1414 bool have_sigma = (sigma ? true : false);
1415 bool have_b = ! b.isempty ();
1416 bool note3 = false;
1417 char bmat = 'I';
1418 F77_INT mode = 1;
1419 int err = 0;
1420 M bt;
1421
1422 if (resid.isempty ())
1423 {
1424 std::string rand_dist = octave::rand::distribution ();
1425 octave::rand::distribution ("uniform");
1426 resid = ColumnVector (octave::rand::vector (n));
1427 octave::rand::distribution (rand_dist);
1428 }
1429 else if (n != resid.numel ())
1430 (*current_liboctave_error_handler) ("eigs: opts.v0 must be n-by-1");
1431
1432 if (n < 3)
1433 (*current_liboctave_error_handler) ("eigs: n must be at least 3");
1434
1435 if (p < 0)
1436 {
1437 p = k * 2;
1438
1439 if (p < 20)
1440 p = 20;
1441
1442 if (p > n)
1443 p = n;
1444 }
1445
1446 if (k <= 0 || k >= n - 1)
1447 (*current_liboctave_error_handler)
1448 ("eigs: Invalid number of eigenvalues to extract"
1449 " (must be 0 < k < n-1).\n"
1450 " Use 'eig (full (A))' instead");
1451
1452 if (p <= k || p > n)
1453 (*current_liboctave_error_handler)
1454 ("eigs: opts.p must be greater than k and less than or equal to n");
1455
1456 if (have_b && cholB && ! permB.isempty ())
1457 {
1458 // Check the we really have a permutation vector
1459 if (permB.numel () != n)
1460 (*current_liboctave_error_handler) ("eigs: permB vector invalid");
1461
1462 Array<bool> checked (dim_vector (n, 1), false);
1463 for (F77_INT i = 0; i < n; i++)
1464 {
1465 octave_idx_type bidx = static_cast<octave_idx_type> (permB(i));
1466
1467 if (checked(bidx) || bidx < 0 || bidx >= n)
1468 (*current_liboctave_error_handler) ("eigs: permB vector invalid");
1469 }
1470 }
1471
1472 if (! have_sigma)
1473 {
1474 if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA"
1475 && typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI"
1476 && typ != "SI")
1477 (*current_liboctave_error_handler) ("eigs: unrecognized sigma value");
1478
1479 if (typ == "LI" || typ == "SI" || typ == "LR" || typ == "SR")
1480 (*current_liboctave_error_handler)
1481 ("eigs: invalid sigma value for real symmetric problem");
1482
1483 if (typ != "SM" && have_b)
1484 note3 = true;
1485
1486 if (typ == "SM")
1487 {
1488 typ = "LM";
1489 sigma = 0.;
1490 mode = 3;
1491 if (have_b)
1492 bmat = 'G';
1493 }
1494 }
1495 else if (! std::abs (sigma))
1496 {
1497 typ = "SM";
1498 if (have_b)
1499 bmat = 'G';
1500 }
1501 else
1502 {
1503 typ = "LM";
1504 mode = 3;
1505 if (have_b)
1506 bmat = 'G';
1507 }
1508
1509 if (mode == 1 && have_b)
1510 {
1511 // See Note 3 dsaupd
1512 note3 = true;
1513 if (cholB)
1514 {
1515 bt = b;
1516 b = b.transpose ();
1517 if (permB.isempty ())
1518 {
1519 permB = ColumnVector (n);
1520 for (F77_INT i = 0; i < n; i++)
1521 permB(i) = i;
1522 }
1523 }
1524 else
1525 {
1526 if (! make_cholb (b, bt, permB))
1527 (*current_liboctave_error_handler)
1528 ("eigs: The matrix B is not positive definite");
1529 }
1530 }
1531
1532 Array<F77_INT> ip (dim_vector (11, 1));
1533 F77_INT *iparam = ip.rwdata ();
1534
1535 ip(0) = 1; //ishift
1536 ip(1) = 0; // ip(1) not referenced
1537 ip(2) = maxit; // mxiter, maximum number of iterations
1538 ip(3) = 1; // NB blocksize in recurrence
1539 ip(4) = 0; // nconv, number of Ritz values that satisfy convergence
1540 ip(5) = 0; //ip(5) not referenced
1541 ip(6) = mode; // mode
1542 ip(7) = 0;
1543 ip(8) = 0;
1544 ip(9) = 0;
1545 ip(10) = 0;
1546 // ip(7) to ip(10) return values
1547
1548 Array<F77_INT> iptr (dim_vector (14, 1));
1549 F77_INT *ipntr = iptr.rwdata ();
1550
1551 F77_INT ido = 0;
1552 int iter = 0;
1553 F77_INT elems;
1554 F77_INT lwork;
1555
1556 if (octave::math::int_multiply_overflow (n, p, &elems))
1557 (*current_liboctave_error_handler)
1558 ("eigs: array too large for Fortran integers");
1559 if (octave::math::int_multiply_overflow (p, p + 8, &lwork))
1560 (*current_liboctave_error_handler)
1561 ("eigs: array too large for Fortran integers");
1562
1563 OCTAVE_LOCAL_BUFFER (double, v, elems);
1564 OCTAVE_LOCAL_BUFFER (double, workl, lwork);
1565 OCTAVE_LOCAL_BUFFER (double, workd, 3 * n);
1566 double *presid = resid.rwdata ();
1567
1568 do
1569 {
1570 F77_INT tmp_info = octave::to_f77_int (info);
1571
1572 F77_FUNC (dsaupd, DSAUPD)
1573 (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
1574 F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2),
1575 k, tol, presid, p, v, n, iparam,
1576 ipntr, workd, workl, lwork, tmp_info
1577 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
1578
1579 info = tmp_info;
1580
1581 if (disp > 0 && ! octave::math::isnan (workl[iptr (5)-1]))
1582 {
1583 if (iter++)
1584 {
1585 os << "Iteration " << iter - 1 <<
1586 ": a few Ritz values of the " << p << "-by-" <<
1587 p << " matrix\n";
1588 if (ido == 99) // convergence
1589 {
1590 for (F77_INT i = 0; i < k; i++)
1591 os << " " << workl[iptr(5)+i-1] << "\n";
1592 }
1593 else
1594 {
1595 // the wanted Ritz estimates are at the end
1596 for (F77_INT i = p - k; i < p; i++)
1597 os << " " << workl[iptr(5)+i-1] << "\n";
1598 }
1599 }
1600
1601 // This is a kludge, as ARPACK doesn't give its
1602 // iteration pointer. But as workl[iptr(5)-1] is
1603 // an output value updated at each iteration, setting
1604 // a value in this array to NaN and testing for it
1605 // is a way of obtaining the iteration counter.
1606 if (ido != 99)
1607 workl[iptr(5)-1] = octave::numeric_limits<double>::NaN ();
1608 }
1609
1610 if (ido == -1 || ido == 1 || ido == 2)
1611 {
1612 if (have_b)
1613 {
1614 if (mode == 1) // regular mode with factorized B
1615 {
1616 Matrix mtmp (n, 1);
1617 for (F77_INT i = 0; i < n; i++)
1618 mtmp(i, 0) = workd[i + iptr(0) - 1];
1619
1620 mtmp = utsolve (bt, permB, mtmp);
1621 ColumnVector y = fcn (mtmp, err);
1622
1623 if (err)
1624 return false;
1625
1626 mtmp = ltsolve (b, permB, y);
1627
1628 for (F77_INT i = 0; i < n; i++)
1629 workd[i+iptr(1)-1] = mtmp(i, 0);
1630 }
1631 else // shift-invert mode
1632 {
1633 if (ido == -1)
1634 {
1635 OCTAVE_LOCAL_BUFFER (double, dtmp, n);
1636
1637 vector_product (b, workd+iptr(0)-1, dtmp);
1638
1639 ColumnVector x(n);
1640
1641 for (F77_INT i = 0; i < n; i++)
1642 x(i) = dtmp[i];
1643
1644 ColumnVector y = fcn (x, err);
1645
1646 if (err)
1647 return false;
1648
1649 double *ip2 = workd + iptr(1) - 1;
1650 for (F77_INT i = 0; i < n; i++)
1651 ip2[i] = y(i);
1652 }
1653 else if (ido == 2)
1654 vector_product (b, workd+iptr(0)-1, workd+iptr(1)-1);
1655 else
1656 {
1657 double *ip2 = workd+iptr(2)-1;
1658 ColumnVector x(n);
1659
1660 for (F77_INT i = 0; i < n; i++)
1661 x(i) = *ip2++;
1662
1663 ColumnVector y = fcn (x, err);
1664
1665 if (err)
1666 return false;
1667
1668 ip2 = workd + iptr(1) - 1;
1669 for (F77_INT i = 0; i < n; i++)
1670 *ip2++ = y(i);
1671 }
1672 }
1673 }
1674 else
1675 {
1676 double *ip2 = workd + iptr(0) - 1;
1677 ColumnVector x(n);
1678
1679 for (F77_INT i = 0; i < n; i++)
1680 x(i) = *ip2++;
1681
1682 ColumnVector y = fcn (x, err);
1683
1684 if (err)
1685 return false;
1686
1687 ip2 = workd + iptr(1) - 1;
1688 for (F77_INT i = 0; i < n; i++)
1689 *ip2++ = y(i);
1690 }
1691 }
1692 else
1693 {
1694 if (info < 0)
1695 (*current_liboctave_error_handler)
1696 ("eigs: error in dsaupd: %s",
1697 arpack_errno2str (info, "dsaupd").c_str ());
1698
1699 break;
1700 }
1701 }
1702 while (1);
1703
1704 F77_INT info2;
1705
1706 // We have a problem in that the size of the C++ bool
1707 // type relative to the fortran logical type. It appears
1708 // that fortran uses 4- or 8-bytes per logical and C++ 1-byte
1709 // per bool, though this might be system dependent. As
1710 // long as the HOWMNY arg is not "S", the logical array
1711 // is just workspace for ARPACK, so use int type to
1712 // avoid problems.
1713 Array<F77_INT> s (dim_vector (p, 1));
1714 F77_INT *sel = s.rwdata ();
1715
1716 eig_vec.resize (n, k);
1717 double *z = eig_vec.rwdata ();
1718
1719 eig_val.resize (k);
1720 double *d = eig_val.rwdata ();
1721
1722 F77_FUNC (dseupd, DSEUPD)
1723 (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma,
1724 F77_CONST_CHAR_ARG2 (&bmat, 1), n,
1725 F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2),
1726 k, tol, presid, p, v, n, iparam, ipntr, workd, workl, lwork, info2
1727 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
1728
1729 if (info2 == 0)
1730 {
1731 for (F77_INT i = ip(4); i < k; i++)
1732 d[i] = octave::numeric_limits<double>::NaN ();
1733 F77_INT k2 = ip(4) / 2;
1734 if (mode == 3 || (mode == 1 && typ != "SM" && typ != "BE"))
1735 {
1736 for (F77_INT i = 0; i < k2; i++)
1737 {
1738 double dtmp = d[i];
1739 d[i] = d[ip(4) - i - 1];
1740 d[ip(4) - i - 1] = dtmp;
1741 }
1742 }
1743
1744 if (rvec)
1745 {
1746 for (F77_INT i = ip(4); i < k; i++)
1747 {
1748 F77_INT off1 = i * n;
1749 for (F77_INT j = 0; j < n; j++)
1750 z[off1 + j] = octave::numeric_limits<double>::NaN ();
1751 }
1752 if (mode == 3 || (mode == 1 && typ != "SM" && typ != "BE"))
1753 {
1754 OCTAVE_LOCAL_BUFFER (double, dtmp, n);
1755
1756 for (F77_INT i = 0; i < k2; i++)
1757 {
1758 F77_INT off1 = i * n;
1759 F77_INT off2 = (ip(4) - i - 1) * n;
1760
1761 if (off1 == off2)
1762 continue;
1763
1764 for (F77_INT j = 0; j < n; j++)
1765 dtmp[j] = z[off1 + j];
1766
1767 for (F77_INT j = 0; j < n; j++)
1768 z[off1 + j] = z[off2 + j];
1769
1770 for (F77_INT j = 0; j < n; j++)
1771 z[off2 + j] = dtmp[j];
1772 }
1773 }
1774 if (note3)
1775 eig_vec = utsolve (bt, permB, eig_vec);
1776 }
1777 }
1778 else
1779 (*current_liboctave_error_handler)
1780 ("eigs: error in dseupd: %s",
1781 arpack_errno2str (info2, "dseupd").c_str ());
1782
1783 return ip(4);
1784}
1785
1786template <typename M>
1787octave_idx_type
1788EigsRealNonSymmetricMatrix (const M& m, const std::string typ,
1789 octave_idx_type k_arg, octave_idx_type p_arg,
1790 octave_idx_type& info, ComplexMatrix& eig_vec,
1791 ComplexColumnVector& eig_val, const M& _b,
1792 ColumnVector& permB, ColumnVector& resid,
1793 std::ostream& os, double tol, bool rvec,
1794 bool cholB, int disp, int maxit)
1795{
1796 F77_INT k = octave::to_f77_int (k_arg);
1797 F77_INT p = octave::to_f77_int (p_arg);
1798 M b(_b);
1799 F77_INT n = octave::to_f77_int (m.cols ());
1800 F77_INT mode = 1;
1801 bool have_b = ! b.isempty ();
1802 bool note3 = false;
1803 char bmat = 'I';
1804 double sigmar = 0.;
1805 double sigmai = 0.;
1806 M bt;
1807
1808 if (m.rows () != m.cols ())
1809 (*current_liboctave_error_handler) ("eigs: A must be square");
1810 if (have_b && (m.rows () != b.rows () || m.rows () != b.cols ()))
1811 (*current_liboctave_error_handler)
1812 ("eigs: B must be square and the same size as A");
1813
1814 if (resid.isempty ())
1815 {
1816 std::string rand_dist = octave::rand::distribution ();
1817 octave::rand::distribution ("uniform");
1818 resid = ColumnVector (octave::rand::vector (n));
1819 octave::rand::distribution (rand_dist);
1820 }
1821 else if (m.cols () != resid.numel ())
1822 (*current_liboctave_error_handler) ("eigs: opts.v0 must be n-by-1");
1823
1824 if (n < 3)
1825 (*current_liboctave_error_handler) ("eigs: n must be at least 3");
1826
1827 if (p < 0)
1828 {
1829 p = k * 2 + 1;
1830
1831 if (p < 20)
1832 p = 20;
1833
1834 if (p > n)
1835 p = n;
1836 }
1837
1838 if (k <= 0 || k >= n - 1)
1839 (*current_liboctave_error_handler)
1840 ("eigs: Invalid number of eigenvalues to extract"
1841 " (must be 0 < k < n-1).\n"
1842 " Use 'eig (full (A))' instead");
1843
1844 if (p <= k || p > n)
1845 (*current_liboctave_error_handler)
1846 ("eigs: opts.p must be greater than k and less than or equal to n");
1847
1848 if (have_b && cholB && ! permB.isempty ())
1849 {
1850 // Check the we really have a permutation vector
1851 if (permB.numel () != n)
1852 (*current_liboctave_error_handler) ("eigs: permB vector invalid");
1853
1854 Array<bool> checked (dim_vector (n, 1), false);
1855 for (F77_INT i = 0; i < n; i++)
1856 {
1857 octave_idx_type bidx = static_cast<octave_idx_type> (permB(i));
1858
1859 if (checked(bidx) || bidx < 0 || bidx >= n)
1860 (*current_liboctave_error_handler) ("eigs: permB vector invalid");
1861 }
1862 }
1863
1864 if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA"
1865 && typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI"
1866 && typ != "SI")
1867 (*current_liboctave_error_handler) ("eigs: unrecognized sigma value");
1868
1869 if (typ == "LA" || typ == "SA" || typ == "BE")
1870 (*current_liboctave_error_handler)
1871 ("eigs: invalid sigma value for unsymmetric problem");
1872
1873 if (have_b)
1874 {
1875 // See Note 3 dsaupd
1876 note3 = true;
1877 if (cholB)
1878 {
1879 bt = b;
1880 b = b.transpose ();
1881 if (permB.isempty ())
1882 {
1883 permB = ColumnVector (n);
1884 for (F77_INT i = 0; i < n; i++)
1885 permB(i) = i;
1886 }
1887 }
1888 else
1889 {
1890 if (! make_cholb (b, bt, permB))
1891 (*current_liboctave_error_handler)
1892 ("eigs: The matrix B is not positive definite");
1893 }
1894 }
1895
1896 Array<F77_INT> ip (dim_vector (11, 1));
1897 F77_INT *iparam = ip.rwdata ();
1898
1899 ip(0) = 1; //ishift
1900 ip(1) = 0; // ip(1) not referenced
1901 ip(2) = maxit; // mxiter, maximum number of iterations
1902 ip(3) = 1; // NB blocksize in recurrence
1903 ip(4) = 0; // nconv, number of Ritz values that satisfy convergence
1904 ip(5) = 0; //ip(5) not referenced
1905 ip(6) = mode; // mode
1906 ip(7) = 0;
1907 ip(8) = 0;
1908 ip(9) = 0;
1909 ip(10) = 0;
1910 // ip(7) to ip(10) return values
1911
1912 Array<F77_INT> iptr (dim_vector (14, 1));
1913 F77_INT *ipntr = iptr.rwdata ();
1914
1915 F77_INT ido = 0;
1916 int iter = 0;
1917 F77_INT elems;
1918 F77_INT lwork;
1919
1920 if (octave::math::int_multiply_overflow (n, p + 1, &elems))
1921 (*current_liboctave_error_handler)
1922 ("eigs: array too large for Fortran integers");
1923 if (octave::math::int_multiply_overflow (3 * p, p + 2, &lwork))
1924 (*current_liboctave_error_handler)
1925 ("eigs: array too large for Fortran integers");
1926
1927 OCTAVE_LOCAL_BUFFER (double, v, elems);
1928 OCTAVE_LOCAL_BUFFER (double, workl, lwork + 1);
1929 OCTAVE_LOCAL_BUFFER (double, workd, 3 * n + 1);
1930 double *presid = resid.rwdata ();
1931
1932 do
1933 {
1934 F77_INT tmp_info = octave::to_f77_int (info);
1935
1936 // On exit, ip(4) <= k + 1 is the number of converged eigenvalues.
1937 // See dnaupd2.
1938 F77_FUNC (dnaupd, DNAUPD)
1939 (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
1940 F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2),
1941 k, tol, presid, p, v, n, iparam,
1942 ipntr, workd, workl, lwork, tmp_info
1943 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
1944 // k is not changed
1945
1946 info = tmp_info;
1947
1948 if (disp > 0 && ! octave::math::isnan(workl[iptr(5)-1]))
1949 {
1950 if (iter++)
1951 {
1952 os << "Iteration " << iter - 1 <<
1953 ": a few Ritz values of the " << p << "-by-" <<
1954 p << " matrix\n";
1955 if (ido == 99) // convergence
1956 {
1957 os << " " << workl[iptr(5)+k] << "\n";
1958 for (F77_INT i = 0; i < k; i++)
1959 os << " " << workl[iptr(5)+i-1] << "\n";
1960 }
1961 else
1962 {
1963 // the wanted Ritz estimates are at the end
1964 for (F77_INT i = p - k - 1; i < p; i++)
1965 os << " " << workl[iptr(5)+i-1] << "\n";
1966 }
1967 }
1968
1969 // This is a kludge, as ARPACK doesn't give its
1970 // iteration pointer. But as workl[iptr(5)-1] is
1971 // an output value updated at each iteration, setting
1972 // a value in this array to NaN and testing for it
1973 // is a way of obtaining the iteration counter.
1974 if (ido != 99)
1975 workl[iptr(5)-1] = octave::numeric_limits<double>::NaN ();
1976 }
1977
1978 if (ido == -1 || ido == 1 || ido == 2)
1979 {
1980 if (have_b)
1981 {
1982 Matrix mtmp (n, 1);
1983 for (F77_INT i = 0; i < n; i++)
1984 mtmp(i, 0) = workd[i + iptr(0) - 1];
1985
1986 mtmp = ltsolve (b, permB, m * utsolve (bt, permB, mtmp));
1987
1988 for (F77_INT i = 0; i < n; i++)
1989 workd[i+iptr(1)-1] = mtmp(i, 0);
1990 }
1991 else if (! vector_product (m, workd + iptr(0) - 1,
1992 workd + iptr(1) - 1))
1993 break;
1994 }
1995 else
1996 {
1997 if (info < 0)
1998 (*current_liboctave_error_handler)
1999 ("eigs: error in dnaupd: %s",
2000 arpack_errno2str (info, "dnaupd").c_str ());
2001
2002 break;
2003 }
2004 }
2005 while (1);
2006
2007 F77_INT info2;
2008
2009 // We have a problem in that the size of the C++ bool
2010 // type relative to the fortran logical type. It appears
2011 // that fortran uses 4- or 8-bytes per logical and C++ 1-byte
2012 // per bool, though this might be system dependent. As
2013 // long as the HOWMNY arg is not "S", the logical array
2014 // is just workspace for ARPACK, so use int type to
2015 // avoid problems.
2016 Array<F77_INT> s (dim_vector (p, 1));
2017 F77_INT *sel = s.rwdata ();
2018
2019 // FIXME: initialize eig_vec2 to zero; apparently dneupd can skip
2020 // the assignment to elements of Z that represent imaginary parts.
2021 // Found with valgrind and
2022 //
2023 // A = [1,0,0,-1;0,1,0,0;0,0,1,0;0,0,2,1];
2024 // [vecs, vals, f] = eigs (A, 1)
2025
2026 Matrix eig_vec2 (n, k + 1, 0.0);
2027 double *z = eig_vec2.rwdata ();
2028
2029 OCTAVE_LOCAL_BUFFER (double, dr, k + 1);
2030 OCTAVE_LOCAL_BUFFER (double, di, k + 1);
2031 OCTAVE_LOCAL_BUFFER (double, workev, 3 * p);
2032 for (F77_INT i = 0; i < k+1; i++)
2033 dr[i] = di[i] = 0.;
2034
2035 F77_INT k0 = k; // original number of eigenvalues required
2036 F77_FUNC (dneupd, DNEUPD)
2037 (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, dr, di, z, n, sigmar,
2038 sigmai, workev, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
2039 F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), k, tol, presid, p, v, n, iparam,
2040 ipntr, workd, workl, lwork, info2 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1)
2041 F77_CHAR_ARG_LEN(2));
2042 // on exit, if (and only if) rvec == true, k may have been increased by one
2043 // and be equal to ip(4), see dngets.
2044
2045 if (! rvec && ip(4) > k)
2046 k = ip(4);
2047
2048 eig_val.resize (k);
2049 Complex *d = eig_val.rwdata ();
2050
2051 if (info2 == 0)
2052 {
2053 bool have_cplx_eig = false;
2054 for (F77_INT i = 0; i < ip(4); i++)
2055 {
2056 if (di[i] == 0)
2057 d[i] = Complex (dr[i], 0.);
2058 else
2059 {
2060 have_cplx_eig = true;
2061 d[i] = Complex (dr[i], di[i]);
2062 }
2063 }
2064 if (have_cplx_eig)
2065 {
2066 for (F77_INT i = ip(4); i < k; i++)
2067 d[i] = Complex (octave::numeric_limits<double>::NaN (),
2068 octave::numeric_limits<double>::NaN ());
2069 }
2070 else
2071 {
2072 for (F77_INT i = ip(4); i < k; i++)
2073 d[i] = Complex (octave::numeric_limits<double>::NaN (), 0.);
2074 }
2075 if (! rvec)
2076 {
2077 // ARPACK seems to give the eigenvalues in reversed order
2078 F77_INT k2 = ip(4) / 2;
2079 for (F77_INT i = 0; i < k2; i++)
2080 {
2081 Complex dtmp = d[i];
2082 d[i] = d[ip(4) - i - 1];
2083 d[ip(4) - i - 1] = dtmp;
2084 }
2085 }
2086 else
2087 {
2088 // When eigenvectors required, ARPACK seems to give the right order
2089 eig_vec.resize (n, k);
2090 F77_INT i = 0;
2091 while (i < ip(4))
2092 {
2093 F77_INT off1 = i * n;
2094 F77_INT off2 = (i+1) * n;
2095 if (std::imag (eig_val(i)) == 0)
2096 {
2097 for (F77_INT j = 0; j < n; j++)
2098 eig_vec(j, i) = Complex (z[j+off1], 0.);
2099 i++;
2100 }
2101 else
2102 {
2103 for (F77_INT j = 0; j < n; j++)
2104 {
2105 eig_vec(j, i) = Complex (z[j+off1], z[j+off2]);
2106 if (i < ip(4) - 1)
2107 eig_vec(j, i+1) = Complex (z[j+off1], -z[j+off2]);
2108 }
2109 i+=2;
2110 }
2111 }
2112 if (have_cplx_eig)
2113 {
2114 for (F77_INT ii = ip(4); ii < k; ii++)
2115 for (F77_INT jj = 0; jj < n; jj++)
2116 eig_vec(jj, ii)
2117 = Complex (octave::numeric_limits<double>::NaN (),
2118 octave::numeric_limits<double>::NaN ());
2119 }
2120 else
2121 {
2122 for (F77_INT ii = ip(4); ii < k; ii++)
2123 for (F77_INT jj = 0; jj < n; jj++)
2124 eig_vec(jj, ii)
2125 = Complex (octave::numeric_limits<double>::NaN (), 0.);
2126 }
2127 if (note3)
2128 eig_vec = utsolve (bt, permB, eig_vec);
2129 }
2130 if (k0 < k)
2131 {
2132 eig_val.resize (k0);
2133 eig_vec.resize (n, k0);
2134 }
2135 }
2136 else
2137 (*current_liboctave_error_handler)
2138 ("eigs: error in dneupd: %s",
2139 arpack_errno2str (info2, "dneupd").c_str ());
2140
2141 return ip(4);
2142}
2143
2144template <typename M>
2145octave_idx_type
2146EigsRealNonSymmetricMatrixShift (const M& m, double sigmar,
2147 octave_idx_type k_arg, octave_idx_type p_arg,
2148 octave_idx_type& info,
2149 ComplexMatrix& eig_vec,
2150 ComplexColumnVector& eig_val, const M& _b,
2151 ColumnVector& permB, ColumnVector& resid,
2152 std::ostream& os, double tol, bool rvec,
2153 bool cholB, int disp, int maxit)
2154{
2155 F77_INT k = octave::to_f77_int (k_arg);
2156 F77_INT p = octave::to_f77_int (p_arg);
2157 M b(_b);
2158 F77_INT n = octave::to_f77_int (m.cols ());
2159 F77_INT mode = 3;
2160 bool have_b = ! b.isempty ();
2161 std::string typ = "LM";
2162 double sigmai = 0.;
2163
2164 if (m.rows () != m.cols ())
2165 (*current_liboctave_error_handler) ("eigs: A must be square");
2166 if (have_b && (m.rows () != b.rows () || m.rows () != b.cols ()))
2167 (*current_liboctave_error_handler)
2168 ("eigs: B must be square and the same size as A");
2169
2170 // FIXME: The "SM" type for mode 1 seems unstable though faster!!
2171 //if (! std::abs (sigmar))
2172 // return EigsRealNonSymmetricMatrix (m, "SM", k, p, info, eig_vec, eig_val,
2173 // _b, permB, resid, os, tol, rvec, cholB,
2174 // disp, maxit);
2175
2176 if (resid.isempty ())
2177 {
2178 std::string rand_dist = octave::rand::distribution ();
2179 octave::rand::distribution ("uniform");
2180 resid = ColumnVector (octave::rand::vector (n));
2181 octave::rand::distribution (rand_dist);
2182 }
2183 else if (m.cols () != resid.numel ())
2184 (*current_liboctave_error_handler) ("eigs: opts.v0 must be n-by-1");
2185
2186 if (n < 3)
2187 (*current_liboctave_error_handler) ("eigs: n must be at least 3");
2188
2189 if (p < 0)
2190 {
2191 p = k * 2 + 1;
2192
2193 if (p < 20)
2194 p = 20;
2195
2196 if (p > n)
2197 p = n;
2198 }
2199
2200 if (k <= 0 || k >= n - 1)
2201 (*current_liboctave_error_handler)
2202 ("eigs: Invalid number of eigenvalues to extract"
2203 " (must be 0 < k < n-1).\n"
2204 " Use 'eig (full (A))' instead");
2205
2206 if (p <= k || p > n)
2207 (*current_liboctave_error_handler)
2208 ("eigs: opts.p must be greater than k and less than or equal to n");
2209
2210 if (have_b && cholB && ! permB.isempty ())
2211 {
2212 // Check that we really have a permutation vector
2213 if (permB.numel () != n)
2214 (*current_liboctave_error_handler) ("eigs: permB vector invalid");
2215
2216 Array<bool> checked (dim_vector (n, 1), false);
2217 for (F77_INT i = 0; i < n; i++)
2218 {
2219 octave_idx_type bidx = static_cast<octave_idx_type> (permB(i));
2220
2221 if (checked(bidx) || bidx < 0 || bidx >= n)
2222 (*current_liboctave_error_handler) ("eigs: permB vector invalid");
2223 }
2224 }
2225
2226 char bmat = 'I';
2227 if (have_b)
2228 bmat = 'G';
2229
2230 Array<F77_INT> ip (dim_vector (11, 1));
2231 F77_INT *iparam = ip.rwdata ();
2232
2233 ip(0) = 1; //ishift
2234 ip(1) = 0; // ip(1) not referenced
2235 ip(2) = maxit; // mxiter, maximum number of iterations
2236 ip(3) = 1; // NB blocksize in recurrence
2237 ip(4) = 0; // nconv, number of Ritz values that satisfy convergence
2238 ip(5) = 0; //ip(5) not referenced
2239 ip(6) = mode; // mode
2240 ip(7) = 0;
2241 ip(8) = 0;
2242 ip(9) = 0;
2243 ip(10) = 0;
2244 // ip(7) to ip(10) return values
2245
2246 Array<F77_INT> iptr (dim_vector (14, 1));
2247 F77_INT *ipntr = iptr.rwdata ();
2248
2249 F77_INT ido = 0;
2250 int iter = 0;
2251 M L, U;
2252 ColumnVector r(n);
2253
2254 OCTAVE_LOCAL_BUFFER (octave_idx_type, P, (have_b ? b.rows () : m.rows ()));
2255 OCTAVE_LOCAL_BUFFER (octave_idx_type, Q, (have_b ? b.cols () : m.cols ()));
2256
2257 if (! LuAminusSigmaB (m, b, cholB, permB, sigmar, L, U, P, Q, r))
2258 return -1;
2259
2260 F77_INT elems;
2261 F77_INT lwork;
2262
2263 if (octave::math::int_multiply_overflow (n, p + 1, &elems))
2264 (*current_liboctave_error_handler)
2265 ("eigs: array too large for Fortran integers");
2266 if (octave::math::int_multiply_overflow (3 * p, p + 2, &lwork))
2267 (*current_liboctave_error_handler)
2268 ("eigs: array too large for Fortran integers");
2269
2270 OCTAVE_LOCAL_BUFFER (double, v, elems);
2271 OCTAVE_LOCAL_BUFFER (double, workl, lwork + 1);
2272 OCTAVE_LOCAL_BUFFER (double, workd, 3 * n + 1);
2273 double *presid = resid.rwdata ();
2274
2275 do
2276 {
2277 F77_INT tmp_info = octave::to_f77_int (info);
2278
2279 // On exit, ip(4) <= k + 1 is the number of converged eigenvalues.
2280 // See dnaupd2.
2281 F77_FUNC (dnaupd, DNAUPD)
2282 (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
2283 F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2),
2284 k, tol, presid, p, v, n, iparam,
2285 ipntr, workd, workl, lwork, tmp_info
2286 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
2287 // k is not changed
2288
2289 info = tmp_info;
2290
2291 if (disp > 0 && ! octave::math::isnan (workl[iptr (5)-1]))
2292 {
2293 if (iter++)
2294 {
2295 os << "Iteration " << iter - 1 <<
2296 ": a few Ritz values of the " << p << "-by-" <<
2297 p << " matrix\n";
2298 if (ido == 99) // convergence
2299 {
2300 os << " " << workl[iptr(5)+k] << "\n";
2301 for (F77_INT i = 0; i < k; i++)
2302 os << " " << workl[iptr(5)+i-1] << "\n";
2303 }
2304 else
2305 {
2306 // the wanted Ritz estimates are at the end
2307 for (F77_INT i = p - k - 1; i < p; i++)
2308 os << " " << workl[iptr(5)+i-1] << "\n";
2309 }
2310 }
2311
2312 // This is a kludge, as ARPACK doesn't give its
2313 // iteration pointer. But as workl[iptr(5)-1] is
2314 // an output value updated at each iteration, setting
2315 // a value in this array to NaN and testing for it
2316 // is a way of obtaining the iteration counter.
2317 if (ido != 99)
2318 workl[iptr(5)-1] = octave::numeric_limits<double>::NaN ();
2319 }
2320
2321 if (ido == -1 || ido == 1 || ido == 2)
2322 {
2323 if (have_b)
2324 {
2325 if (ido == -1)
2326 {
2327 OCTAVE_LOCAL_BUFFER (double, dtmp, n);
2328
2329 vector_product (b, workd+iptr(0)-1, dtmp);
2330
2331 Matrix tmp (n, 1);
2332
2333 for (F77_INT i = 0; i < n; i++)
2334 tmp(i, 0) = dtmp[P[i]] / r(P[i]);
2335
2336 lusolve (L, U, tmp);
2337
2338 double *ip2 = workd+iptr(1)-1;
2339 for (F77_INT i = 0; i < n; i++)
2340 ip2[Q[i]] = tmp(i, 0);
2341 }
2342 else if (ido == 2)
2343 vector_product (b, workd+iptr(0)-1, workd+iptr(1)-1);
2344 else
2345 {
2346 double *ip2 = workd+iptr(2)-1;
2347 Matrix tmp (n, 1);
2348
2349 for (F77_INT i = 0; i < n; i++)
2350 tmp(i, 0) = ip2[P[i]] / r(P[i]);
2351
2352 lusolve (L, U, tmp);
2353
2354 ip2 = workd+iptr(1)-1;
2355 for (F77_INT i = 0; i < n; i++)
2356 ip2[Q[i]] = tmp(i, 0);
2357 }
2358 }
2359 else
2360 {
2361 // ido cannot be 2 for non-generalized problems (see dnaupd2).
2362 double *ip2 = workd+iptr(0)-1;
2363 Matrix tmp (n, 1);
2364
2365 for (F77_INT i = 0; i < n; i++)
2366 tmp(i, 0) = ip2[P[i]] / r(P[i]);
2367
2368 lusolve (L, U, tmp);
2369
2370 ip2 = workd+iptr(1)-1;
2371 for (F77_INT i = 0; i < n; i++)
2372 ip2[Q[i]] = tmp(i, 0);
2373 }
2374 }
2375 else
2376 {
2377 if (info < 0)
2378 (*current_liboctave_error_handler)
2379 ("eigs: error in dnaupd: %s",
2380 arpack_errno2str (info, "dnaupd").c_str ());
2381
2382 break;
2383 }
2384 }
2385 while (1);
2386
2387 F77_INT info2;
2388
2389 // We have a problem in that the size of the C++ bool
2390 // type relative to the fortran logical type. It appears
2391 // that fortran uses 4- or 8-bytes per logical and C++ 1-byte
2392 // per bool, though this might be system dependent. As
2393 // long as the HOWMNY arg is not "S", the logical array
2394 // is just workspace for ARPACK, so use int type to
2395 // avoid problems.
2396 Array<F77_INT> s (dim_vector (p, 1));
2397 F77_INT *sel = s.rwdata ();
2398
2399 // FIXME: initialize eig_vec2 to zero; apparently dneupd can skip
2400 // the assignment to elements of Z that represent imaginary parts.
2401 // Found with valgrind and
2402 //
2403 // A = [1,0,0,-1;0,1,0,0;0,0,1,0;0,0,2,1];
2404 // [vecs, vals, f] = eigs (A, 1)
2405
2406 Matrix eig_vec2 (n, k + 1, 0.0);
2407 double *z = eig_vec2.rwdata ();
2408
2409 OCTAVE_LOCAL_BUFFER (double, dr, k + 1);
2410 OCTAVE_LOCAL_BUFFER (double, di, k + 1);
2411 OCTAVE_LOCAL_BUFFER (double, workev, 3 * p);
2412 for (F77_INT i = 0; i < k+1; i++)
2413 dr[i] = di[i] = 0.;
2414
2415 F77_INT k0 = k; // original number of eigenvalues required
2416 F77_FUNC (dneupd, DNEUPD)
2417 (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, dr, di, z, n, sigmar,
2418 sigmai, workev, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
2419 F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), k, tol, presid, p, v, n, iparam,
2420 ipntr, workd, workl, lwork, info2 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1)
2421 F77_CHAR_ARG_LEN(2));
2422 // On exit, if (and only if) rvec == true, k may have been increased by one
2423 // and be equal to ip(4), see dngets.
2424
2425 if (! rvec && ip(4) > k)
2426 k = ip(4);
2427
2428 eig_val.resize (k);
2429 Complex *d = eig_val.rwdata ();
2430
2431 if (info2 == 0)
2432 {
2433 bool have_cplx_eig = false;
2434 for (F77_INT i = 0; i < ip(4); i++)
2435 {
2436 if (di[i] == 0.)
2437 d[i] = Complex (dr[i], 0.);
2438 else
2439 {
2440 have_cplx_eig = true;
2441 d[i] = Complex (dr[i], di[i]);
2442 }
2443 }
2444 if (have_cplx_eig)
2445 {
2446 for (F77_INT i = ip(4); i < k; i++)
2447 d[i] = Complex (octave::numeric_limits<double>::NaN (),
2448 octave::numeric_limits<double>::NaN ());
2449 }
2450 else
2451 {
2452 for (F77_INT i = ip(4); i < k; i++)
2453 d[i] = Complex (octave::numeric_limits<double>::NaN (), 0.);
2454 }
2455
2456 if (! rvec)
2457 {
2458 // ARPACK seems to give the eigenvalues in reversed order
2459 F77_INT k2 = ip(4) / 2;
2460 for (F77_INT i = 0; i < k2; i++)
2461 {
2462 Complex dtmp = d[i];
2463 d[i] = d[ip(4) - i - 1];
2464 d[ip(4) - i - 1] = dtmp;
2465 }
2466 }
2467 else
2468 {
2469 // When eigenvectors required, ARPACK seems to give the right order
2470 eig_vec.resize (n, k);
2471 F77_INT i = 0;
2472 while (i < ip(4))
2473 {
2474 F77_INT off1 = i * n;
2475 F77_INT off2 = (i+1) * n;
2476 if (std::imag (eig_val(i)) == 0)
2477 {
2478 for (F77_INT j = 0; j < n; j++)
2479 eig_vec(j, i) = Complex (z[j+off1], 0.);
2480 i++;
2481 }
2482 else
2483 {
2484 for (F77_INT j = 0; j < n; j++)
2485 {
2486 eig_vec(j, i) = Complex (z[j+off1], z[j+off2]);
2487 if (i < ip(4) - 1)
2488 eig_vec(j, i+1) = Complex (z[j+off1], -z[j+off2]);
2489 }
2490 i+=2;
2491 }
2492 }
2493 if (have_cplx_eig)
2494 {
2495 for (F77_INT ii = ip(4); ii < k; ii++)
2496 for (F77_INT jj = 0; jj < n; jj++)
2497 eig_vec(jj, ii)
2498 = Complex (octave::numeric_limits<double>::NaN (),
2499 octave::numeric_limits<double>::NaN ());
2500 }
2501 else
2502 {
2503 for (F77_INT ii = ip(4); ii < k; ii++)
2504 for (F77_INT jj = 0; jj < n; jj++)
2505 eig_vec(jj, ii)
2506 = Complex (octave::numeric_limits<double>::NaN (), 0.);
2507 }
2508 }
2509 if (k0 < k)
2510 {
2511 eig_val.resize (k0);
2512 eig_vec.resize (n, k0);
2513 }
2514 }
2515 else
2516 (*current_liboctave_error_handler)
2517 ("eigs: error in dneupd: %s",
2518 arpack_errno2str (info2, "dneupd").c_str ());
2519
2520 return ip(4);
2521}
2522
2523template <typename M>
2524octave_idx_type
2525EigsRealNonSymmetricFunc (EigsFunc fcn, octave_idx_type n_arg,
2526 const std::string& _typ, double sigmar,
2527 octave_idx_type k_arg, octave_idx_type p_arg,
2528 octave_idx_type& info, ComplexMatrix& eig_vec,
2529 ComplexColumnVector& eig_val, const M& _b,
2530 ColumnVector& permB, ColumnVector& resid,
2531 std::ostream& os, double tol, bool rvec,
2532 bool cholB, int disp, int maxit)
2533{
2534 F77_INT n = octave::to_f77_int (n_arg);
2535 F77_INT k = octave::to_f77_int (k_arg);
2536 F77_INT p = octave::to_f77_int (p_arg);
2537 M b(_b);
2538 std::string typ (_typ);
2539 bool have_sigma = (sigmar ? true : false);
2540 double sigmai = 0.;
2541 F77_INT mode = 1;
2542 bool have_b = ! b.isempty ();
2543 bool note3 = false;
2544 char bmat = 'I';
2545 int err = 0;
2546 M bt;
2547
2548 if (resid.isempty ())
2549 {
2550 std::string rand_dist = octave::rand::distribution ();
2551 octave::rand::distribution ("uniform");
2552 resid = ColumnVector (octave::rand::vector (n));
2553 octave::rand::distribution (rand_dist);
2554 }
2555 else if (n != resid.numel ())
2556 (*current_liboctave_error_handler) ("eigs: opts.v0 must be n-by-1");
2557
2558 if (n < 3)
2559 (*current_liboctave_error_handler) ("eigs: n must be at least 3");
2560
2561 if (p < 0)
2562 {
2563 p = k * 2 + 1;
2564
2565 if (p < 20)
2566 p = 20;
2567
2568 if (p > n)
2569 p = n;
2570 }
2571
2572 if (k <= 0 || k >= n - 1)
2573 (*current_liboctave_error_handler)
2574 ("eigs: Invalid number of eigenvalues to extract"
2575 " (must be 0 < k < n-1).\n"
2576 " Use 'eig (full (A))' instead");
2577
2578 if (p <= k || p > n)
2579 (*current_liboctave_error_handler)
2580 ("eigs: opts.p must be greater than k and less than or equal to n");
2581
2582 if (have_b && cholB && ! permB.isempty ())
2583 {
2584 // Check the we really have a permutation vector
2585 if (permB.numel () != n)
2586 (*current_liboctave_error_handler) ("eigs: permB vector invalid");
2587
2588 Array<bool> checked (dim_vector (n, 1), false);
2589 for (F77_INT i = 0; i < n; i++)
2590 {
2591 octave_idx_type bidx = static_cast<octave_idx_type> (permB(i));
2592
2593 if (checked(bidx) || bidx < 0 || bidx >= n)
2594 (*current_liboctave_error_handler) ("eigs: permB vector invalid");
2595 }
2596 }
2597
2598 if (! have_sigma)
2599 {
2600 if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA"
2601 && typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI"
2602 && typ != "SI")
2603 (*current_liboctave_error_handler) ("eigs: unrecognized sigma value");
2604
2605 if (typ == "LA" || typ == "SA" || typ == "BE")
2606 (*current_liboctave_error_handler)
2607 ("eigs: invalid sigma value for unsymmetric problem");
2608
2609 if (typ != "SM" && have_b)
2610 note3 = true;
2611
2612 if (typ == "SM")
2613 {
2614 typ = "LM";
2615 sigmar = 0.;
2616 mode = 3;
2617 if (have_b)
2618 bmat = 'G';
2619 }
2620 }
2621 else if (! std::abs (sigmar))
2622 {
2623 typ = "SM";
2624 if (have_b)
2625 bmat = 'G';
2626 }
2627 else
2628 {
2629 typ = "LM";
2630 mode = 3;
2631 if (have_b)
2632 bmat = 'G';
2633 }
2634
2635 if (mode == 1 && have_b)
2636 {
2637 // See Note 3 dsaupd
2638 note3 = true;
2639 if (cholB)
2640 {
2641 bt = b;
2642 b = b.transpose ();
2643 if (permB.isempty ())
2644 {
2645 permB = ColumnVector (n);
2646 for (F77_INT i = 0; i < n; i++)
2647 permB(i) = i;
2648 }
2649 }
2650 else
2651 {
2652 if (! make_cholb (b, bt, permB))
2653 (*current_liboctave_error_handler)
2654 ("eigs: The matrix B is not positive definite");
2655 }
2656 }
2657
2658 Array<F77_INT> ip (dim_vector (11, 1));
2659 F77_INT *iparam = ip.rwdata ();
2660
2661 ip(0) = 1; //ishift
2662 ip(1) = 0; // ip(1) not referenced
2663 ip(2) = maxit; // mxiter, maximum number of iterations
2664 ip(3) = 1; // NB blocksize in recurrence
2665 ip(4) = 0; // nconv, number of Ritz values that satisfy convergence
2666 ip(5) = 0; //ip(5) not referenced
2667 ip(6) = mode; // mode
2668 ip(7) = 0;
2669 ip(8) = 0;
2670 ip(9) = 0;
2671 ip(10) = 0;
2672 // ip(7) to ip(10) return values
2673
2674 Array<F77_INT> iptr (dim_vector (14, 1));
2675 F77_INT *ipntr = iptr.rwdata ();
2676
2677 F77_INT ido = 0;
2678 int iter = 0;
2679 F77_INT elems;
2680 F77_INT lwork;
2681
2682 if (octave::math::int_multiply_overflow (n, p + 1, &elems))
2683 (*current_liboctave_error_handler)
2684 ("eigs: array too large for Fortran integers");
2685 if (octave::math::int_multiply_overflow (3 * p, p + 2, &lwork))
2686 (*current_liboctave_error_handler)
2687 ("eigs: array too large for Fortran integers");
2688
2689 OCTAVE_LOCAL_BUFFER (double, v, elems);
2690 OCTAVE_LOCAL_BUFFER (double, workl, lwork + 1);
2691 OCTAVE_LOCAL_BUFFER (double, workd, 3 * n + 1);
2692
2693 double *presid = resid.rwdata ();
2694
2695 do
2696 {
2697 F77_INT tmp_info = octave::to_f77_int (info);
2698
2699 // On exit, ip(4) <= k + 1 is the number of converged eigenvalues
2700 // see dnaupd2.
2701 F77_FUNC (dnaupd, DNAUPD)
2702 (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
2703 F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2),
2704 k, tol, presid, p, v, n, iparam,
2705 ipntr, workd, workl, lwork, tmp_info
2706 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
2707 // k is not changed
2708
2709 info = tmp_info;
2710
2711 if (disp > 0 && ! octave::math::isnan(workl[iptr(5)-1]))
2712 {
2713 if (iter++)
2714 {
2715 os << "Iteration " << iter - 1 <<
2716 ": a few Ritz values of the " << p << "-by-" <<
2717 p << " matrix\n";
2718 if (ido == 99) // convergence
2719 {
2720 os << " " << workl[iptr(5)+k] << "\n";
2721 for (F77_INT i = 0; i < k; i++)
2722 os << " " << workl[iptr(5)+i-1] << "\n";
2723 }
2724 else
2725 {
2726 // the wanted Ritz estimates are at the end
2727 for (F77_INT i = p - k - 1; i < p; i++)
2728 os << " " << workl[iptr(5)+i-1] << "\n";
2729 }
2730 }
2731
2732 // This is a kludge, as ARPACK doesn't give its
2733 // iteration pointer. But as workl[iptr(5)-1] is
2734 // an output value updated at each iteration, setting
2735 // a value in this array to NaN and testing for it
2736 // is a way of obtaining the iteration counter.
2737 if (ido != 99)
2738 workl[iptr(5)-1] = octave::numeric_limits<double>::NaN ();
2739 }
2740
2741 if (ido == -1 || ido == 1 || ido == 2)
2742 {
2743 if (have_b)
2744 {
2745 if (mode == 1) // regular mode with factorized B
2746 {
2747 Matrix mtmp (n, 1);
2748 for (F77_INT i = 0; i < n; i++)
2749 mtmp(i, 0) = workd[i + iptr(0) - 1];
2750
2751 mtmp = utsolve (bt, permB, mtmp);
2752 ColumnVector y = fcn (mtmp, err);
2753
2754 if (err)
2755 return false;
2756
2757 mtmp = ltsolve (b, permB, y);
2758
2759 for (F77_INT i = 0; i < n; i++)
2760 workd[i+iptr(1)-1] = mtmp(i, 0);
2761 }
2762 else // shift-invert mode
2763 {
2764 if (ido == -1)
2765 {
2766 OCTAVE_LOCAL_BUFFER (double, dtmp, n);
2767
2768 vector_product (b, workd+iptr(0)-1, dtmp);
2769
2770 ColumnVector x(n);
2771
2772 for (F77_INT i = 0; i < n; i++)
2773 x(i) = dtmp[i];
2774
2775 ColumnVector y = fcn (x, err);
2776
2777 if (err)
2778 return false;
2779
2780 double *ip2 = workd + iptr(1) - 1;
2781 for (F77_INT i = 0; i < n; i++)
2782 ip2[i] = y(i);
2783 }
2784 else if (ido == 2)
2785 vector_product (b, workd+iptr(0)-1, workd+iptr(1)-1);
2786 else
2787 {
2788 double *ip2 = workd+iptr(2)-1;
2789 ColumnVector x(n);
2790
2791 for (F77_INT i = 0; i < n; i++)
2792 x(i) = *ip2++;
2793
2794 ColumnVector y = fcn (x, err);
2795
2796 if (err)
2797 return false;
2798
2799 ip2 = workd + iptr(1) - 1;
2800 for (F77_INT i = 0; i < n; i++)
2801 *ip2++ = y(i);
2802 }
2803 }
2804 }
2805 else
2806 {
2807 double *ip2 = workd + iptr(0) - 1;
2808 ColumnVector x(n);
2809
2810 for (F77_INT i = 0; i < n; i++)
2811 x(i) = *ip2++;
2812
2813 ColumnVector y = fcn (x, err);
2814
2815 if (err)
2816 return false;
2817
2818 ip2 = workd + iptr(1) - 1;
2819 for (F77_INT i = 0; i < n; i++)
2820 *ip2++ = y(i);
2821 }
2822 }
2823 else
2824 {
2825 if (info < 0)
2826 (*current_liboctave_error_handler)
2827 ("eigs: error in dnaupd: %s",
2828 arpack_errno2str (info, "dnaupd").c_str ());
2829
2830 break;
2831 }
2832 }
2833 while (1);
2834
2835 F77_INT info2;
2836
2837 // We have a problem in that the size of the C++ bool
2838 // type relative to the fortran logical type. It appears
2839 // that fortran uses 4- or 8-bytes per logical and C++ 1-byte
2840 // per bool, though this might be system dependent. As
2841 // long as the HOWMNY arg is not "S", the logical array
2842 // is just workspace for ARPACK, so use int type to
2843 // avoid problems.
2844 Array<F77_INT> s (dim_vector (p, 1));
2845 F77_INT *sel = s.rwdata ();
2846
2847 // FIXME: initialize eig_vec2 to zero; apparently dneupd can skip
2848 // the assignment to elements of Z that represent imaginary parts.
2849 // Found with valgrind and
2850 //
2851 // A = [1,0,0,-1;0,1,0,0;0,0,1,0;0,0,2,1];
2852 // [vecs, vals, f] = eigs (A, 1)
2853
2854 Matrix eig_vec2 (n, k + 1, 0.0);
2855 double *z = eig_vec2.rwdata ();
2856
2857 OCTAVE_LOCAL_BUFFER (double, dr, k + 1);
2858 OCTAVE_LOCAL_BUFFER (double, di, k + 1);
2859 OCTAVE_LOCAL_BUFFER (double, workev, 3 * p);
2860 for (F77_INT i = 0; i < k+1; i++)
2861 dr[i] = di[i] = 0.;
2862
2863 F77_INT k0 = k; // original number of eigenvalues required
2864 F77_FUNC (dneupd, DNEUPD)
2865 (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, dr, di, z, n, sigmar,
2866 sigmai, workev, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
2867 F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), k, tol, presid, p, v, n, iparam,
2868 ipntr, workd, workl, lwork, info2 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1)
2869 F77_CHAR_ARG_LEN(2));
2870 // On exit, if (and only if) rvec == true, k may have been increased by one
2871 // and be equal to ip(4), see dngets.
2872
2873 if (! rvec && ip(4) > k)
2874 k = ip(4);
2875
2876 eig_val.resize (k);
2877 Complex *d = eig_val.rwdata ();
2878
2879 if (info2 == 0)
2880 {
2881 bool have_cplx_eig = false;
2882 for (F77_INT i = 0; i < ip(4); i++)
2883 {
2884 if (di[i] == 0.)
2885 d[i] = Complex (dr[i], 0.);
2886 else
2887 {
2888 have_cplx_eig = true;
2889 d[i] = Complex (dr[i], di[i]);
2890 }
2891 }
2892 if (have_cplx_eig)
2893 {
2894 for (F77_INT i = ip(4); i < k; i++)
2895 d[i] = Complex (octave::numeric_limits<double>::NaN (),
2896 octave::numeric_limits<double>::NaN ());
2897 }
2898 else
2899 {
2900 for (F77_INT i = ip(4); i < k; i++)
2901 d[i] = Complex (octave::numeric_limits<double>::NaN (), 0.);
2902 }
2903
2904 if (! rvec)
2905 {
2906 // ARPACK seems to give the eigenvalues in reversed order
2907 octave_idx_type k2 = ip(4) / 2;
2908 for (F77_INT i = 0; i < k2; i++)
2909 {
2910 Complex dtmp = d[i];
2911 d[i] = d[ip(4) - i - 1];
2912 d[ip(4) - i - 1] = dtmp;
2913 }
2914 }
2915 else
2916 {
2917 // ARPACK seems to give the eigenvalues in reversed order
2918 eig_vec.resize (n, k);
2919 F77_INT i = 0;
2920 while (i < ip(4))
2921 {
2922 F77_INT off1 = i * n;
2923 F77_INT off2 = (i+1) * n;
2924 if (std::imag (eig_val(i)) == 0)
2925 {
2926 for (F77_INT j = 0; j < n; j++)
2927 eig_vec(j, i) = Complex (z[j+off1], 0.);
2928 i++;
2929 }
2930 else
2931 {
2932 for (F77_INT j = 0; j < n; j++)
2933 {
2934 eig_vec(j, i) = Complex (z[j+off1], z[j+off2]);
2935 if (i < ip(4) - 1)
2936 eig_vec(j, i+1) = Complex (z[j+off1], -z[j+off2]);
2937 }
2938 i+=2;
2939 }
2940 }
2941 if (have_cplx_eig)
2942 {
2943 for (F77_INT ii = ip(4); ii < k; ii++)
2944 for (F77_INT jj = 0; jj < n; jj++)
2945 eig_vec(jj, ii)
2946 = Complex (octave::numeric_limits<double>::NaN (),
2947 octave::numeric_limits<double>::NaN ());
2948 }
2949 else
2950 {
2951 for (F77_INT ii = ip(4); ii < k; ii++)
2952 for (F77_INT jj = 0; jj < n; jj++)
2953 eig_vec(jj, ii)
2954 = Complex (octave::numeric_limits<double>::NaN (), 0.);
2955 }
2956 if (note3)
2957 eig_vec = utsolve (bt, permB, eig_vec);
2958 }
2959 if (k0 < k)
2960 {
2961 eig_val.resize (k0);
2962 eig_vec.resize (n, k0);
2963 }
2964 }
2965 else
2966 (*current_liboctave_error_handler)
2967 ("eigs: error in dneupd: %s",
2968 arpack_errno2str (info2, "dneupd").c_str ());
2969
2970 return ip(4);
2971}
2972
2973template <typename M>
2974octave_idx_type
2975EigsComplexNonSymmetricMatrix (const M& m, const std::string typ,
2976 octave_idx_type k_arg, octave_idx_type p_arg,
2977 octave_idx_type& info, ComplexMatrix& eig_vec,
2978 ComplexColumnVector& eig_val, const M& _b,
2979 ColumnVector& permB,
2980 ComplexColumnVector& cresid,
2981 std::ostream& os, double tol, bool rvec,
2982 bool cholB, int disp, int maxit)
2983{
2984 F77_INT k = octave::to_f77_int (k_arg);
2985 F77_INT p = octave::to_f77_int (p_arg);
2986 M b(_b);
2987 F77_INT n = octave::to_f77_int (m.cols ());
2988 F77_INT mode = 1;
2989 bool have_b = ! b.isempty ();
2990 bool note3 = false;
2991 char bmat = 'I';
2992 Complex sigma = 0.;
2993 M bt;
2994
2995 if (m.rows () != m.cols ())
2996 (*current_liboctave_error_handler) ("eigs: A must be square");
2997 if (have_b && (m.rows () != b.rows () || m.rows () != b.cols ()))
2998 (*current_liboctave_error_handler)
2999 ("eigs: B must be square and the same size as A");
3000
3001 if (cresid.isempty ())
3002 {
3003 std::string rand_dist = octave::rand::distribution ();
3004 octave::rand::distribution ("uniform");
3005 Array<double> rr (octave::rand::vector (n));
3006 Array<double> ri (octave::rand::vector (n));
3007 cresid = ComplexColumnVector (n);
3008 for (F77_INT i = 0; i < n; i++)
3009 cresid(i) = Complex (rr(i), ri(i));
3010 octave::rand::distribution (rand_dist);
3011 }
3012 else if (m.cols () != cresid.numel ())
3013 (*current_liboctave_error_handler) ("eigs: opts.v0 must be n-by-1");
3014
3015 if (n < 3)
3016 (*current_liboctave_error_handler) ("eigs: n must be at least 3");
3017
3018 if (p < 0)
3019 {
3020 p = k * 2 + 1;
3021
3022 if (p < 20)
3023 p = 20;
3024
3025 if (p > n)
3026 p = n;
3027 }
3028
3029 if (k <= 0 || k >= n - 1)
3030 (*current_liboctave_error_handler)
3031 ("eigs: Invalid number of eigenvalues to extract"
3032 " (must be 0 < k < n-1).\n"
3033 " Use 'eig (full (A))' instead");
3034
3035 if (p <= k || p > n)
3036 (*current_liboctave_error_handler)
3037 ("eigs: opts.p must be greater than k and less than or equal to n");
3038
3039 if (have_b && cholB && ! permB.isempty ())
3040 {
3041 // Check the we really have a permutation vector
3042 if (permB.numel () != n)
3043 (*current_liboctave_error_handler) ("eigs: permB vector invalid");
3044
3045 Array<bool> checked (dim_vector (n, 1), false);
3046 for (F77_INT i = 0; i < n; i++)
3047 {
3048 octave_idx_type bidx = static_cast<octave_idx_type> (permB(i));
3049
3050 if (checked(bidx) || bidx < 0 || bidx >= n)
3051 (*current_liboctave_error_handler) ("eigs: permB vector invalid");
3052 }
3053 }
3054
3055 if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA"
3056 && typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI"
3057 && typ != "SI")
3058 (*current_liboctave_error_handler) ("eigs: unrecognized sigma value");
3059
3060 if (typ == "LA" || typ == "SA" || typ == "BE")
3061 (*current_liboctave_error_handler)
3062 ("eigs: invalid sigma value for complex problem");
3063
3064 if (have_b)
3065 {
3066 // See Note 3 dsaupd
3067 note3 = true;
3068 if (cholB)
3069 {
3070 bt = b;
3071 b = b.hermitian ();
3072 if (permB.isempty ())
3073 {
3074 permB = ColumnVector (n);
3075 for (F77_INT i = 0; i < n; i++)
3076 permB(i) = i;
3077 }
3078 }
3079 else
3080 {
3081 if (! make_cholb (b, bt, permB))
3082 (*current_liboctave_error_handler)
3083 ("eigs: The matrix B is not positive definite");
3084 }
3085 }
3086
3087 Array<F77_INT> ip (dim_vector (11, 1));
3088 F77_INT *iparam = ip.rwdata ();
3089
3090 ip(0) = 1; //ishift
3091 ip(1) = 0; // ip(1) not referenced
3092 ip(2) = maxit; // mxiter, maximum number of iterations
3093 ip(3) = 1; // NB blocksize in recurrence
3094 ip(4) = 0; // nconv, number of Ritz values that satisfy convergence
3095 ip(5) = 0; //ip(5) not referenced
3096 ip(6) = mode; // mode
3097 ip(7) = 0;
3098 ip(8) = 0;
3099 ip(9) = 0;
3100 ip(10) = 0;
3101 // ip(7) to ip(10) return values
3102
3103 Array<F77_INT> iptr (dim_vector (14, 1));
3104 F77_INT *ipntr = iptr.rwdata ();
3105
3106 F77_INT ido = 0;
3107 int iter = 0;
3108 F77_INT elems;
3109 F77_INT lwork;
3110
3111 if (octave::math::int_multiply_overflow (n, p, &elems))
3112 (*current_liboctave_error_handler)
3113 ("eigs: array too large for Fortran integers");
3114 if (octave::math::int_multiply_overflow (p, 3 * p + 5, &lwork))
3115 (*current_liboctave_error_handler)
3116 ("eigs: array too large for Fortran integers");
3117
3118 OCTAVE_LOCAL_BUFFER (Complex, v, elems);
3119 OCTAVE_LOCAL_BUFFER (Complex, workl, lwork);
3120 OCTAVE_LOCAL_BUFFER (Complex, workd, 3 * n);
3121 OCTAVE_LOCAL_BUFFER (double, rwork, p);
3122 Complex *presid = cresid.rwdata ();
3123
3124 do
3125 {
3126 F77_INT tmp_info = octave::to_f77_int (info);
3127
3128 F77_FUNC (znaupd, ZNAUPD)
3129 (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
3130 F77_CONST_CHAR_ARG2 (typ.c_str (), 2),
3131 k, tol, F77_DBLE_CMPLX_ARG (presid), p, F77_DBLE_CMPLX_ARG (v), n,
3132 iparam, ipntr,
3133 F77_DBLE_CMPLX_ARG (workd), F77_DBLE_CMPLX_ARG (workl), lwork, rwork,
3134 tmp_info F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
3135
3136 info = tmp_info;
3137
3138 if (disp > 0 && ! octave::math::isnan (workl[iptr (5)-1]))
3139 {
3140 if (iter++)
3141 {
3142 os << "Iteration " << iter - 1 <<
3143 ": a few Ritz values of the " << p << "-by-" <<
3144 p << " matrix\n";
3145 if (ido == 99) // convergence
3146 {
3147 for (F77_INT i = 0; i < k; i++)
3148 os << " " << workl[iptr(5)+i-1] << "\n";
3149 }
3150 else
3151 {
3152 // the wanted Ritz estimates are at the end
3153 for (F77_INT i = p - k; i < p; i++)
3154 os << " " << workl[iptr(5)+i-1] << "\n";
3155 }
3156 }
3157
3158 // This is a kludge, as ARPACK doesn't give its
3159 // iteration pointer. But as workl[iptr(5)-1] is
3160 // an output value updated at each iteration, setting
3161 // a value in this array to NaN and testing for it
3162 // is a way of obtaining the iteration counter.
3163 if (ido != 99)
3164 workl[iptr(5)-1] = octave::numeric_limits<double>::NaN ();
3165 }
3166
3167 if (ido == -1 || ido == 1 || ido == 2)
3168 {
3169 if (have_b)
3170 {
3171 ComplexMatrix mtmp (n, 1);
3172 for (F77_INT i = 0; i < n; i++)
3173 mtmp(i, 0) = workd[i + iptr(0) - 1];
3174 mtmp = ltsolve (b, permB, m * utsolve (bt, permB, mtmp));
3175 for (F77_INT i = 0; i < n; i++)
3176 workd[i+iptr(1)-1] = mtmp(i, 0);
3177
3178 }
3179 else if (! vector_product (m, workd + iptr(0) - 1,
3180 workd + iptr(1) - 1))
3181 break;
3182 }
3183 else
3184 {
3185 if (info < 0)
3186 (*current_liboctave_error_handler)
3187 ("eigs: error %" OCTAVE_IDX_TYPE_FORMAT " in znaupd: %s",
3188 info, arpack_errno2str (info, "znaupd").c_str ());
3189
3190 break;
3191 }
3192 }
3193 while (1);
3194
3195 F77_INT info2;
3196
3197 // We have a problem in that the size of the C++ bool
3198 // type relative to the fortran logical type. It appears
3199 // that fortran uses 4- or 8-bytes per logical and C++ 1-byte
3200 // per bool, though this might be system dependent. As
3201 // long as the HOWMNY arg is not "S", the logical array
3202 // is just workspace for ARPACK, so use int type to
3203 // avoid problems.
3204 Array<F77_INT> s (dim_vector (p, 1));
3205 F77_INT *sel = s.rwdata ();
3206
3207 eig_vec.resize (n, k);
3208 Complex *z = eig_vec.rwdata ();
3209
3210 eig_val.resize (k+1);
3211 Complex *d = eig_val.rwdata ();
3212
3213 OCTAVE_LOCAL_BUFFER (Complex, workev, 2 * p);
3214
3215 F77_FUNC (zneupd, ZNEUPD)
3216 (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, F77_DBLE_CMPLX_ARG (d),
3217 F77_DBLE_CMPLX_ARG (z), n, F77_CONST_DBLE_CMPLX_ARG (&sigma),
3218 F77_DBLE_CMPLX_ARG (workev),
3219 F77_CONST_CHAR_ARG2 (&bmat, 1), n,
3220 F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2),
3221 k, tol, F77_DBLE_CMPLX_ARG (presid), p, F77_DBLE_CMPLX_ARG (v), n,
3222 iparam, ipntr,
3223 F77_DBLE_CMPLX_ARG (workd), F77_DBLE_CMPLX_ARG (workl), lwork, rwork,
3224 info2 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
3225
3226 if (info2 == 0)
3227 {
3228 for (F77_INT i = ip(4); i < k; i++)
3229 d[i] = Complex (octave::numeric_limits<double>::NaN (),
3230 octave::numeric_limits<double>::NaN ());
3231
3232 F77_INT k2 = ip(4) / 2;
3233 for (F77_INT i = 0; i < k2; i++)
3234 {
3235 Complex ctmp = d[i];
3236 d[i] = d[ip(4) - i - 1];
3237 d[ip(4) - i - 1] = ctmp;
3238 }
3239 eig_val.resize (k);
3240
3241 if (rvec)
3242 {
3243 OCTAVE_LOCAL_BUFFER (Complex, ctmp, n);
3244
3245 for (F77_INT i = ip(4); i < k; i++)
3246 {
3247 F77_INT off1 = i * n;
3248 for (F77_INT j = 0; j < n; j++)
3249 z[off1 + j] = Complex (octave::numeric_limits<double>::NaN (),
3250 octave::numeric_limits<double>::NaN ());
3251 }
3252
3253 for (F77_INT i = 0; i < k2; i++)
3254 {
3255 F77_INT off1 = i * n;
3256 F77_INT off2 = (ip(4) - i - 1) * n;
3257
3258 if (off1 == off2)
3259 continue;
3260
3261 for (F77_INT j = 0; j < n; j++)
3262 ctmp[j] = z[off1 + j];
3263
3264 for (F77_INT j = 0; j < n; j++)
3265 z[off1 + j] = z[off2 + j];
3266
3267 for (F77_INT j = 0; j < n; j++)
3268 z[off2 + j] = ctmp[j];
3269 }
3270
3271 if (note3)
3272 eig_vec = utsolve (bt, permB, eig_vec);
3273 }
3274 }
3275 else
3276 (*current_liboctave_error_handler)
3277 ("eigs: error in zneupd: %s",
3278 arpack_errno2str (info2, "zneupd").c_str ());
3279
3280 return ip(4);
3281}
3282
3283template <typename M>
3284octave_idx_type
3285EigsComplexNonSymmetricMatrixShift (const M& m, Complex sigma,
3286 octave_idx_type k_arg, octave_idx_type p_arg,
3287 octave_idx_type& info,
3288 ComplexMatrix& eig_vec,
3289 ComplexColumnVector& eig_val, const M& _b,
3290 ColumnVector& permB,
3291 ComplexColumnVector& cresid,
3292 std::ostream& os, double tol, bool rvec,
3293 bool cholB, int disp, int maxit)
3294{
3295 F77_INT k = octave::to_f77_int (k_arg);
3296 F77_INT p = octave::to_f77_int (p_arg);
3297 M b(_b);
3298 F77_INT n = octave::to_f77_int (m.cols ());
3299 F77_INT mode = 3;
3300 bool have_b = ! b.isempty ();
3301 std::string typ = "LM";
3302
3303 if (m.rows () != m.cols ())
3304 (*current_liboctave_error_handler) ("eigs: A must be square");
3305 if (have_b && (m.rows () != b.rows () || m.rows () != b.cols ()))
3306 (*current_liboctave_error_handler)
3307 ("eigs: B must be square and the same size as A");
3308
3309 // FIXME: The "SM" type for mode 1 seems unstable though faster!!
3310 //if (! std::abs (sigma))
3311 // return EigsComplexNonSymmetricMatrix (m, "SM", k, p, info, eig_vec,
3312 // eig_val, _b, permB, cresid, os, tol,
3313 // rvec, cholB, disp, maxit);
3314
3315 if (cresid.isempty ())
3316 {
3317 std::string rand_dist = octave::rand::distribution ();
3318 octave::rand::distribution ("uniform");
3319 Array<double> rr (octave::rand::vector (n));
3320 Array<double> ri (octave::rand::vector (n));
3321 cresid = ComplexColumnVector (n);
3322 for (F77_INT i = 0; i < n; i++)
3323 cresid(i) = Complex (rr(i), ri(i));
3324 octave::rand::distribution (rand_dist);
3325 }
3326 else if (m.cols () != cresid.numel ())
3327 (*current_liboctave_error_handler) ("eigs: opts.v0 must be n-by-1");
3328
3329 if (n < 3)
3330 (*current_liboctave_error_handler) ("eigs: n must be at least 3");
3331
3332 if (p < 0)
3333 {
3334 p = k * 2 + 1;
3335
3336 if (p < 20)
3337 p = 20;
3338
3339 if (p > n)
3340 p = n;
3341 }
3342
3343 if (k <= 0 || k >= n - 1)
3344 (*current_liboctave_error_handler)
3345 ("eigs: Invalid number of eigenvalues to extract"
3346 " (must be 0 < k < n-1).\n"
3347 " Use 'eig (full (A))' instead");
3348
3349 if (p <= k || p > n)
3350 (*current_liboctave_error_handler)
3351 ("eigs: opts.p must be greater than k and less than or equal to n");
3352
3353 if (have_b && cholB && ! permB.isempty ())
3354 {
3355 // Check that we really have a permutation vector
3356 if (permB.numel () != n)
3357 (*current_liboctave_error_handler) ("eigs: permB vector invalid");
3358
3359 Array<bool> checked (dim_vector (n, 1), false);
3360 for (F77_INT i = 0; i < n; i++)
3361 {
3362 octave_idx_type bidx = static_cast<octave_idx_type> (permB(i));
3363
3364 if (checked(bidx) || bidx < 0 || bidx >= n)
3365 (*current_liboctave_error_handler) ("eigs: permB vector invalid");
3366 }
3367 }
3368
3369 char bmat = 'I';
3370 if (have_b)
3371 bmat = 'G';
3372
3373 Array<F77_INT> ip (dim_vector (11, 1));
3374 F77_INT *iparam = ip.rwdata ();
3375
3376 ip(0) = 1; //ishift
3377 ip(1) = 0; // ip(1) not referenced
3378 ip(2) = maxit; // mxiter, maximum number of iterations
3379 ip(3) = 1; // NB blocksize in recurrence
3380 ip(4) = 0; // nconv, number of Ritz values that satisfy convergence
3381 ip(5) = 0; //ip(5) not referenced
3382 ip(6) = mode; // mode
3383 ip(7) = 0;
3384 ip(8) = 0;
3385 ip(9) = 0;
3386 ip(10) = 0;
3387 // ip(7) to ip(10) return values
3388
3389 Array<F77_INT> iptr (dim_vector (14, 1));
3390 F77_INT *ipntr = iptr.rwdata ();
3391
3392 F77_INT ido = 0;
3393 int iter = 0;
3394 M L, U;
3395 ColumnVector r(n);
3396
3397 OCTAVE_LOCAL_BUFFER (octave_idx_type, P, (have_b ? b.rows () : m.rows ()));
3398 OCTAVE_LOCAL_BUFFER (octave_idx_type, Q, (have_b ? b.cols () : m.cols ()));
3399
3400 if (! LuAminusSigmaB (m, b, cholB, permB, sigma, L, U, P, Q, r))
3401 return -1;
3402
3403 F77_INT elems;
3404 F77_INT lwork;
3405
3406 if (octave::math::int_multiply_overflow (n, p, &elems))
3407 (*current_liboctave_error_handler)
3408 ("eigs: array too large for Fortran integers");
3409 if (octave::math::int_multiply_overflow (p, 3 * p + 5, &lwork))
3410 (*current_liboctave_error_handler)
3411 ("eigs: array too large for Fortran integers");
3412
3413 OCTAVE_LOCAL_BUFFER (Complex, v, elems);
3414 OCTAVE_LOCAL_BUFFER (Complex, workl, lwork);
3415 OCTAVE_LOCAL_BUFFER (Complex, workd, 3 * n);
3416 OCTAVE_LOCAL_BUFFER (double, rwork, p);
3417 Complex *presid = cresid.rwdata ();
3418
3419 do
3420 {
3421 F77_INT tmp_info = octave::to_f77_int (info);
3422
3423 F77_FUNC (znaupd, ZNAUPD)
3424 (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
3425 F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2),
3426 k, tol, F77_DBLE_CMPLX_ARG (presid), p, F77_DBLE_CMPLX_ARG (v), n,
3427 iparam, ipntr,
3428 F77_DBLE_CMPLX_ARG (workd), F77_DBLE_CMPLX_ARG (workl), lwork, rwork,
3429 tmp_info F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
3430
3431 info = tmp_info;
3432
3433 if (disp > 0 && ! octave::math::isnan(workl[iptr(5)-1]))
3434 {
3435 if (iter++)
3436 {
3437 os << "Iteration " << iter - 1 <<
3438 ": a few Ritz values of the " << p << "-by-" <<
3439 p << " matrix\n";
3440 if (ido == 99) // convergence
3441 {
3442 for (F77_INT i = 0; i < k; i++)
3443 os << " " << workl[iptr(5)+i-1] << "\n";
3444 }
3445 else
3446 {
3447 // the wanted Ritz estimates are at the end
3448 for (F77_INT i = p - k; i < p; i++)
3449 os << " " << workl[iptr(5)+i-1] << "\n";
3450 }
3451 }
3452
3453 // This is a kludge, as ARPACK doesn't give its
3454 // iteration pointer. But as workl[iptr(5)-1] is
3455 // an output value updated at each iteration, setting
3456 // a value in this array to NaN and testing for it
3457 // is a way of obtaining the iteration counter.
3458 if (ido != 99)
3459 workl[iptr(5)-1] = octave::numeric_limits<double>::NaN ();
3460 }
3461
3462 if (ido == -1 || ido == 1 || ido == 2)
3463 {
3464 if (have_b)
3465 {
3466 if (ido == -1)
3467 {
3468 OCTAVE_LOCAL_BUFFER (Complex, ctmp, n);
3469
3470 vector_product (b, workd+iptr(0)-1, ctmp);
3471
3472 ComplexMatrix tmp (n, 1);
3473
3474 for (F77_INT i = 0; i < n; i++)
3475 tmp(i, 0) = ctmp[P[i]] / r(P[i]);
3476
3477 lusolve (L, U, tmp);
3478
3479 Complex *ip2 = workd+iptr(1)-1;
3480 for (F77_INT i = 0; i < n; i++)
3481 ip2[Q[i]] = tmp(i, 0);
3482 }
3483 else if (ido == 2)
3484 vector_product (b, workd + iptr(0) - 1, workd + iptr(1) - 1);
3485 else
3486 {
3487 Complex *ip2 = workd+iptr(2)-1;
3488 ComplexMatrix tmp (n, 1);
3489
3490 for (F77_INT i = 0; i < n; i++)
3491 tmp(i, 0) = ip2[P[i]] / r(P[i]);
3492
3493 lusolve (L, U, tmp);
3494
3495 ip2 = workd+iptr(1)-1;
3496 for (F77_INT i = 0; i < n; i++)
3497 ip2[Q[i]] = tmp(i, 0);
3498 }
3499 }
3500 else
3501 {
3502 // ido cannot be 2 for non-generalized problems (see znaup2).
3503 Complex *ip2 = workd+iptr(0)-1;
3504 ComplexMatrix tmp (n, 1);
3505
3506 for (F77_INT i = 0; i < n; i++)
3507 tmp(i, 0) = ip2[P[i]] / r(P[i]);
3508
3509 lusolve (L, U, tmp);
3510
3511 ip2 = workd+iptr(1)-1;
3512 for (F77_INT i = 0; i < n; i++)
3513 ip2[Q[i]] = tmp(i, 0);
3514 }
3515 }
3516 else
3517 {
3518 if (info < 0)
3519 (*current_liboctave_error_handler)
3520 ("eigs: error %" OCTAVE_IDX_TYPE_FORMAT " in znaupd: %s",
3521 info, arpack_errno2str (info, "znaupd").c_str ());
3522
3523 break;
3524 }
3525 }
3526 while (1);
3527
3528 F77_INT info2;
3529
3530 // We have a problem in that the size of the C++ bool
3531 // type relative to the fortran logical type. It appears
3532 // that fortran uses 4- or 8-bytes per logical and C++ 1-byte
3533 // per bool, though this might be system dependent. As
3534 // long as the HOWMNY arg is not "S", the logical array
3535 // is just workspace for ARPACK, so use int type to
3536 // avoid problems.
3537 Array<F77_INT> s (dim_vector (p, 1));
3538 F77_INT *sel = s.rwdata ();
3539
3540 eig_vec.resize (n, k);
3541 Complex *z = eig_vec.rwdata ();
3542
3543 eig_val.resize (k+1);
3544 Complex *d = eig_val.rwdata ();
3545
3546 OCTAVE_LOCAL_BUFFER (Complex, workev, 2 * p);
3547
3548 F77_FUNC (zneupd, ZNEUPD)
3549 (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, F77_DBLE_CMPLX_ARG (d),
3550 F77_DBLE_CMPLX_ARG (z), n, F77_CONST_DBLE_CMPLX_ARG (&sigma),
3551 F77_DBLE_CMPLX_ARG (workev),
3552 F77_CONST_CHAR_ARG2 (&bmat, 1), n,
3553 F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2),
3554 k, tol, F77_DBLE_CMPLX_ARG (presid), p, F77_DBLE_CMPLX_ARG (v), n,
3555 iparam, ipntr,
3556 F77_DBLE_CMPLX_ARG (workd), F77_DBLE_CMPLX_ARG (workl), lwork, rwork,
3557 info2 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
3558
3559 if (info2 == 0)
3560 {
3561 for (F77_INT i = ip(4); i < k; i++)
3562 d[i] = Complex (octave::numeric_limits<double>::NaN (),
3563 octave::numeric_limits<double>::NaN ());
3564
3565 F77_INT k2 = ip(4) / 2;
3566 for (F77_INT i = 0; i < k2; i++)
3567 {
3568 Complex ctmp = d[i];
3569 d[i] = d[ip(4) - i - 1];
3570 d[ip(4) - i - 1] = ctmp;
3571 }
3572 eig_val.resize (k);
3573
3574 if (rvec)
3575 {
3576 OCTAVE_LOCAL_BUFFER (Complex, ctmp, n);
3577
3578 for (F77_INT i = ip(4); i < k; i++)
3579 {
3580 F77_INT off1 = i * n;
3581 for (F77_INT j = 0; j < n; j++)
3582 z[off1 + j] = Complex (octave::numeric_limits<double>::NaN (),
3583 octave::numeric_limits<double>::NaN ());
3584 }
3585
3586 for (F77_INT i = 0; i < k2; i++)
3587 {
3588 F77_INT off1 = i * n;
3589 F77_INT off2 = (ip(4) - i - 1) * n;
3590
3591 if (off1 == off2)
3592 continue;
3593
3594 for (F77_INT j = 0; j < n; j++)
3595 ctmp[j] = z[off1 + j];
3596
3597 for (F77_INT j = 0; j < n; j++)
3598 z[off1 + j] = z[off2 + j];
3599
3600 for (F77_INT j = 0; j < n; j++)
3601 z[off2 + j] = ctmp[j];
3602 }
3603 }
3604 }
3605 else
3606 (*current_liboctave_error_handler)
3607 ("eigs: error in zneupd: %s",
3608 arpack_errno2str (info2, "zneupd").c_str ());
3609
3610 return ip(4);
3611}
3612
3613template <typename M>
3614octave_idx_type
3615EigsComplexNonSymmetricFunc (EigsComplexFunc fcn, octave_idx_type n_arg,
3616 const std::string& _typ, Complex sigma,
3617 octave_idx_type k_arg, octave_idx_type p_arg,
3618 octave_idx_type& info, ComplexMatrix& eig_vec,
3619 ComplexColumnVector& eig_val, const M& _b,
3620 ColumnVector& permB, ComplexColumnVector& cresid,
3621 std::ostream& os, double tol, bool rvec,
3622 bool cholB, int disp, int maxit)
3623{
3624 F77_INT n = octave::to_f77_int (n_arg);
3625 F77_INT k = octave::to_f77_int (k_arg);
3626 F77_INT p = octave::to_f77_int (p_arg);
3627 M b(_b);
3628 std::string typ (_typ);
3629 bool have_sigma = (std::abs (sigma) ? true : false);
3630 F77_INT mode = 1;
3631 bool have_b = ! b.isempty ();
3632 bool note3 = false;
3633 char bmat = 'I';
3634 int err = 0;
3635 M bt;
3636
3637 if (cresid.isempty ())
3638 {
3639 std::string rand_dist = octave::rand::distribution ();
3640 octave::rand::distribution ("uniform");
3641 Array<double> rr (octave::rand::vector (n));
3642 Array<double> ri (octave::rand::vector (n));
3643 cresid = ComplexColumnVector (n);
3644 for (F77_INT i = 0; i < n; i++)
3645 cresid(i) = Complex (rr(i), ri(i));
3646 octave::rand::distribution (rand_dist);
3647 }
3648 else if (n != cresid.numel ())
3649 (*current_liboctave_error_handler) ("eigs: opts.v0 must be n-by-1");
3650
3651 if (n < 3)
3652 (*current_liboctave_error_handler) ("eigs: n must be at least 3");
3653
3654 if (p < 0)
3655 {
3656 p = k * 2 + 1;
3657
3658 if (p < 20)
3659 p = 20;
3660
3661 if (p > n)
3662 p = n;
3663 }
3664
3665 if (k <= 0 || k >= n - 1)
3666 (*current_liboctave_error_handler)
3667 ("eigs: Invalid number of eigenvalues to extract"
3668 " (must be 0 < k < n-1).\n"
3669 " Use 'eig (full (A))' instead");
3670
3671 if (p <= k || p > n)
3672 (*current_liboctave_error_handler)
3673 ("eigs: opts.p must be greater than k and less than or equal to n");
3674
3675 if (have_b && cholB && ! permB.isempty ())
3676 {
3677 // Check the we really have a permutation vector
3678 if (permB.numel () != n)
3679 (*current_liboctave_error_handler) ("eigs: permB vector invalid");
3680
3681 Array<bool> checked (dim_vector (n, 1), false);
3682 for (F77_INT i = 0; i < n; i++)
3683 {
3684 octave_idx_type bidx = static_cast<octave_idx_type> (permB(i));
3685
3686 if (checked(bidx) || bidx < 0 || bidx >= n)
3687 (*current_liboctave_error_handler) ("eigs: permB vector invalid");
3688 }
3689 }
3690
3691 if (! have_sigma)
3692 {
3693 if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA"
3694 && typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI"
3695 && typ != "SI")
3696 (*current_liboctave_error_handler) ("eigs: unrecognized sigma value");
3697
3698 if (typ == "LA" || typ == "SA" || typ == "BE")
3699 (*current_liboctave_error_handler)
3700 ("eigs: invalid sigma value for complex problem");
3701
3702 if (typ != "SM" && have_b)
3703 note3 = true;
3704
3705 if (typ == "SM")
3706 {
3707 typ = "LM";
3708 sigma = 0.;
3709 mode = 3;
3710 if (have_b)
3711 bmat ='G';
3712 }
3713 }
3714 else if (! std::abs (sigma))
3715 {
3716 typ = "SM";
3717 if (have_b)
3718 bmat = 'G';
3719 }
3720 else
3721 {
3722 typ = "LM";
3723 mode = 3;
3724 if (have_b)
3725 bmat = 'G';
3726 }
3727
3728 if (mode == 1 && have_b)
3729 {
3730 // See Note 3 dsaupd
3731 note3 = true;
3732 if (cholB)
3733 {
3734 bt = b;
3735 b = b.hermitian ();
3736 if (permB.isempty ())
3737 {
3738 permB = ColumnVector (n);
3739 for (F77_INT i = 0; i < n; i++)
3740 permB(i) = i;
3741 }
3742 }
3743 else
3744 {
3745 if (! make_cholb (b, bt, permB))
3746 (*current_liboctave_error_handler)
3747 ("eigs: The matrix B is not positive definite");
3748 }
3749 }
3750
3751 Array<F77_INT> ip (dim_vector (11, 1));
3752 F77_INT *iparam = ip.rwdata ();
3753
3754 ip(0) = 1; //ishift
3755 ip(1) = 0; // ip(1) not referenced
3756 ip(2) = maxit; // mxiter, maximum number of iterations
3757 ip(3) = 1; // NB blocksize in recurrence
3758 ip(4) = 0; // nconv, number of Ritz values that satisfy convergence
3759 ip(5) = 0; //ip(5) not referenced
3760 ip(6) = mode; // mode
3761 ip(7) = 0;
3762 ip(8) = 0;
3763 ip(9) = 0;
3764 ip(10) = 0;
3765 // ip(7) to ip(10) return values
3766
3767 Array<F77_INT> iptr (dim_vector (14, 1));
3768 F77_INT *ipntr = iptr.rwdata ();
3769
3770 F77_INT ido = 0;
3771 int iter = 0;
3772 F77_INT elems;
3773 F77_INT lwork;
3774
3775 if (octave::math::int_multiply_overflow (n, p, &elems))
3776 (*current_liboctave_error_handler)
3777 ("eigs: array too large for Fortran integers");
3778 if (octave::math::int_multiply_overflow (p, 3 * p + 5, &lwork))
3779 (*current_liboctave_error_handler)
3780 ("eigs: array too large for Fortran integers");
3781
3782 OCTAVE_LOCAL_BUFFER (Complex, v, elems);
3783 OCTAVE_LOCAL_BUFFER (Complex, workl, lwork);
3784 OCTAVE_LOCAL_BUFFER (Complex, workd, 3 * n);
3785 OCTAVE_LOCAL_BUFFER (double, rwork, p);
3786 Complex *presid = cresid.rwdata ();
3787
3788 do
3789 {
3790 F77_INT tmp_info = octave::to_f77_int (info);
3791
3792 F77_FUNC (znaupd, ZNAUPD)
3793 (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
3794 F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2),
3795 k, tol, F77_DBLE_CMPLX_ARG (presid), p, F77_DBLE_CMPLX_ARG (v), n,
3796 iparam, ipntr,
3797 F77_DBLE_CMPLX_ARG (workd), F77_DBLE_CMPLX_ARG (workl), lwork, rwork,
3798 tmp_info F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
3799
3800 info = tmp_info;
3801
3802 if (disp > 0 && ! octave::math::isnan(workl[iptr(5)-1]))
3803 {
3804 if (iter++)
3805 {
3806 os << "Iteration " << iter - 1 <<
3807 ": a few Ritz values of the " << p << "-by-" <<
3808 p << " matrix\n";
3809 if (ido == 99) // convergence
3810 {
3811 for (F77_INT i = 0; i < k; i++)
3812 os << " " << workl[iptr(5)+i-1] << "\n";
3813 }
3814 else
3815 {
3816 // the wanted Ritz estimates are at the end
3817 for (F77_INT i = p - k; i < p; i++)
3818 os << " " << workl[iptr(5)+i-1] << "\n";
3819 }
3820 }
3821
3822 // This is a kludge, as ARPACK doesn't give its
3823 // iteration pointer. But as workl[iptr(5)-1] is
3824 // an output value updated at each iteration, setting
3825 // a value in this array to NaN and testing for it
3826 // is a way of obtaining the iteration counter.
3827 if (ido != 99)
3828 workl[iptr(5)-1] = octave::numeric_limits<double>::NaN ();
3829 }
3830
3831 if (ido == -1 || ido == 1 || ido == 2)
3832 {
3833 if (have_b)
3834 {
3835 if (mode == 1) // regular mode with factorized B
3836 {
3837 ComplexMatrix mtmp (n, 1);
3838 for (F77_INT i = 0; i < n; i++)
3839 mtmp(i, 0) = workd[i + iptr(0) - 1];
3840
3841 mtmp = utsolve (bt, permB, mtmp);
3842 ComplexColumnVector y = fcn (mtmp, err);
3843
3844 if (err)
3845 return false;
3846
3847 mtmp = ltsolve (b, permB, y);
3848
3849 for (F77_INT i = 0; i < n; i++)
3850 workd[i+iptr(1)-1] = mtmp(i, 0);
3851 }
3852 else // shift-invert mode
3853 {
3854 if (ido == -1)
3855 {
3856 OCTAVE_LOCAL_BUFFER (Complex, ctmp, n);
3857
3858 vector_product (b, workd+iptr(0)-1, ctmp);
3859
3860 ComplexColumnVector x(n);
3861
3862 for (F77_INT i = 0; i < n; i++)
3863 x(i) = ctmp[i];
3864
3865 ComplexColumnVector y = fcn (x, err);
3866
3867 if (err)
3868 return false;
3869
3870 Complex *ip2 = workd+iptr(1)-1;
3871 for (F77_INT i = 0; i < n; i++)
3872 ip2[i] = y(i);
3873 }
3874 else if (ido == 2)
3875 vector_product (b, workd+iptr(0)-1, workd+iptr(1)-1);
3876 else
3877 {
3878 Complex *ip2 = workd+iptr(2)-1;
3879 ComplexColumnVector x(n);
3880
3881 for (F77_INT i = 0; i < n; i++)
3882 x(i) = *ip2++;
3883
3884 ComplexColumnVector y = fcn (x, err);
3885
3886 if (err)
3887 return false;
3888
3889 ip2 = workd + iptr(1) - 1;
3890 for (F77_INT i = 0; i < n; i++)
3891 *ip2++ = y(i);
3892 }
3893 }
3894 }
3895 else
3896 {
3897 Complex *ip2 = workd + iptr(0) - 1;
3898 ComplexColumnVector x(n);
3899
3900 for (F77_INT i = 0; i < n; i++)
3901 x(i) = *ip2++;
3902
3903 ComplexColumnVector y = fcn (x, err);
3904
3905 if (err)
3906 return false;
3907
3908 ip2 = workd + iptr(1) - 1;
3909 for (F77_INT i = 0; i < n; i++)
3910 *ip2++ = y(i);
3911 }
3912 }
3913 else
3914 {
3915 if (info < 0)
3916 (*current_liboctave_error_handler)
3917 ("eigs: error %" OCTAVE_IDX_TYPE_FORMAT " in znaupd: %s",
3918 info, arpack_errno2str (info, "znaupd").c_str ());
3919
3920 break;
3921 }
3922 }
3923 while (1);
3924
3925 F77_INT info2;
3926
3927 // We have a problem in that the size of the C++ bool
3928 // type relative to the fortran logical type. It appears
3929 // that fortran uses 4- or 8-bytes per logical and C++ 1-byte
3930 // per bool, though this might be system dependent. As
3931 // long as the HOWMNY arg is not "S", the logical array
3932 // is just workspace for ARPACK, so use int type to
3933 // avoid problems.
3934 Array<F77_INT> s (dim_vector (p, 1));
3935 F77_INT *sel = s.rwdata ();
3936
3937 eig_vec.resize (n, k);
3938 Complex *z = eig_vec.rwdata ();
3939
3940 eig_val.resize (k+1);
3941 Complex *d = eig_val.rwdata ();
3942
3943 OCTAVE_LOCAL_BUFFER (Complex, workev, 2 * p);
3944
3945 F77_FUNC (zneupd, ZNEUPD)
3946 (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, F77_DBLE_CMPLX_ARG (d),
3947 F77_DBLE_CMPLX_ARG (z), n, F77_DBLE_CMPLX_ARG (&sigma),
3948 F77_DBLE_CMPLX_ARG (workev),
3949 F77_CONST_CHAR_ARG2 (&bmat, 1), n,
3950 F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2),
3951 k, tol, F77_DBLE_CMPLX_ARG (presid), p, F77_DBLE_CMPLX_ARG (v), n,
3952 iparam, ipntr,
3953 F77_DBLE_CMPLX_ARG (workd), F77_DBLE_CMPLX_ARG (workl), lwork, rwork,
3954 info2 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
3955
3956 if (info2 == 0)
3957 {
3958 for (F77_INT i = ip(4); i < k; i++)
3959 d[i] = Complex (octave::numeric_limits<double>::NaN (),
3960 octave::numeric_limits<double>::NaN ());
3961
3962 F77_INT k2 = ip(4) / 2;
3963 for (F77_INT i = 0; i < k2; i++)
3964 {
3965 Complex ctmp = d[i];
3966 d[i] = d[ip(4) - i - 1];
3967 d[ip(4) - i - 1] = ctmp;
3968 }
3969 eig_val.resize (k);
3970
3971 if (rvec)
3972 {
3973 OCTAVE_LOCAL_BUFFER (Complex, ctmp, n);
3974
3975 for (F77_INT i = ip(4); i < k; i++)
3976 {
3977 F77_INT off1 = i * n;
3978 for (F77_INT j = 0; j < n; j++)
3979 z[off1 + j] = Complex (octave::numeric_limits<double>::NaN (),
3980 octave::numeric_limits<double>::NaN ());
3981 }
3982
3983 for (F77_INT i = 0; i < k2; i++)
3984 {
3985 F77_INT off1 = i * n;
3986 F77_INT off2 = (ip(4) - i - 1) * n;
3987
3988 if (off1 == off2)
3989 continue;
3990
3991 for (F77_INT j = 0; j < n; j++)
3992 ctmp[j] = z[off1 + j];
3993
3994 for (F77_INT j = 0; j < n; j++)
3995 z[off1 + j] = z[off2 + j];
3996
3997 for (F77_INT j = 0; j < n; j++)
3998 z[off2 + j] = ctmp[j];
3999 }
4000 if (note3)
4001 eig_vec = utsolve (bt, permB, eig_vec);
4002 }
4003 }
4004 else
4005 (*current_liboctave_error_handler)
4006 ("eigs: error in zneupd: %s",
4007 arpack_errno2str (info2, "zneupd").c_str ());
4008
4009 return ip(4);
4010}
4011
4012// Instantiations for the types we need.
4013
4014// Matrix
4015
4016template
4017OCTAVE_API octave_idx_type
4018EigsRealSymmetricMatrix<Matrix>
4019 (const Matrix& m, const std::string typ, octave_idx_type k,
4020 octave_idx_type p, octave_idx_type& info, Matrix& eig_vec,
4021 ColumnVector& eig_val, const Matrix& _b, ColumnVector& permB,
4022 ColumnVector& resid, std::ostream& os, double tol, bool rvec,
4023 bool cholB, int disp, int maxit);
4024
4025template
4026OCTAVE_API octave_idx_type
4027EigsRealSymmetricMatrixShift<Matrix>
4028 (const Matrix& m, double sigma, octave_idx_type k,
4029 octave_idx_type p, octave_idx_type& info, Matrix& eig_vec,
4030 ColumnVector& eig_val, const Matrix& _b, ColumnVector& permB,
4031 ColumnVector& resid, std::ostream& os, double tol, bool rvec,
4032 bool cholB, int disp, int maxit);
4033
4034template
4035OCTAVE_API octave_idx_type
4036EigsRealSymmetricFunc<Matrix>
4037 (EigsFunc fcn, octave_idx_type n, const std::string& _typ, double sigma,
4038 octave_idx_type k, octave_idx_type p, octave_idx_type& info,
4039 Matrix& eig_vec, ColumnVector& eig_val, const Matrix& _b,
4040 ColumnVector& permB, ColumnVector& resid, std::ostream& os, double tol,
4041 bool rvec, bool cholB, int disp, int maxit);
4042
4043template
4044OCTAVE_API octave_idx_type
4045EigsRealNonSymmetricMatrix<Matrix>
4046 (const Matrix& m, const std::string typ, octave_idx_type k,
4047 octave_idx_type p, octave_idx_type& info, ComplexMatrix& eig_vec,
4048 ComplexColumnVector& eig_val, const Matrix& _b, ColumnVector& permB,
4049 ColumnVector& resid, std::ostream& os, double tol, bool rvec,
4050 bool cholB, int disp, int maxit);
4051
4052template
4053OCTAVE_API octave_idx_type
4054EigsRealNonSymmetricMatrixShift<Matrix>
4055 (const Matrix& m, double sigmar, octave_idx_type k,
4056 octave_idx_type p, octave_idx_type& info, ComplexMatrix& eig_vec,
4057 ComplexColumnVector& eig_val, const Matrix& _b, ColumnVector& permB,
4058 ColumnVector& resid, std::ostream& os, double tol, bool rvec,
4059 bool cholB, int disp, int maxit);
4060
4061template
4062OCTAVE_API octave_idx_type
4063EigsRealNonSymmetricFunc<Matrix>
4064 (EigsFunc fcn, octave_idx_type n, const std::string& _typ, double sigmar,
4065 octave_idx_type k, octave_idx_type p, octave_idx_type& info,
4066 ComplexMatrix& eig_vec, ComplexColumnVector& eig_val, const Matrix& _b,
4067 ColumnVector& permB, ColumnVector& resid, std::ostream& os, double tol,
4068 bool rvec, bool cholB, int disp, int maxit);
4069
4070// SparseMatrix
4071
4072template
4073OCTAVE_API octave_idx_type
4074EigsRealSymmetricMatrix<SparseMatrix>
4075 (const SparseMatrix& m, const std::string typ, octave_idx_type k,
4076 octave_idx_type p, octave_idx_type& info, Matrix& eig_vec,
4077 ColumnVector& eig_val, const SparseMatrix& _b, ColumnVector& permB,
4078 ColumnVector& resid, std::ostream& os, double tol, bool rvec,
4079 bool cholB, int disp, int maxit);
4080
4081template
4082OCTAVE_API octave_idx_type
4083EigsRealSymmetricMatrixShift<SparseMatrix>
4084 (const SparseMatrix& m, double sigma, octave_idx_type k,
4085 octave_idx_type p, octave_idx_type& info, Matrix& eig_vec,
4086 ColumnVector& eig_val, const SparseMatrix& _b, ColumnVector& permB,
4087 ColumnVector& resid, std::ostream& os, double tol, bool rvec,
4088 bool cholB, int disp, int maxit);
4089
4090template
4091OCTAVE_API octave_idx_type
4092EigsRealSymmetricFunc<SparseMatrix>
4093 (EigsFunc fcn, octave_idx_type n, const std::string& _typ, double sigma,
4094 octave_idx_type k, octave_idx_type p, octave_idx_type& info,
4095 Matrix& eig_vec, ColumnVector& eig_val, const SparseMatrix& _b,
4096 ColumnVector& permB, ColumnVector& resid, std::ostream& os, double tol,
4097 bool rvec, bool cholB, int disp, int maxit);
4098
4099template
4100OCTAVE_API octave_idx_type
4101EigsRealNonSymmetricMatrix<SparseMatrix>
4102 (const SparseMatrix& m, const std::string typ, octave_idx_type k,
4103 octave_idx_type p, octave_idx_type& info, ComplexMatrix& eig_vec,
4104 ComplexColumnVector& eig_val, const SparseMatrix& _b, ColumnVector& permB,
4105 ColumnVector& resid, std::ostream& os, double tol, bool rvec,
4106 bool cholB, int disp, int maxit);
4107
4108template
4109OCTAVE_API octave_idx_type
4110EigsRealNonSymmetricMatrixShift<SparseMatrix>
4111 (const SparseMatrix& m, double sigmar, octave_idx_type k,
4112 octave_idx_type p, octave_idx_type& info, ComplexMatrix& eig_vec,
4113 ComplexColumnVector& eig_val, const SparseMatrix& _b, ColumnVector& permB,
4114 ColumnVector& resid, std::ostream& os, double tol, bool rvec,
4115 bool cholB, int disp, int maxit);
4116
4117template
4118OCTAVE_API octave_idx_type
4119EigsRealNonSymmetricFunc<SparseMatrix>
4120 (EigsFunc fcn, octave_idx_type n, const std::string& _typ, double sigmar,
4121 octave_idx_type k, octave_idx_type p, octave_idx_type& info,
4122 ComplexMatrix& eig_vec, ComplexColumnVector& eig_val,
4123 const SparseMatrix& _b, ColumnVector& permB, ColumnVector& resid,
4124 std::ostream& os, double tol, bool rvec, bool cholB, int disp, int maxit);
4125
4126// ComplexMatrix
4127
4128template
4129OCTAVE_API octave_idx_type
4130EigsComplexNonSymmetricMatrix<ComplexMatrix>
4131 (const ComplexMatrix& m, const std::string typ, octave_idx_type k,
4132 octave_idx_type p, octave_idx_type& info, ComplexMatrix& eig_vec,
4133 ComplexColumnVector& eig_val, const ComplexMatrix& _b, ColumnVector& permB,
4134 ComplexColumnVector& cresid, std::ostream& os, double tol,
4135 bool rvec, bool cholB, int disp, int maxit);
4136
4137template
4138OCTAVE_API octave_idx_type
4139EigsComplexNonSymmetricMatrixShift<ComplexMatrix>
4140 (const ComplexMatrix& m, Complex sigma, octave_idx_type k,
4141 octave_idx_type p, octave_idx_type& info, ComplexMatrix& eig_vec,
4142 ComplexColumnVector& eig_val, const ComplexMatrix& _b, ColumnVector& permB,
4143 ComplexColumnVector& cresid, std::ostream& os, double tol,
4144 bool rvec, bool cholB, int disp, int maxit);
4145
4146template
4147OCTAVE_API octave_idx_type
4148EigsComplexNonSymmetricFunc<ComplexMatrix>
4149 (EigsComplexFunc fcn, octave_idx_type n, const std::string& _typ, Complex sigma,
4150 octave_idx_type k, octave_idx_type p, octave_idx_type& info,
4151 ComplexMatrix& eig_vec, ComplexColumnVector& eig_val,
4152 const ComplexMatrix& _b, ColumnVector& permB, ComplexColumnVector& cresid,
4153 std::ostream& os, double tol, bool rvec, bool cholB, int disp, int maxit);
4154
4155// SparseComplexMatrix
4156
4157template
4158OCTAVE_API octave_idx_type
4159EigsComplexNonSymmetricMatrix<SparseComplexMatrix>
4160 (const SparseComplexMatrix& m, const std::string typ, octave_idx_type k,
4161 octave_idx_type p, octave_idx_type& info, ComplexMatrix& eig_vec,
4162 ComplexColumnVector& eig_val, const SparseComplexMatrix& _b,
4163 ColumnVector& permB, ComplexColumnVector& cresid, std::ostream& os,
4164 double tol, bool rvec, bool cholB, int disp, int maxit);
4165
4166template
4167OCTAVE_API octave_idx_type
4168EigsComplexNonSymmetricMatrixShift<SparseComplexMatrix>
4169 (const SparseComplexMatrix& m, Complex sigma, octave_idx_type k,
4170 octave_idx_type p, octave_idx_type& info, ComplexMatrix& eig_vec,
4171 ComplexColumnVector& eig_val, const SparseComplexMatrix& _b,
4172 ColumnVector& permB, ComplexColumnVector& cresid, std::ostream& os,
4173 double tol, bool rvec, bool cholB, int disp, int maxit);
4174
4175template
4176OCTAVE_API octave_idx_type
4177EigsComplexNonSymmetricFunc<SparseComplexMatrix>
4178 (EigsComplexFunc fcn, octave_idx_type n, const std::string& _typ, Complex sigma,
4179 octave_idx_type k, octave_idx_type p, octave_idx_type& info,
4180 ComplexMatrix& eig_vec, ComplexColumnVector& eig_val,
4181 const SparseComplexMatrix& _b, ColumnVector& permB,
4182 ComplexColumnVector& cresid, std::ostream& os, double tol, bool rvec,
4183 bool cholB, int disp, int maxit);
4184#endif
Array
N Dimensional Array with copy-on-write semantics.
Definition Array-base.h:130
Array::xelem
T & xelem(octave_idx_type n)
Size of the specified dimension.
Definition Array-base.h:547
Array::rows
octave_idx_type rows() const
Definition Array-base.h:485
Array::cols
octave_idx_type cols() const
Definition Array-base.h:495
Array::isempty
bool isempty() const
Size of the specified dimension.
Definition Array-base.h:674
Array::data
const T * data() const
Size of the specified dimension.
Definition Array-base.h:687
Array::rwdata
T * rwdata()
Size of the specified dimension.
Definition Array-base.cc:1790
Array::numel
octave_idx_type numel() const
Number of elements in the array.
Definition Array-base.h:440
ColumnVector::resize
void resize(octave_idx_type n, const double &rfv=0)
Definition dColVector.h:112
ComplexColumnVector::resize
void resize(octave_idx_type n, const Complex &rfv=Complex(0))
Definition CColVector.h:146
ComplexMatrix::resize
void resize(octave_idx_type nr, octave_idx_type nc, const Complex &rfv=Complex(0))
Definition CMatrix.h:191
ComplexMatrix::hermitian
ComplexMatrix hermitian() const
Definition CMatrix.h:168
Matrix::transpose
Matrix transpose() const
Definition dMatrix.h:138
Matrix::resize
void resize(octave_idx_type nr, octave_idx_type nc, double rfv=0)
Definition dMatrix.h:156
SparseComplexMatrix::hermitian
SparseComplexMatrix hermitian() const
Definition CSparse.cc:853
SparseMatrix::transpose
SparseMatrix transpose() const
Definition dSparse.h:131
Sparse::cols
octave_idx_type cols() const
Definition Sparse.h:351
Sparse::cidx
octave_idx_type * cidx()
Definition Sparse.h:595
Sparse::data
T * data()
Definition Sparse.h:573
Sparse::xdata
T * xdata()
Definition Sparse.h:575
Sparse::ridx
octave_idx_type * ridx()
Definition Sparse.h:582
Sparse::rows
octave_idx_type rows() const
Definition Sparse.h:350
Sparse::xcidx
octave_idx_type * xcidx()
Definition Sparse.h:601
Sparse::xridx
octave_idx_type * xridx()
Definition Sparse.h:588
Sparse::isempty
bool isempty() const
Definition Sparse.h:569
dim_vector
Vector representing the dimensions (size) of an Array.
Definition dim-vector.h:92
octave_idx_type
EigsComplexNonSymmetricFunc< ComplexMatrix >
template octave_idx_type EigsComplexNonSymmetricFunc< ComplexMatrix >(EigsComplexFunc fcn, octave_idx_type n, const std::string &_typ, Complex sigma, octave_idx_type k, octave_idx_type p, octave_idx_type &info, ComplexMatrix &eig_vec, ComplexColumnVector &eig_val, const ComplexMatrix &_b, ColumnVector &permB, ComplexColumnVector &cresid, std::ostream &os, double tol, bool rvec, bool cholB, int disp, int maxit)
EigsComplexNonSymmetricMatrix< SparseComplexMatrix >
template octave_idx_type EigsComplexNonSymmetricMatrix< SparseComplexMatrix >(const SparseComplexMatrix &m, const std::string typ, octave_idx_type k, octave_idx_type p, octave_idx_type &info, ComplexMatrix &eig_vec, ComplexColumnVector &eig_val, const SparseComplexMatrix &_b, ColumnVector &permB, ComplexColumnVector &cresid, std::ostream &os, double tol, bool rvec, bool cholB, int disp, int maxit)
EigsComplexNonSymmetricMatrixShift< SparseComplexMatrix >
template octave_idx_type EigsComplexNonSymmetricMatrixShift< SparseComplexMatrix >(const SparseComplexMatrix &m, Complex sigma, octave_idx_type k, octave_idx_type p, octave_idx_type &info, ComplexMatrix &eig_vec, ComplexColumnVector &eig_val, const SparseComplexMatrix &_b, ColumnVector &permB, ComplexColumnVector &cresid, std::ostream &os, double tol, bool rvec, bool cholB, int disp, int maxit)
EigsRealSymmetricFunc< Matrix >
template octave_idx_type EigsRealSymmetricFunc< Matrix >(EigsFunc fcn, octave_idx_type n, const std::string &_typ, double sigma, octave_idx_type k, octave_idx_type p, octave_idx_type &info, Matrix &eig_vec, ColumnVector &eig_val, const Matrix &_b, ColumnVector &permB, ColumnVector &resid, std::ostream &os, double tol, bool rvec, bool cholB, int disp, int maxit)
EigsRealNonSymmetricMatrixShift
octave_idx_type EigsRealNonSymmetricMatrixShift(const M &m, double sigmar, octave_idx_type k_arg, octave_idx_type p_arg, octave_idx_type &info, ComplexMatrix &eig_vec, ComplexColumnVector &eig_val, const M &_b, ColumnVector &permB, ColumnVector &resid, std::ostream &os, double tol, bool rvec, bool cholB, int disp, int maxit)
Definition eigs-base.cc:2146
EigsRealSymmetricMatrixShift
octave_idx_type EigsRealSymmetricMatrixShift(const M &m, double sigma, octave_idx_type k_arg, octave_idx_type p_arg, octave_idx_type &info, Matrix &eig_vec, ColumnVector &eig_val, const M &_b, ColumnVector &permB, ColumnVector &resid, std::ostream &os, double tol, bool rvec, bool cholB, int disp, int maxit)
Definition eigs-base.cc:1090
EigsRealNonSymmetricMatrixShift< Matrix >
template octave_idx_type EigsRealNonSymmetricMatrixShift< Matrix >(const Matrix &m, double sigmar, octave_idx_type k, octave_idx_type p, octave_idx_type &info, ComplexMatrix &eig_vec, ComplexColumnVector &eig_val, const Matrix &_b, ColumnVector &permB, ColumnVector &resid, std::ostream &os, double tol, bool rvec, bool cholB, int disp, int maxit)
EigsRealNonSymmetricMatrix
octave_idx_type EigsRealNonSymmetricMatrix(const M &m, const std::string typ, octave_idx_type k_arg, octave_idx_type p_arg, octave_idx_type &info, ComplexMatrix &eig_vec, ComplexColumnVector &eig_val, const M &_b, ColumnVector &permB, ColumnVector &resid, std::ostream &os, double tol, bool rvec, bool cholB, int disp, int maxit)
Definition eigs-base.cc:1788
EigsComplexNonSymmetricMatrix
octave_idx_type EigsComplexNonSymmetricMatrix(const M &m, const std::string typ, octave_idx_type k_arg, octave_idx_type p_arg, octave_idx_type &info, ComplexMatrix &eig_vec, ComplexColumnVector &eig_val, const M &_b, ColumnVector &permB, ComplexColumnVector &cresid, std::ostream &os, double tol, bool rvec, bool cholB, int disp, int maxit)
Definition eigs-base.cc:2975
EigsRealSymmetricFunc
octave_idx_type EigsRealSymmetricFunc(EigsFunc fcn, octave_idx_type n_arg, const std::string &_typ, double sigma, octave_idx_type k_arg, octave_idx_type p_arg, octave_idx_type &info, Matrix &eig_vec, ColumnVector &eig_val, const M &_b, ColumnVector &permB, ColumnVector &resid, std::ostream &os, double tol, bool rvec, bool cholB, int disp, int maxit)
Definition eigs-base.cc:1400
arpack_errno2str
std::string arpack_errno2str(const octave_idx_type &errnum, const std::string &fcn_name)
Definition eigs-base.cc:66
EigsRealSymmetricMatrixShift< SparseMatrix >
template octave_idx_type EigsRealSymmetricMatrixShift< SparseMatrix >(const SparseMatrix &m, double sigma, octave_idx_type k, octave_idx_type p, octave_idx_type &info, Matrix &eig_vec, ColumnVector &eig_val, const SparseMatrix &_b, ColumnVector &permB, ColumnVector &resid, std::ostream &os, double tol, bool rvec, bool cholB, int disp, int maxit)
EigsRealNonSymmetricMatrixShift< SparseMatrix >
template octave_idx_type EigsRealNonSymmetricMatrixShift< SparseMatrix >(const SparseMatrix &m, double sigmar, octave_idx_type k, octave_idx_type p, octave_idx_type &info, ComplexMatrix &eig_vec, ComplexColumnVector &eig_val, const SparseMatrix &_b, ColumnVector &permB, ColumnVector &resid, std::ostream &os, double tol, bool rvec, bool cholB, int disp, int maxit)
EigsComplexNonSymmetricMatrixShift< ComplexMatrix >
template octave_idx_type EigsComplexNonSymmetricMatrixShift< ComplexMatrix >(const ComplexMatrix &m, Complex sigma, octave_idx_type k, octave_idx_type p, octave_idx_type &info, ComplexMatrix &eig_vec, ComplexColumnVector &eig_val, const ComplexMatrix &_b, ColumnVector &permB, ComplexColumnVector &cresid, std::ostream &os, double tol, bool rvec, bool cholB, int disp, int maxit)
EigsRealNonSymmetricFunc< SparseMatrix >
template octave_idx_type EigsRealNonSymmetricFunc< SparseMatrix >(EigsFunc fcn, octave_idx_type n, const std::string &_typ, double sigmar, octave_idx_type k, octave_idx_type p, octave_idx_type &info, ComplexMatrix &eig_vec, ComplexColumnVector &eig_val, const SparseMatrix &_b, ColumnVector &permB, ColumnVector &resid, std::ostream &os, double tol, bool rvec, bool cholB, int disp, int maxit)
EigsRealNonSymmetricMatrix< Matrix >
template octave_idx_type EigsRealNonSymmetricMatrix< Matrix >(const Matrix &m, const std::string typ, octave_idx_type k, octave_idx_type p, octave_idx_type &info, ComplexMatrix &eig_vec, ComplexColumnVector &eig_val, const Matrix &_b, ColumnVector &permB, ColumnVector &resid, std::ostream &os, double tol, bool rvec, bool cholB, int disp, int maxit)
EigsRealSymmetricMatrixShift< Matrix >
template octave_idx_type EigsRealSymmetricMatrixShift< Matrix >(const Matrix &m, double sigma, octave_idx_type k, octave_idx_type p, octave_idx_type &info, Matrix &eig_vec, ColumnVector &eig_val, const Matrix &_b, ColumnVector &permB, ColumnVector &resid, std::ostream &os, double tol, bool rvec, bool cholB, int disp, int maxit)
EigsRealSymmetricMatrix< Matrix >
template octave_idx_type EigsRealSymmetricMatrix< Matrix >(const Matrix &m, const std::string typ, octave_idx_type k, octave_idx_type p, octave_idx_type &info, Matrix &eig_vec, ColumnVector &eig_val, const Matrix &_b, ColumnVector &permB, ColumnVector &resid, std::ostream &os, double tol, bool rvec, bool cholB, int disp, int maxit)
EigsRealNonSymmetricFunc
octave_idx_type EigsRealNonSymmetricFunc(EigsFunc fcn, octave_idx_type n_arg, const std::string &_typ, double sigmar, octave_idx_type k_arg, octave_idx_type p_arg, octave_idx_type &info, ComplexMatrix &eig_vec, ComplexColumnVector &eig_val, const M &_b, ColumnVector &permB, ColumnVector &resid, std::ostream &os, double tol, bool rvec, bool cholB, int disp, int maxit)
Definition eigs-base.cc:2525
EigsComplexNonSymmetricMatrixShift
octave_idx_type EigsComplexNonSymmetricMatrixShift(const M &m, Complex sigma, octave_idx_type k_arg, octave_idx_type p_arg, octave_idx_type &info, ComplexMatrix &eig_vec, ComplexColumnVector &eig_val, const M &_b, ColumnVector &permB, ComplexColumnVector &cresid, std::ostream &os, double tol, bool rvec, bool cholB, int disp, int maxit)
Definition eigs-base.cc:3285
EigsRealNonSymmetricMatrix< SparseMatrix >
template octave_idx_type EigsRealNonSymmetricMatrix< SparseMatrix >(const SparseMatrix &m, const std::string typ, octave_idx_type k, octave_idx_type p, octave_idx_type &info, ComplexMatrix &eig_vec, ComplexColumnVector &eig_val, const SparseMatrix &_b, ColumnVector &permB, ColumnVector &resid, std::ostream &os, double tol, bool rvec, bool cholB, int disp, int maxit)
EigsRealSymmetricFunc< SparseMatrix >
template octave_idx_type EigsRealSymmetricFunc< SparseMatrix >(EigsFunc fcn, octave_idx_type n, const std::string &_typ, double sigma, octave_idx_type k, octave_idx_type p, octave_idx_type &info, Matrix &eig_vec, ColumnVector &eig_val, const SparseMatrix &_b, ColumnVector &permB, ColumnVector &resid, std::ostream &os, double tol, bool rvec, bool cholB, int disp, int maxit)
EigsRealSymmetricMatrix
octave_idx_type EigsRealSymmetricMatrix(const M &m, const std::string typ, octave_idx_type k_arg, octave_idx_type p_arg, octave_idx_type &info, Matrix &eig_vec, ColumnVector &eig_val, const M &_b, ColumnVector &permB, ColumnVector &resid, std::ostream &os, double tol, bool rvec, bool cholB, int disp, int maxit)
Definition eigs-base.cc:791
EigsComplexNonSymmetricFunc
octave_idx_type EigsComplexNonSymmetricFunc(EigsComplexFunc fcn, octave_idx_type n_arg, const std::string &_typ, Complex sigma, octave_idx_type k_arg, octave_idx_type p_arg, octave_idx_type &info, ComplexMatrix &eig_vec, ComplexColumnVector &eig_val, const M &_b, ColumnVector &permB, ComplexColumnVector &cresid, std::ostream &os, double tol, bool rvec, bool cholB, int disp, int maxit)
Definition eigs-base.cc:3615
EigsComplexNonSymmetricFunc< SparseComplexMatrix >
template octave_idx_type EigsComplexNonSymmetricFunc< SparseComplexMatrix >(EigsComplexFunc fcn, octave_idx_type n, const std::string &_typ, Complex sigma, octave_idx_type k, octave_idx_type p, octave_idx_type &info, ComplexMatrix &eig_vec, ComplexColumnVector &eig_val, const SparseComplexMatrix &_b, ColumnVector &permB, ComplexColumnVector &cresid, std::ostream &os, double tol, bool rvec, bool cholB, int disp, int maxit)
EigsComplexNonSymmetricMatrix< ComplexMatrix >
template octave_idx_type EigsComplexNonSymmetricMatrix< ComplexMatrix >(const ComplexMatrix &m, const std::string typ, octave_idx_type k, octave_idx_type p, octave_idx_type &info, ComplexMatrix &eig_vec, ComplexColumnVector &eig_val, const ComplexMatrix &_b, ColumnVector &permB, ComplexColumnVector &cresid, std::ostream &os, double tol, bool rvec, bool cholB, int disp, int maxit)
EigsRealSymmetricMatrix< SparseMatrix >
template octave_idx_type EigsRealSymmetricMatrix< SparseMatrix >(const SparseMatrix &m, const std::string typ, octave_idx_type k, octave_idx_type p, octave_idx_type &info, Matrix &eig_vec, ColumnVector &eig_val, const SparseMatrix &_b, ColumnVector &permB, ColumnVector &resid, std::ostream &os, double tol, bool rvec, bool cholB, int disp, int maxit)
EigsRealNonSymmetricFunc< Matrix >
template octave_idx_type EigsRealNonSymmetricFunc< Matrix >(EigsFunc fcn, octave_idx_type n, const std::string &_typ, double sigmar, octave_idx_type k, octave_idx_type p, octave_idx_type &info, ComplexMatrix &eig_vec, ComplexColumnVector &eig_val, const Matrix &_b, ColumnVector &permB, ColumnVector &resid, std::ostream &os, double tol, bool rvec, bool cholB, int disp, int maxit)
EigsFunc
std::function< ColumnVector(const ColumnVector &x, int &eigs_error)> EigsFunc
Definition eigs-base.h:39
EigsComplexFunc
std::function< ComplexColumnVector(const ComplexColumnVector &x, int &eigs_error)> EigsComplexFunc
Definition eigs-base.h:43
F77_DBLE_CMPLX_ARG
#define F77_DBLE_CMPLX_ARG(x)
Definition f77-fcn.h:316
F77_XFCN
#define F77_XFCN(f, F, args)
Definition f77-fcn.h:45
F77_INT
octave_f77_int_type F77_INT
Definition f77-fcn.h:306
F77_CONST_DBLE_CMPLX_ARG
#define F77_CONST_DBLE_CMPLX_ARG(x)
Definition f77-fcn.h:319
Q
F77_RET_T const F77_INT const F77_INT const F77_INT F77_DBLE const F77_INT F77_DBLE const F77_INT F77_DBLE * Q
Definition lapack-proto.h:995
M
F77_RET_T const F77_INT const F77_INT const F77_INT const F77_DBLE const F77_DBLE F77_INT & M
Definition lapack-proto.h:904
current_liboctave_error_handler
OCTAVE_NORETURN liboctave_error_handler current_liboctave_error_handler
Definition lo-error.c:41
OCTAVE_API
#define OCTAVE_API
Definition main.in.cc:55
Complex
std::complex< double > Complex
Definition oct-cmplx.h:33
OCTAVE_LOCAL_BUFFER
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition oct-locbuf.h:44
d
F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE * d
Definition slatec-proto.h:39
x
F77_RET_T const F77_DBLE * x
Definition slatec-proto.h:38
F77_FUNC
F77_RET_T F77_FUNC(xerbla, XERBLA)(F77_CONST_CHAR_ARG_DEF(s_arg