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