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