GNU Octave  3.8.0
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__.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2005-2013 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 #ifdef HAVE_CONFIG_H
24 #include <config.h>
25 #endif
26 
27 #include "ov.h"
28 #include "defun-dld.h"
29 #include "error.h"
30 #include "gripes.h"
31 #include "quit.h"
32 #include "variables.h"
33 #include "ov-re-sparse.h"
34 #include "ov-cx-sparse.h"
35 #include "oct-map.h"
36 #include "pager.h"
37 #include "unwind-prot.h"
38 
39 #include "eigs-base.cc"
40 
41 // Global pointer for user defined function.
43 
44 // Have we warned about imaginary values returned from user function?
45 static bool warned_imaginary = false;
46 
47 // Is this a recursive call?
48 static int call_depth = 0;
49 
51 eigs_func (const ColumnVector &x, int &eigs_error)
52 {
53  ColumnVector retval;
54  octave_value_list args;
55  args(0) = x;
56 
57  if (eigs_fcn)
58  {
59  octave_value_list tmp = eigs_fcn->do_multi_index_op (1, args);
60 
61  if (error_state)
62  {
63  eigs_error = 1;
64  gripe_user_supplied_eval ("eigs");
65  return retval;
66  }
67 
68  if (tmp.length () && tmp(0).is_defined ())
69  {
70  if (! warned_imaginary && tmp(0).is_complex_type ())
71  {
72  warning ("eigs: ignoring imaginary part returned from user-supplied function");
73  warned_imaginary = true;
74  }
75 
76  retval = ColumnVector (tmp(0).vector_value ());
77 
78  if (error_state)
79  {
80  eigs_error = 1;
81  gripe_user_supplied_eval ("eigs");
82  }
83  }
84  else
85  {
86  eigs_error = 1;
87  gripe_user_supplied_eval ("eigs");
88  }
89  }
90 
91  return retval;
92 }
93 
95 eigs_complex_func (const ComplexColumnVector &x, int &eigs_error)
96 {
97  ComplexColumnVector retval;
98  octave_value_list args;
99  args(0) = x;
100 
101  if (eigs_fcn)
102  {
103  octave_value_list tmp = eigs_fcn->do_multi_index_op (1, args);
104 
105  if (error_state)
106  {
107  eigs_error = 1;
108  gripe_user_supplied_eval ("eigs");
109  return retval;
110  }
111 
112  if (tmp.length () && tmp(0).is_defined ())
113  {
114  retval = ComplexColumnVector (tmp(0).complex_vector_value ());
115 
116  if (error_state)
117  {
118  eigs_error = 1;
119  gripe_user_supplied_eval ("eigs");
120  }
121  }
122  else
123  {
124  eigs_error = 1;
125  gripe_user_supplied_eval ("eigs");
126  }
127  }
128 
129  return retval;
130 }
131 
132 DEFUN_DLD (__eigs__, args, nargout,
133  "-*- texinfo -*-\n\
134 @deftypefn {Loadable Function} {@var{d} =} __eigs__ (@var{A})\n\
135 @deftypefnx {Loadable Function} {@var{d} =} __eigs__ (@var{A}, @var{k})\n\
136 @deftypefnx {Loadable Function} {@var{d} =} __eigs__ (@var{A}, @var{k}, @var{sigma})\n\
137 @deftypefnx {Loadable Function} {@var{d} =} __eigs__ (@var{A}, @var{k}, @var{sigma}, @var{opts})\n\
138 @deftypefnx {Loadable Function} {@var{d} =} __eigs__ (@var{A}, @var{B})\n\
139 @deftypefnx {Loadable Function} {@var{d} =} __eigs__ (@var{A}, @var{B}, @var{k})\n\
140 @deftypefnx {Loadable Function} {@var{d} =} __eigs__ (@var{A}, @var{B}, @var{k}, @var{sigma})\n\
141 @deftypefnx {Loadable Function} {@var{d} =} __eigs__ (@var{A}, @var{B}, @var{k}, @var{sigma}, @var{opts})\n\
142 @deftypefnx {Loadable Function} {@var{d} =} __eigs__ (@var{af}, @var{n})\n\
143 @deftypefnx {Loadable Function} {@var{d} =} __eigs__ (@var{af}, @var{n}, @var{B})\n\
144 @deftypefnx {Loadable Function} {@var{d} =} __eigs__ (@var{af}, @var{n}, @var{k})\n\
145 @deftypefnx {Loadable Function} {@var{d} =} __eigs__ (@var{af}, @var{n}, @var{B}, @var{k})\n\
146 @deftypefnx {Loadable Function} {@var{d} =} __eigs__ (@var{af}, @var{n}, @var{k}, @var{sigma})\n\
147 @deftypefnx {Loadable Function} {@var{d} =} __eigs__ (@var{af}, @var{n}, @var{B}, @var{k}, @var{sigma})\n\
148 @deftypefnx {Loadable Function} {@var{d} =} __eigs__ (@var{af}, @var{n}, @var{k}, @var{sigma}, @var{opts})\n\
149 @deftypefnx {Loadable Function} {@var{d} =} __eigs__ (@var{af}, @var{n}, @var{B}, @var{k}, @var{sigma}, @var{opts})\n\
150 @deftypefnx {Loadable Function} {[@var{V}, @var{d}] =} __eigs__ (@var{A}, @dots{})\n\
151 @deftypefnx {Loadable Function} {[@var{V}, @var{d}] =} __eigs__ (@var{af}, @var{n}, @dots{})\n\
152 @deftypefnx {Loadable Function} {[@var{V}, @var{d}, @var{flag}] =} __eigs__ (@var{A}, @dots{})\n\
153 @deftypefnx {Loadable Function} {[@var{V}, @var{d}, @var{flag}] =} __eigs__ (@var{af}, @var{n}, @dots{})\n\
154 Undocumented internal function.\n\
155 @end deftypefn")
156 {
157  octave_value_list retval;
158 #ifdef HAVE_ARPACK
159  int nargin = args.length ();
160  std::string fcn_name;
161  octave_idx_type n = 0;
162  octave_idx_type k = 6;
163  Complex sigma = 0.;
164  double sigmar, sigmai;
165  bool have_sigma = false;
166  std::string typ = "LM";
167  Matrix amm, bmm, bmt;
168  ComplexMatrix acm, bcm, bct;
169  SparseMatrix asmm, bsmm, bsmt;
170  SparseComplexMatrix ascm, bscm, bsct;
171  int b_arg = 0;
172  bool have_b = false;
173  bool have_a_fun = false;
174  bool a_is_complex = false;
175  bool b_is_complex = false;
176  bool symmetric = false;
177  bool sym_tested = false;
178  bool cholB = false;
179  bool a_is_sparse = false;
180  ColumnVector permB;
181  int arg_offset = 0;
182  double tol = std::numeric_limits<double>::epsilon ();
183  int maxit = 300;
184  int disp = 0;
185  octave_idx_type p = -1;
186  ColumnVector resid;
187  ComplexColumnVector cresid;
188  octave_idx_type info = 1;
189 
190  warned_imaginary = false;
191 
192  unwind_protect frame;
193 
194  frame.protect_var (call_depth);
195  call_depth++;
196 
197  if (call_depth > 1)
198  {
199  error ("eigs: invalid recursive call");
200  if (fcn_name.length ())
201  clear_function (fcn_name);
202  return retval;
203  }
204 
205  if (nargin == 0)
206  print_usage ();
207  else if (args(0).is_function_handle () || args(0).is_inline_function ()
208  || args(0).is_string ())
209  {
210  if (args(0).is_string ())
211  {
212  std::string name = args(0).string_value ();
213  std::string fname = "function y = ";
214  fcn_name = unique_symbol_name ("__eigs_fcn_");
215  fname.append (fcn_name);
216  fname.append ("(x) y = ");
217  eigs_fcn = extract_function (args(0), "eigs", fcn_name, fname,
218  "; endfunction");
219  }
220  else
221  eigs_fcn = args(0).function_value ();
222 
223  if (!eigs_fcn)
224  {
225  error ("eigs: unknown function");
226  return retval;
227  }
228 
229  if (nargin < 2)
230  {
231  error ("eigs: incorrect number of arguments");
232  return retval;
233  }
234  else
235  {
236  n = args(1).nint_value ();
237  arg_offset = 1;
238  have_a_fun = true;
239  }
240  }
241  else
242  {
243  if (args(0).is_complex_type ())
244  {
245  if (args(0).is_sparse_type ())
246  {
247  ascm = (args(0).sparse_complex_matrix_value ());
248  a_is_sparse = true;
249  }
250  else
251  acm = (args(0).complex_matrix_value ());
252  a_is_complex = true;
253  symmetric = false; // ARPACK doesn't special case complex symmetric
254  sym_tested = true;
255  }
256  else
257  {
258  if (args(0).is_sparse_type ())
259  {
260  asmm = (args(0).sparse_matrix_value ());
261  a_is_sparse = true;
262  }
263  else
264  {
265  amm = (args(0).matrix_value ());
266  }
267  }
268 
269  }
270 
271  // Note hold off reading B till later to avoid issues of double
272  // copies of the matrix if B is full/real while A is complex.
273  if (!error_state && nargin > 1 + arg_offset &&
274  !(args(1 + arg_offset).is_real_scalar ()))
275  {
276  if (args(1+arg_offset).is_complex_type ())
277  {
278  b_arg = 1+arg_offset;
279  have_b = true;
280  b_is_complex = true;
281  arg_offset++;
282  }
283  else
284  {
285  b_arg = 1+arg_offset;
286  have_b = true;
287  arg_offset++;
288  }
289  }
290 
291  if (!error_state && nargin > (1+arg_offset))
292  k = args(1+arg_offset).nint_value ();
293 
294  if (!error_state && nargin > (2+arg_offset))
295  {
296  if (args(2+arg_offset).is_string ())
297  {
298  typ = args(2+arg_offset).string_value ();
299 
300  // Use STL function to convert to upper case
301  transform (typ.begin (), typ.end (), typ.begin (), toupper);
302 
303  sigma = 0.;
304  }
305  else
306  {
307  sigma = args(2+arg_offset).complex_value ();
308 
309  if (! error_state)
310  have_sigma = true;
311  else
312  {
313  error ("eigs: SIGMA must be a scalar or a string");
314  return retval;
315  }
316  }
317  }
318 
319  sigmar = std::real (sigma);
320  sigmai = std::imag (sigma);
321 
322  if (!error_state && nargin > (3+arg_offset))
323  {
324  if (args(3+arg_offset).is_map ())
325  {
326  octave_scalar_map map = args(3+arg_offset).scalar_map_value ();
327 
328  if (! error_state)
329  {
330  octave_value tmp;
331 
332  // issym is ignored for complex matrix inputs
333  tmp = map.getfield ("issym");
334  if (tmp.is_defined () && !sym_tested)
335  {
336  symmetric = tmp.double_value () != 0.;
337  sym_tested = true;
338  }
339 
340  // isreal is ignored if A is not a function
341  tmp = map.getfield ("isreal");
342  if (tmp.is_defined () && have_a_fun)
343  a_is_complex = ! (tmp.double_value () != 0.);
344 
345  tmp = map.getfield ("tol");
346  if (tmp.is_defined ())
347  tol = tmp.double_value ();
348 
349  tmp = map.getfield ("maxit");
350  if (tmp.is_defined ())
351  maxit = tmp.nint_value ();
352 
353  tmp = map.getfield ("p");
354  if (tmp.is_defined ())
355  p = tmp.nint_value ();
356 
357  tmp = map.getfield ("v0");
358  if (tmp.is_defined ())
359  {
360  if (a_is_complex || b_is_complex)
361  cresid = ComplexColumnVector (tmp.complex_vector_value ());
362  else
363  resid = ColumnVector (tmp.vector_value ());
364  }
365 
366  tmp = map.getfield ("disp");
367  if (tmp.is_defined ())
368  disp = tmp.nint_value ();
369 
370  tmp = map.getfield ("cholB");
371  if (tmp.is_defined ())
372  cholB = tmp.double_value () != 0.;
373 
374  tmp = map.getfield ("permB");
375  if (tmp.is_defined ())
376  permB = ColumnVector (tmp.vector_value ()) - 1.0;
377  }
378  else
379  {
380  error ("eigs: OPTS argument must be a scalar structure");
381  return retval;
382  }
383  }
384  else
385  {
386  error ("eigs: OPTS argument must be a structure");
387  return retval;
388  }
389  }
390 
391  if (nargin > (4+arg_offset))
392  {
393  error ("eigs: incorrect number of arguments");
394  return retval;
395  }
396 
397  // Test undeclared (no issym) matrix inputs for symmetry
398  if (!sym_tested && !have_a_fun)
399  {
400  if (a_is_sparse)
401  symmetric = asmm.is_symmetric ();
402  else
403  symmetric = amm.is_symmetric ();
404  }
405 
406  if (have_b)
407  {
408  if (a_is_complex || b_is_complex)
409  {
410  if (a_is_sparse)
411  bscm = args(b_arg).sparse_complex_matrix_value ();
412  else
413  bcm = args(b_arg).complex_matrix_value ();
414  }
415  else
416  {
417  if (a_is_sparse)
418  bsmm = args(b_arg).sparse_matrix_value ();
419  else
420  bmm = args(b_arg).matrix_value ();
421  }
422  }
423 
424  // Mode 1 for SM mode seems unstable for some reason.
425  // Use Mode 3 instead, with sigma = 0.
426  if (!error_state && !have_sigma && typ == "SM")
427  have_sigma = true;
428 
429  if (!error_state)
430  {
431  octave_idx_type nconv;
432  if (a_is_complex || b_is_complex)
433  {
434  ComplexMatrix eig_vec;
435  ComplexColumnVector eig_val;
436 
437 
438  if (have_a_fun)
440  (eigs_complex_func, n, typ, sigma, k, p, info, eig_vec,
441  eig_val, cresid, octave_stdout, tol, (nargout > 1), cholB,
442  disp, maxit);
443  else if (have_sigma)
444  {
445  if (a_is_sparse)
447  (ascm, sigma, k, p, info, eig_vec, eig_val, bscm, permB,
448  cresid, octave_stdout, tol, (nargout > 1), cholB, disp,
449  maxit);
450  else
452  (acm, sigma, k, p, info, eig_vec, eig_val, bcm, permB,
453  cresid, octave_stdout, tol, (nargout > 1), cholB, disp,
454  maxit);
455  }
456  else
457  {
458  if (a_is_sparse)
460  (ascm, typ, k, p, info, eig_vec, eig_val, bscm, permB,
461  cresid, octave_stdout, tol, (nargout > 1), cholB, disp,
462  maxit);
463  else
465  (acm, typ, k, p, info, eig_vec, eig_val, bcm, permB,
466  cresid, octave_stdout, tol, (nargout > 1), cholB, disp,
467  maxit);
468  }
469 
470  if (nargout < 2)
471  retval(0) = eig_val;
472  else
473  {
474  retval(2) = double (info);
475  retval(1) = ComplexDiagMatrix (eig_val);
476  retval(0) = eig_vec;
477  }
478  }
479  else if (sigmai != 0.)
480  {
481  // Promote real problem to a complex one.
482  ComplexMatrix eig_vec;
483  ComplexColumnVector eig_val;
484 
485  if (have_a_fun)
487  (eigs_complex_func, n, typ, sigma, k, p, info, eig_vec,
488  eig_val, cresid, octave_stdout, tol, (nargout > 1), cholB,
489  disp, maxit);
490  else
491  {
492  if (a_is_sparse)
494  (SparseComplexMatrix (asmm), sigma, k, p, info, eig_vec,
495  eig_val, SparseComplexMatrix (bsmm), permB, cresid,
496  octave_stdout, tol, (nargout > 1), cholB, disp, maxit);
497  else
499  (ComplexMatrix (amm), sigma, k, p, info, eig_vec,
500  eig_val, ComplexMatrix (bmm), permB, cresid,
501  octave_stdout, tol, (nargout > 1), cholB, disp, maxit);
502  }
503 
504  if (nargout < 2)
505  retval(0) = eig_val;
506  else
507  {
508  retval(2) = double (info);
509  retval(1) = ComplexDiagMatrix (eig_val);
510  retval(0) = eig_vec;
511  }
512  }
513  else
514  {
515  if (symmetric)
516  {
517  Matrix eig_vec;
518  ColumnVector eig_val;
519 
520  if (have_a_fun)
521  nconv = EigsRealSymmetricFunc
522  (eigs_func, n, typ, sigmar, k, p, info, eig_vec,
523  eig_val, resid, octave_stdout, tol, (nargout > 1),
524  cholB, disp, maxit);
525  else if (have_sigma)
526  {
527  if (a_is_sparse)
529  (asmm, sigmar, k, p, info, eig_vec, eig_val, bsmm,
530  permB, resid, octave_stdout, tol, (nargout > 1),
531  cholB, disp, maxit);
532  else
534  (amm, sigmar, k, p, info, eig_vec, eig_val, bmm,
535  permB, resid, octave_stdout, tol, (nargout > 1),
536  cholB, disp, maxit);
537  }
538  else
539  {
540  if (a_is_sparse)
542  (asmm, typ, k, p, info, eig_vec, eig_val, bsmm,
543  permB, resid, octave_stdout, tol, (nargout > 1),
544  cholB, disp, maxit);
545  else
547  (amm, typ, k, p, info, eig_vec, eig_val, bmm, permB,
548  resid, octave_stdout, tol, (nargout > 1), cholB,
549  disp, maxit);
550  }
551 
552  if (nargout < 2)
553  retval(0) = eig_val;
554  else
555  {
556  retval(2) = double (info);
557  retval(1) = DiagMatrix (eig_val);
558  retval(0) = eig_vec;
559  }
560  }
561  else
562  {
563  ComplexMatrix eig_vec;
564  ComplexColumnVector eig_val;
565 
566  if (have_a_fun)
568  (eigs_func, n, typ, sigmar, k, p, info, eig_vec,
569  eig_val, resid, octave_stdout, tol, (nargout > 1),
570  cholB, disp, maxit);
571  else if (have_sigma)
572  {
573  if (a_is_sparse)
575  (asmm, sigmar, k, p, info, eig_vec, eig_val, bsmm,
576  permB, resid, octave_stdout, tol, (nargout > 1),
577  cholB, disp, maxit);
578  else
580  (amm, sigmar, k, p, info, eig_vec, eig_val, bmm,
581  permB, resid, octave_stdout, tol, (nargout > 1),
582  cholB, disp, maxit);
583  }
584  else
585  {
586  if (a_is_sparse)
588  (asmm, typ, k, p, info, eig_vec, eig_val, bsmm,
589  permB, resid, octave_stdout, tol, (nargout > 1),
590  cholB, disp, maxit);
591  else
593  (amm, typ, k, p, info, eig_vec, eig_val, bmm, permB,
594  resid, octave_stdout, tol, (nargout > 1), cholB,
595  disp, maxit);
596  }
597 
598  if (nargout < 2)
599  retval(0) = eig_val;
600  else
601  {
602  retval(2) = double (info);
603  retval(1) = ComplexDiagMatrix (eig_val);
604  retval(0) = eig_vec;
605  }
606  }
607  }
608 
609  if (nconv <= 0)
610  warning ("eigs: None of the %d requested eigenvalues converged", k);
611  else if (nconv < k)
612  warning ("eigs: Only %d of the %d requested eigenvalues converged",
613  nconv, k);
614  }
615 
616  if (! fcn_name.empty ())
617  clear_function (fcn_name);
618 #else
619  error ("eigs: not available in this version of Octave");
620 #endif
621 
622  return retval;
623 }