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