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