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