Mercurial > octave-nkf

comparison liboctave/eigs-base.cc @ 11586:12df7854fa7c

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
strip trailing whitespace from source files
author John W. Eaton <jwe@octave.org>
date Thu, 20 Jan 2011 17:24:59 -0500
parents 57632dea2446
children d25dfa9ed18b
comparison
equal deleted inserted replaced
11585:1473d0cf86d2 11586:12df7854fa7c
44 #include "dbleLU.h" 44 #include "dbleLU.h"
45 #include "CmplxLU.h" 45 #include "CmplxLU.h"
46 46
47 #ifdef HAVE_ARPACK 47 #ifdef HAVE_ARPACK
48 typedef ColumnVector (*EigsFunc) (const ColumnVector &x, int &eigs_error); 48 typedef ColumnVector (*EigsFunc) (const ColumnVector &x, int &eigs_error);
49 typedef ComplexColumnVector (*EigsComplexFunc) 49 typedef ComplexColumnVector (*EigsComplexFunc)
50 (const ComplexColumnVector &x, int &eigs_error); 50 (const ComplexColumnVector &x, int &eigs_error);
51 51
52 // Arpack and blas fortran functions we call. 52 // Arpack and blas fortran functions we call.
53 extern "C" 53 extern "C"
54 { 54 {
55 F77_RET_T 55 F77_RET_T
56 F77_FUNC (dsaupd, DSAUPD) (octave_idx_type&, 56 F77_FUNC (dsaupd, DSAUPD) (octave_idx_type&,
57 F77_CONST_CHAR_ARG_DECL, 57 F77_CONST_CHAR_ARG_DECL,
58 const octave_idx_type&, 58 const octave_idx_type&,
59 F77_CONST_CHAR_ARG_DECL, 59 F77_CONST_CHAR_ARG_DECL,
60 const octave_idx_type&, const double&, 60 const octave_idx_type&, const double&,
61 double*, const octave_idx_type&, double*, 61 double*, const octave_idx_type&, double*,
62 const octave_idx_type&, octave_idx_type*, 62 const octave_idx_type&, octave_idx_type*,
63 octave_idx_type*, double*, double*, 63 octave_idx_type*, double*, double*,
64 const octave_idx_type&, octave_idx_type& 64 const octave_idx_type&, octave_idx_type&
65 F77_CHAR_ARG_LEN_DECL 65 F77_CHAR_ARG_LEN_DECL
66 F77_CHAR_ARG_LEN_DECL); 66 F77_CHAR_ARG_LEN_DECL);
67 67
68 F77_RET_T 68 F77_RET_T
69 F77_FUNC (dseupd, DSEUPD) (const octave_idx_type&, 69 F77_FUNC (dseupd, DSEUPD) (const octave_idx_type&,
70 F77_CONST_CHAR_ARG_DECL, 70 F77_CONST_CHAR_ARG_DECL,
71 octave_idx_type*, double*, double*, 71 octave_idx_type*, double*, double*,
72 const octave_idx_type&, const double&, 72 const octave_idx_type&, const double&,
73 F77_CONST_CHAR_ARG_DECL, 73 F77_CONST_CHAR_ARG_DECL,
74 const octave_idx_type&, 74 const octave_idx_type&,
75 F77_CONST_CHAR_ARG_DECL, 75 F77_CONST_CHAR_ARG_DECL,
76 const octave_idx_type&, const double&, double*, 76 const octave_idx_type&, const double&, double*,
77 const octave_idx_type&, double*, 77 const octave_idx_type&, double*,
78 const octave_idx_type&, octave_idx_type*, 78 const octave_idx_type&, octave_idx_type*,
79 octave_idx_type*, double*, double*, 79 octave_idx_type*, double*, double*,
80 const octave_idx_type&, octave_idx_type& 80 const octave_idx_type&, octave_idx_type&
81 F77_CHAR_ARG_LEN_DECL 81 F77_CHAR_ARG_LEN_DECL
82 F77_CHAR_ARG_LEN_DECL 82 F77_CHAR_ARG_LEN_DECL
83 F77_CHAR_ARG_LEN_DECL); 83 F77_CHAR_ARG_LEN_DECL);
84 84
85 F77_RET_T 85 F77_RET_T
86 F77_FUNC (dnaupd, DNAUPD) (octave_idx_type&, 86 F77_FUNC (dnaupd, DNAUPD) (octave_idx_type&,
87 F77_CONST_CHAR_ARG_DECL, 87 F77_CONST_CHAR_ARG_DECL,
88 const octave_idx_type&, 88 const octave_idx_type&,
89 F77_CONST_CHAR_ARG_DECL, 89 F77_CONST_CHAR_ARG_DECL,
90 octave_idx_type&, const double&, 90 octave_idx_type&, const double&,
91 double*, const octave_idx_type&, double*, 91 double*, const octave_idx_type&, double*,
92 const octave_idx_type&, octave_idx_type*, 92 const octave_idx_type&, octave_idx_type*,
93 octave_idx_type*, double*, double*, 93 octave_idx_type*, double*, double*,
94 const octave_idx_type&, octave_idx_type& 94 const octave_idx_type&, octave_idx_type&
95 F77_CHAR_ARG_LEN_DECL 95 F77_CHAR_ARG_LEN_DECL
96 F77_CHAR_ARG_LEN_DECL); 96 F77_CHAR_ARG_LEN_DECL);
97 97
98 F77_RET_T 98 F77_RET_T
99 F77_FUNC (dneupd, DNEUPD) (const octave_idx_type&, 99 F77_FUNC (dneupd, DNEUPD) (const octave_idx_type&,
100 F77_CONST_CHAR_ARG_DECL, 100 F77_CONST_CHAR_ARG_DECL,
101 octave_idx_type*, double*, double*, 101 octave_idx_type*, double*, double*,
102 double*, const octave_idx_type&, const double&, 102 double*, const octave_idx_type&, const double&,
103 const double&, double*, 103 const double&, double*,
104 F77_CONST_CHAR_ARG_DECL, 104 F77_CONST_CHAR_ARG_DECL,
105 const octave_idx_type&, 105 const octave_idx_type&,
106 F77_CONST_CHAR_ARG_DECL, 106 F77_CONST_CHAR_ARG_DECL,
107 octave_idx_type&, const double&, double*, 107 octave_idx_type&, const double&, double*,
108 const octave_idx_type&, double*, 108 const octave_idx_type&, double*,
109 const octave_idx_type&, octave_idx_type*, 109 const octave_idx_type&, octave_idx_type*,
110 octave_idx_type*, double*, double*, 110 octave_idx_type*, double*, double*,
111 const octave_idx_type&, octave_idx_type& 111 const octave_idx_type&, octave_idx_type&
112 F77_CHAR_ARG_LEN_DECL 112 F77_CHAR_ARG_LEN_DECL
113 F77_CHAR_ARG_LEN_DECL 113 F77_CHAR_ARG_LEN_DECL
114 F77_CHAR_ARG_LEN_DECL); 114 F77_CHAR_ARG_LEN_DECL);
115 115
119 const octave_idx_type&, 119 const octave_idx_type&,
120 F77_CONST_CHAR_ARG_DECL, 120 F77_CONST_CHAR_ARG_DECL,
121 const octave_idx_type&, const double&, 121 const octave_idx_type&, const double&,
122 Complex*, const octave_idx_type&, Complex*, 122 Complex*, const octave_idx_type&, Complex*,
123 const octave_idx_type&, octave_idx_type*, 123 const octave_idx_type&, octave_idx_type*,
124 octave_idx_type*, Complex*, Complex*, 124 octave_idx_type*, Complex*, Complex*,
125 const octave_idx_type&, double *, octave_idx_type& 125 const octave_idx_type&, double *, octave_idx_type&
126 F77_CHAR_ARG_LEN_DECL 126 F77_CHAR_ARG_LEN_DECL
127 F77_CHAR_ARG_LEN_DECL); 127 F77_CHAR_ARG_LEN_DECL);
128 128
129 F77_RET_T 129 F77_RET_T
130 F77_FUNC (zneupd, ZNEUPD) (const octave_idx_type&, 130 F77_FUNC (zneupd, ZNEUPD) (const octave_idx_type&,
131 F77_CONST_CHAR_ARG_DECL, 131 F77_CONST_CHAR_ARG_DECL,
132 octave_idx_type*, Complex*, Complex*, 132 octave_idx_type*, Complex*, Complex*,
133 const octave_idx_type&, const Complex&, Complex*, 133 const octave_idx_type&, const Complex&, Complex*,
134 F77_CONST_CHAR_ARG_DECL, 134 F77_CONST_CHAR_ARG_DECL,
135 const octave_idx_type&, 135 const octave_idx_type&,
136 F77_CONST_CHAR_ARG_DECL, 136 F77_CONST_CHAR_ARG_DECL,
137 const octave_idx_type&, const double&, 137 const octave_idx_type&, const double&,
138 Complex*, const octave_idx_type&, Complex*, 138 Complex*, const octave_idx_type&, Complex*,
139 const octave_idx_type&, octave_idx_type*, 139 const octave_idx_type&, octave_idx_type*,
140 octave_idx_type*, Complex*, Complex*, 140 octave_idx_type*, Complex*, Complex*,
141 const octave_idx_type&, double *, octave_idx_type& 141 const octave_idx_type&, double *, octave_idx_type&
142 F77_CHAR_ARG_LEN_DECL 142 F77_CHAR_ARG_LEN_DECL
143 F77_CHAR_ARG_LEN_DECL 143 F77_CHAR_ARG_LEN_DECL
144 F77_CHAR_ARG_LEN_DECL); 144 F77_CHAR_ARG_LEN_DECL);
145 145
166 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL) 166 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
167 static octave_idx_type 167 static octave_idx_type
168 lusolve (const SparseMatrix&, const SparseMatrix&, Matrix&); 168 lusolve (const SparseMatrix&, const SparseMatrix&, Matrix&);
169 169
170 static octave_idx_type 170 static octave_idx_type
171 lusolve (const SparseComplexMatrix&, const SparseComplexMatrix&, 171 lusolve (const SparseComplexMatrix&, const SparseComplexMatrix&,
172 ComplexMatrix&); 172 ComplexMatrix&);
173 173
174 static octave_idx_type 174 static octave_idx_type
175 lusolve (const Matrix&, const Matrix&, Matrix&); 175 lusolve (const Matrix&, const Matrix&, Matrix&);
176 176
177 static octave_idx_type 177 static octave_idx_type
178 lusolve (const ComplexMatrix&, const ComplexMatrix&, ComplexMatrix&); 178 lusolve (const ComplexMatrix&, const ComplexMatrix&, ComplexMatrix&);
179 179
180 static ComplexMatrix 180 static ComplexMatrix
181 ltsolve (const SparseComplexMatrix&, const ColumnVector&, 181 ltsolve (const SparseComplexMatrix&, const ColumnVector&,
182 const ComplexMatrix&); 182 const ComplexMatrix&);
183 183
184 static Matrix 184 static Matrix
185 ltsolve (const SparseMatrix&, const ColumnVector&, const Matrix&,); 185 ltsolve (const SparseMatrix&, const ColumnVector&, const Matrix&,);
186 186
239 { 239 {
240 retval.resize (n, b_nc); 240 retval.resize (n, b_nc);
241 for (octave_idx_type j = 0; j < b_nc; j++) 241 for (octave_idx_type j = 0; j < b_nc; j++)
242 { 242 {
243 for (octave_idx_type i = 0; i < n; i++) 243 for (octave_idx_type i = 0; i < n; i++)
244 retval.elem(static_cast<octave_idx_type>(qv[i]), j) = 244 retval.elem(static_cast<octave_idx_type>(qv[i]), j) =
245 tmp.elem(i,j); 245 tmp.elem(i,j);
246 } 246 }
247 } 247 }
248 248
249 return retval; 249 return retval;
295 x, 1, 0.0, y, 1 295 x, 1, 0.0, y, 1
296 F77_CHAR_ARG_LEN (1))); 296 F77_CHAR_ARG_LEN (1)));
297 297
298 if (f77_exception_encountered) 298 if (f77_exception_encountered)
299 { 299 {
300 (*current_liboctave_error_handler) 300 (*current_liboctave_error_handler)
301 ("eigs: unrecoverable error in dgemv"); 301 ("eigs: unrecoverable error in dgemv");
302 return false; 302 return false;
303 } 303 }
304 else 304 else
305 return true; 305 return true;
306 } 306 }
307 307
308 static bool 308 static bool
309 vector_product (const SparseComplexMatrix& m, const Complex* x, 309 vector_product (const SparseComplexMatrix& m, const Complex* x,
310 Complex* y) 310 Complex* y)
311 { 311 {
312 octave_idx_type nc = m.cols (); 312 octave_idx_type nc = m.cols ();
313 313
314 for (octave_idx_type j = 0; j < nc; j++) 314 for (octave_idx_type j = 0; j < nc; j++)
332 x, 1, 0.0, y, 1 332 x, 1, 0.0, y, 1
333 F77_CHAR_ARG_LEN (1))); 333 F77_CHAR_ARG_LEN (1)));
334 334
335 if (f77_exception_encountered) 335 if (f77_exception_encountered)
336 { 336 {
337 (*current_liboctave_error_handler) 337 (*current_liboctave_error_handler)
338 ("eigs: unrecoverable error in zgemv"); 338 ("eigs: unrecoverable error in zgemv");
339 return false; 339 return false;
340 } 340 }
341 else 341 else
342 return true; 342 return true;
398 return true; 398 return true;
399 } 399 }
400 } 400 }
401 401
402 static bool 402 static bool
403 make_cholb (SparseComplexMatrix& b, SparseComplexMatrix& bt, 403 make_cholb (SparseComplexMatrix& b, SparseComplexMatrix& bt,
404 ColumnVector& permB) 404 ColumnVector& permB)
405 { 405 {
406 octave_idx_type info; 406 octave_idx_type info;
407 SparseComplexCHOL fact (b, info, false); 407 SparseComplexCHOL fact (b, info, false);
408 408
416 return true; 416 return true;
417 } 417 }
418 } 418 }
419 419
420 static bool 420 static bool
421 LuAminusSigmaB (const SparseMatrix &m, const SparseMatrix &b, 421 LuAminusSigmaB (const SparseMatrix &m, const SparseMatrix &b,
422 bool cholB, const ColumnVector& permB, double sigma, 422 bool cholB, const ColumnVector& permB, double sigma,
423 SparseMatrix &L, SparseMatrix &U, octave_idx_type *P, 423 SparseMatrix &L, SparseMatrix &U, octave_idx_type *P,
424 octave_idx_type *Q) 424 octave_idx_type *Q)
425 { 425 {
426 bool have_b = ! b.is_empty (); 426 bool have_b = ! b.is_empty ();
427 octave_idx_type n = m.rows(); 427 octave_idx_type n = m.rows();
428 428
437 { 437 {
438 SparseMatrix tmp(n,n,n); 438 SparseMatrix tmp(n,n,n);
439 for (octave_idx_type i = 0; i < n; i++) 439 for (octave_idx_type i = 0; i < n; i++)
440 { 440 {
441 tmp.xcidx(i) = i; 441 tmp.xcidx(i) = i;
442 tmp.xridx(i) = 442 tmp.xridx(i) =
443 static_cast<octave_idx_type>(permB(i)); 443 static_cast<octave_idx_type>(permB(i));
444 tmp.xdata(i) = 1; 444 tmp.xdata(i) = 1;
445 } 445 }
446 tmp.xcidx(n) = n; 446 tmp.xcidx(n) = n;
447 447
514 514
515 return true; 515 return true;
516 } 516 }
517 517
518 static bool 518 static bool
519 LuAminusSigmaB (const Matrix &m, const Matrix &b, 519 LuAminusSigmaB (const Matrix &m, const Matrix &b,
520 bool cholB, const ColumnVector& permB, double sigma, 520 bool cholB, const ColumnVector& permB, double sigma,
521 Matrix &L, Matrix &U, octave_idx_type *P, 521 Matrix &L, Matrix &U, octave_idx_type *P,
522 octave_idx_type *Q) 522 octave_idx_type *Q)
523 { 523 {
524 bool have_b = ! b.is_empty (); 524 bool have_b = ! b.is_empty ();
525 octave_idx_type n = m.cols(); 525 octave_idx_type n = m.cols();
526 526
535 const double *pB = permB.fortran_vec(); 535 const double *pB = permB.fortran_vec();
536 double *p = AminusSigmaB.fortran_vec(); 536 double *p = AminusSigmaB.fortran_vec();
537 537
538 if (permB.length()) 538 if (permB.length())
539 { 539 {
540 for (octave_idx_type j = 0; 540 for (octave_idx_type j = 0;
541 j < b.cols(); j++) 541 j < b.cols(); j++)
542 for (octave_idx_type i = 0; 542 for (octave_idx_type i = 0;
543 i < b.rows(); i++) 543 i < b.rows(); i++)
544 *p++ -= tmp.xelem (static_cast<octave_idx_type>(pB[i]), 544 *p++ -= tmp.xelem (static_cast<octave_idx_type>(pB[i]),
545 static_cast<octave_idx_type>(pB[j])); 545 static_cast<octave_idx_type>(pB[j]));
546 } 546 }
547 else 547 else
561 LU fact (AminusSigmaB); 561 LU fact (AminusSigmaB);
562 562
563 L = fact.P().transpose() * fact.L (); 563 L = fact.P().transpose() * fact.L ();
564 U = fact.U (); 564 U = fact.U ();
565 for (octave_idx_type j = 0; j < n; j++) 565 for (octave_idx_type j = 0; j < n; j++)
566 P[j] = Q[j] = j; 566 P[j] = Q[j] = j;
567 567
568 // Test condition number of LU decomposition 568 // Test condition number of LU decomposition
569 double minU = octave_NaN; 569 double minU = octave_NaN;
570 double maxU = octave_NaN; 570 double maxU = octave_NaN;
571 for (octave_idx_type j = 0; j < n; j++) 571 for (octave_idx_type j = 0; j < n; j++)
581 double rcond = (minU / maxU); 581 double rcond = (minU / maxU);
582 volatile double rcond_plus_one = rcond + 1.0; 582 volatile double rcond_plus_one = rcond + 1.0;
583 583
584 if (rcond_plus_one == 1.0 || xisnan (rcond)) 584 if (rcond_plus_one == 1.0 || xisnan (rcond))
585 { 585 {
586 (*current_liboctave_warning_handler) 586 (*current_liboctave_warning_handler)
587 ("eigs: 'A - sigma*B' is singular, indicating sigma is exactly"); 587 ("eigs: 'A - sigma*B' is singular, indicating sigma is exactly");
588 (*current_liboctave_warning_handler) 588 (*current_liboctave_warning_handler)
589 (" an eigenvalue. Convergence is not guaranteed"); 589 (" an eigenvalue. Convergence is not guaranteed");
590 } 590 }
591 591
592 return true; 592 return true;
593 } 593 }
594 594
595 static bool 595 static bool
596 LuAminusSigmaB (const SparseComplexMatrix &m, const SparseComplexMatrix &b, 596 LuAminusSigmaB (const SparseComplexMatrix &m, const SparseComplexMatrix &b,
597 bool cholB, const ColumnVector& permB, Complex sigma, 597 bool cholB, const ColumnVector& permB, Complex sigma,
598 SparseComplexMatrix &L, SparseComplexMatrix &U, 598 SparseComplexMatrix &L, SparseComplexMatrix &U,
599 octave_idx_type *P, octave_idx_type *Q) 599 octave_idx_type *P, octave_idx_type *Q)
600 { 600 {
601 bool have_b = ! b.is_empty (); 601 bool have_b = ! b.is_empty ();
612 { 612 {
613 SparseMatrix tmp(n,n,n); 613 SparseMatrix tmp(n,n,n);
614 for (octave_idx_type i = 0; i < n; i++) 614 for (octave_idx_type i = 0; i < n; i++)
615 { 615 {
616 tmp.xcidx(i) = i; 616 tmp.xcidx(i) = i;
617 tmp.xridx(i) = 617 tmp.xridx(i) =
618 static_cast<octave_idx_type>(permB(i)); 618 static_cast<octave_idx_type>(permB(i));
619 tmp.xdata(i) = 1; 619 tmp.xdata(i) = 1;
620 } 620 }
621 tmp.xcidx(n) = n; 621 tmp.xcidx(n) = n;
622 622
623 AminusSigmaB = AminusSigmaB - tmp * b.hermitian() * b * 623 AminusSigmaB = AminusSigmaB - tmp * b.hermitian() * b *
624 tmp.transpose() * sigma; 624 tmp.transpose() * sigma;
625 } 625 }
626 else 626 else
627 AminusSigmaB = AminusSigmaB - sigma * b.hermitian() * b; 627 AminusSigmaB = AminusSigmaB - sigma * b.hermitian() * b;
628 } 628 }
688 688
689 return true; 689 return true;
690 } 690 }
691 691
692 static bool 692 static bool
693 LuAminusSigmaB (const ComplexMatrix &m, const ComplexMatrix &b, 693 LuAminusSigmaB (const ComplexMatrix &m, const ComplexMatrix &b,
694 bool cholB, const ColumnVector& permB, Complex sigma, 694 bool cholB, const ColumnVector& permB, Complex sigma,
695 ComplexMatrix &L, ComplexMatrix &U, octave_idx_type *P, 695 ComplexMatrix &L, ComplexMatrix &U, octave_idx_type *P,
696 octave_idx_type *Q) 696 octave_idx_type *Q)
697 { 697 {
698 bool have_b = ! b.is_empty (); 698 bool have_b = ! b.is_empty ();
699 octave_idx_type n = m.cols(); 699 octave_idx_type n = m.cols();
700 700
709 const double *pB = permB.fortran_vec(); 709 const double *pB = permB.fortran_vec();
710 Complex *p = AminusSigmaB.fortran_vec(); 710 Complex *p = AminusSigmaB.fortran_vec();
711 711
712 if (permB.length()) 712 if (permB.length())
713 { 713 {
714 for (octave_idx_type j = 0; 714 for (octave_idx_type j = 0;
715 j < b.cols(); j++) 715 j < b.cols(); j++)
716 for (octave_idx_type i = 0; 716 for (octave_idx_type i = 0;
717 i < b.rows(); i++) 717 i < b.rows(); i++)
718 *p++ -= tmp.xelem (static_cast<octave_idx_type>(pB[i]), 718 *p++ -= tmp.xelem (static_cast<octave_idx_type>(pB[i]),
719 static_cast<octave_idx_type>(pB[j])); 719 static_cast<octave_idx_type>(pB[j]));
720 } 720 }
721 else 721 else
735 ComplexLU fact (AminusSigmaB); 735 ComplexLU fact (AminusSigmaB);
736 736
737 L = fact.P().transpose() * fact.L (); 737 L = fact.P().transpose() * fact.L ();
738 U = fact.U (); 738 U = fact.U ();
739 for (octave_idx_type j = 0; j < n; j++) 739 for (octave_idx_type j = 0; j < n; j++)
740 P[j] = Q[j] = j; 740 P[j] = Q[j] = j;
741 741
742 // Test condition number of LU decomposition 742 // Test condition number of LU decomposition
743 double minU = octave_NaN; 743 double minU = octave_NaN;
744 double maxU = octave_NaN; 744 double maxU = octave_NaN;
745 for (octave_idx_type j = 0; j < n; j++) 745 for (octave_idx_type j = 0; j < n; j++)
755 double rcond = (minU / maxU); 755 double rcond = (minU / maxU);
756 volatile double rcond_plus_one = rcond + 1.0; 756 volatile double rcond_plus_one = rcond + 1.0;
757 757
758 if (rcond_plus_one == 1.0 || xisnan (rcond)) 758 if (rcond_plus_one == 1.0 || xisnan (rcond))
759 { 759 {
760 (*current_liboctave_warning_handler) 760 (*current_liboctave_warning_handler)
761 ("eigs: 'A - sigma*B' is singular, indicating sigma is exactly"); 761 ("eigs: 'A - sigma*B' is singular, indicating sigma is exactly");
762 (*current_liboctave_warning_handler) 762 (*current_liboctave_warning_handler)
763 (" an eigenvalue. Convergence is not guaranteed"); 763 (" an eigenvalue. Convergence is not guaranteed");
764 } 764 }
765 765
766 return true; 766 return true;
767 } 767 }
768 768
769 template <class M> 769 template <class M>
770 octave_idx_type 770 octave_idx_type
771 EigsRealSymmetricMatrix (const M& m, const std::string typ, 771 EigsRealSymmetricMatrix (const M& m, const std::string typ,
772 octave_idx_type k, octave_idx_type p, 772 octave_idx_type k, octave_idx_type p,
773 octave_idx_type &info, Matrix &eig_vec, 773 octave_idx_type &info, Matrix &eig_vec,
774 ColumnVector &eig_val, const M& _b, 774 ColumnVector &eig_val, const M& _b,
775 ColumnVector &permB, ColumnVector &resid, 775 ColumnVector &permB, ColumnVector &resid,
776 std::ostream& os, double tol, bool rvec, 776 std::ostream& os, double tol, bool rvec,
777 bool cholB, int disp, int maxit) 777 bool cholB, int disp, int maxit)
778 { 778 {
779 M b(_b); 779 M b(_b);
780 octave_idx_type n = m.cols (); 780 octave_idx_type n = m.cols ();
781 octave_idx_type mode = 1; 781 octave_idx_type mode = 1;
790 (*current_liboctave_error_handler) ("eigs: A must be square"); 790 (*current_liboctave_error_handler) ("eigs: A must be square");
791 return -1; 791 return -1;
792 } 792 }
793 if (have_b && (m.rows() != b.rows() || m.rows() != b.cols())) 793 if (have_b && (m.rows() != b.rows() || m.rows() != b.cols()))
794 { 794 {
795 (*current_liboctave_error_handler) 795 (*current_liboctave_error_handler)
796 ("eigs: B must be square and the same size as A"); 796 ("eigs: B must be square and the same size as A");
797 return -1; 797 return -1;
798 } 798 }
799 799
800 if (resid.is_empty()) 800 if (resid.is_empty())
816 { 816 {
817 p = k * 2; 817 p = k * 2;
818 818
819 if (p < 20) 819 if (p < 20)
820 p = 20; 820 p = 20;
821 821
822 if (p > n - 1) 822 if (p > n - 1)
823 p = n - 1 ; 823 p = n - 1 ;
824 } 824 }
825 825
826 if (k < 1 || k > n - 2) 826 if (k < 1 || k > n - 2)
827 { 827 {
828 (*current_liboctave_error_handler) 828 (*current_liboctave_error_handler)
829 ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1-1).\n" 829 ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1-1).\n"
830 " Use 'eig(full(A))' instead"); 830 " Use 'eig(full(A))' instead");
831 return -1; 831 return -1;
832 } 832 }
833 833
834 if (p <= k || p >= n) 834 if (p <= k || p >= n)
835 { 835 {
836 (*current_liboctave_error_handler) 836 (*current_liboctave_error_handler)
837 ("eigs: opts.p must be greater than k and less than n"); 837 ("eigs: opts.p must be greater than k and less than n");
838 return -1; 838 return -1;
839 } 839 }
840 840
841 if (have_b && cholB && permB.length() != 0) 841 if (have_b && cholB && permB.length() != 0)
842 { 842 {
843 // Check the we really have a permutation vector 843 // Check the we really have a permutation vector
844 if (permB.length() != n) 844 if (permB.length() != n)
845 { 845 {
846 (*current_liboctave_error_handler) 846 (*current_liboctave_error_handler)
847 ("eigs: permB vector invalid"); 847 ("eigs: permB vector invalid");
848 return -1; 848 return -1;
849 } 849 }
850 else 850 else
851 { 851 {
852 Array<bool> checked (dim_vector (n, 1), false); 852 Array<bool> checked (dim_vector (n, 1), false);
853 for (octave_idx_type i = 0; i < n; i++) 853 for (octave_idx_type i = 0; i < n; i++)
854 { 854 {
855 octave_idx_type bidx = 855 octave_idx_type bidx =
856 static_cast<octave_idx_type> (permB(i)); 856 static_cast<octave_idx_type> (permB(i));
857 if (checked(bidx) || bidx < 0 || 857 if (checked(bidx) || bidx < 0 ||
858 bidx >= n || D_NINT (bidx) != bidx) 858 bidx >= n || D_NINT (bidx) != bidx)
859 { 859 {
860 (*current_liboctave_error_handler) 860 (*current_liboctave_error_handler)
861 ("eigs: permB vector invalid"); 861 ("eigs: permB vector invalid");
862 return -1; 862 return -1;
863 } 863 }
864 } 864 }
865 } 865 }
866 } 866 }
867 867
868 if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" && 868 if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" &&
869 typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" && 869 typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" &&
870 typ != "SI") 870 typ != "SI")
871 { 871 {
872 (*current_liboctave_error_handler) 872 (*current_liboctave_error_handler)
873 ("eigs: unrecognized sigma value"); 873 ("eigs: unrecognized sigma value");
874 return -1; 874 return -1;
875 } 875 }
876 876
877 if (typ == "LI" || typ == "SI" || typ == "LR" || typ == "SR") 877 if (typ == "LI" || typ == "SI" || typ == "LR" || typ == "SR")
878 { 878 {
879 (*current_liboctave_error_handler) 879 (*current_liboctave_error_handler)
880 ("eigs: invalid sigma value for real symmetric problem"); 880 ("eigs: invalid sigma value for real symmetric problem");
881 return -1; 881 return -1;
882 } 882 }
883 883
884 if (have_b) 884 if (have_b)
898 } 898 }
899 else 899 else
900 { 900 {
901 if (! make_cholb(b, bt, permB)) 901 if (! make_cholb(b, bt, permB))
902 { 902 {
903 (*current_liboctave_error_handler) 903 (*current_liboctave_error_handler)
904 ("eigs: The matrix B is not positive definite"); 904 ("eigs: The matrix B is not positive definite");
905 return -1; 905 return -1;
906 } 906 }
907 } 907 }
908 } 908 }
920 ip(7) = 0; 920 ip(7) = 0;
921 ip(8) = 0; 921 ip(8) = 0;
922 ip(9) = 0; 922 ip(9) = 0;
923 ip(10) = 0; 923 ip(10) = 0;
924 // ip(7) to ip(10) return values 924 // ip(7) to ip(10) return values
925 925
926 Array<octave_idx_type> iptr (dim_vector (14, 1)); 926 Array<octave_idx_type> iptr (dim_vector (14, 1));
927 octave_idx_type *ipntr = iptr.fortran_vec (); 927 octave_idx_type *ipntr = iptr.fortran_vec ();
928 928
929 octave_idx_type ido = 0; 929 octave_idx_type ido = 0;
930 int iter = 0; 930 int iter = 0;
933 OCTAVE_LOCAL_BUFFER (double, v, n * p); 933 OCTAVE_LOCAL_BUFFER (double, v, n * p);
934 OCTAVE_LOCAL_BUFFER (double, workl, lwork); 934 OCTAVE_LOCAL_BUFFER (double, workl, lwork);
935 OCTAVE_LOCAL_BUFFER (double, workd, 3 * n); 935 OCTAVE_LOCAL_BUFFER (double, workd, 3 * n);
936 double *presid = resid.fortran_vec (); 936 double *presid = resid.fortran_vec ();
937 937
938 do 938 do
939 { 939 {
940 F77_FUNC (dsaupd, DSAUPD) 940 F77_FUNC (dsaupd, DSAUPD)
941 (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, 941 (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
942 F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), 942 F77_CONST_CHAR_ARG2 ((typ.c_str()), 2),
943 k, tol, presid, p, v, n, iparam, 943 k, tol, presid, p, v, n, iparam,
944 ipntr, workd, workl, lwork, info 944 ipntr, workd, workl, lwork, info
945 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); 945 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
946 946
947 if (f77_exception_encountered) 947 if (f77_exception_encountered)
948 { 948 {
949 (*current_liboctave_error_handler) 949 (*current_liboctave_error_handler)
950 ("eigs: unrecoverable exception encountered in dsaupd"); 950 ("eigs: unrecoverable exception encountered in dsaupd");
951 return -1; 951 return -1;
952 } 952 }
953 953
954 if (disp > 0 && !xisnan(workl[iptr(5)-1])) 954 if (disp > 0 && !xisnan(workl[iptr(5)-1]))
955 { 955 {
956 if (iter++) 956 if (iter++)
957 { 957 {
958 os << "Iteration " << iter - 1 << 958 os << "Iteration " << iter - 1 <<
959 ": a few Ritz values of the " << p << "-by-" << 959 ": a few Ritz values of the " << p << "-by-" <<
960 p << " matrix\n"; 960 p << " matrix\n";
961 for (int i = 0 ; i < k; i++) 961 for (int i = 0 ; i < k; i++)
962 os << " " << workl[iptr(5)+i-1] << "\n"; 962 os << " " << workl[iptr(5)+i-1] << "\n";
963 } 963 }
966 // iteration pointer. But as workl[iptr(5)-1] is 966 // iteration pointer. But as workl[iptr(5)-1] is
967 // an output value updated at each iteration, setting 967 // an output value updated at each iteration, setting
968 // a value in this array to NaN and testing for it 968 // a value in this array to NaN and testing for it
969 // is a way of obtaining the iteration counter. 969 // is a way of obtaining the iteration counter.
970 if (ido != 99) 970 if (ido != 99)
971 workl[iptr(5)-1] = octave_NaN; 971 workl[iptr(5)-1] = octave_NaN;
972 } 972 }
973 973
974 if (ido == -1 || ido == 1 || ido == 2) 974 if (ido == -1 || ido == 1 || ido == 2)
975 { 975 {
976 if (have_b) 976 if (have_b)
977 { 977 {
978 Matrix mtmp (n,1); 978 Matrix mtmp (n,1);
979 for (octave_idx_type i = 0; i < n; i++) 979 for (octave_idx_type i = 0; i < n; i++)
980 mtmp(i,0) = workd[i + iptr(0) - 1]; 980 mtmp(i,0) = workd[i + iptr(0) - 1];
981 981
982 mtmp = utsolve(bt, permB, m * ltsolve(b, permB, mtmp)); 982 mtmp = utsolve(bt, permB, m * ltsolve(b, permB, mtmp));
983 983
984 for (octave_idx_type i = 0; i < n; i++) 984 for (octave_idx_type i = 0; i < n; i++)
985 workd[i+iptr(1)-1] = mtmp(i,0); 985 workd[i+iptr(1)-1] = mtmp(i,0);
986 } 986 }
987 else if (!vector_product (m, workd + iptr(0) - 1, 987 else if (!vector_product (m, workd + iptr(0) - 1,
988 workd + iptr(1) - 1)) 988 workd + iptr(1) - 1))
989 break; 989 break;
990 } 990 }
991 else 991 else
992 { 992 {
993 if (info < 0) 993 if (info < 0)
994 { 994 {
995 (*current_liboctave_error_handler) 995 (*current_liboctave_error_handler)
996 ("eigs: error %d in dsaupd", info); 996 ("eigs: error %d in dsaupd", info);
997 return -1; 997 return -1;
998 } 998 }
999 break; 999 break;
1000 } 1000 }
1001 } 1001 }
1002 while (1); 1002 while (1);
1003 1003
1004 octave_idx_type info2; 1004 octave_idx_type info2;
1005 1005
1006 // We have a problem in that the size of the C++ bool 1006 // We have a problem in that the size of the C++ bool
1007 // type relative to the fortran logical type. It appears 1007 // type relative to the fortran logical type. It appears
1008 // that fortran uses 4- or 8-bytes per logical and C++ 1-byte 1008 // that fortran uses 4- or 8-bytes per logical and C++ 1-byte
1009 // per bool, though this might be system dependent. As 1009 // per bool, though this might be system dependent. As
1010 // long as the HOWMNY arg is not "S", the logical array 1010 // long as the HOWMNY arg is not "S", the logical array
1011 // is just workspace for ARPACK, so use int type to 1011 // is just workspace for ARPACK, so use int type to
1012 // avoid problems. 1012 // avoid problems.
1013 Array<octave_idx_type> s (dim_vector (p, 1)); 1013 Array<octave_idx_type> s (dim_vector (p, 1));
1014 octave_idx_type *sel = s.fortran_vec (); 1014 octave_idx_type *sel = s.fortran_vec ();
1015 1015
1016 eig_vec.resize (n, k); 1016 eig_vec.resize (n, k);
1017 double *z = eig_vec.fortran_vec (); 1017 double *z = eig_vec.fortran_vec ();
1018 1018
1019 eig_val.resize (k); 1019 eig_val.resize (k);
1020 double *d = eig_val.fortran_vec (); 1020 double *d = eig_val.fortran_vec ();
1021 1021
1022 F77_FUNC (dseupd, DSEUPD) 1022 F77_FUNC (dseupd, DSEUPD)
1023 (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma, 1023 (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma,
1024 F77_CONST_CHAR_ARG2 (&bmat, 1), n, 1024 F77_CONST_CHAR_ARG2 (&bmat, 1), n,
1025 F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), k, tol, presid, p, v, n, iparam, 1025 F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), k, tol, presid, p, v, n, iparam,
1026 ipntr, workd, workl, lwork, info2 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) 1026 ipntr, workd, workl, lwork, info2 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1)
1027 F77_CHAR_ARG_LEN(2)); 1027 F77_CHAR_ARG_LEN(2));
1028 1028
1029 if (f77_exception_encountered) 1029 if (f77_exception_encountered)
1030 { 1030 {
1031 (*current_liboctave_error_handler) 1031 (*current_liboctave_error_handler)
1076 eig_vec = ltsolve(b, permB, eig_vec); 1076 eig_vec = ltsolve(b, permB, eig_vec);
1077 } 1077 }
1078 } 1078 }
1079 else 1079 else
1080 { 1080 {
1081 (*current_liboctave_error_handler) 1081 (*current_liboctave_error_handler)
1082 ("eigs: error %d in dseupd", info2); 1082 ("eigs: error %d in dseupd", info2);
1083 return -1; 1083 return -1;
1084 } 1084 }
1085 } 1085 }
1086 1086
1088 } 1088 }
1089 1089
1090 template <class M> 1090 template <class M>
1091 octave_idx_type 1091 octave_idx_type
1092 EigsRealSymmetricMatrixShift (const M& m, double sigma, 1092 EigsRealSymmetricMatrixShift (const M& m, double sigma,
1093 octave_idx_type k, octave_idx_type p, 1093 octave_idx_type k, octave_idx_type p,
1094 octave_idx_type &info, Matrix &eig_vec, 1094 octave_idx_type &info, Matrix &eig_vec,
1095 ColumnVector &eig_val, const M& _b, 1095 ColumnVector &eig_val, const M& _b,
1096 ColumnVector &permB, ColumnVector &resid, 1096 ColumnVector &permB, ColumnVector &resid,
1097 std::ostream& os, double tol, bool rvec, 1097 std::ostream& os, double tol, bool rvec,
1098 bool cholB, int disp, int maxit) 1098 bool cholB, int disp, int maxit)
1099 { 1099 {
1100 M b(_b); 1100 M b(_b);
1101 octave_idx_type n = m.cols (); 1101 octave_idx_type n = m.cols ();
1102 octave_idx_type mode = 3; 1102 octave_idx_type mode = 3;
1108 (*current_liboctave_error_handler) ("eigs: A must be square"); 1108 (*current_liboctave_error_handler) ("eigs: A must be square");
1109 return -1; 1109 return -1;
1110 } 1110 }
1111 if (have_b && (m.rows() != b.rows() || m.rows() != b.cols())) 1111 if (have_b && (m.rows() != b.rows() || m.rows() != b.cols()))
1112 { 1112 {
1113 (*current_liboctave_error_handler) 1113 (*current_liboctave_error_handler)
1114 ("eigs: B must be square and the same size as A"); 1114 ("eigs: B must be square and the same size as A");
1115 return -1; 1115 return -1;
1116 } 1116 }
1117 1117
1118 // FIXME: The "SM" type for mode 1 seems unstable though faster!! 1118 // FIXME: The "SM" type for mode 1 seems unstable though faster!!
1136 return -1; 1136 return -1;
1137 } 1137 }
1138 1138
1139 if (k <= 0 || k >= n - 1) 1139 if (k <= 0 || k >= n - 1)
1140 { 1140 {
1141 (*current_liboctave_error_handler) 1141 (*current_liboctave_error_handler)
1142 ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1-1).\n" 1142 ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1-1).\n"
1143 " Use 'eig(full(A))' instead"); 1143 " Use 'eig(full(A))' instead");
1144 return -1; 1144 return -1;
1145 } 1145 }
1146 1146
1148 { 1148 {
1149 p = k * 2; 1149 p = k * 2;
1150 1150
1151 if (p < 20) 1151 if (p < 20)
1152 p = 20; 1152 p = 20;
1153 1153
1154 if (p > n - 1) 1154 if (p > n - 1)
1155 p = n - 1 ; 1155 p = n - 1 ;
1156 } 1156 }
1157 1157
1158 if (p <= k || p >= n) 1158 if (p <= k || p >= n)
1159 { 1159 {
1160 (*current_liboctave_error_handler) 1160 (*current_liboctave_error_handler)
1161 ("eigs: opts.p must be greater than k and less than n"); 1161 ("eigs: opts.p must be greater than k and less than n");
1162 return -1; 1162 return -1;
1163 } 1163 }
1164 1164
1165 if (have_b && cholB && permB.length() != 0) 1165 if (have_b && cholB && permB.length() != 0)
1166 { 1166 {
1167 // Check the we really have a permutation vector 1167 // Check the we really have a permutation vector
1168 if (permB.length() != n) 1168 if (permB.length() != n)
1169 { 1169 {
1170 (*current_liboctave_error_handler) ("eigs: permB vector invalid"); 1170 (*current_liboctave_error_handler) ("eigs: permB vector invalid");
1173 else 1173 else
1174 { 1174 {
1175 Array<bool> checked (dim_vector (n, 1), false); 1175 Array<bool> checked (dim_vector (n, 1), false);
1176 for (octave_idx_type i = 0; i < n; i++) 1176 for (octave_idx_type i = 0; i < n; i++)
1177 { 1177 {
1178 octave_idx_type bidx = 1178 octave_idx_type bidx =
1179 static_cast<octave_idx_type> (permB(i)); 1179 static_cast<octave_idx_type> (permB(i));
1180 if (checked(bidx) || bidx < 0 || 1180 if (checked(bidx) || bidx < 0 ||
1181 bidx >= n || D_NINT (bidx) != bidx) 1181 bidx >= n || D_NINT (bidx) != bidx)
1182 { 1182 {
1183 (*current_liboctave_error_handler) 1183 (*current_liboctave_error_handler)
1184 ("eigs: permB vector invalid"); 1184 ("eigs: permB vector invalid");
1185 return -1; 1185 return -1;
1186 } 1186 }
1187 } 1187 }
1188 } 1188 }
1226 OCTAVE_LOCAL_BUFFER (double, v, n * p); 1226 OCTAVE_LOCAL_BUFFER (double, v, n * p);
1227 OCTAVE_LOCAL_BUFFER (double, workl, lwork); 1227 OCTAVE_LOCAL_BUFFER (double, workl, lwork);
1228 OCTAVE_LOCAL_BUFFER (double, workd, 3 * n); 1228 OCTAVE_LOCAL_BUFFER (double, workd, 3 * n);
1229 double *presid = resid.fortran_vec (); 1229 double *presid = resid.fortran_vec ();
1230 1230
1231 do 1231 do
1232 { 1232 {
1233 F77_FUNC (dsaupd, DSAUPD) 1233 F77_FUNC (dsaupd, DSAUPD)
1234 (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, 1234 (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
1235 F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), 1235 F77_CONST_CHAR_ARG2 ((typ.c_str()), 2),
1236 k, tol, presid, p, v, n, iparam, 1236 k, tol, presid, p, v, n, iparam,
1237 ipntr, workd, workl, lwork, info 1237 ipntr, workd, workl, lwork, info
1238 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); 1238 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
1239 1239
1240 if (f77_exception_encountered) 1240 if (f77_exception_encountered)
1241 { 1241 {
1242 (*current_liboctave_error_handler) 1242 (*current_liboctave_error_handler)
1243 ("eigs: unrecoverable exception encountered in dsaupd"); 1243 ("eigs: unrecoverable exception encountered in dsaupd");
1244 return -1; 1244 return -1;
1245 } 1245 }
1246 1246
1247 if (disp > 0 && !xisnan(workl[iptr(5)-1])) 1247 if (disp > 0 && !xisnan(workl[iptr(5)-1]))
1248 { 1248 {
1249 if (iter++) 1249 if (iter++)
1250 { 1250 {
1251 os << "Iteration " << iter - 1 << 1251 os << "Iteration " << iter - 1 <<
1252 ": a few Ritz values of the " << p << "-by-" << 1252 ": a few Ritz values of the " << p << "-by-" <<
1253 p << " matrix\n"; 1253 p << " matrix\n";
1254 for (int i = 0 ; i < k; i++) 1254 for (int i = 0 ; i < k; i++)
1255 os << " " << workl[iptr(5)+i-1] << "\n"; 1255 os << " " << workl[iptr(5)+i-1] << "\n";
1256 } 1256 }
1259 // iteration pointer. But as workl[iptr(5)-1] is 1259 // iteration pointer. But as workl[iptr(5)-1] is
1260 // an output value updated at each iteration, setting 1260 // an output value updated at each iteration, setting
1261 // a value in this array to NaN and testing for it 1261 // a value in this array to NaN and testing for it
1262 // is a way of obtaining the iteration counter. 1262 // is a way of obtaining the iteration counter.
1263 if (ido != 99) 1263 if (ido != 99)
1264 workl[iptr(5)-1] = octave_NaN; 1264 workl[iptr(5)-1] = octave_NaN;
1265 } 1265 }
1266 1266
1267 if (ido == -1 || ido == 1 || ido == 2) 1267 if (ido == -1 || ido == 1 || ido == 2)
1268 { 1268 {
1269 if (have_b) 1269 if (have_b)
1276 1276
1277 Matrix tmp(n, 1); 1277 Matrix tmp(n, 1);
1278 1278
1279 for (octave_idx_type i = 0; i < n; i++) 1279 for (octave_idx_type i = 0; i < n; i++)
1280 tmp(i,0) = dtmp[P[i]]; 1280 tmp(i,0) = dtmp[P[i]];
1281 1281
1282 lusolve (L, U, tmp); 1282 lusolve (L, U, tmp);
1283 1283
1284 double *ip2 = workd+iptr(1)-1; 1284 double *ip2 = workd+iptr(1)-1;
1285 for (octave_idx_type i = 0; i < n; i++) 1285 for (octave_idx_type i = 0; i < n; i++)
1286 ip2[Q[i]] = tmp(i,0); 1286 ip2[Q[i]] = tmp(i,0);
1292 double *ip2 = workd+iptr(2)-1; 1292 double *ip2 = workd+iptr(2)-1;
1293 Matrix tmp(n, 1); 1293 Matrix tmp(n, 1);
1294 1294
1295 for (octave_idx_type i = 0; i < n; i++) 1295 for (octave_idx_type i = 0; i < n; i++)
1296 tmp(i,0) = ip2[P[i]]; 1296 tmp(i,0) = ip2[P[i]];
1297 1297
1298 lusolve (L, U, tmp); 1298 lusolve (L, U, tmp);
1299 1299
1300 ip2 = workd+iptr(1)-1; 1300 ip2 = workd+iptr(1)-1;
1301 for (octave_idx_type i = 0; i < n; i++) 1301 for (octave_idx_type i = 0; i < n; i++)
1302 ip2[Q[i]] = tmp(i,0); 1302 ip2[Q[i]] = tmp(i,0);
1314 double *ip2 = workd+iptr(0)-1; 1314 double *ip2 = workd+iptr(0)-1;
1315 Matrix tmp(n, 1); 1315 Matrix tmp(n, 1);
1316 1316
1317 for (octave_idx_type i = 0; i < n; i++) 1317 for (octave_idx_type i = 0; i < n; i++)
1318 tmp(i,0) = ip2[P[i]]; 1318 tmp(i,0) = ip2[P[i]];
1319 1319
1320 lusolve (L, U, tmp); 1320 lusolve (L, U, tmp);
1321 1321
1322 ip2 = workd+iptr(1)-1; 1322 ip2 = workd+iptr(1)-1;
1323 for (octave_idx_type i = 0; i < n; i++) 1323 for (octave_idx_type i = 0; i < n; i++)
1324 ip2[Q[i]] = tmp(i,0); 1324 ip2[Q[i]] = tmp(i,0);
1327 } 1327 }
1328 else 1328 else
1329 { 1329 {
1330 if (info < 0) 1330 if (info < 0)
1331 { 1331 {
1332 (*current_liboctave_error_handler) 1332 (*current_liboctave_error_handler)
1333 ("eigs: error %d in dsaupd", info); 1333 ("eigs: error %d in dsaupd", info);
1334 return -1; 1334 return -1;
1335 } 1335 }
1336 break; 1336 break;
1337 } 1337 }
1338 } 1338 }
1339 while (1); 1339 while (1);
1340 1340
1341 octave_idx_type info2; 1341 octave_idx_type info2;
1342 1342
1343 // We have a problem in that the size of the C++ bool 1343 // We have a problem in that the size of the C++ bool
1344 // type relative to the fortran logical type. It appears 1344 // type relative to the fortran logical type. It appears
1345 // that fortran uses 4- or 8-bytes per logical and C++ 1-byte 1345 // that fortran uses 4- or 8-bytes per logical and C++ 1-byte
1346 // per bool, though this might be system dependent. As 1346 // per bool, though this might be system dependent. As
1347 // long as the HOWMNY arg is not "S", the logical array 1347 // long as the HOWMNY arg is not "S", the logical array
1348 // is just workspace for ARPACK, so use int type to 1348 // is just workspace for ARPACK, so use int type to
1349 // avoid problems. 1349 // avoid problems.
1350 Array<octave_idx_type> s (dim_vector (p, 1)); 1350 Array<octave_idx_type> s (dim_vector (p, 1));
1351 octave_idx_type *sel = s.fortran_vec (); 1351 octave_idx_type *sel = s.fortran_vec ();
1352 1352
1353 eig_vec.resize (n, k); 1353 eig_vec.resize (n, k);
1354 double *z = eig_vec.fortran_vec (); 1354 double *z = eig_vec.fortran_vec ();
1355 1355
1356 eig_val.resize (k); 1356 eig_val.resize (k);
1357 double *d = eig_val.fortran_vec (); 1357 double *d = eig_val.fortran_vec ();
1358 1358
1359 F77_FUNC (dseupd, DSEUPD) 1359 F77_FUNC (dseupd, DSEUPD)
1360 (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma, 1360 (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma,
1361 F77_CONST_CHAR_ARG2 (&bmat, 1), n, 1361 F77_CONST_CHAR_ARG2 (&bmat, 1), n,
1362 F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), 1362 F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2),
1363 k, tol, presid, p, v, n, iparam, ipntr, workd, workl, lwork, info2 1363 k, tol, presid, p, v, n, iparam, ipntr, workd, workl, lwork, info2
1364 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); 1364 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
1365 1365
1416 } 1416 }
1417 1417
1418 octave_idx_type 1418 octave_idx_type
1419 EigsRealSymmetricFunc (EigsFunc fun, octave_idx_type n, 1419 EigsRealSymmetricFunc (EigsFunc fun, octave_idx_type n,
1420 const std::string &_typ, double sigma, 1420 const std::string &_typ, double sigma,
1421 octave_idx_type k, octave_idx_type p, 1421 octave_idx_type k, octave_idx_type p,
1422 octave_idx_type &info, Matrix &eig_vec, 1422 octave_idx_type &info, Matrix &eig_vec,
1423 ColumnVector &eig_val, ColumnVector &resid, 1423 ColumnVector &eig_val, ColumnVector &resid,
1424 std::ostream& os, double tol, bool rvec, 1424 std::ostream& os, double tol, bool rvec,
1425 bool /* cholB */, int disp, int maxit) 1425 bool /* cholB */, int disp, int maxit)
1426 { 1426 {
1427 std::string typ (_typ); 1427 std::string typ (_typ);
1428 bool have_sigma = (sigma ? true : false); 1428 bool have_sigma = (sigma ? true : false);
1449 { 1449 {
1450 p = k * 2; 1450 p = k * 2;
1451 1451
1452 if (p < 20) 1452 if (p < 20)
1453 p = 20; 1453 p = 20;
1454 1454
1455 if (p > n - 1) 1455 if (p > n - 1)
1456 p = n - 1 ; 1456 p = n - 1 ;
1457 } 1457 }
1458 1458
1459 if (k <= 0 || k >= n - 1) 1459 if (k <= 0 || k >= n - 1)
1460 { 1460 {
1461 (*current_liboctave_error_handler) 1461 (*current_liboctave_error_handler)
1462 ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1).\n" 1462 ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1).\n"
1463 " Use 'eig(full(A))' instead"); 1463 " Use 'eig(full(A))' instead");
1471 return -1; 1471 return -1;
1472 } 1472 }
1473 1473
1474 if (! have_sigma) 1474 if (! have_sigma)
1475 { 1475 {
1476 if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" && 1476 if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" &&
1477 typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" && 1477 typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" &&
1478 typ != "SI") 1478 typ != "SI")
1479 (*current_liboctave_error_handler) 1479 (*current_liboctave_error_handler)
1480 ("eigs: unrecognized sigma value"); 1480 ("eigs: unrecognized sigma value");
1481 1481
1482 if (typ == "LI" || typ == "SI" || typ == "LR" || typ == "SR") 1482 if (typ == "LI" || typ == "SI" || typ == "LR" || typ == "SR")
1483 { 1483 {
1484 (*current_liboctave_error_handler) 1484 (*current_liboctave_error_handler)
1485 ("eigs: invalid sigma value for real symmetric problem"); 1485 ("eigs: invalid sigma value for real symmetric problem");
1486 return -1; 1486 return -1;
1487 } 1487 }
1488 1488
1489 if (typ == "SM") 1489 if (typ == "SM")
1514 ip(7) = 0; 1514 ip(7) = 0;
1515 ip(8) = 0; 1515 ip(8) = 0;
1516 ip(9) = 0; 1516 ip(9) = 0;
1517 ip(10) = 0; 1517 ip(10) = 0;
1518 // ip(7) to ip(10) return values 1518 // ip(7) to ip(10) return values
1519 1519
1520 Array<octave_idx_type> iptr (dim_vector (14, 1)); 1520 Array<octave_idx_type> iptr (dim_vector (14, 1));
1521 octave_idx_type *ipntr = iptr.fortran_vec (); 1521 octave_idx_type *ipntr = iptr.fortran_vec ();
1522 1522
1523 octave_idx_type ido = 0; 1523 octave_idx_type ido = 0;
1524 int iter = 0; 1524 int iter = 0;
1527 OCTAVE_LOCAL_BUFFER (double, v, n * p); 1527 OCTAVE_LOCAL_BUFFER (double, v, n * p);
1528 OCTAVE_LOCAL_BUFFER (double, workl, lwork); 1528 OCTAVE_LOCAL_BUFFER (double, workl, lwork);
1529 OCTAVE_LOCAL_BUFFER (double, workd, 3 * n); 1529 OCTAVE_LOCAL_BUFFER (double, workd, 3 * n);
1530 double *presid = resid.fortran_vec (); 1530 double *presid = resid.fortran_vec ();
1531 1531
1532 do 1532 do
1533 { 1533 {
1534 F77_FUNC (dsaupd, DSAUPD) 1534 F77_FUNC (dsaupd, DSAUPD)
1535 (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, 1535 (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
1536 F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), 1536 F77_CONST_CHAR_ARG2 ((typ.c_str()), 2),
1537 k, tol, presid, p, v, n, iparam, 1537 k, tol, presid, p, v, n, iparam,
1538 ipntr, workd, workl, lwork, info 1538 ipntr, workd, workl, lwork, info
1539 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); 1539 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
1540 1540
1541 if (f77_exception_encountered) 1541 if (f77_exception_encountered)
1542 { 1542 {
1543 (*current_liboctave_error_handler) 1543 (*current_liboctave_error_handler)
1544 ("eigs: unrecoverable exception encountered in dsaupd"); 1544 ("eigs: unrecoverable exception encountered in dsaupd");
1545 return -1; 1545 return -1;
1546 } 1546 }
1547 1547
1548 if (disp > 0 && !xisnan(workl[iptr(5)-1])) 1548 if (disp > 0 && !xisnan(workl[iptr(5)-1]))
1549 { 1549 {
1550 if (iter++) 1550 if (iter++)
1551 { 1551 {
1552 os << "Iteration " << iter - 1 << 1552 os << "Iteration " << iter - 1 <<
1553 ": a few Ritz values of the " << p << "-by-" << 1553 ": a few Ritz values of the " << p << "-by-" <<
1554 p << " matrix\n"; 1554 p << " matrix\n";
1555 for (int i = 0 ; i < k; i++) 1555 for (int i = 0 ; i < k; i++)
1556 os << " " << workl[iptr(5)+i-1] << "\n"; 1556 os << " " << workl[iptr(5)+i-1] << "\n";
1557 } 1557 }
1560 // iteration pointer. But as workl[iptr(5)-1] is 1560 // iteration pointer. But as workl[iptr(5)-1] is
1561 // an output value updated at each iteration, setting 1561 // an output value updated at each iteration, setting
1562 // a value in this array to NaN and testing for it 1562 // a value in this array to NaN and testing for it
1563 // is a way of obtaining the iteration counter. 1563 // is a way of obtaining the iteration counter.
1564 if (ido != 99) 1564 if (ido != 99)
1565 workl[iptr(5)-1] = octave_NaN; 1565 workl[iptr(5)-1] = octave_NaN;
1566 } 1566 }
1567 1567
1568 1568
1569 if (ido == -1 || ido == 1 || ido == 2) 1569 if (ido == -1 || ido == 1 || ido == 2)
1570 { 1570 {
1585 } 1585 }
1586 else 1586 else
1587 { 1587 {
1588 if (info < 0) 1588 if (info < 0)
1589 { 1589 {
1590 (*current_liboctave_error_handler) 1590 (*current_liboctave_error_handler)
1591 ("eigs: error %d in dsaupd", info); 1591 ("eigs: error %d in dsaupd", info);
1592 return -1; 1592 return -1;
1593 } 1593 }
1594 break; 1594 break;
1595 } 1595 }
1596 } 1596 }
1597 while (1); 1597 while (1);
1598 1598
1599 octave_idx_type info2; 1599 octave_idx_type info2;
1600 1600
1601 // We have a problem in that the size of the C++ bool 1601 // We have a problem in that the size of the C++ bool
1602 // type relative to the fortran logical type. It appears 1602 // type relative to the fortran logical type. It appears
1603 // that fortran uses 4- or 8-bytes per logical and C++ 1-byte 1603 // that fortran uses 4- or 8-bytes per logical and C++ 1-byte
1604 // per bool, though this might be system dependent. As 1604 // per bool, though this might be system dependent. As
1605 // long as the HOWMNY arg is not "S", the logical array 1605 // long as the HOWMNY arg is not "S", the logical array
1606 // is just workspace for ARPACK, so use int type to 1606 // is just workspace for ARPACK, so use int type to
1607 // avoid problems. 1607 // avoid problems.
1608 Array<octave_idx_type> s (dim_vector (p, 1)); 1608 Array<octave_idx_type> s (dim_vector (p, 1));
1609 octave_idx_type *sel = s.fortran_vec (); 1609 octave_idx_type *sel = s.fortran_vec ();
1610 1610
1611 eig_vec.resize (n, k); 1611 eig_vec.resize (n, k);
1612 double *z = eig_vec.fortran_vec (); 1612 double *z = eig_vec.fortran_vec ();
1613 1613
1614 eig_val.resize (k); 1614 eig_val.resize (k);
1615 double *d = eig_val.fortran_vec (); 1615 double *d = eig_val.fortran_vec ();
1616 1616
1617 F77_FUNC (dseupd, DSEUPD) 1617 F77_FUNC (dseupd, DSEUPD)
1618 (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma, 1618 (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma,
1619 F77_CONST_CHAR_ARG2 (&bmat, 1), n, 1619 F77_CONST_CHAR_ARG2 (&bmat, 1), n,
1620 F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), 1620 F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2),
1621 k, tol, presid, p, v, n, iparam, ipntr, workd, workl, lwork, info2 1621 k, tol, presid, p, v, n, iparam, ipntr, workd, workl, lwork, info2
1622 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); 1622 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
1623 1623
1679 return ip(4); 1679 return ip(4);
1680 } 1680 }
1681 1681
1682 template <class M> 1682 template <class M>
1683 octave_idx_type 1683 octave_idx_type
1684 EigsRealNonSymmetricMatrix (const M& m, const std::string typ, 1684 EigsRealNonSymmetricMatrix (const M& m, const std::string typ,
1685 octave_idx_type k, octave_idx_type p, 1685 octave_idx_type k, octave_idx_type p,
1686 octave_idx_type &info, ComplexMatrix &eig_vec, 1686 octave_idx_type &info, ComplexMatrix &eig_vec,
1687 ComplexColumnVector &eig_val, const M& _b, 1687 ComplexColumnVector &eig_val, const M& _b,
1688 ColumnVector &permB, ColumnVector &resid, 1688 ColumnVector &permB, ColumnVector &resid,
1689 std::ostream& os, double tol, bool rvec, 1689 std::ostream& os, double tol, bool rvec,
1690 bool cholB, int disp, int maxit) 1690 bool cholB, int disp, int maxit)
1691 { 1691 {
1692 M b(_b); 1692 M b(_b);
1693 octave_idx_type n = m.cols (); 1693 octave_idx_type n = m.cols ();
1694 octave_idx_type mode = 1; 1694 octave_idx_type mode = 1;
1704 (*current_liboctave_error_handler) ("eigs: A must be square"); 1704 (*current_liboctave_error_handler) ("eigs: A must be square");
1705 return -1; 1705 return -1;
1706 } 1706 }
1707 if (have_b && (m.rows() != b.rows() || m.rows() != b.cols())) 1707 if (have_b && (m.rows() != b.rows() || m.rows() != b.cols()))
1708 { 1708 {
1709 (*current_liboctave_error_handler) 1709 (*current_liboctave_error_handler)
1710 ("eigs: B must be square and the same size as A"); 1710 ("eigs: B must be square and the same size as A");
1711 return -1; 1711 return -1;
1712 } 1712 }
1713 1713
1714 if (resid.is_empty()) 1714 if (resid.is_empty())
1730 { 1730 {
1731 p = k * 2 + 1; 1731 p = k * 2 + 1;
1732 1732
1733 if (p < 20) 1733 if (p < 20)
1734 p = 20; 1734 p = 20;
1735 1735
1736 if (p > n - 1) 1736 if (p > n - 1)
1737 p = n - 1 ; 1737 p = n - 1 ;
1738 } 1738 }
1739 1739
1740 if (k <= 0 || k >= n - 1) 1740 if (k <= 0 || k >= n - 1)
1741 { 1741 {
1742 (*current_liboctave_error_handler) 1742 (*current_liboctave_error_handler)
1743 ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1).\n" 1743 ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1).\n"
1744 " Use 'eig(full(A))' instead"); 1744 " Use 'eig(full(A))' instead");
1745 return -1; 1745 return -1;
1746 } 1746 }
1747 1747
1748 if (p <= k || p >= n) 1748 if (p <= k || p >= n)
1749 { 1749 {
1750 (*current_liboctave_error_handler) 1750 (*current_liboctave_error_handler)
1751 ("eigs: opts.p must be greater than k and less than n"); 1751 ("eigs: opts.p must be greater than k and less than n");
1752 return -1; 1752 return -1;
1753 } 1753 }
1754 1754
1755 if (have_b && cholB && permB.length() != 0) 1755 if (have_b && cholB && permB.length() != 0)
1756 { 1756 {
1757 // Check the we really have a permutation vector 1757 // Check the we really have a permutation vector
1758 if (permB.length() != n) 1758 if (permB.length() != n)
1759 { 1759 {
1760 (*current_liboctave_error_handler) 1760 (*current_liboctave_error_handler)
1761 ("eigs: permB vector invalid"); 1761 ("eigs: permB vector invalid");
1762 return -1; 1762 return -1;
1763 } 1763 }
1764 else 1764 else
1765 { 1765 {
1766 Array<bool> checked (dim_vector (n, 1), false); 1766 Array<bool> checked (dim_vector (n, 1), false);
1767 for (octave_idx_type i = 0; i < n; i++) 1767 for (octave_idx_type i = 0; i < n; i++)
1768 { 1768 {
1769 octave_idx_type bidx = 1769 octave_idx_type bidx =
1770 static_cast<octave_idx_type> (permB(i)); 1770 static_cast<octave_idx_type> (permB(i));
1771 if (checked(bidx) || bidx < 0 || 1771 if (checked(bidx) || bidx < 0 ||
1772 bidx >= n || D_NINT (bidx) != bidx) 1772 bidx >= n || D_NINT (bidx) != bidx)
1773 { 1773 {
1774 (*current_liboctave_error_handler) 1774 (*current_liboctave_error_handler)
1775 ("eigs: permB vector invalid"); 1775 ("eigs: permB vector invalid");
1776 return -1; 1776 return -1;
1777 } 1777 }
1778 } 1778 }
1779 } 1779 }
1780 } 1780 }
1781 1781
1782 if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" && 1782 if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" &&
1783 typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" && 1783 typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" &&
1784 typ != "SI") 1784 typ != "SI")
1785 { 1785 {
1786 (*current_liboctave_error_handler) 1786 (*current_liboctave_error_handler)
1787 ("eigs: unrecognized sigma value"); 1787 ("eigs: unrecognized sigma value");
1788 return -1; 1788 return -1;
1789 } 1789 }
1790 1790
1791 if (typ == "LA" || typ == "SA" || typ == "BE") 1791 if (typ == "LA" || typ == "SA" || typ == "BE")
1792 { 1792 {
1793 (*current_liboctave_error_handler) 1793 (*current_liboctave_error_handler)
1794 ("eigs: invalid sigma value for unsymmetric problem"); 1794 ("eigs: invalid sigma value for unsymmetric problem");
1795 return -1; 1795 return -1;
1796 } 1796 }
1797 1797
1798 if (have_b) 1798 if (have_b)
1812 } 1812 }
1813 else 1813 else
1814 { 1814 {
1815 if (! make_cholb(b, bt, permB)) 1815 if (! make_cholb(b, bt, permB))
1816 { 1816 {
1817 (*current_liboctave_error_handler) 1817 (*current_liboctave_error_handler)
1818 ("eigs: The matrix B is not positive definite"); 1818 ("eigs: The matrix B is not positive definite");
1819 return -1; 1819 return -1;
1820 } 1820 }
1821 } 1821 }
1822 } 1822 }
1834 ip(7) = 0; 1834 ip(7) = 0;
1835 ip(8) = 0; 1835 ip(8) = 0;
1836 ip(9) = 0; 1836 ip(9) = 0;
1837 ip(10) = 0; 1837 ip(10) = 0;
1838 // ip(7) to ip(10) return values 1838 // ip(7) to ip(10) return values
1839 1839
1840 Array<octave_idx_type> iptr (dim_vector (14, 1)); 1840 Array<octave_idx_type> iptr (dim_vector (14, 1));
1841 octave_idx_type *ipntr = iptr.fortran_vec (); 1841 octave_idx_type *ipntr = iptr.fortran_vec ();
1842 1842
1843 octave_idx_type ido = 0; 1843 octave_idx_type ido = 0;
1844 int iter = 0; 1844 int iter = 0;
1847 OCTAVE_LOCAL_BUFFER (double, v, n * (p + 1)); 1847 OCTAVE_LOCAL_BUFFER (double, v, n * (p + 1));
1848 OCTAVE_LOCAL_BUFFER (double, workl, lwork + 1); 1848 OCTAVE_LOCAL_BUFFER (double, workl, lwork + 1);
1849 OCTAVE_LOCAL_BUFFER (double, workd, 3 * n + 1); 1849 OCTAVE_LOCAL_BUFFER (double, workd, 3 * n + 1);
1850 double *presid = resid.fortran_vec (); 1850 double *presid = resid.fortran_vec ();
1851 1851
1852 do 1852 do
1853 { 1853 {
1854 F77_FUNC (dnaupd, DNAUPD) 1854 F77_FUNC (dnaupd, DNAUPD)
1855 (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, 1855 (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
1856 F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), 1856 F77_CONST_CHAR_ARG2 ((typ.c_str()), 2),
1857 k, tol, presid, p, v, n, iparam, 1857 k, tol, presid, p, v, n, iparam,
1858 ipntr, workd, workl, lwork, info 1858 ipntr, workd, workl, lwork, info
1859 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); 1859 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
1860 1860
1861 if (f77_exception_encountered) 1861 if (f77_exception_encountered)
1862 { 1862 {
1863 (*current_liboctave_error_handler) 1863 (*current_liboctave_error_handler)
1864 ("eigs: unrecoverable exception encountered in dnaupd"); 1864 ("eigs: unrecoverable exception encountered in dnaupd");
1865 return -1; 1865 return -1;
1866 } 1866 }
1867 1867
1868 if (disp > 0 && !xisnan(workl[iptr(5)-1])) 1868 if (disp > 0 && !xisnan(workl[iptr(5)-1]))
1869 { 1869 {
1870 if (iter++) 1870 if (iter++)
1871 { 1871 {
1872 os << "Iteration " << iter - 1 << 1872 os << "Iteration " << iter - 1 <<
1873 ": a few Ritz values of the " << p << "-by-" << 1873 ": a few Ritz values of the " << p << "-by-" <<
1874 p << " matrix\n"; 1874 p << " matrix\n";
1875 for (int i = 0 ; i < k; i++) 1875 for (int i = 0 ; i < k; i++)
1876 os << " " << workl[iptr(5)+i-1] << "\n"; 1876 os << " " << workl[iptr(5)+i-1] << "\n";
1877 } 1877 }
1880 // iteration pointer. But as workl[iptr(5)-1] is 1880 // iteration pointer. But as workl[iptr(5)-1] is
1881 // an output value updated at each iteration, setting 1881 // an output value updated at each iteration, setting
1882 // a value in this array to NaN and testing for it 1882 // a value in this array to NaN and testing for it
1883 // is a way of obtaining the iteration counter. 1883 // is a way of obtaining the iteration counter.
1884 if (ido != 99) 1884 if (ido != 99)
1885 workl[iptr(5)-1] = octave_NaN; 1885 workl[iptr(5)-1] = octave_NaN;
1886 } 1886 }
1887 1887
1888 if (ido == -1 || ido == 1 || ido == 2) 1888 if (ido == -1 || ido == 1 || ido == 2)
1889 { 1889 {
1890 if (have_b) 1890 if (have_b)
1891 { 1891 {
1892 Matrix mtmp (n,1); 1892 Matrix mtmp (n,1);
1893 for (octave_idx_type i = 0; i < n; i++) 1893 for (octave_idx_type i = 0; i < n; i++)
1894 mtmp(i,0) = workd[i + iptr(0) - 1]; 1894 mtmp(i,0) = workd[i + iptr(0) - 1];
1895 1895
1896 mtmp = utsolve(bt, permB, m * ltsolve(b, permB, mtmp)); 1896 mtmp = utsolve(bt, permB, m * ltsolve(b, permB, mtmp));
1897 1897
1898 for (octave_idx_type i = 0; i < n; i++) 1898 for (octave_idx_type i = 0; i < n; i++)
1899 workd[i+iptr(1)-1] = mtmp(i,0); 1899 workd[i+iptr(1)-1] = mtmp(i,0);
1900 } 1900 }
1901 else if (!vector_product (m, workd + iptr(0) - 1, 1901 else if (!vector_product (m, workd + iptr(0) - 1,
1902 workd + iptr(1) - 1)) 1902 workd + iptr(1) - 1))
1903 break; 1903 break;
1904 } 1904 }
1905 else 1905 else
1906 { 1906 {
1907 if (info < 0) 1907 if (info < 0)
1908 { 1908 {
1909 (*current_liboctave_error_handler) 1909 (*current_liboctave_error_handler)
1910 ("eigs: error %d in dnaupd", info); 1910 ("eigs: error %d in dnaupd", info);
1911 return -1; 1911 return -1;
1912 } 1912 }
1913 break; 1913 break;
1914 } 1914 }
1915 } 1915 }
1916 while (1); 1916 while (1);
1917 1917
1918 octave_idx_type info2; 1918 octave_idx_type info2;
1919 1919
1920 // We have a problem in that the size of the C++ bool 1920 // We have a problem in that the size of the C++ bool
1921 // type relative to the fortran logical type. It appears 1921 // type relative to the fortran logical type. It appears
1922 // that fortran uses 4- or 8-bytes per logical and C++ 1-byte 1922 // that fortran uses 4- or 8-bytes per logical and C++ 1-byte
1923 // per bool, though this might be system dependent. As 1923 // per bool, though this might be system dependent. As
1924 // long as the HOWMNY arg is not "S", the logical array 1924 // long as the HOWMNY arg is not "S", the logical array
1925 // is just workspace for ARPACK, so use int type to 1925 // is just workspace for ARPACK, so use int type to
1926 // avoid problems. 1926 // avoid problems.
1927 Array<octave_idx_type> s (dim_vector (p, 1)); 1927 Array<octave_idx_type> s (dim_vector (p, 1));
1928 octave_idx_type *sel = s.fortran_vec (); 1928 octave_idx_type *sel = s.fortran_vec ();
1929 1929
1930 Matrix eig_vec2 (n, k + 1); 1930 Matrix eig_vec2 (n, k + 1);
1934 OCTAVE_LOCAL_BUFFER (double, di, k + 1); 1934 OCTAVE_LOCAL_BUFFER (double, di, k + 1);
1935 OCTAVE_LOCAL_BUFFER (double, workev, 3 * p); 1935 OCTAVE_LOCAL_BUFFER (double, workev, 3 * p);
1936 for (octave_idx_type i = 0; i < k+1; i++) 1936 for (octave_idx_type i = 0; i < k+1; i++)
1937 dr[i] = di[i] = 0.; 1937 dr[i] = di[i] = 0.;
1938 1938
1939 F77_FUNC (dneupd, DNEUPD) 1939 F77_FUNC (dneupd, DNEUPD)
1940 (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, dr, di, z, n, sigmar, 1940 (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, dr, di, z, n, sigmar,
1941 sigmai, workev, F77_CONST_CHAR_ARG2 (&bmat, 1), n, 1941 sigmai, workev, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
1942 F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), k, tol, presid, p, v, n, iparam, 1942 F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), k, tol, presid, p, v, n, iparam,
1943 ipntr, workd, workl, lwork, info2 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) 1943 ipntr, workd, workl, lwork, info2 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1)
1944 F77_CHAR_ARG_LEN(2)); 1944 F77_CHAR_ARG_LEN(2));
1945 1945
1946 if (f77_exception_encountered) 1946 if (f77_exception_encountered)
1947 { 1947 {
1948 (*current_liboctave_error_handler) 1948 (*current_liboctave_error_handler)
2006 octave_idx_type off1 = i * n; 2006 octave_idx_type off1 = i * n;
2007 octave_idx_type off2 = (i+1) * n; 2007 octave_idx_type off2 = (i+1) * n;
2008 if (std::imag(eig_val(i)) == 0) 2008 if (std::imag(eig_val(i)) == 0)
2009 { 2009 {
2010 for (octave_idx_type j = 0; j < n; j++) 2010 for (octave_idx_type j = 0; j < n; j++)
2011 eig_vec(j,i) = 2011 eig_vec(j,i) =
2012 Complex(z[j+off1],0.); 2012 Complex(z[j+off1],0.);
2013 i++; 2013 i++;
2014 } 2014 }
2015 else 2015 else
2016 { 2016 {
2017 for (octave_idx_type j = 0; j < n; j++) 2017 for (octave_idx_type j = 0; j < n; j++)
2018 { 2018 {
2019 eig_vec(j,i) = 2019 eig_vec(j,i) =
2020 Complex(z[j+off1],z[j+off2]); 2020 Complex(z[j+off1],z[j+off2]);
2021 if (i < k - 1) 2021 if (i < k - 1)
2022 eig_vec(j,i+1) = 2022 eig_vec(j,i+1) =
2023 Complex(z[j+off1],-z[j+off2]); 2023 Complex(z[j+off1],-z[j+off2]);
2024 } 2024 }
2025 i+=2; 2025 i+=2;
2026 } 2026 }
2027 } 2027 }
2030 eig_vec = ltsolve(M (b), permB, eig_vec); 2030 eig_vec = ltsolve(M (b), permB, eig_vec);
2031 } 2031 }
2032 } 2032 }
2033 else 2033 else
2034 { 2034 {
2035 (*current_liboctave_error_handler) 2035 (*current_liboctave_error_handler)
2036 ("eigs: error %d in dneupd", info2); 2036 ("eigs: error %d in dneupd", info2);
2037 return -1; 2037 return -1;
2038 } 2038 }
2039 } 2039 }
2040 2040
2042 } 2042 }
2043 2043
2044 template <class M> 2044 template <class M>
2045 octave_idx_type 2045 octave_idx_type
2046 EigsRealNonSymmetricMatrixShift (const M& m, double sigmar, 2046 EigsRealNonSymmetricMatrixShift (const M& m, double sigmar,
2047 octave_idx_type k, octave_idx_type p, 2047 octave_idx_type k, octave_idx_type p,
2048 octave_idx_type &info, 2048 octave_idx_type &info,
2049 ComplexMatrix &eig_vec, 2049 ComplexMatrix &eig_vec,
2050 ComplexColumnVector &eig_val, const M& _b, 2050 ComplexColumnVector &eig_val, const M& _b,
2051 ColumnVector &permB, ColumnVector &resid, 2051 ColumnVector &permB, ColumnVector &resid,
2052 std::ostream& os, double tol, bool rvec, 2052 std::ostream& os, double tol, bool rvec,
2053 bool cholB, int disp, int maxit) 2053 bool cholB, int disp, int maxit)
2054 { 2054 {
2055 M b(_b); 2055 M b(_b);
2056 octave_idx_type n = m.cols (); 2056 octave_idx_type n = m.cols ();
2057 octave_idx_type mode = 3; 2057 octave_idx_type mode = 3;
2064 (*current_liboctave_error_handler) ("eigs: A must be square"); 2064 (*current_liboctave_error_handler) ("eigs: A must be square");
2065 return -1; 2065 return -1;
2066 } 2066 }
2067 if (have_b && (m.rows() != b.rows() || m.rows() != b.cols())) 2067 if (have_b && (m.rows() != b.rows() || m.rows() != b.cols()))
2068 { 2068 {
2069 (*current_liboctave_error_handler) 2069 (*current_liboctave_error_handler)
2070 ("eigs: B must be square and the same size as A"); 2070 ("eigs: B must be square and the same size as A");
2071 return -1; 2071 return -1;
2072 } 2072 }
2073 2073
2074 // FIXME: The "SM" type for mode 1 seems unstable though faster!! 2074 // FIXME: The "SM" type for mode 1 seems unstable though faster!!
2096 { 2096 {
2097 p = k * 2 + 1; 2097 p = k * 2 + 1;
2098 2098
2099 if (p < 20) 2099 if (p < 20)
2100 p = 20; 2100 p = 20;
2101 2101
2102 if (p > n - 1) 2102 if (p > n - 1)
2103 p = n - 1 ; 2103 p = n - 1 ;
2104 } 2104 }
2105 2105
2106 if (k <= 0 || k >= n - 1) 2106 if (k <= 0 || k >= n - 1)
2107 { 2107 {
2108 (*current_liboctave_error_handler) 2108 (*current_liboctave_error_handler)
2109 ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1).\n" 2109 ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1).\n"
2110 " Use 'eig(full(A))' instead"); 2110 " Use 'eig(full(A))' instead");
2111 return -1; 2111 return -1;
2112 } 2112 }
2113 2113
2114 if (p <= k || p >= n) 2114 if (p <= k || p >= n)
2115 { 2115 {
2116 (*current_liboctave_error_handler) 2116 (*current_liboctave_error_handler)
2117 ("eigs: opts.p must be greater than k and less than n"); 2117 ("eigs: opts.p must be greater than k and less than n");
2118 return -1; 2118 return -1;
2119 } 2119 }
2120 2120
2121 if (have_b && cholB && permB.length() != 0) 2121 if (have_b && cholB && permB.length() != 0)
2122 { 2122 {
2123 // Check that we really have a permutation vector 2123 // Check that we really have a permutation vector
2124 if (permB.length() != n) 2124 if (permB.length() != n)
2125 { 2125 {
2126 (*current_liboctave_error_handler) ("eigs: permB vector invalid"); 2126 (*current_liboctave_error_handler) ("eigs: permB vector invalid");
2129 else 2129 else
2130 { 2130 {
2131 Array<bool> checked (dim_vector (n, 1), false); 2131 Array<bool> checked (dim_vector (n, 1), false);
2132 for (octave_idx_type i = 0; i < n; i++) 2132 for (octave_idx_type i = 0; i < n; i++)
2133 { 2133 {
2134 octave_idx_type bidx = 2134 octave_idx_type bidx =
2135 static_cast<octave_idx_type> (permB(i)); 2135 static_cast<octave_idx_type> (permB(i));
2136 if (checked(bidx) || bidx < 0 || 2136 if (checked(bidx) || bidx < 0 ||
2137 bidx >= n || D_NINT (bidx) != bidx) 2137 bidx >= n || D_NINT (bidx) != bidx)
2138 { 2138 {
2139 (*current_liboctave_error_handler) 2139 (*current_liboctave_error_handler)
2140 ("eigs: permB vector invalid"); 2140 ("eigs: permB vector invalid");
2141 return -1; 2141 return -1;
2142 } 2142 }
2143 } 2143 }
2144 } 2144 }
2182 OCTAVE_LOCAL_BUFFER (double, v, n * (p + 1)); 2182 OCTAVE_LOCAL_BUFFER (double, v, n * (p + 1));
2183 OCTAVE_LOCAL_BUFFER (double, workl, lwork + 1); 2183 OCTAVE_LOCAL_BUFFER (double, workl, lwork + 1);
2184 OCTAVE_LOCAL_BUFFER (double, workd, 3 * n + 1); 2184 OCTAVE_LOCAL_BUFFER (double, workd, 3 * n + 1);
2185 double *presid = resid.fortran_vec (); 2185 double *presid = resid.fortran_vec ();
2186 2186
2187 do 2187 do
2188 { 2188 {
2189 F77_FUNC (dnaupd, DNAUPD) 2189 F77_FUNC (dnaupd, DNAUPD)
2190 (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, 2190 (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
2191 F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), 2191 F77_CONST_CHAR_ARG2 ((typ.c_str()), 2),
2192 k, tol, presid, p, v, n, iparam, 2192 k, tol, presid, p, v, n, iparam,
2193 ipntr, workd, workl, lwork, info 2193 ipntr, workd, workl, lwork, info
2194 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); 2194 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
2195 2195
2196 if (f77_exception_encountered) 2196 if (f77_exception_encountered)
2197 { 2197 {
2198 (*current_liboctave_error_handler) 2198 (*current_liboctave_error_handler)
2199 ("eigs: unrecoverable exception encountered in dsaupd"); 2199 ("eigs: unrecoverable exception encountered in dsaupd");
2200 return -1; 2200 return -1;
2201 } 2201 }
2202 2202
2203 if (disp > 0 && !xisnan(workl[iptr(5)-1])) 2203 if (disp > 0 && !xisnan(workl[iptr(5)-1]))
2204 { 2204 {
2205 if (iter++) 2205 if (iter++)
2206 { 2206 {
2207 os << "Iteration " << iter - 1 << 2207 os << "Iteration " << iter - 1 <<
2208 ": a few Ritz values of the " << p << "-by-" << 2208 ": a few Ritz values of the " << p << "-by-" <<
2209 p << " matrix\n"; 2209 p << " matrix\n";
2210 for (int i = 0 ; i < k; i++) 2210 for (int i = 0 ; i < k; i++)
2211 os << " " << workl[iptr(5)+i-1] << "\n"; 2211 os << " " << workl[iptr(5)+i-1] << "\n";
2212 } 2212 }
2215 // iteration pointer. But as workl[iptr(5)-1] is 2215 // iteration pointer. But as workl[iptr(5)-1] is
2216 // an output value updated at each iteration, setting 2216 // an output value updated at each iteration, setting
2217 // a value in this array to NaN and testing for it 2217 // a value in this array to NaN and testing for it
2218 // is a way of obtaining the iteration counter. 2218 // is a way of obtaining the iteration counter.
2219 if (ido != 99) 2219 if (ido != 99)
2220 workl[iptr(5)-1] = octave_NaN; 2220 workl[iptr(5)-1] = octave_NaN;
2221 } 2221 }
2222 2222
2223 if (ido == -1 || ido == 1 || ido == 2) 2223 if (ido == -1 || ido == 1 || ido == 2)
2224 { 2224 {
2225 if (have_b) 2225 if (have_b)
2232 2232
2233 Matrix tmp(n, 1); 2233 Matrix tmp(n, 1);
2234 2234
2235 for (octave_idx_type i = 0; i < n; i++) 2235 for (octave_idx_type i = 0; i < n; i++)
2236 tmp(i,0) = dtmp[P[i]]; 2236 tmp(i,0) = dtmp[P[i]];
2237 2237
2238 lusolve (L, U, tmp); 2238 lusolve (L, U, tmp);
2239 2239
2240 double *ip2 = workd+iptr(1)-1; 2240 double *ip2 = workd+iptr(1)-1;
2241 for (octave_idx_type i = 0; i < n; i++) 2241 for (octave_idx_type i = 0; i < n; i++)
2242 ip2[Q[i]] = tmp(i,0); 2242 ip2[Q[i]] = tmp(i,0);
2248 double *ip2 = workd+iptr(2)-1; 2248 double *ip2 = workd+iptr(2)-1;
2249 Matrix tmp(n, 1); 2249 Matrix tmp(n, 1);
2250 2250
2251 for (octave_idx_type i = 0; i < n; i++) 2251 for (octave_idx_type i = 0; i < n; i++)
2252 tmp(i,0) = ip2[P[i]]; 2252 tmp(i,0) = ip2[P[i]];
2253 2253
2254 lusolve (L, U, tmp); 2254 lusolve (L, U, tmp);
2255 2255
2256 ip2 = workd+iptr(1)-1; 2256 ip2 = workd+iptr(1)-1;
2257 for (octave_idx_type i = 0; i < n; i++) 2257 for (octave_idx_type i = 0; i < n; i++)
2258 ip2[Q[i]] = tmp(i,0); 2258 ip2[Q[i]] = tmp(i,0);
2270 double *ip2 = workd+iptr(0)-1; 2270 double *ip2 = workd+iptr(0)-1;
2271 Matrix tmp(n, 1); 2271 Matrix tmp(n, 1);
2272 2272
2273 for (octave_idx_type i = 0; i < n; i++) 2273 for (octave_idx_type i = 0; i < n; i++)
2274 tmp(i,0) = ip2[P[i]]; 2274 tmp(i,0) = ip2[P[i]];
2275 2275
2276 lusolve (L, U, tmp); 2276 lusolve (L, U, tmp);
2277 2277
2278 ip2 = workd+iptr(1)-1; 2278 ip2 = workd+iptr(1)-1;
2279 for (octave_idx_type i = 0; i < n; i++) 2279 for (octave_idx_type i = 0; i < n; i++)
2280 ip2[Q[i]] = tmp(i,0); 2280 ip2[Q[i]] = tmp(i,0);
2283 } 2283 }
2284 else 2284 else
2285 { 2285 {
2286 if (info < 0) 2286 if (info < 0)
2287 { 2287 {
2288 (*current_liboctave_error_handler) 2288 (*current_liboctave_error_handler)
2289 ("eigs: error %d in dsaupd", info); 2289 ("eigs: error %d in dsaupd", info);
2290 return -1; 2290 return -1;
2291 } 2291 }
2292 break; 2292 break;
2293 } 2293 }
2294 } 2294 }
2295 while (1); 2295 while (1);
2296 2296
2297 octave_idx_type info2; 2297 octave_idx_type info2;
2298 2298
2299 // We have a problem in that the size of the C++ bool 2299 // We have a problem in that the size of the C++ bool
2300 // type relative to the fortran logical type. It appears 2300 // type relative to the fortran logical type. It appears
2301 // that fortran uses 4- or 8-bytes per logical and C++ 1-byte 2301 // that fortran uses 4- or 8-bytes per logical and C++ 1-byte
2302 // per bool, though this might be system dependent. As 2302 // per bool, though this might be system dependent. As
2303 // long as the HOWMNY arg is not "S", the logical array 2303 // long as the HOWMNY arg is not "S", the logical array
2304 // is just workspace for ARPACK, so use int type to 2304 // is just workspace for ARPACK, so use int type to
2305 // avoid problems. 2305 // avoid problems.
2306 Array<octave_idx_type> s (dim_vector (p, 1)); 2306 Array<octave_idx_type> s (dim_vector (p, 1));
2307 octave_idx_type *sel = s.fortran_vec (); 2307 octave_idx_type *sel = s.fortran_vec ();
2308 2308
2309 Matrix eig_vec2 (n, k + 1); 2309 Matrix eig_vec2 (n, k + 1);
2310 double *z = eig_vec2.fortran_vec (); 2310 double *z = eig_vec2.fortran_vec ();
2311 2311
2312 OCTAVE_LOCAL_BUFFER (double, dr, k + 1); 2312 OCTAVE_LOCAL_BUFFER (double, dr, k + 1);
2313 OCTAVE_LOCAL_BUFFER (double, di, k + 1); 2313 OCTAVE_LOCAL_BUFFER (double, di, k + 1);
2314 OCTAVE_LOCAL_BUFFER (double, workev, 3 * p); 2314 OCTAVE_LOCAL_BUFFER (double, workev, 3 * p);
2315 for (octave_idx_type i = 0; i < k+1; i++) 2315 for (octave_idx_type i = 0; i < k+1; i++)
2316 dr[i] = di[i] = 0.; 2316 dr[i] = di[i] = 0.;
2317 2317
2318 F77_FUNC (dneupd, DNEUPD) 2318 F77_FUNC (dneupd, DNEUPD)
2319 (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, dr, di, z, n, sigmar, 2319 (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, dr, di, z, n, sigmar,
2320 sigmai, workev, F77_CONST_CHAR_ARG2 (&bmat, 1), n, 2320 sigmai, workev, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
2321 F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), k, tol, presid, p, v, n, iparam, 2321 F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), k, tol, presid, p, v, n, iparam,
2322 ipntr, workd, workl, lwork, info2 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) 2322 ipntr, workd, workl, lwork, info2 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1)
2323 F77_CHAR_ARG_LEN(2)); 2323 F77_CHAR_ARG_LEN(2));
2324 2324
2325 if (f77_exception_encountered) 2325 if (f77_exception_encountered)
2326 { 2326 {
2327 (*current_liboctave_error_handler) 2327 (*current_liboctave_error_handler)
2385 octave_idx_type off1 = i * n; 2385 octave_idx_type off1 = i * n;
2386 octave_idx_type off2 = (i+1) * n; 2386 octave_idx_type off2 = (i+1) * n;
2387 if (std::imag(eig_val(i)) == 0) 2387 if (std::imag(eig_val(i)) == 0)
2388 { 2388 {
2389 for (octave_idx_type j = 0; j < n; j++) 2389 for (octave_idx_type j = 0; j < n; j++)
2390 eig_vec(j,i) = 2390 eig_vec(j,i) =
2391 Complex(z[j+off1],0.); 2391 Complex(z[j+off1],0.);
2392 i++; 2392 i++;
2393 } 2393 }
2394 else 2394 else
2395 { 2395 {
2396 for (octave_idx_type j = 0; j < n; j++) 2396 for (octave_idx_type j = 0; j < n; j++)
2397 { 2397 {
2398 eig_vec(j,i) = 2398 eig_vec(j,i) =
2399 Complex(z[j+off1],z[j+off2]); 2399 Complex(z[j+off1],z[j+off2]);
2400 if (i < k - 1) 2400 if (i < k - 1)
2401 eig_vec(j,i+1) = 2401 eig_vec(j,i+1) =
2402 Complex(z[j+off1],-z[j+off2]); 2402 Complex(z[j+off1],-z[j+off2]);
2403 } 2403 }
2404 i+=2; 2404 i+=2;
2405 } 2405 }
2406 } 2406 }
2407 } 2407 }
2408 } 2408 }
2409 else 2409 else
2410 { 2410 {
2411 (*current_liboctave_error_handler) 2411 (*current_liboctave_error_handler)
2412 ("eigs: error %d in dneupd", info2); 2412 ("eigs: error %d in dneupd", info2);
2413 return -1; 2413 return -1;
2414 } 2414 }
2415 } 2415 }
2416 2416
2418 } 2418 }
2419 2419
2420 octave_idx_type 2420 octave_idx_type
2421 EigsRealNonSymmetricFunc (EigsFunc fun, octave_idx_type n, 2421 EigsRealNonSymmetricFunc (EigsFunc fun, octave_idx_type n,
2422 const std::string &_typ, double sigmar, 2422 const std::string &_typ, double sigmar,
2423 octave_idx_type k, octave_idx_type p, 2423 octave_idx_type k, octave_idx_type p,
2424 octave_idx_type &info, ComplexMatrix &eig_vec, 2424 octave_idx_type &info, ComplexMatrix &eig_vec,
2425 ComplexColumnVector &eig_val, ColumnVector &resid, 2425 ComplexColumnVector &eig_val, ColumnVector &resid,
2426 std::ostream& os, double tol, bool rvec, 2426 std::ostream& os, double tol, bool rvec,
2427 bool /* cholB */, int disp, int maxit) 2427 bool /* cholB */, int disp, int maxit)
2428 { 2428 {
2429 std::string typ (_typ); 2429 std::string typ (_typ);
2430 bool have_sigma = (sigmar ? true : false); 2430 bool have_sigma = (sigmar ? true : false);
2452 { 2452 {
2453 p = k * 2 + 1; 2453 p = k * 2 + 1;
2454 2454
2455 if (p < 20) 2455 if (p < 20)
2456 p = 20; 2456 p = 20;
2457 2457
2458 if (p > n - 1) 2458 if (p > n - 1)
2459 p = n - 1 ; 2459 p = n - 1 ;
2460 } 2460 }
2461 2461
2462 if (k <= 0 || k >= n - 1) 2462 if (k <= 0 || k >= n - 1)
2475 } 2475 }
2476 2476
2477 2477
2478 if (! have_sigma) 2478 if (! have_sigma)
2479 { 2479 {
2480 if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" && 2480 if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" &&
2481 typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" && 2481 typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" &&
2482 typ != "SI") 2482 typ != "SI")
2483 (*current_liboctave_error_handler) 2483 (*current_liboctave_error_handler)
2484 ("eigs: unrecognized sigma value"); 2484 ("eigs: unrecognized sigma value");
2485 2485
2486 if (typ == "LA" || typ == "SA" || typ == "BE") 2486 if (typ == "LA" || typ == "SA" || typ == "BE")
2487 { 2487 {
2488 (*current_liboctave_error_handler) 2488 (*current_liboctave_error_handler)
2489 ("eigs: invalid sigma value for unsymmetric problem"); 2489 ("eigs: invalid sigma value for unsymmetric problem");
2490 return -1; 2490 return -1;
2491 } 2491 }
2492 2492
2493 if (typ == "SM") 2493 if (typ == "SM")
2518 ip(7) = 0; 2518 ip(7) = 0;
2519 ip(8) = 0; 2519 ip(8) = 0;
2520 ip(9) = 0; 2520 ip(9) = 0;
2521 ip(10) = 0; 2521 ip(10) = 0;
2522 // ip(7) to ip(10) return values 2522 // ip(7) to ip(10) return values
2523 2523
2524 Array<octave_idx_type> iptr (dim_vector (14, 1)); 2524 Array<octave_idx_type> iptr (dim_vector (14, 1));
2525 octave_idx_type *ipntr = iptr.fortran_vec (); 2525 octave_idx_type *ipntr = iptr.fortran_vec ();
2526 2526
2527 octave_idx_type ido = 0; 2527 octave_idx_type ido = 0;
2528 int iter = 0; 2528 int iter = 0;
2531 OCTAVE_LOCAL_BUFFER (double, v, n * (p + 1)); 2531 OCTAVE_LOCAL_BUFFER (double, v, n * (p + 1));
2532 OCTAVE_LOCAL_BUFFER (double, workl, lwork + 1); 2532 OCTAVE_LOCAL_BUFFER (double, workl, lwork + 1);
2533 OCTAVE_LOCAL_BUFFER (double, workd, 3 * n + 1); 2533 OCTAVE_LOCAL_BUFFER (double, workd, 3 * n + 1);
2534 double *presid = resid.fortran_vec (); 2534 double *presid = resid.fortran_vec ();
2535 2535
2536 do 2536 do
2537 { 2537 {
2538 F77_FUNC (dnaupd, DNAUPD) 2538 F77_FUNC (dnaupd, DNAUPD)
2539 (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, 2539 (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
2540 F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), 2540 F77_CONST_CHAR_ARG2 ((typ.c_str()), 2),
2541 k, tol, presid, p, v, n, iparam, 2541 k, tol, presid, p, v, n, iparam,
2542 ipntr, workd, workl, lwork, info 2542 ipntr, workd, workl, lwork, info
2543 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); 2543 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
2544 2544
2545 if (f77_exception_encountered) 2545 if (f77_exception_encountered)
2546 { 2546 {
2547 (*current_liboctave_error_handler) 2547 (*current_liboctave_error_handler)
2548 ("eigs: unrecoverable exception encountered in dnaupd"); 2548 ("eigs: unrecoverable exception encountered in dnaupd");
2549 return -1; 2549 return -1;
2550 } 2550 }
2551 2551
2552 if (disp > 0 && !xisnan(workl[iptr(5)-1])) 2552 if (disp > 0 && !xisnan(workl[iptr(5)-1]))
2553 { 2553 {
2554 if (iter++) 2554 if (iter++)
2555 { 2555 {
2556 os << "Iteration " << iter - 1 << 2556 os << "Iteration " << iter - 1 <<
2557 ": a few Ritz values of the " << p << "-by-" << 2557 ": a few Ritz values of the " << p << "-by-" <<
2558 p << " matrix\n"; 2558 p << " matrix\n";
2559 for (int i = 0 ; i < k; i++) 2559 for (int i = 0 ; i < k; i++)
2560 os << " " << workl[iptr(5)+i-1] << "\n"; 2560 os << " " << workl[iptr(5)+i-1] << "\n";
2561 } 2561 }
2564 // iteration pointer. But as workl[iptr(5)-1] is 2564 // iteration pointer. But as workl[iptr(5)-1] is
2565 // an output value updated at each iteration, setting 2565 // an output value updated at each iteration, setting
2566 // a value in this array to NaN and testing for it 2566 // a value in this array to NaN and testing for it
2567 // is a way of obtaining the iteration counter. 2567 // is a way of obtaining the iteration counter.
2568 if (ido != 99) 2568 if (ido != 99)
2569 workl[iptr(5)-1] = octave_NaN; 2569 workl[iptr(5)-1] = octave_NaN;
2570 } 2570 }
2571 2571
2572 if (ido == -1 || ido == 1 || ido == 2) 2572 if (ido == -1 || ido == 1 || ido == 2)
2573 { 2573 {
2574 double *ip2 = workd + iptr(0) - 1; 2574 double *ip2 = workd + iptr(0) - 1;
2588 } 2588 }
2589 else 2589 else
2590 { 2590 {
2591 if (info < 0) 2591 if (info < 0)
2592 { 2592 {
2593 (*current_liboctave_error_handler) 2593 (*current_liboctave_error_handler)
2594 ("eigs: error %d in dsaupd", info); 2594 ("eigs: error %d in dsaupd", info);
2595 return -1; 2595 return -1;
2596 } 2596 }
2597 break; 2597 break;
2598 } 2598 }
2599 } 2599 }
2600 while (1); 2600 while (1);
2601 2601
2602 octave_idx_type info2; 2602 octave_idx_type info2;
2603 2603
2604 // We have a problem in that the size of the C++ bool 2604 // We have a problem in that the size of the C++ bool
2605 // type relative to the fortran logical type. It appears 2605 // type relative to the fortran logical type. It appears
2606 // that fortran uses 4- or 8-bytes per logical and C++ 1-byte 2606 // that fortran uses 4- or 8-bytes per logical and C++ 1-byte
2607 // per bool, though this might be system dependent. As 2607 // per bool, though this might be system dependent. As
2608 // long as the HOWMNY arg is not "S", the logical array 2608 // long as the HOWMNY arg is not "S", the logical array
2609 // is just workspace for ARPACK, so use int type to 2609 // is just workspace for ARPACK, so use int type to
2610 // avoid problems. 2610 // avoid problems.
2611 Array<octave_idx_type> s (dim_vector (p, 1)); 2611 Array<octave_idx_type> s (dim_vector (p, 1));
2612 octave_idx_type *sel = s.fortran_vec (); 2612 octave_idx_type *sel = s.fortran_vec ();
2613 2613
2614 Matrix eig_vec2 (n, k + 1); 2614 Matrix eig_vec2 (n, k + 1);
2618 OCTAVE_LOCAL_BUFFER (double, di, k + 1); 2618 OCTAVE_LOCAL_BUFFER (double, di, k + 1);
2619 OCTAVE_LOCAL_BUFFER (double, workev, 3 * p); 2619 OCTAVE_LOCAL_BUFFER (double, workev, 3 * p);
2620 for (octave_idx_type i = 0; i < k+1; i++) 2620 for (octave_idx_type i = 0; i < k+1; i++)
2621 dr[i] = di[i] = 0.; 2621 dr[i] = di[i] = 0.;
2622 2622
2623 F77_FUNC (dneupd, DNEUPD) 2623 F77_FUNC (dneupd, DNEUPD)
2624 (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, dr, di, z, n, sigmar, 2624 (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, dr, di, z, n, sigmar,
2625 sigmai, workev, F77_CONST_CHAR_ARG2 (&bmat, 1), n, 2625 sigmai, workev, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
2626 F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), k, tol, presid, p, v, n, iparam, 2626 F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), k, tol, presid, p, v, n, iparam,
2627 ipntr, workd, workl, lwork, info2 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) 2627 ipntr, workd, workl, lwork, info2 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1)
2628 F77_CHAR_ARG_LEN(2)); 2628 F77_CHAR_ARG_LEN(2));
2629 2629
2630 if (f77_exception_encountered) 2630 if (f77_exception_encountered)
2631 { 2631 {
2632 (*current_liboctave_error_handler) 2632 (*current_liboctave_error_handler)
2690 octave_idx_type off1 = i * n; 2690 octave_idx_type off1 = i * n;
2691 octave_idx_type off2 = (i+1) * n; 2691 octave_idx_type off2 = (i+1) * n;
2692 if (std::imag(eig_val(i)) == 0) 2692 if (std::imag(eig_val(i)) == 0)
2693 { 2693 {
2694 for (octave_idx_type j = 0; j < n; j++) 2694 for (octave_idx_type j = 0; j < n; j++)
2695 eig_vec(j,i) = 2695 eig_vec(j,i) =
2696 Complex(z[j+off1],0.); 2696 Complex(z[j+off1],0.);
2697 i++; 2697 i++;
2698 } 2698 }
2699 else 2699 else
2700 { 2700 {
2701 for (octave_idx_type j = 0; j < n; j++) 2701 for (octave_idx_type j = 0; j < n; j++)
2702 { 2702 {
2703 eig_vec(j,i) = 2703 eig_vec(j,i) =
2704 Complex(z[j+off1],z[j+off2]); 2704 Complex(z[j+off1],z[j+off2]);
2705 if (i < k - 1) 2705 if (i < k - 1)
2706 eig_vec(j,i+1) = 2706 eig_vec(j,i+1) =
2707 Complex(z[j+off1],-z[j+off2]); 2707 Complex(z[j+off1],-z[j+off2]);
2708 } 2708 }
2709 i+=2; 2709 i+=2;
2710 } 2710 }
2711 } 2711 }
2712 } 2712 }
2713 } 2713 }
2714 else 2714 else
2715 { 2715 {
2716 (*current_liboctave_error_handler) 2716 (*current_liboctave_error_handler)
2717 ("eigs: error %d in dneupd", info2); 2717 ("eigs: error %d in dneupd", info2);
2718 return -1; 2718 return -1;
2719 } 2719 }
2720 } 2720 }
2721 2721
2722 return ip(4); 2722 return ip(4);
2723 } 2723 }
2724 2724
2725 template <class M> 2725 template <class M>
2726 octave_idx_type 2726 octave_idx_type
2727 EigsComplexNonSymmetricMatrix (const M& m, const std::string typ, 2727 EigsComplexNonSymmetricMatrix (const M& m, const std::string typ,
2728 octave_idx_type k, octave_idx_type p, 2728 octave_idx_type k, octave_idx_type p,
2729 octave_idx_type &info, ComplexMatrix &eig_vec, 2729 octave_idx_type &info, ComplexMatrix &eig_vec,
2730 ComplexColumnVector &eig_val, const M& _b, 2730 ComplexColumnVector &eig_val, const M& _b,
2731 ColumnVector &permB, 2731 ColumnVector &permB,
2732 ComplexColumnVector &cresid, 2732 ComplexColumnVector &cresid,
2733 std::ostream& os, double tol, bool rvec, 2733 std::ostream& os, double tol, bool rvec,
2734 bool cholB, int disp, int maxit) 2734 bool cholB, int disp, int maxit)
2735 { 2735 {
2736 M b(_b); 2736 M b(_b);
2737 octave_idx_type n = m.cols (); 2737 octave_idx_type n = m.cols ();
2738 octave_idx_type mode = 1; 2738 octave_idx_type mode = 1;
2747 (*current_liboctave_error_handler) ("eigs: A must be square"); 2747 (*current_liboctave_error_handler) ("eigs: A must be square");
2748 return -1; 2748 return -1;
2749 } 2749 }
2750 if (have_b && (m.rows() != b.rows() || m.rows() != b.cols())) 2750 if (have_b && (m.rows() != b.rows() || m.rows() != b.cols()))
2751 { 2751 {
2752 (*current_liboctave_error_handler) 2752 (*current_liboctave_error_handler)
2753 ("eigs: B must be square and the same size as A"); 2753 ("eigs: B must be square and the same size as A");
2754 return -1; 2754 return -1;
2755 } 2755 }
2756 2756
2757 if (cresid.is_empty()) 2757 if (cresid.is_empty())
2777 { 2777 {
2778 p = k * 2 + 1; 2778 p = k * 2 + 1;
2779 2779
2780 if (p < 20) 2780 if (p < 20)
2781 p = 20; 2781 p = 20;
2782 2782
2783 if (p > n - 1) 2783 if (p > n - 1)
2784 p = n - 1 ; 2784 p = n - 1 ;
2785 } 2785 }
2786 2786
2787 if (k <= 0 || k >= n - 1) 2787 if (k <= 0 || k >= n - 1)
2788 { 2788 {
2789 (*current_liboctave_error_handler) 2789 (*current_liboctave_error_handler)
2790 ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1).\n" 2790 ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1).\n"
2791 " Use 'eig(full(A))' instead"); 2791 " Use 'eig(full(A))' instead");
2792 return -1; 2792 return -1;
2793 } 2793 }
2794 2794
2795 if (p <= k || p >= n) 2795 if (p <= k || p >= n)
2796 { 2796 {
2797 (*current_liboctave_error_handler) 2797 (*current_liboctave_error_handler)
2798 ("eigs: opts.p must be greater than k and less than n"); 2798 ("eigs: opts.p must be greater than k and less than n");
2799 return -1; 2799 return -1;
2800 } 2800 }
2801 2801
2802 if (have_b && cholB && permB.length() != 0) 2802 if (have_b && cholB && permB.length() != 0)
2803 { 2803 {
2804 // Check the we really have a permutation vector 2804 // Check the we really have a permutation vector
2805 if (permB.length() != n) 2805 if (permB.length() != n)
2806 { 2806 {
2807 (*current_liboctave_error_handler) 2807 (*current_liboctave_error_handler)
2808 ("eigs: permB vector invalid"); 2808 ("eigs: permB vector invalid");
2809 return -1; 2809 return -1;
2810 } 2810 }
2811 else 2811 else
2812 { 2812 {
2813 Array<bool> checked (dim_vector (n, 1), false); 2813 Array<bool> checked (dim_vector (n, 1), false);
2814 for (octave_idx_type i = 0; i < n; i++) 2814 for (octave_idx_type i = 0; i < n; i++)
2815 { 2815 {
2816 octave_idx_type bidx = 2816 octave_idx_type bidx =
2817 static_cast<octave_idx_type> (permB(i)); 2817 static_cast<octave_idx_type> (permB(i));
2818 if (checked(bidx) || bidx < 0 || 2818 if (checked(bidx) || bidx < 0 ||
2819 bidx >= n || D_NINT (bidx) != bidx) 2819 bidx >= n || D_NINT (bidx) != bidx)
2820 { 2820 {
2821 (*current_liboctave_error_handler) 2821 (*current_liboctave_error_handler)
2822 ("eigs: permB vector invalid"); 2822 ("eigs: permB vector invalid");
2823 return -1; 2823 return -1;
2824 } 2824 }
2825 } 2825 }
2826 } 2826 }
2827 } 2827 }
2828 2828
2829 if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" && 2829 if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" &&
2830 typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" && 2830 typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" &&
2831 typ != "SI") 2831 typ != "SI")
2832 { 2832 {
2833 (*current_liboctave_error_handler) 2833 (*current_liboctave_error_handler)
2834 ("eigs: unrecognized sigma value"); 2834 ("eigs: unrecognized sigma value");
2835 return -1; 2835 return -1;
2836 } 2836 }
2837 2837
2838 if (typ == "LA" || typ == "SA" || typ == "BE") 2838 if (typ == "LA" || typ == "SA" || typ == "BE")
2839 { 2839 {
2840 (*current_liboctave_error_handler) 2840 (*current_liboctave_error_handler)
2841 ("eigs: invalid sigma value for complex problem"); 2841 ("eigs: invalid sigma value for complex problem");
2842 return -1; 2842 return -1;
2843 } 2843 }
2844 2844
2845 if (have_b) 2845 if (have_b)
2859 } 2859 }
2860 else 2860 else
2861 { 2861 {
2862 if (! make_cholb(b, bt, permB)) 2862 if (! make_cholb(b, bt, permB))
2863 { 2863 {
2864 (*current_liboctave_error_handler) 2864 (*current_liboctave_error_handler)
2865 ("eigs: The matrix B is not positive definite"); 2865 ("eigs: The matrix B is not positive definite");
2866 return -1; 2866 return -1;
2867 } 2867 }
2868 } 2868 }
2869 } 2869 }
2881 ip(7) = 0; 2881 ip(7) = 0;
2882 ip(8) = 0; 2882 ip(8) = 0;
2883 ip(9) = 0; 2883 ip(9) = 0;
2884 ip(10) = 0; 2884 ip(10) = 0;
2885 // ip(7) to ip(10) return values 2885 // ip(7) to ip(10) return values
2886 2886
2887 Array<octave_idx_type> iptr (dim_vector (14, 1)); 2887 Array<octave_idx_type> iptr (dim_vector (14, 1));
2888 octave_idx_type *ipntr = iptr.fortran_vec (); 2888 octave_idx_type *ipntr = iptr.fortran_vec ();
2889 2889
2890 octave_idx_type ido = 0; 2890 octave_idx_type ido = 0;
2891 int iter = 0; 2891 int iter = 0;
2892 octave_idx_type lwork = p * (3 * p + 5); 2892 octave_idx_type lwork = p * (3 * p + 5);
2893 2893
2894 OCTAVE_LOCAL_BUFFER (Complex, v, n * p); 2894 OCTAVE_LOCAL_BUFFER (Complex, v, n * p);
2895 OCTAVE_LOCAL_BUFFER (Complex, workl, lwork); 2895 OCTAVE_LOCAL_BUFFER (Complex, workl, lwork);
2896 OCTAVE_LOCAL_BUFFER (Complex, workd, 3 * n); 2896 OCTAVE_LOCAL_BUFFER (Complex, workd, 3 * n);
2897 OCTAVE_LOCAL_BUFFER (double, rwork, p); 2897 OCTAVE_LOCAL_BUFFER (double, rwork, p);
2898 Complex *presid = cresid.fortran_vec (); 2898 Complex *presid = cresid.fortran_vec ();
2899 2899
2900 do 2900 do
2901 { 2901 {
2902 F77_FUNC (znaupd, ZNAUPD) 2902 F77_FUNC (znaupd, ZNAUPD)
2903 (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, 2903 (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
2904 F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), 2904 F77_CONST_CHAR_ARG2 ((typ.c_str()), 2),
2905 k, tol, presid, p, v, n, iparam, 2905 k, tol, presid, p, v, n, iparam,
2906 ipntr, workd, workl, lwork, rwork, info 2906 ipntr, workd, workl, lwork, rwork, info
2907 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); 2907 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
2908 2908
2909 if (f77_exception_encountered) 2909 if (f77_exception_encountered)
2910 { 2910 {
2911 (*current_liboctave_error_handler) 2911 (*current_liboctave_error_handler)
2912 ("eigs: unrecoverable exception encountered in znaupd"); 2912 ("eigs: unrecoverable exception encountered in znaupd");
2913 return -1; 2913 return -1;
2914 } 2914 }
2915 2915
2916 if (disp > 0 && !xisnan(workl[iptr(5)-1])) 2916 if (disp > 0 && !xisnan(workl[iptr(5)-1]))
2917 { 2917 {
2918 if (iter++) 2918 if (iter++)
2919 { 2919 {
2920 os << "Iteration " << iter - 1 << 2920 os << "Iteration " << iter - 1 <<
2921 ": a few Ritz values of the " << p << "-by-" << 2921 ": a few Ritz values of the " << p << "-by-" <<
2922 p << " matrix\n"; 2922 p << " matrix\n";
2923 for (int i = 0 ; i < k; i++) 2923 for (int i = 0 ; i < k; i++)
2924 os << " " << workl[iptr(5)+i-1] << "\n"; 2924 os << " " << workl[iptr(5)+i-1] << "\n";
2925 } 2925 }
2926 2926
2927 // This is a kludge, as ARPACK doesn't give its 2927 // This is a kludge, as ARPACK doesn't give its
2928 // iteration pointer. But as workl[iptr(5)-1] is 2928 // iteration pointer. But as workl[iptr(5)-1] is
2929 // an output value updated at each iteration, setting 2929 // an output value updated at each iteration, setting
2930 // a value in this array to NaN and testing for it 2930 // a value in this array to NaN and testing for it
2931 // is a way of obtaining the iteration counter. 2931 // is a way of obtaining the iteration counter.
2932 if (ido != 99) 2932 if (ido != 99)
2933 workl[iptr(5)-1] = octave_NaN; 2933 workl[iptr(5)-1] = octave_NaN;
2934 } 2934 }
2935 2935
2936 if (ido == -1 || ido == 1 || ido == 2) 2936 if (ido == -1 || ido == 1 || ido == 2)
2937 { 2937 {
2938 if (have_b) 2938 if (have_b)
2943 mtmp = utsolve(bt, permB, m * ltsolve(b, permB, mtmp)); 2943 mtmp = utsolve(bt, permB, m * ltsolve(b, permB, mtmp));
2944 for (octave_idx_type i = 0; i < n; i++) 2944 for (octave_idx_type i = 0; i < n; i++)
2945 workd[i+iptr(1)-1] = mtmp(i,0); 2945 workd[i+iptr(1)-1] = mtmp(i,0);
2946 2946
2947 } 2947 }
2948 else if (!vector_product (m, workd + iptr(0) - 1, 2948 else if (!vector_product (m, workd + iptr(0) - 1,
2949 workd + iptr(1) - 1)) 2949 workd + iptr(1) - 1))
2950 break; 2950 break;
2951 } 2951 }
2952 else 2952 else
2953 { 2953 {
2954 if (info < 0) 2954 if (info < 0)
2955 { 2955 {
2956 (*current_liboctave_error_handler) 2956 (*current_liboctave_error_handler)
2957 ("eigs: error %d in znaupd", info); 2957 ("eigs: error %d in znaupd", info);
2958 return -1; 2958 return -1;
2959 } 2959 }
2960 break; 2960 break;
2961 } 2961 }
2962 } 2962 }
2963 while (1); 2963 while (1);
2964 2964
2965 octave_idx_type info2; 2965 octave_idx_type info2;
2966 2966
2967 // We have a problem in that the size of the C++ bool 2967 // We have a problem in that the size of the C++ bool
2968 // type relative to the fortran logical type. It appears 2968 // type relative to the fortran logical type. It appears
2969 // that fortran uses 4- or 8-bytes per logical and C++ 1-byte 2969 // that fortran uses 4- or 8-bytes per logical and C++ 1-byte
2970 // per bool, though this might be system dependent. As 2970 // per bool, though this might be system dependent. As
2971 // long as the HOWMNY arg is not "S", the logical array 2971 // long as the HOWMNY arg is not "S", the logical array
2972 // is just workspace for ARPACK, so use int type to 2972 // is just workspace for ARPACK, so use int type to
2973 // avoid problems. 2973 // avoid problems.
2974 Array<octave_idx_type> s (dim_vector (p, 1)); 2974 Array<octave_idx_type> s (dim_vector (p, 1));
2975 octave_idx_type *sel = s.fortran_vec (); 2975 octave_idx_type *sel = s.fortran_vec ();
2976 2976
2977 eig_vec.resize (n, k); 2977 eig_vec.resize (n, k);
2980 eig_val.resize (k+1); 2980 eig_val.resize (k+1);
2981 Complex *d = eig_val.fortran_vec (); 2981 Complex *d = eig_val.fortran_vec ();
2982 2982
2983 OCTAVE_LOCAL_BUFFER (Complex, workev, 2 * p); 2983 OCTAVE_LOCAL_BUFFER (Complex, workev, 2 * p);
2984 2984
2985 F77_FUNC (zneupd, ZNEUPD) 2985 F77_FUNC (zneupd, ZNEUPD)
2986 (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma, workev, 2986 (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma, workev,
2987 F77_CONST_CHAR_ARG2 (&bmat, 1), n, 2987 F77_CONST_CHAR_ARG2 (&bmat, 1), n,
2988 F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), 2988 F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2),
2989 k, tol, presid, p, v, n, iparam, ipntr, workd, workl, lwork, rwork, info2 2989 k, tol, presid, p, v, n, iparam, ipntr, workd, workl, lwork, rwork, info2
2990 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); 2990 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
2991 2991
2992 if (f77_exception_encountered) 2992 if (f77_exception_encountered)
2993 { 2993 {
2994 (*current_liboctave_error_handler) 2994 (*current_liboctave_error_handler)
2995 ("eigs: unrecoverable exception encountered in zneupd"); 2995 ("eigs: unrecoverable exception encountered in zneupd");
2996 return -1; 2996 return -1;
2997 } 2997 }
2998 2998
2999 if (info2 == 0) 2999 if (info2 == 0)
3033 eig_vec = ltsolve(b, permB, eig_vec); 3033 eig_vec = ltsolve(b, permB, eig_vec);
3034 } 3034 }
3035 } 3035 }
3036 else 3036 else
3037 { 3037 {
3038 (*current_liboctave_error_handler) 3038 (*current_liboctave_error_handler)
3039 ("eigs: error %d in zneupd", info2); 3039 ("eigs: error %d in zneupd", info2);
3040 return -1; 3040 return -1;
3041 } 3041 }
3042 3042
3043 return ip(4); 3043 return ip(4);
3044 } 3044 }
3045 3045
3046 template <class M> 3046 template <class M>
3047 octave_idx_type 3047 octave_idx_type
3048 EigsComplexNonSymmetricMatrixShift (const M& m, Complex sigma, 3048 EigsComplexNonSymmetricMatrixShift (const M& m, Complex sigma,
3049 octave_idx_type k, octave_idx_type p, 3049 octave_idx_type k, octave_idx_type p,
3050 octave_idx_type &info, 3050 octave_idx_type &info,
3051 ComplexMatrix &eig_vec, 3051 ComplexMatrix &eig_vec,
3052 ComplexColumnVector &eig_val, const M& _b, 3052 ComplexColumnVector &eig_val, const M& _b,
3053 ColumnVector &permB, 3053 ColumnVector &permB,
3054 ComplexColumnVector &cresid, 3054 ComplexColumnVector &cresid,
3055 std::ostream& os, double tol, bool rvec, 3055 std::ostream& os, double tol, bool rvec,
3056 bool cholB, int disp, int maxit) 3056 bool cholB, int disp, int maxit)
3057 { 3057 {
3058 M b(_b); 3058 M b(_b);
3059 octave_idx_type n = m.cols (); 3059 octave_idx_type n = m.cols ();
3060 octave_idx_type mode = 3; 3060 octave_idx_type mode = 3;
3066 (*current_liboctave_error_handler) ("eigs: A must be square"); 3066 (*current_liboctave_error_handler) ("eigs: A must be square");
3067 return -1; 3067 return -1;
3068 } 3068 }
3069 if (have_b && (m.rows() != b.rows() || m.rows() != b.cols())) 3069 if (have_b && (m.rows() != b.rows() || m.rows() != b.cols()))
3070 { 3070 {
3071 (*current_liboctave_error_handler) 3071 (*current_liboctave_error_handler)
3072 ("eigs: B must be square and the same size as A"); 3072 ("eigs: B must be square and the same size as A");
3073 return -1; 3073 return -1;
3074 } 3074 }
3075 3075
3076 // FIXME: The "SM" type for mode 1 seems unstable though faster!! 3076 // FIXME: The "SM" type for mode 1 seems unstable though faster!!
3102 { 3102 {
3103 p = k * 2 + 1; 3103 p = k * 2 + 1;
3104 3104
3105 if (p < 20) 3105 if (p < 20)
3106 p = 20; 3106 p = 20;
3107 3107
3108 if (p > n - 1) 3108 if (p > n - 1)
3109 p = n - 1 ; 3109 p = n - 1 ;
3110 } 3110 }
3111 3111
3112 if (k <= 0 || k >= n - 1) 3112 if (k <= 0 || k >= n - 1)
3113 { 3113 {
3114 (*current_liboctave_error_handler) 3114 (*current_liboctave_error_handler)
3115 ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1).\n" 3115 ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1).\n"
3116 " Use 'eig(full(A))' instead"); 3116 " Use 'eig(full(A))' instead");
3117 return -1; 3117 return -1;
3118 } 3118 }
3119 3119
3120 if (p <= k || p >= n) 3120 if (p <= k || p >= n)
3121 { 3121 {
3122 (*current_liboctave_error_handler) 3122 (*current_liboctave_error_handler)
3123 ("eigs: opts.p must be greater than k and less than n"); 3123 ("eigs: opts.p must be greater than k and less than n");
3124 return -1; 3124 return -1;
3125 } 3125 }
3126 3126
3127 if (have_b && cholB && permB.length() != 0) 3127 if (have_b && cholB && permB.length() != 0)
3128 { 3128 {
3129 // Check that we really have a permutation vector 3129 // Check that we really have a permutation vector
3130 if (permB.length() != n) 3130 if (permB.length() != n)
3131 { 3131 {
3132 (*current_liboctave_error_handler) ("eigs: permB vector invalid"); 3132 (*current_liboctave_error_handler) ("eigs: permB vector invalid");
3135 else 3135 else
3136 { 3136 {
3137 Array<bool> checked (dim_vector (n, 1), false); 3137 Array<bool> checked (dim_vector (n, 1), false);
3138 for (octave_idx_type i = 0; i < n; i++) 3138 for (octave_idx_type i = 0; i < n; i++)
3139 { 3139 {
3140 octave_idx_type bidx = 3140 octave_idx_type bidx =
3141 static_cast<octave_idx_type> (permB(i)); 3141 static_cast<octave_idx_type> (permB(i));
3142 if (checked(bidx) || bidx < 0 || 3142 if (checked(bidx) || bidx < 0 ||
3143 bidx >= n || D_NINT (bidx) != bidx) 3143 bidx >= n || D_NINT (bidx) != bidx)
3144 { 3144 {
3145 (*current_liboctave_error_handler) 3145 (*current_liboctave_error_handler)
3146 ("eigs: permB vector invalid"); 3146 ("eigs: permB vector invalid");
3147 return -1; 3147 return -1;
3148 } 3148 }
3149 } 3149 }
3150 } 3150 }
3182 3182
3183 if (! LuAminusSigmaB(m, b, cholB, permB, sigma, L, U, P, Q)) 3183 if (! LuAminusSigmaB(m, b, cholB, permB, sigma, L, U, P, Q))
3184 return -1; 3184 return -1;
3185 3185
3186 octave_idx_type lwork = p * (3 * p + 5); 3186 octave_idx_type lwork = p * (3 * p + 5);
3187 3187
3188 OCTAVE_LOCAL_BUFFER (Complex, v, n * p); 3188 OCTAVE_LOCAL_BUFFER (Complex, v, n * p);
3189 OCTAVE_LOCAL_BUFFER (Complex, workl, lwork); 3189 OCTAVE_LOCAL_BUFFER (Complex, workl, lwork);
3190 OCTAVE_LOCAL_BUFFER (Complex, workd, 3 * n); 3190 OCTAVE_LOCAL_BUFFER (Complex, workd, 3 * n);
3191 OCTAVE_LOCAL_BUFFER (double, rwork, p); 3191 OCTAVE_LOCAL_BUFFER (double, rwork, p);
3192 Complex *presid = cresid.fortran_vec (); 3192 Complex *presid = cresid.fortran_vec ();
3193 3193
3194 do 3194 do
3195 { 3195 {
3196 F77_FUNC (znaupd, ZNAUPD) 3196 F77_FUNC (znaupd, ZNAUPD)
3197 (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, 3197 (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
3198 F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), 3198 F77_CONST_CHAR_ARG2 ((typ.c_str()), 2),
3199 k, tol, presid, p, v, n, iparam, 3199 k, tol, presid, p, v, n, iparam,
3200 ipntr, workd, workl, lwork, rwork, info 3200 ipntr, workd, workl, lwork, rwork, info
3201 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); 3201 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
3202 3202
3203 if (f77_exception_encountered) 3203 if (f77_exception_encountered)
3204 { 3204 {
3205 (*current_liboctave_error_handler) 3205 (*current_liboctave_error_handler)
3206 ("eigs: unrecoverable exception encountered in znaupd"); 3206 ("eigs: unrecoverable exception encountered in znaupd");
3207 return -1; 3207 return -1;
3208 } 3208 }
3209 3209
3210 if (disp > 0 && !xisnan(workl[iptr(5)-1])) 3210 if (disp > 0 && !xisnan(workl[iptr(5)-1]))
3211 { 3211 {
3212 if (iter++) 3212 if (iter++)
3213 { 3213 {
3214 os << "Iteration " << iter - 1 << 3214 os << "Iteration " << iter - 1 <<
3215 ": a few Ritz values of the " << p << "-by-" << 3215 ": a few Ritz values of the " << p << "-by-" <<
3216 p << " matrix\n"; 3216 p << " matrix\n";
3217 for (int i = 0 ; i < k; i++) 3217 for (int i = 0 ; i < k; i++)
3218 os << " " << workl[iptr(5)+i-1] << "\n"; 3218 os << " " << workl[iptr(5)+i-1] << "\n";
3219 } 3219 }
3220 3220
3221 // This is a kludge, as ARPACK doesn't give its 3221 // This is a kludge, as ARPACK doesn't give its
3222 // iteration pointer. But as workl[iptr(5)-1] is 3222 // iteration pointer. But as workl[iptr(5)-1] is
3223 // an output value updated at each iteration, setting 3223 // an output value updated at each iteration, setting
3224 // a value in this array to NaN and testing for it 3224 // a value in this array to NaN and testing for it
3225 // is a way of obtaining the iteration counter. 3225 // is a way of obtaining the iteration counter.
3226 if (ido != 99) 3226 if (ido != 99)
3227 workl[iptr(5)-1] = octave_NaN; 3227 workl[iptr(5)-1] = octave_NaN;
3228 } 3228 }
3229 3229
3230 if (ido == -1 || ido == 1 || ido == 2) 3230 if (ido == -1 || ido == 1 || ido == 2)
3231 { 3231 {
3232 if (have_b) 3232 if (have_b)
3239 3239
3240 ComplexMatrix tmp(n, 1); 3240 ComplexMatrix tmp(n, 1);
3241 3241
3242 for (octave_idx_type i = 0; i < n; i++) 3242 for (octave_idx_type i = 0; i < n; i++)
3243 tmp(i,0) = ctmp[P[i]]; 3243 tmp(i,0) = ctmp[P[i]];
3244 3244
3245 lusolve (L, U, tmp); 3245 lusolve (L, U, tmp);
3246 3246
3247 Complex *ip2 = workd+iptr(1)-1; 3247 Complex *ip2 = workd+iptr(1)-1;
3248 for (octave_idx_type i = 0; i < n; i++) 3248 for (octave_idx_type i = 0; i < n; i++)
3249 ip2[Q[i]] = tmp(i,0); 3249 ip2[Q[i]] = tmp(i,0);
3255 Complex *ip2 = workd+iptr(2)-1; 3255 Complex *ip2 = workd+iptr(2)-1;
3256 ComplexMatrix tmp(n, 1); 3256 ComplexMatrix tmp(n, 1);
3257 3257
3258 for (octave_idx_type i = 0; i < n; i++) 3258 for (octave_idx_type i = 0; i < n; i++)
3259 tmp(i,0) = ip2[P[i]]; 3259 tmp(i,0) = ip2[P[i]];
3260 3260
3261 lusolve (L, U, tmp); 3261 lusolve (L, U, tmp);
3262 3262
3263 ip2 = workd+iptr(1)-1; 3263 ip2 = workd+iptr(1)-1;
3264 for (octave_idx_type i = 0; i < n; i++) 3264 for (octave_idx_type i = 0; i < n; i++)
3265 ip2[Q[i]] = tmp(i,0); 3265 ip2[Q[i]] = tmp(i,0);
3278 Complex *ip2 = workd+iptr(0)-1; 3278 Complex *ip2 = workd+iptr(0)-1;
3279 ComplexMatrix tmp(n, 1); 3279 ComplexMatrix tmp(n, 1);
3280 3280
3281 for (octave_idx_type i = 0; i < n; i++) 3281 for (octave_idx_type i = 0; i < n; i++)
3282 tmp(i,0) = ip2[P[i]]; 3282 tmp(i,0) = ip2[P[i]];
3283 3283
3284 lusolve (L, U, tmp); 3284 lusolve (L, U, tmp);
3285 3285
3286 ip2 = workd+iptr(1)-1; 3286 ip2 = workd+iptr(1)-1;
3287 for (octave_idx_type i = 0; i < n; i++) 3287 for (octave_idx_type i = 0; i < n; i++)
3288 ip2[Q[i]] = tmp(i,0); 3288 ip2[Q[i]] = tmp(i,0);
3291 } 3291 }
3292 else 3292 else
3293 { 3293 {
3294 if (info < 0) 3294 if (info < 0)
3295 { 3295 {
3296 (*current_liboctave_error_handler) 3296 (*current_liboctave_error_handler)
3297 ("eigs: error %d in dsaupd", info); 3297 ("eigs: error %d in dsaupd", info);
3298 return -1; 3298 return -1;
3299 } 3299 }
3300 break; 3300 break;
3301 } 3301 }
3302 } 3302 }
3303 while (1); 3303 while (1);
3304 3304
3305 octave_idx_type info2; 3305 octave_idx_type info2;
3306 3306
3307 // We have a problem in that the size of the C++ bool 3307 // We have a problem in that the size of the C++ bool
3308 // type relative to the fortran logical type. It appears 3308 // type relative to the fortran logical type. It appears
3309 // that fortran uses 4- or 8-bytes per logical and C++ 1-byte 3309 // that fortran uses 4- or 8-bytes per logical and C++ 1-byte
3310 // per bool, though this might be system dependent. As 3310 // per bool, though this might be system dependent. As
3311 // long as the HOWMNY arg is not "S", the logical array 3311 // long as the HOWMNY arg is not "S", the logical array
3312 // is just workspace for ARPACK, so use int type to 3312 // is just workspace for ARPACK, so use int type to
3313 // avoid problems. 3313 // avoid problems.
3314 Array<octave_idx_type> s (dim_vector (p, 1)); 3314 Array<octave_idx_type> s (dim_vector (p, 1));
3315 octave_idx_type *sel = s.fortran_vec (); 3315 octave_idx_type *sel = s.fortran_vec ();
3316 3316
3317 eig_vec.resize (n, k); 3317 eig_vec.resize (n, k);
3320 eig_val.resize (k+1); 3320 eig_val.resize (k+1);
3321 Complex *d = eig_val.fortran_vec (); 3321 Complex *d = eig_val.fortran_vec ();
3322 3322
3323 OCTAVE_LOCAL_BUFFER (Complex, workev, 2 * p); 3323 OCTAVE_LOCAL_BUFFER (Complex, workev, 2 * p);
3324 3324
3325 F77_FUNC (zneupd, ZNEUPD) 3325 F77_FUNC (zneupd, ZNEUPD)
3326 (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma, workev, 3326 (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma, workev,
3327 F77_CONST_CHAR_ARG2 (&bmat, 1), n, 3327 F77_CONST_CHAR_ARG2 (&bmat, 1), n,
3328 F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), 3328 F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2),
3329 k, tol, presid, p, v, n, iparam, ipntr, workd, workl, lwork, rwork, info2 3329 k, tol, presid, p, v, n, iparam, ipntr, workd, workl, lwork, rwork, info2
3330 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); 3330 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
3331 3331
3332 if (f77_exception_encountered) 3332 if (f77_exception_encountered)
3333 { 3333 {
3334 (*current_liboctave_error_handler) 3334 (*current_liboctave_error_handler)
3335 ("eigs: unrecoverable exception encountered in zneupd"); 3335 ("eigs: unrecoverable exception encountered in zneupd");
3336 return -1; 3336 return -1;
3337 } 3337 }
3338 3338
3339 if (info2 == 0) 3339 if (info2 == 0)
3370 } 3370 }
3371 } 3371 }
3372 } 3372 }
3373 else 3373 else
3374 { 3374 {
3375 (*current_liboctave_error_handler) 3375 (*current_liboctave_error_handler)
3376 ("eigs: error %d in zneupd", info2); 3376 ("eigs: error %d in zneupd", info2);
3377 return -1; 3377 return -1;
3378 } 3378 }
3379 3379
3380 return ip(4); 3380 return ip(4);
3381 } 3381 }
3382 3382
3383 octave_idx_type 3383 octave_idx_type
3384 EigsComplexNonSymmetricFunc (EigsComplexFunc fun, octave_idx_type n, 3384 EigsComplexNonSymmetricFunc (EigsComplexFunc fun, octave_idx_type n,
3385 const std::string &_typ, Complex sigma, 3385 const std::string &_typ, Complex sigma,
3386 octave_idx_type k, octave_idx_type p, 3386 octave_idx_type k, octave_idx_type p,
3387 octave_idx_type &info, ComplexMatrix &eig_vec, 3387 octave_idx_type &info, ComplexMatrix &eig_vec,
3388 ComplexColumnVector &eig_val, 3388 ComplexColumnVector &eig_val,
3389 ComplexColumnVector &cresid, std::ostream& os, 3389 ComplexColumnVector &cresid, std::ostream& os,
3390 double tol, bool rvec, bool /* cholB */, 3390 double tol, bool rvec, bool /* cholB */,
3391 int disp, int maxit) 3391 int disp, int maxit)
3392 { 3392 {
3393 std::string typ (_typ); 3393 std::string typ (_typ);
3394 bool have_sigma = (std::abs(sigma) ? true : false); 3394 bool have_sigma = (std::abs(sigma) ? true : false);
3419 { 3419 {
3420 p = k * 2 + 1; 3420 p = k * 2 + 1;
3421 3421
3422 if (p < 20) 3422 if (p < 20)
3423 p = 20; 3423 p = 20;
3424 3424
3425 if (p > n - 1) 3425 if (p > n - 1)
3426 p = n - 1 ; 3426 p = n - 1 ;
3427 } 3427 }
3428 3428
3429 if (k <= 0 || k >= n - 1) 3429 if (k <= 0 || k >= n - 1)
3441 return -1; 3441 return -1;
3442 } 3442 }
3443 3443
3444 if (! have_sigma) 3444 if (! have_sigma)
3445 { 3445 {
3446 if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" && 3446 if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" &&
3447 typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" && 3447 typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" &&
3448 typ != "SI") 3448 typ != "SI")
3449 (*current_liboctave_error_handler) 3449 (*current_liboctave_error_handler)
3450 ("eigs: unrecognized sigma value"); 3450 ("eigs: unrecognized sigma value");
3451 3451
3452 if (typ == "LA" || typ == "SA" || typ == "BE") 3452 if (typ == "LA" || typ == "SA" || typ == "BE")
3453 { 3453 {
3454 (*current_liboctave_error_handler) 3454 (*current_liboctave_error_handler)
3455 ("eigs: invalid sigma value for complex problem"); 3455 ("eigs: invalid sigma value for complex problem");
3456 return -1; 3456 return -1;
3457 } 3457 }
3458 3458
3459 if (typ == "SM") 3459 if (typ == "SM")
3484 ip(7) = 0; 3484 ip(7) = 0;
3485 ip(8) = 0; 3485 ip(8) = 0;
3486 ip(9) = 0; 3486 ip(9) = 0;
3487 ip(10) = 0; 3487 ip(10) = 0;
3488 // ip(7) to ip(10) return values 3488 // ip(7) to ip(10) return values
3489 3489
3490 Array<octave_idx_type> iptr (dim_vector (14, 1)); 3490 Array<octave_idx_type> iptr (dim_vector (14, 1));
3491 octave_idx_type *ipntr = iptr.fortran_vec (); 3491 octave_idx_type *ipntr = iptr.fortran_vec ();
3492 3492
3493 octave_idx_type ido = 0; 3493 octave_idx_type ido = 0;
3494 int iter = 0; 3494 int iter = 0;
3495 octave_idx_type lwork = p * (3 * p + 5); 3495 octave_idx_type lwork = p * (3 * p + 5);
3496 3496
3497 OCTAVE_LOCAL_BUFFER (Complex, v, n * p); 3497 OCTAVE_LOCAL_BUFFER (Complex, v, n * p);
3498 OCTAVE_LOCAL_BUFFER (Complex, workl, lwork); 3498 OCTAVE_LOCAL_BUFFER (Complex, workl, lwork);
3499 OCTAVE_LOCAL_BUFFER (Complex, workd, 3 * n); 3499 OCTAVE_LOCAL_BUFFER (Complex, workd, 3 * n);
3500 OCTAVE_LOCAL_BUFFER (double, rwork, p); 3500 OCTAVE_LOCAL_BUFFER (double, rwork, p);
3501 Complex *presid = cresid.fortran_vec (); 3501 Complex *presid = cresid.fortran_vec ();
3502 3502
3503 do 3503 do
3504 { 3504 {
3505 F77_FUNC (znaupd, ZNAUPD) 3505 F77_FUNC (znaupd, ZNAUPD)
3506 (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, 3506 (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
3507 F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), 3507 F77_CONST_CHAR_ARG2 ((typ.c_str()), 2),
3508 k, tol, presid, p, v, n, iparam, 3508 k, tol, presid, p, v, n, iparam,
3509 ipntr, workd, workl, lwork, rwork, info 3509 ipntr, workd, workl, lwork, rwork, info
3510 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); 3510 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
3511 3511
3512 if (f77_exception_encountered) 3512 if (f77_exception_encountered)
3513 { 3513 {
3514 (*current_liboctave_error_handler) 3514 (*current_liboctave_error_handler)
3515 ("eigs: unrecoverable exception encountered in znaupd"); 3515 ("eigs: unrecoverable exception encountered in znaupd");
3516 return -1; 3516 return -1;
3517 } 3517 }
3518 3518
3519 if (disp > 0 && !xisnan(workl[iptr(5)-1])) 3519 if (disp > 0 && !xisnan(workl[iptr(5)-1]))
3520 { 3520 {
3521 if (iter++) 3521 if (iter++)
3522 { 3522 {
3523 os << "Iteration " << iter - 1 << 3523 os << "Iteration " << iter - 1 <<
3524 ": a few Ritz values of the " << p << "-by-" << 3524 ": a few Ritz values of the " << p << "-by-" <<
3525 p << " matrix\n"; 3525 p << " matrix\n";
3526 for (int i = 0 ; i < k; i++) 3526 for (int i = 0 ; i < k; i++)
3527 os << " " << workl[iptr(5)+i-1] << "\n"; 3527 os << " " << workl[iptr(5)+i-1] << "\n";
3528 } 3528 }
3529 3529
3530 // This is a kludge, as ARPACK doesn't give its 3530 // This is a kludge, as ARPACK doesn't give its
3531 // iteration pointer. But as workl[iptr(5)-1] is 3531 // iteration pointer. But as workl[iptr(5)-1] is
3532 // an output value updated at each iteration, setting 3532 // an output value updated at each iteration, setting
3533 // a value in this array to NaN and testing for it 3533 // a value in this array to NaN and testing for it
3534 // is a way of obtaining the iteration counter. 3534 // is a way of obtaining the iteration counter.
3535 if (ido != 99) 3535 if (ido != 99)
3536 workl[iptr(5)-1] = octave_NaN; 3536 workl[iptr(5)-1] = octave_NaN;
3537 } 3537 }
3538 3538
3539 if (ido == -1 || ido == 1 || ido == 2) 3539 if (ido == -1 || ido == 1 || ido == 2)
3540 { 3540 {
3541 Complex *ip2 = workd + iptr(0) - 1; 3541 Complex *ip2 = workd + iptr(0) - 1;
3555 } 3555 }
3556 else 3556 else
3557 { 3557 {
3558 if (info < 0) 3558 if (info < 0)
3559 { 3559 {
3560 (*current_liboctave_error_handler) 3560 (*current_liboctave_error_handler)
3561 ("eigs: error %d in dsaupd", info); 3561 ("eigs: error %d in dsaupd", info);
3562 return -1; 3562 return -1;
3563 } 3563 }
3564 break; 3564 break;
3565 } 3565 }
3566 } 3566 }
3567 while (1); 3567 while (1);
3568 3568
3569 octave_idx_type info2; 3569 octave_idx_type info2;
3570 3570
3571 // We have a problem in that the size of the C++ bool 3571 // We have a problem in that the size of the C++ bool
3572 // type relative to the fortran logical type. It appears 3572 // type relative to the fortran logical type. It appears
3573 // that fortran uses 4- or 8-bytes per logical and C++ 1-byte 3573 // that fortran uses 4- or 8-bytes per logical and C++ 1-byte
3574 // per bool, though this might be system dependent. As 3574 // per bool, though this might be system dependent. As
3575 // long as the HOWMNY arg is not "S", the logical array 3575 // long as the HOWMNY arg is not "S", the logical array
3576 // is just workspace for ARPACK, so use int type to 3576 // is just workspace for ARPACK, so use int type to
3577 // avoid problems. 3577 // avoid problems.
3578 Array<octave_idx_type> s (dim_vector (p, 1)); 3578 Array<octave_idx_type> s (dim_vector (p, 1));
3579 octave_idx_type *sel = s.fortran_vec (); 3579 octave_idx_type *sel = s.fortran_vec ();
3580 3580
3581 eig_vec.resize (n, k); 3581 eig_vec.resize (n, k);
3584 eig_val.resize (k+1); 3584 eig_val.resize (k+1);
3585 Complex *d = eig_val.fortran_vec (); 3585 Complex *d = eig_val.fortran_vec ();
3586 3586
3587 OCTAVE_LOCAL_BUFFER (Complex, workev, 2 * p); 3587 OCTAVE_LOCAL_BUFFER (Complex, workev, 2 * p);
3588 3588
3589 F77_FUNC (zneupd, ZNEUPD) 3589 F77_FUNC (zneupd, ZNEUPD)
3590 (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma, workev, 3590 (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma, workev,
3591 F77_CONST_CHAR_ARG2 (&bmat, 1), n, 3591 F77_CONST_CHAR_ARG2 (&bmat, 1), n,
3592 F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), 3592 F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2),
3593 k, tol, presid, p, v, n, iparam, ipntr, workd, workl, lwork, rwork, info2 3593 k, tol, presid, p, v, n, iparam, ipntr, workd, workl, lwork, rwork, info2
3594 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); 3594 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
3595 3595
3596 if (f77_exception_encountered) 3596 if (f77_exception_encountered)
3597 { 3597 {
3598 (*current_liboctave_error_handler) 3598 (*current_liboctave_error_handler)
3599 ("eigs: unrecoverable exception encountered in zneupd"); 3599 ("eigs: unrecoverable exception encountered in zneupd");
3600 return -1; 3600 return -1;
3601 } 3601 }
3602 3602
3603 if (info2 == 0) 3603 if (info2 == 0)
3634 } 3634 }
3635 } 3635 }
3636 } 3636 }
3637 else 3637 else
3638 { 3638 {
3639 (*current_liboctave_error_handler) 3639 (*current_liboctave_error_handler)
3640 ("eigs: error %d in zneupd", info2); 3640 ("eigs: error %d in zneupd", info2);
3641 return -1; 3641 return -1;
3642 } 3642 }
3643 3643
3644 return ip(4); 3644 return ip(4);
3645 } 3645 }
3646 3646
3647 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL) 3647 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
3648 extern octave_idx_type 3648 extern octave_idx_type
3649 EigsRealSymmetricMatrix (const Matrix& m, const std::string typ, 3649 EigsRealSymmetricMatrix (const Matrix& m, const std::string typ,
3650 octave_idx_type k, octave_idx_type p, 3650 octave_idx_type k, octave_idx_type p,
3651 octave_idx_type &info, Matrix &eig_vec, 3651 octave_idx_type &info, Matrix &eig_vec,
3652 ColumnVector &eig_val, const Matrix& b, 3652 ColumnVector &eig_val, const Matrix& b,
3653 ColumnVector &permB, ColumnVector &resid, 3653 ColumnVector &permB, ColumnVector &resid,
3654 std::ostream &os, double tol = DBL_EPSILON, 3654 std::ostream &os, double tol = DBL_EPSILON,
3655 bool rvec = false, bool cholB = 0, int disp = 0, 3655 bool rvec = false, bool cholB = 0, int disp = 0,
3656 int maxit = 300); 3656 int maxit = 300);
3657 3657
3658 extern octave_idx_type 3658 extern octave_idx_type
3659 EigsRealSymmetricMatrix (const SparseMatrix& m, const std::string typ, 3659 EigsRealSymmetricMatrix (const SparseMatrix& m, const std::string typ,
3660 octave_idx_type k, octave_idx_type p, 3660 octave_idx_type k, octave_idx_type p,
3661 octave_idx_type &info, Matrix &eig_vec, 3661 octave_idx_type &info, Matrix &eig_vec,
3662 ColumnVector &eig_val, const SparseMatrix& b, 3662 ColumnVector &eig_val, const SparseMatrix& b,
3663 ColumnVector &permB, ColumnVector &resid, 3663 ColumnVector &permB, ColumnVector &resid,
3664 std::ostream& os, double tol = DBL_EPSILON, 3664 std::ostream& os, double tol = DBL_EPSILON,
3665 bool rvec = false, bool cholB = 0, int disp = 0, 3665 bool rvec = false, bool cholB = 0, int disp = 0,
3666 int maxit = 300); 3666 int maxit = 300);
3667 3667
3668 extern octave_idx_type 3668 extern octave_idx_type
3669 EigsRealSymmetricMatrixShift (const Matrix& m, double sigma, 3669 EigsRealSymmetricMatrixShift (const Matrix& m, double sigma,
3670 octave_idx_type k, octave_idx_type p, 3670 octave_idx_type k, octave_idx_type p,
3671 octave_idx_type &info, Matrix &eig_vec, 3671 octave_idx_type &info, Matrix &eig_vec,
3672 ColumnVector &eig_val, const Matrix& b, 3672 ColumnVector &eig_val, const Matrix& b,
3673 ColumnVector &permB, ColumnVector &resid, 3673 ColumnVector &permB, ColumnVector &resid,
3674 std::ostream &os, double tol = DBL_EPSILON, 3674 std::ostream &os, double tol = DBL_EPSILON,
3675 bool rvec = false, bool cholB = 0, int disp = 0, 3675 bool rvec = false, bool cholB = 0, int disp = 0,
3676 int maxit = 300); 3676 int maxit = 300);
3677 3677
3678 extern octave_idx_type 3678 extern octave_idx_type
3679 EigsRealSymmetricMatrixShift (const SparseMatrix& m, double sigma, 3679 EigsRealSymmetricMatrixShift (const SparseMatrix& m, double sigma,
3680 octave_idx_type k, octave_idx_type p, 3680 octave_idx_type k, octave_idx_type p,
3681 octave_idx_type &info, Matrix &eig_vec, 3681 octave_idx_type &info, Matrix &eig_vec,
3682 ColumnVector &eig_val, const SparseMatrix& b, 3682 ColumnVector &eig_val, const SparseMatrix& b,
3683 ColumnVector &permB, ColumnVector &resid, 3683 ColumnVector &permB, ColumnVector &resid,
3684 std::ostream &os, double tol = DBL_EPSILON, 3684 std::ostream &os, double tol = DBL_EPSILON,
3685 bool rvec = false, bool cholB = 0, int disp = 0, 3685 bool rvec = false, bool cholB = 0, int disp = 0,
3686 int maxit = 300); 3686 int maxit = 300);
3687 3687
3688 extern octave_idx_type 3688 extern octave_idx_type
3689 EigsRealSymmetricFunc (EigsFunc fun, octave_idx_type n, 3689 EigsRealSymmetricFunc (EigsFunc fun, octave_idx_type n,
3690 const std::string &typ, double sigma, 3690 const std::string &typ, double sigma,
3691 octave_idx_type k, octave_idx_type p, 3691 octave_idx_type k, octave_idx_type p,
3692 octave_idx_type &info, 3692 octave_idx_type &info,
3693 Matrix &eig_vec, ColumnVector &eig_val, 3693 Matrix &eig_vec, ColumnVector &eig_val,
3694 ColumnVector &resid, std::ostream &os, 3694 ColumnVector &resid, std::ostream &os,
3695 double tol = DBL_EPSILON, bool rvec = false, 3695 double tol = DBL_EPSILON, bool rvec = false,
3696 bool cholB = 0, int disp = 0, int maxit = 300); 3696 bool cholB = 0, int disp = 0, int maxit = 300);
3697 3697
3698 extern octave_idx_type 3698 extern octave_idx_type
3699 EigsRealNonSymmetricMatrix (const Matrix& m, const std::string typ, 3699 EigsRealNonSymmetricMatrix (const Matrix& m, const std::string typ,
3700 octave_idx_type k, octave_idx_type p, 3700 octave_idx_type k, octave_idx_type p,
3701 octave_idx_type &info, ComplexMatrix &eig_vec, 3701 octave_idx_type &info, ComplexMatrix &eig_vec,
3702 ComplexColumnVector &eig_val, const Matrix& b, 3702 ComplexColumnVector &eig_val, const Matrix& b,
3703 ColumnVector &permB, ColumnVector &resid, 3703 ColumnVector &permB, ColumnVector &resid,
3704 std::ostream &os, double tol = DBL_EPSILON, 3704 std::ostream &os, double tol = DBL_EPSILON,
3705 bool rvec = false, bool cholB = 0, int disp = 0, 3705 bool rvec = false, bool cholB = 0, int disp = 0,
3706 int maxit = 300); 3706 int maxit = 300);
3707 3707
3708 extern octave_idx_type 3708 extern octave_idx_type
3709 EigsRealNonSymmetricMatrix (const SparseMatrix& m, const std::string typ, 3709 EigsRealNonSymmetricMatrix (const SparseMatrix& m, const std::string typ,
3710 octave_idx_type k, octave_idx_type p, 3710 octave_idx_type k, octave_idx_type p,
3711 octave_idx_type &info, ComplexMatrix &eig_vec, 3711 octave_idx_type &info, ComplexMatrix &eig_vec,
3712 ComplexColumnVector &eig_val, 3712 ComplexColumnVector &eig_val,
3713 const SparseMatrix& b, 3713 const SparseMatrix& b,
3714 ColumnVector &permB, ColumnVector &resid, 3714 ColumnVector &permB, ColumnVector &resid,
3715 std::ostream &os, double tol = DBL_EPSILON, 3715 std::ostream &os, double tol = DBL_EPSILON,
3716 bool rvec = false, bool cholB = 0, int disp = 0, 3716 bool rvec = false, bool cholB = 0, int disp = 0,
3717 int maxit = 300); 3717 int maxit = 300);
3718 3718
3719 extern octave_idx_type 3719 extern octave_idx_type
3720 EigsRealNonSymmetricMatrixShift (const Matrix& m, double sigma, 3720 EigsRealNonSymmetricMatrixShift (const Matrix& m, double sigma,
3721 octave_idx_type k, octave_idx_type p, 3721 octave_idx_type k, octave_idx_type p,
3722 octave_idx_type &info, 3722 octave_idx_type &info,
3723 ComplexMatrix &eig_vec, 3723 ComplexMatrix &eig_vec,
3724 ComplexColumnVector &eig_val, const Matrix& b, 3724 ComplexColumnVector &eig_val, const Matrix& b,
3725 ColumnVector &permB, ColumnVector &resid, 3725 ColumnVector &permB, ColumnVector &resid,
3726 std::ostream &os, double tol = DBL_EPSILON, 3726 std::ostream &os, double tol = DBL_EPSILON,
3727 bool rvec = false, bool cholB = 0, 3727 bool rvec = false, bool cholB = 0,
3728 int disp = 0, int maxit = 300); 3728 int disp = 0, int maxit = 300);
3729 3729
3730 extern octave_idx_type 3730 extern octave_idx_type
3731 EigsRealNonSymmetricMatrixShift (const SparseMatrix& m, double sigma, 3731 EigsRealNonSymmetricMatrixShift (const SparseMatrix& m, double sigma,
3732 octave_idx_type k, octave_idx_type p, 3732 octave_idx_type k, octave_idx_type p,
3733 octave_idx_type &info, 3733 octave_idx_type &info,
3734 ComplexMatrix &eig_vec, 3734 ComplexMatrix &eig_vec,
3735 ComplexColumnVector &eig_val, 3735 ComplexColumnVector &eig_val,
3736 const SparseMatrix& b, 3736 const SparseMatrix& b,
3737 ColumnVector &permB, ColumnVector &resid, 3737 ColumnVector &permB, ColumnVector &resid,
3738 std::ostream &os, double tol = DBL_EPSILON, 3738 std::ostream &os, double tol = DBL_EPSILON,
3739 bool rvec = false, bool cholB = 0, 3739 bool rvec = false, bool cholB = 0,
3740 int disp = 0, int maxit = 300); 3740 int disp = 0, int maxit = 300);
3741 3741
3742 extern octave_idx_type 3742 extern octave_idx_type
3743 EigsRealNonSymmetricFunc (EigsFunc fun, octave_idx_type n, 3743 EigsRealNonSymmetricFunc (EigsFunc fun, octave_idx_type n,
3744 const std::string &_typ, double sigma, 3744 const std::string &_typ, double sigma,
3745 octave_idx_type k, octave_idx_type p, 3745 octave_idx_type k, octave_idx_type p,
3746 octave_idx_type &info, ComplexMatrix &eig_vec, 3746 octave_idx_type &info, ComplexMatrix &eig_vec,
3747 ComplexColumnVector &eig_val, 3747 ComplexColumnVector &eig_val,
3748 ColumnVector &resid, std::ostream& os, 3748 ColumnVector &resid, std::ostream& os,
3749 double tol = DBL_EPSILON, bool rvec = false, 3749 double tol = DBL_EPSILON, bool rvec = false,
3750 bool cholB = 0, int disp = 0, int maxit = 300); 3750 bool cholB = 0, int disp = 0, int maxit = 300);
3751 3751
3752 extern octave_idx_type 3752 extern octave_idx_type
3753 EigsComplexNonSymmetricMatrix (const ComplexMatrix& m, const std::string typ, 3753 EigsComplexNonSymmetricMatrix (const ComplexMatrix& m, const std::string typ,
3754 octave_idx_type k, octave_idx_type p, 3754 octave_idx_type k, octave_idx_type p,
3755 octave_idx_type &info, ComplexMatrix &eig_vec, 3755 octave_idx_type &info, ComplexMatrix &eig_vec,
3756 ComplexColumnVector &eig_val, 3756 ComplexColumnVector &eig_val,
3757 const ComplexMatrix& b, ColumnVector &permB, 3757 const ComplexMatrix& b, ColumnVector &permB,
3758 ComplexColumnVector &resid, 3758 ComplexColumnVector &resid,
3759 std::ostream &os, double tol = DBL_EPSILON, 3759 std::ostream &os, double tol = DBL_EPSILON,
3760 bool rvec = false, bool cholB = 0, int disp = 0, 3760 bool rvec = false, bool cholB = 0, int disp = 0,
3761 int maxit = 300); 3761 int maxit = 300);
3762 3762
3763 extern octave_idx_type 3763 extern octave_idx_type
3764 EigsComplexNonSymmetricMatrix (const SparseComplexMatrix& m, 3764 EigsComplexNonSymmetricMatrix (const SparseComplexMatrix& m,
3765 const std::string typ, octave_idx_type k, 3765 const std::string typ, octave_idx_type k,
3766 octave_idx_type p, octave_idx_type &info, 3766 octave_idx_type p, octave_idx_type &info,
3767 ComplexMatrix &eig_vec, 3767 ComplexMatrix &eig_vec,
3768 ComplexColumnVector &eig_val, 3768 ComplexColumnVector &eig_val,
3769 const SparseComplexMatrix& b, 3769 const SparseComplexMatrix& b,
3770 ColumnVector &permB, 3770 ColumnVector &permB,
3771 ComplexColumnVector &resid, 3771 ComplexColumnVector &resid,
3772 std::ostream &os, double tol = DBL_EPSILON, 3772 std::ostream &os, double tol = DBL_EPSILON,
3773 bool rvec = false, bool cholB = 0, int disp = 0, 3773 bool rvec = false, bool cholB = 0, int disp = 0,
3774 int maxit = 300); 3774 int maxit = 300);
3775 3775
3776 extern octave_idx_type 3776 extern octave_idx_type
3777 EigsComplexNonSymmetricMatrixShift (const ComplexMatrix& m, Complex sigma, 3777 EigsComplexNonSymmetricMatrixShift (const ComplexMatrix& m, Complex sigma,
3778 octave_idx_type k, octave_idx_type p, 3778 octave_idx_type k, octave_idx_type p,
3779 octave_idx_type &info, 3779 octave_idx_type &info,
3780 ComplexMatrix &eig_vec, 3780 ComplexMatrix &eig_vec,
3781 ComplexColumnVector &eig_val, 3781 ComplexColumnVector &eig_val,
3782 const ComplexMatrix& b, 3782 const ComplexMatrix& b,
3783 ColumnVector &permB, 3783 ColumnVector &permB,
3784 ComplexColumnVector &resid, 3784 ComplexColumnVector &resid,
3785 std::ostream &os, double tol = DBL_EPSILON, 3785 std::ostream &os, double tol = DBL_EPSILON,
3786 bool rvec = false, bool cholB = 0, 3786 bool rvec = false, bool cholB = 0,
3787 int disp = 0, int maxit = 300); 3787 int disp = 0, int maxit = 300);
3788 3788
3789 extern octave_idx_type 3789 extern octave_idx_type
3790 EigsComplexNonSymmetricMatrixShift (const SparseComplexMatrix& m, 3790 EigsComplexNonSymmetricMatrixShift (const SparseComplexMatrix& m,
3791 Complex sigma, 3791 Complex sigma,
3792 octave_idx_type k, octave_idx_type p, 3792 octave_idx_type k, octave_idx_type p,
3793 octave_idx_type &info, 3793 octave_idx_type &info,
3794 ComplexMatrix &eig_vec, 3794 ComplexMatrix &eig_vec,
3795 ComplexColumnVector &eig_val, 3795 ComplexColumnVector &eig_val,
3796 const SparseComplexMatrix& b, 3796 const SparseComplexMatrix& b,
3797 ColumnVector &permB, 3797 ColumnVector &permB,
3798 ComplexColumnVector &resid, 3798 ComplexColumnVector &resid,
3799 std::ostream &os, double tol = DBL_EPSILON, 3799 std::ostream &os, double tol = DBL_EPSILON,
3800 bool rvec = false, bool cholB = 0, 3800 bool rvec = false, bool cholB = 0,
3801 int disp = 0, int maxit = 300); 3801 int disp = 0, int maxit = 300);
3802 3802
3803 extern octave_idx_type 3803 extern octave_idx_type
3804 EigsComplexNonSymmetricFunc (EigsComplexFunc fun, octave_idx_type n, 3804 EigsComplexNonSymmetricFunc (EigsComplexFunc fun, octave_idx_type n,
3805 const std::string &_typ, Complex sigma, 3805 const std::string &_typ, Complex sigma,
3806 octave_idx_type k, octave_idx_type p, 3806 octave_idx_type k, octave_idx_type p,
3807 octave_idx_type &info, ComplexMatrix &eig_vec, 3807 octave_idx_type &info, ComplexMatrix &eig_vec,
3808 ComplexColumnVector &eig_val, 3808 ComplexColumnVector &eig_val,
3809 ComplexColumnVector &resid, std::ostream& os, 3809 ComplexColumnVector &resid, std::ostream& os,
3810 double tol = DBL_EPSILON, bool rvec = false, 3810 double tol = DBL_EPSILON, bool rvec = false,
3811 bool cholB = 0, int disp = 0, int maxit = 300); 3811 bool cholB = 0, int disp = 0, int maxit = 300);
3812 #endif 3812 #endif
3813 3813
3814 #ifndef _MSC_VER 3814 #ifndef _MSC_VER
3815 template static octave_idx_type 3815 template static octave_idx_type
3816 lusolve (const SparseMatrix&, const SparseMatrix&, Matrix&); 3816 lusolve (const SparseMatrix&, const SparseMatrix&, Matrix&);
3817 3817
3818 template static octave_idx_type 3818 template static octave_idx_type
3819 lusolve (const SparseComplexMatrix&, const SparseComplexMatrix&, 3819 lusolve (const SparseComplexMatrix&, const SparseComplexMatrix&,
3820 ComplexMatrix&); 3820 ComplexMatrix&);
3821 3821
3822 template static octave_idx_type 3822 template static octave_idx_type
3823 lusolve (const Matrix&, const Matrix&, Matrix&); 3823 lusolve (const Matrix&, const Matrix&, Matrix&);
3824 3824
3825 template static octave_idx_type 3825 template static octave_idx_type
3826 lusolve (const ComplexMatrix&, const ComplexMatrix&, ComplexMatrix&); 3826 lusolve (const ComplexMatrix&, const ComplexMatrix&, ComplexMatrix&);
3827 3827
3828 template static ComplexMatrix 3828 template static ComplexMatrix
3829 ltsolve (const SparseComplexMatrix&, const ColumnVector&, 3829 ltsolve (const SparseComplexMatrix&, const ColumnVector&,
3830 const ComplexMatrix&); 3830 const ComplexMatrix&);
3831 3831
3832 template static Matrix 3832 template static Matrix
3833 ltsolve (const SparseMatrix&, const ColumnVector&, const Matrix&); 3833 ltsolve (const SparseMatrix&, const ColumnVector&, const Matrix&);
3834 3834