GNU Octave: liboctave/eigs-base.cc Source File

Go to the documentation of this file.
00001 /*
00002 
00003 Copyright (C) 2005-2012 David Bateman
00004 
00005 This file is part of Octave.
00006 
00007 Octave is free software; you can redistribute it and/or modify it
00008 under the terms of the GNU General Public License as published by the
00009 Free Software Foundation; either version 3 of the License, or (at your
00010 option) any later version.
00011 
00012 Octave is distributed in the hope that it will be useful, but WITHOUT
00013 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00014 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00015 for more details.
00016 
00017 You should have received a copy of the GNU General Public License
00018 along with Octave; see the file COPYING.  If not, see
00019 <http://www.gnu.org/licenses/>.
00020 
00021 */
00022 
00023 #ifdef HAVE_CONFIG_H
00024 #include <config.h>
00025 #endif
00026 
00027 #include <cfloat>
00028 #include <cmath>
00029 #include <vector>
00030 #include <iostream>
00031 
00032 #include "f77-fcn.h"
00033 #include "quit.h"
00034 #include "SparsedbleLU.h"
00035 #include "SparseCmplxLU.h"
00036 #include "dSparse.h"
00037 #include "CSparse.h"
00038 #include "MatrixType.h"
00039 #include "SparsedbleCHOL.h"
00040 #include "SparseCmplxCHOL.h"
00041 #include "oct-rand.h"
00042 #include "dbleCHOL.h"
00043 #include "CmplxCHOL.h"
00044 #include "dbleLU.h"
00045 #include "CmplxLU.h"
00046 
00047 #ifdef HAVE_ARPACK
00048 typedef ColumnVector (*EigsFunc) (const ColumnVector &x, int &eigs_error);
00049 typedef ComplexColumnVector (*EigsComplexFunc)
00050   (const ComplexColumnVector &x, int &eigs_error);
00051 
00052 // Arpack and blas fortran functions we call.
00053 extern "C"
00054 {
00055   F77_RET_T
00056   F77_FUNC (dsaupd, DSAUPD) (octave_idx_type&,
00057                              F77_CONST_CHAR_ARG_DECL,
00058                              const octave_idx_type&,
00059                              F77_CONST_CHAR_ARG_DECL,
00060                              const octave_idx_type&, const double&,
00061                              double*, const octave_idx_type&, double*,
00062                              const octave_idx_type&, octave_idx_type*,
00063                              octave_idx_type*, double*, double*,
00064                              const octave_idx_type&, octave_idx_type&
00065                              F77_CHAR_ARG_LEN_DECL
00066                              F77_CHAR_ARG_LEN_DECL);
00067 
00068   F77_RET_T
00069   F77_FUNC (dseupd, DSEUPD) (const octave_idx_type&,
00070                              F77_CONST_CHAR_ARG_DECL,
00071                              octave_idx_type*, double*, double*,
00072                              const octave_idx_type&, const double&,
00073                              F77_CONST_CHAR_ARG_DECL,
00074                              const octave_idx_type&,
00075                              F77_CONST_CHAR_ARG_DECL,
00076                              const octave_idx_type&, const double&, double*,
00077                              const octave_idx_type&, double*,
00078                              const octave_idx_type&, octave_idx_type*,
00079                              octave_idx_type*, double*, double*,
00080                              const octave_idx_type&, octave_idx_type&
00081                              F77_CHAR_ARG_LEN_DECL
00082                              F77_CHAR_ARG_LEN_DECL
00083                              F77_CHAR_ARG_LEN_DECL);
00084 
00085   F77_RET_T
00086   F77_FUNC (dnaupd, DNAUPD) (octave_idx_type&,
00087                              F77_CONST_CHAR_ARG_DECL,
00088                              const octave_idx_type&,
00089                              F77_CONST_CHAR_ARG_DECL,
00090                              octave_idx_type&, const double&,
00091                              double*, const octave_idx_type&, double*,
00092                              const octave_idx_type&, octave_idx_type*,
00093                              octave_idx_type*, double*, double*,
00094                              const octave_idx_type&, octave_idx_type&
00095                              F77_CHAR_ARG_LEN_DECL
00096                              F77_CHAR_ARG_LEN_DECL);
00097 
00098   F77_RET_T
00099   F77_FUNC (dneupd, DNEUPD) (const octave_idx_type&,
00100                              F77_CONST_CHAR_ARG_DECL,
00101                              octave_idx_type*, double*, double*,
00102                              double*, const octave_idx_type&, const double&,
00103                              const double&, double*,
00104                              F77_CONST_CHAR_ARG_DECL,
00105                              const octave_idx_type&,
00106                              F77_CONST_CHAR_ARG_DECL,
00107                              octave_idx_type&, const double&, double*,
00108                              const octave_idx_type&, double*,
00109                              const octave_idx_type&, octave_idx_type*,
00110                              octave_idx_type*, double*, double*,
00111                              const octave_idx_type&, octave_idx_type&
00112                              F77_CHAR_ARG_LEN_DECL
00113                              F77_CHAR_ARG_LEN_DECL
00114                              F77_CHAR_ARG_LEN_DECL);
00115 
00116   F77_RET_T
00117   F77_FUNC (znaupd, ZNAUPD) (octave_idx_type&,
00118                              F77_CONST_CHAR_ARG_DECL,
00119                              const octave_idx_type&,
00120                              F77_CONST_CHAR_ARG_DECL,
00121                              const octave_idx_type&, const double&,
00122                              Complex*, const octave_idx_type&, Complex*,
00123                              const octave_idx_type&, octave_idx_type*,
00124                              octave_idx_type*, Complex*, Complex*,
00125                              const octave_idx_type&, double *, octave_idx_type&
00126                              F77_CHAR_ARG_LEN_DECL
00127                              F77_CHAR_ARG_LEN_DECL);
00128 
00129   F77_RET_T
00130   F77_FUNC (zneupd, ZNEUPD) (const octave_idx_type&,
00131                              F77_CONST_CHAR_ARG_DECL,
00132                              octave_idx_type*, Complex*, Complex*,
00133                              const octave_idx_type&, const Complex&, Complex*,
00134                              F77_CONST_CHAR_ARG_DECL,
00135                              const octave_idx_type&,
00136                              F77_CONST_CHAR_ARG_DECL,
00137                              const octave_idx_type&, const double&,
00138                              Complex*, const octave_idx_type&, Complex*,
00139                              const octave_idx_type&, octave_idx_type*,
00140                              octave_idx_type*, Complex*, Complex*,
00141                              const octave_idx_type&, double *, octave_idx_type&
00142                              F77_CHAR_ARG_LEN_DECL
00143                              F77_CHAR_ARG_LEN_DECL
00144                              F77_CHAR_ARG_LEN_DECL);
00145 
00146   F77_RET_T
00147   F77_FUNC (dgemv, DGEMV) (F77_CONST_CHAR_ARG_DECL,
00148                            const octave_idx_type&, const octave_idx_type&,
00149                            const double&, const double*,
00150                            const octave_idx_type&, const double*,
00151                            const octave_idx_type&, const double&, double*,
00152                            const octave_idx_type&
00153                            F77_CHAR_ARG_LEN_DECL);
00154 
00155 
00156   F77_RET_T
00157   F77_FUNC (zgemv, ZGEMV) (F77_CONST_CHAR_ARG_DECL,
00158                            const octave_idx_type&, const octave_idx_type&, const Complex&,
00159                            const Complex*, const octave_idx_type&, const Complex*,
00160                            const octave_idx_type&, const Complex&, Complex*, const octave_idx_type&
00161                            F77_CHAR_ARG_LEN_DECL);
00162 
00163 }
00164 
00165 
00166 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
00167 static octave_idx_type
00168 lusolve (const SparseMatrix&, const SparseMatrix&, Matrix&);
00169 
00170 static octave_idx_type
00171 lusolve (const SparseComplexMatrix&, const SparseComplexMatrix&,
00172          ComplexMatrix&);
00173 
00174 static octave_idx_type
00175 lusolve (const Matrix&, const Matrix&, Matrix&);
00176 
00177 static octave_idx_type
00178 lusolve (const ComplexMatrix&, const ComplexMatrix&, ComplexMatrix&);
00179 
00180 static ComplexMatrix
00181 ltsolve (const SparseComplexMatrix&, const ColumnVector&,
00182                 const ComplexMatrix&);
00183 
00184 static Matrix
00185 ltsolve (const SparseMatrix&, const ColumnVector&, const Matrix&,);
00186 
00187 static ComplexMatrix
00188 ltsolve (const ComplexMatrix&, const ColumnVector&, const ComplexMatrix&);
00189 
00190 static Matrix
00191 ltsolve (const Matrix&, const ColumnVector&, const Matrix&,);
00192 
00193 static ComplexMatrix
00194 utsolve (const SparseComplexMatrix&, const ColumnVector&, const ComplexMatrix&);
00195 
00196 static Matrix
00197 utsolve (const SparseMatrix&, const ColumnVector&, const Matrix&);
00198 
00199 static ComplexMatrix
00200 utsolve (const ComplexMatrix&, const ColumnVector&, const ComplexMatrix&);
00201 
00202 static Matrix
00203 utsolve (const Matrix&, const ColumnVector&, const Matrix&);
00204 
00205 #endif
00206 
00207 template <class M, class SM>
00208 static octave_idx_type
00209 lusolve (const SM& L, const SM& U, M& m)
00210 {
00211   octave_idx_type err = 0;
00212   double rcond;
00213   MatrixType utyp (MatrixType::Upper);
00214 
00215   // Sparse L is lower triangular, Dense L is permuted lower triangular!!!
00216   m = L.solve (m, err, rcond, 0);
00217   if (err)
00218     return err;
00219 
00220   m = U.solve (utyp, m, err, rcond, 0);
00221 
00222   return err;
00223 }
00224 
00225 template <class SM, class M>
00226 static M
00227 ltsolve (const SM& L, const ColumnVector& Q, const M& m)
00228 {
00229   octave_idx_type n = L.cols();
00230   octave_idx_type b_nc = m.cols();
00231   octave_idx_type err = 0;
00232   double rcond;
00233   MatrixType ltyp (MatrixType::Lower);
00234   M tmp = L.solve (ltyp, m, err, rcond, 0);
00235   M retval;
00236   const double* qv = Q.fortran_vec();
00237 
00238   if (!err)
00239     {
00240       retval.resize (n, b_nc);
00241       for (octave_idx_type j = 0; j < b_nc; j++)
00242         {
00243           for (octave_idx_type i = 0; i < n; i++)
00244             retval.elem(static_cast<octave_idx_type>(qv[i]), j)  =
00245               tmp.elem(i,j);
00246         }
00247     }
00248 
00249   return retval;
00250 }
00251 
00252 template <class SM, class M>
00253 static M
00254 utsolve (const SM& U, const ColumnVector& Q, const M& m)
00255 {
00256   octave_idx_type n = U.cols();
00257   octave_idx_type b_nc = m.cols();
00258   octave_idx_type err = 0;
00259   double rcond;
00260   MatrixType utyp (MatrixType::Upper);
00261 
00262   M retval (n, b_nc);
00263   const double* qv = Q.fortran_vec();
00264   for (octave_idx_type j = 0; j < b_nc; j++)
00265     {
00266       for (octave_idx_type i = 0; i < n; i++)
00267         retval.elem(i,j) = m.elem(static_cast<octave_idx_type>(qv[i]), j);
00268     }
00269   return U.solve (utyp, retval, err, rcond, 0);
00270 }
00271 
00272 static bool
00273 vector_product (const SparseMatrix& m, const double* x, double* y)
00274 {
00275   octave_idx_type nc = m.cols ();
00276 
00277   for (octave_idx_type j = 0; j < nc; j++)
00278     y[j] = 0.;
00279 
00280   for (octave_idx_type j = 0; j < nc; j++)
00281     for (octave_idx_type i = m.cidx(j); i < m.cidx(j+1); i++)
00282       y[m.ridx(i)] += m.data(i) * x[j];
00283 
00284   return true;
00285 }
00286 
00287 static bool
00288 vector_product (const Matrix& m, const double *x, double *y)
00289 {
00290   octave_idx_type nr = m.rows ();
00291   octave_idx_type nc = m.cols ();
00292 
00293   F77_XFCN (dgemv, DGEMV, (F77_CONST_CHAR_ARG2 ("N", 1),
00294                            nr, nc, 1.0,  m.data (), nr,
00295                            x, 1, 0.0, y, 1
00296                            F77_CHAR_ARG_LEN (1)));
00297 
00298   if (f77_exception_encountered)
00299     {
00300       (*current_liboctave_error_handler)
00301         ("eigs: unrecoverable error in dgemv");
00302       return false;
00303     }
00304   else
00305     return true;
00306 }
00307 
00308 static bool
00309 vector_product (const SparseComplexMatrix& m, const Complex* x,
00310                         Complex* y)
00311 {
00312   octave_idx_type nc = m.cols ();
00313 
00314   for (octave_idx_type j = 0; j < nc; j++)
00315     y[j] = 0.;
00316 
00317   for (octave_idx_type j = 0; j < nc; j++)
00318     for (octave_idx_type i = m.cidx(j); i < m.cidx(j+1); i++)
00319       y[m.ridx(i)] += m.data(i) * x[j];
00320 
00321   return true;
00322 }
00323 
00324 static bool
00325 vector_product (const ComplexMatrix& m, const Complex *x, Complex *y)
00326 {
00327   octave_idx_type nr = m.rows ();
00328   octave_idx_type nc = m.cols ();
00329 
00330   F77_XFCN (zgemv, ZGEMV, (F77_CONST_CHAR_ARG2 ("N", 1),
00331                            nr, nc, 1.0,  m.data (), nr,
00332                            x, 1, 0.0, y, 1
00333                            F77_CHAR_ARG_LEN (1)));
00334 
00335   if (f77_exception_encountered)
00336     {
00337       (*current_liboctave_error_handler)
00338         ("eigs: unrecoverable error in zgemv");
00339       return false;
00340     }
00341   else
00342     return true;
00343 }
00344 
00345 static bool
00346 make_cholb (Matrix& b, Matrix& bt, ColumnVector& permB)
00347 {
00348   octave_idx_type info;
00349   CHOL fact (b, info);
00350   octave_idx_type n = b.cols();
00351 
00352   if (info != 0)
00353     return false;
00354   else
00355     {
00356       bt = fact.chol_matrix ();
00357       b =  bt.transpose();
00358       permB = ColumnVector(n);
00359       for (octave_idx_type i = 0; i < n; i++)
00360         permB(i) = i;
00361       return true;
00362     }
00363 }
00364 
00365 static bool
00366 make_cholb (SparseMatrix& b, SparseMatrix& bt, ColumnVector& permB)
00367 {
00368   octave_idx_type info;
00369   SparseCHOL fact (b, info, false);
00370 
00371   if (fact.P() != 0)
00372     return false;
00373   else
00374     {
00375       b = fact.L();
00376       bt = b.transpose();
00377       permB = fact.perm() - 1.0;
00378       return true;
00379     }
00380 }
00381 
00382 static bool
00383 make_cholb (ComplexMatrix& b, ComplexMatrix& bt, ColumnVector& permB)
00384 {
00385   octave_idx_type info;
00386   ComplexCHOL fact (b, info);
00387   octave_idx_type n = b.cols();
00388 
00389   if (info != 0)
00390     return false;
00391   else
00392     {
00393       bt = fact.chol_matrix ();
00394       b =  bt.hermitian();
00395       permB = ColumnVector(n);
00396       for (octave_idx_type i = 0; i < n; i++)
00397         permB(i) = i;
00398       return true;
00399     }
00400 }
00401 
00402 static bool
00403 make_cholb (SparseComplexMatrix& b, SparseComplexMatrix& bt,
00404             ColumnVector& permB)
00405 {
00406   octave_idx_type info;
00407   SparseComplexCHOL fact (b, info, false);
00408 
00409   if (fact.P() != 0)
00410     return false;
00411   else
00412     {
00413       b = fact.L();
00414       bt = b.hermitian();
00415       permB = fact.perm() - 1.0;
00416       return true;
00417     }
00418 }
00419 
00420 static bool
00421 LuAminusSigmaB (const SparseMatrix &m, const SparseMatrix &b,
00422                 bool cholB, const ColumnVector& permB, double sigma,
00423                 SparseMatrix &L, SparseMatrix &U, octave_idx_type *P,
00424                 octave_idx_type *Q)
00425 {
00426   bool have_b = ! b.is_empty ();
00427   octave_idx_type n = m.rows();
00428 
00429   // Caclulate LU decomposition of 'A - sigma * B'
00430   SparseMatrix AminusSigmaB (m);
00431 
00432   if (have_b)
00433     {
00434       if (cholB)
00435         {
00436           if (permB.length())
00437             {
00438               SparseMatrix tmp(n,n,n);
00439               for (octave_idx_type i = 0; i < n; i++)
00440                 {
00441                   tmp.xcidx(i) = i;
00442                   tmp.xridx(i) =
00443                     static_cast<octave_idx_type>(permB(i));
00444                   tmp.xdata(i) = 1;
00445                 }
00446               tmp.xcidx(n) = n;
00447 
00448               AminusSigmaB = AminusSigmaB - sigma * tmp *
00449                 b.transpose() * b * tmp.transpose();
00450             }
00451           else
00452             AminusSigmaB = AminusSigmaB - sigma *
00453               b.transpose() * b;
00454         }
00455       else
00456         AminusSigmaB = AminusSigmaB - sigma * b;
00457     }
00458   else
00459     {
00460       SparseMatrix sigmat (n, n, n);
00461 
00462           // Create sigma * speye(n,n)
00463           sigmat.xcidx (0) = 0;
00464           for (octave_idx_type i = 0; i < n; i++)
00465             {
00466               sigmat.xdata(i) = sigma;
00467               sigmat.xridx(i) = i;
00468               sigmat.xcidx(i+1) = i + 1;
00469             }
00470 
00471           AminusSigmaB = AminusSigmaB - sigmat;
00472         }
00473 
00474   SparseLU fact (AminusSigmaB);
00475 
00476   L = fact.L ();
00477   U = fact.U ();
00478   const octave_idx_type *P2 = fact.row_perm ();
00479   const octave_idx_type *Q2 = fact.col_perm ();
00480 
00481   for (octave_idx_type j = 0; j < n; j++)
00482     {
00483       P[j] = P2[j];
00484       Q[j] = Q2[j];
00485     }
00486 
00487   // Test condition number of LU decomposition
00488   double minU = octave_NaN;
00489   double maxU = octave_NaN;
00490   for (octave_idx_type j = 0; j < n; j++)
00491     {
00492       double d = 0.;
00493       if (U.xcidx(j+1) > U.xcidx(j) &&
00494           U.xridx (U.xcidx(j+1)-1) == j)
00495         d = std::abs (U.xdata (U.xcidx(j+1)-1));
00496 
00497       if (xisnan (minU) || d < minU)
00498         minU = d;
00499 
00500       if (xisnan (maxU) || d > maxU)
00501         maxU = d;
00502     }
00503 
00504   double rcond = (minU / maxU);
00505   volatile double rcond_plus_one = rcond + 1.0;
00506 
00507   if (rcond_plus_one == 1.0 || xisnan (rcond))
00508     {
00509       (*current_liboctave_warning_handler)
00510         ("eigs: 'A - sigma*B' is singular, indicating sigma is exactly");
00511       (*current_liboctave_warning_handler)
00512         ("       an eigenvalue. Convergence is not guaranteed");
00513     }
00514 
00515   return true;
00516 }
00517 
00518 static bool
00519 LuAminusSigmaB (const Matrix &m, const Matrix &b,
00520                 bool cholB, const ColumnVector& permB, double sigma,
00521                 Matrix &L, Matrix &U, octave_idx_type *P,
00522                 octave_idx_type *Q)
00523 {
00524   bool have_b = ! b.is_empty ();
00525   octave_idx_type n = m.cols();
00526 
00527   // Caclulate LU decomposition of 'A - sigma * B'
00528   Matrix AminusSigmaB (m);
00529 
00530   if (have_b)
00531     {
00532       if (cholB)
00533         {
00534           Matrix tmp = sigma * b.transpose() * b;
00535           const double *pB = permB.fortran_vec();
00536           double *p = AminusSigmaB.fortran_vec();
00537 
00538           if (permB.length())
00539             {
00540               for (octave_idx_type j = 0;
00541                    j < b.cols(); j++)
00542                 for (octave_idx_type i = 0;
00543                      i < b.rows(); i++)
00544                   *p++ -=  tmp.xelem (static_cast<octave_idx_type>(pB[i]),
00545                                       static_cast<octave_idx_type>(pB[j]));
00546             }
00547           else
00548             AminusSigmaB = AminusSigmaB - tmp;
00549         }
00550       else
00551         AminusSigmaB = AminusSigmaB - sigma * b;
00552     }
00553   else
00554     {
00555       double *p = AminusSigmaB.fortran_vec();
00556 
00557       for (octave_idx_type i = 0; i < n; i++)
00558         p[i*(n+1)] -= sigma;
00559     }
00560 
00561   LU fact (AminusSigmaB);
00562 
00563   L = fact.P().transpose() * fact.L ();
00564   U = fact.U ();
00565   for (octave_idx_type j = 0; j < n; j++)
00566     P[j] = Q[j] = j;
00567 
00568   // Test condition number of LU decomposition
00569   double minU = octave_NaN;
00570   double maxU = octave_NaN;
00571   for (octave_idx_type j = 0; j < n; j++)
00572     {
00573       double d = std::abs (U.xelem(j,j));
00574       if (xisnan (minU) || d < minU)
00575         minU = d;
00576 
00577       if (xisnan (maxU) || d > maxU)
00578         maxU = d;
00579     }
00580 
00581   double rcond = (minU / maxU);
00582   volatile double rcond_plus_one = rcond + 1.0;
00583 
00584   if (rcond_plus_one == 1.0 || xisnan (rcond))
00585     {
00586       (*current_liboctave_warning_handler)
00587         ("eigs: 'A - sigma*B' is singular, indicating sigma is exactly");
00588       (*current_liboctave_warning_handler)
00589         ("       an eigenvalue. Convergence is not guaranteed");
00590     }
00591 
00592   return true;
00593 }
00594 
00595 static bool
00596 LuAminusSigmaB (const SparseComplexMatrix &m, const SparseComplexMatrix &b,
00597                 bool cholB, const ColumnVector& permB, Complex sigma,
00598                 SparseComplexMatrix &L, SparseComplexMatrix &U,
00599                 octave_idx_type *P, octave_idx_type *Q)
00600 {
00601   bool have_b = ! b.is_empty ();
00602   octave_idx_type n = m.rows();
00603 
00604   // Caclulate LU decomposition of 'A - sigma * B'
00605   SparseComplexMatrix AminusSigmaB (m);
00606 
00607   if (have_b)
00608     {
00609       if (cholB)
00610         {
00611           if (permB.length())
00612             {
00613               SparseMatrix tmp(n,n,n);
00614               for (octave_idx_type i = 0; i < n; i++)
00615                 {
00616                   tmp.xcidx(i) = i;
00617                   tmp.xridx(i) =
00618                     static_cast<octave_idx_type>(permB(i));
00619                   tmp.xdata(i) = 1;
00620                 }
00621               tmp.xcidx(n) = n;
00622 
00623               AminusSigmaB = AminusSigmaB - tmp * b.hermitian() * b *
00624                 tmp.transpose() * sigma;
00625             }
00626           else
00627             AminusSigmaB = AminusSigmaB - sigma * b.hermitian() * b;
00628         }
00629       else
00630         AminusSigmaB = AminusSigmaB - sigma * b;
00631     }
00632   else
00633     {
00634       SparseComplexMatrix sigmat (n, n, n);
00635 
00636       // Create sigma * speye(n,n)
00637       sigmat.xcidx (0) = 0;
00638       for (octave_idx_type i = 0; i < n; i++)
00639         {
00640           sigmat.xdata(i) = sigma;
00641           sigmat.xridx(i) = i;
00642           sigmat.xcidx(i+1) = i + 1;
00643         }
00644 
00645       AminusSigmaB = AminusSigmaB - sigmat;
00646     }
00647 
00648   SparseComplexLU fact (AminusSigmaB);
00649 
00650   L = fact.L ();
00651   U = fact.U ();
00652   const octave_idx_type *P2 = fact.row_perm ();
00653   const octave_idx_type *Q2 = fact.col_perm ();
00654 
00655   for (octave_idx_type j = 0; j < n; j++)
00656     {
00657       P[j] = P2[j];
00658       Q[j] = Q2[j];
00659     }
00660 
00661   // Test condition number of LU decomposition
00662   double minU = octave_NaN;
00663   double maxU = octave_NaN;
00664   for (octave_idx_type j = 0; j < n; j++)
00665     {
00666       double d = 0.;
00667       if (U.xcidx(j+1) > U.xcidx(j) &&
00668           U.xridx (U.xcidx(j+1)-1) == j)
00669         d = std::abs (U.xdata (U.xcidx(j+1)-1));
00670 
00671       if (xisnan (minU) || d < minU)
00672         minU = d;
00673 
00674       if (xisnan (maxU) || d > maxU)
00675         maxU = d;
00676     }
00677 
00678   double rcond = (minU / maxU);
00679   volatile double rcond_plus_one = rcond + 1.0;
00680 
00681   if (rcond_plus_one == 1.0 || xisnan (rcond))
00682     {
00683       (*current_liboctave_warning_handler)
00684         ("eigs: 'A - sigma*B' is singular, indicating sigma is exactly");
00685       (*current_liboctave_warning_handler)
00686         ("       an eigenvalue. Convergence is not guaranteed");
00687     }
00688 
00689   return true;
00690 }
00691 
00692 static bool
00693 LuAminusSigmaB (const ComplexMatrix &m, const ComplexMatrix &b,
00694                 bool cholB, const ColumnVector& permB, Complex sigma,
00695                 ComplexMatrix &L, ComplexMatrix &U, octave_idx_type *P,
00696                 octave_idx_type *Q)
00697 {
00698   bool have_b = ! b.is_empty ();
00699   octave_idx_type n = m.cols();
00700 
00701   // Caclulate LU decomposition of 'A - sigma * B'
00702   ComplexMatrix AminusSigmaB (m);
00703 
00704   if (have_b)
00705     {
00706       if (cholB)
00707         {
00708           ComplexMatrix tmp = sigma * b.hermitian() * b;
00709           const double *pB = permB.fortran_vec();
00710           Complex *p = AminusSigmaB.fortran_vec();
00711 
00712           if (permB.length())
00713             {
00714               for (octave_idx_type j = 0;
00715                    j < b.cols(); j++)
00716                 for (octave_idx_type i = 0;
00717                      i < b.rows(); i++)
00718                   *p++ -=  tmp.xelem (static_cast<octave_idx_type>(pB[i]),
00719                                       static_cast<octave_idx_type>(pB[j]));
00720             }
00721           else
00722             AminusSigmaB = AminusSigmaB - tmp;
00723         }
00724       else
00725         AminusSigmaB = AminusSigmaB - sigma * b;
00726     }
00727   else
00728     {
00729       Complex *p = AminusSigmaB.fortran_vec();
00730 
00731       for (octave_idx_type i = 0; i < n; i++)
00732         p[i*(n+1)] -= sigma;
00733     }
00734 
00735   ComplexLU fact (AminusSigmaB);
00736 
00737   L = fact.P().transpose() * fact.L ();
00738   U = fact.U ();
00739   for (octave_idx_type j = 0; j < n; j++)
00740     P[j] = Q[j] = j;
00741 
00742   // Test condition number of LU decomposition
00743   double minU = octave_NaN;
00744   double maxU = octave_NaN;
00745   for (octave_idx_type j = 0; j < n; j++)
00746     {
00747       double d = std::abs (U.xelem(j,j));
00748       if (xisnan (minU) || d < minU)
00749         minU = d;
00750 
00751       if (xisnan (maxU) || d > maxU)
00752         maxU = d;
00753     }
00754 
00755   double rcond = (minU / maxU);
00756   volatile double rcond_plus_one = rcond + 1.0;
00757 
00758   if (rcond_plus_one == 1.0 || xisnan (rcond))
00759     {
00760       (*current_liboctave_warning_handler)
00761         ("eigs: 'A - sigma*B' is singular, indicating sigma is exactly");
00762       (*current_liboctave_warning_handler)
00763         ("       an eigenvalue. Convergence is not guaranteed");
00764     }
00765 
00766   return true;
00767 }
00768 
00769 template <class M>
00770 octave_idx_type
00771 EigsRealSymmetricMatrix (const M& m, const std::string typ,
00772                          octave_idx_type k, octave_idx_type p,
00773                          octave_idx_type &info, Matrix &eig_vec,
00774                          ColumnVector &eig_val, const M& _b,
00775                          ColumnVector &permB, ColumnVector &resid,
00776                          std::ostream& os, double tol, bool rvec,
00777                          bool cholB, int disp, int maxit)
00778 {
00779   M b(_b);
00780   octave_idx_type n = m.cols ();
00781   octave_idx_type mode = 1;
00782   bool have_b = ! b.is_empty();
00783   bool note3 = false;
00784   char bmat = 'I';
00785   double sigma = 0.;
00786   M bt;
00787 
00788   if (m.rows() != m.cols())
00789     {
00790       (*current_liboctave_error_handler) ("eigs: A must be square");
00791       return -1;
00792     }
00793   if (have_b && (m.rows() != b.rows() || m.rows() != b.cols()))
00794     {
00795       (*current_liboctave_error_handler)
00796         ("eigs: B must be square and the same size as A");
00797       return -1;
00798     }
00799 
00800   if (resid.is_empty())
00801     {
00802       std::string rand_dist = octave_rand::distribution();
00803       octave_rand::distribution("uniform");
00804       resid = ColumnVector (octave_rand::vector(n));
00805       octave_rand::distribution(rand_dist);
00806     }
00807 
00808   if (n < 3)
00809     {
00810       (*current_liboctave_error_handler)
00811         ("eigs: n must be at least 3");
00812       return -1;
00813     }
00814 
00815   if (p < 0)
00816     {
00817       p = k * 2;
00818 
00819       if (p < 20)
00820         p = 20;
00821 
00822       if (p > n - 1)
00823         p = n - 1 ;
00824     }
00825 
00826   if (k < 1 || k > n - 2)
00827     {
00828       (*current_liboctave_error_handler)
00829         ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1-1).\n"
00830          "      Use 'eig(full(A))' instead");
00831       return -1;
00832     }
00833 
00834   if (p <= k || p >= n)
00835     {
00836       (*current_liboctave_error_handler)
00837         ("eigs: opts.p must be greater than k and less than n");
00838       return -1;
00839     }
00840 
00841   if (have_b && cholB && permB.length() != 0)
00842     {
00843       // Check the we really have a permutation vector
00844       if (permB.length() != n)
00845         {
00846           (*current_liboctave_error_handler)
00847             ("eigs: permB vector invalid");
00848           return -1;
00849         }
00850       else
00851         {
00852           Array<bool> checked (dim_vector (n, 1), false);
00853           for (octave_idx_type i = 0; i < n; i++)
00854             {
00855               octave_idx_type bidx =
00856                 static_cast<octave_idx_type> (permB(i));
00857               if (checked(bidx) || bidx < 0 ||
00858                   bidx >= n || D_NINT (bidx) != bidx)
00859                 {
00860                   (*current_liboctave_error_handler)
00861                     ("eigs: permB vector invalid");
00862                   return -1;
00863                 }
00864             }
00865         }
00866     }
00867 
00868   if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" &&
00869       typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" &&
00870       typ != "SI")
00871     {
00872       (*current_liboctave_error_handler)
00873         ("eigs: unrecognized sigma value");
00874       return -1;
00875     }
00876 
00877   if (typ == "LI" || typ == "SI" || typ == "LR" || typ == "SR")
00878     {
00879       (*current_liboctave_error_handler)
00880         ("eigs: invalid sigma value for real symmetric problem");
00881       return -1;
00882     }
00883 
00884   if (have_b)
00885     {
00886       // See Note 3 dsaupd
00887       note3 = true;
00888       if (cholB)
00889         {
00890           bt = b;
00891           b = b.transpose();
00892           if (permB.length() == 0)
00893             {
00894               permB = ColumnVector(n);
00895               for (octave_idx_type i = 0; i < n; i++)
00896                 permB(i) = i;
00897             }
00898         }
00899       else
00900         {
00901           if (! make_cholb(b, bt, permB))
00902             {
00903               (*current_liboctave_error_handler)
00904                 ("eigs: The matrix B is not positive definite");
00905               return -1;
00906             }
00907         }
00908     }
00909 
00910   Array<octave_idx_type> ip (dim_vector (11, 1));
00911   octave_idx_type *iparam = ip.fortran_vec ();
00912 
00913   ip(0) = 1; //ishift
00914   ip(1) = 0;   // ip(1) not referenced
00915   ip(2) = maxit; // mxiter, maximum number of iterations
00916   ip(3) = 1; // NB blocksize in recurrence
00917   ip(4) = 0; // nconv, number of Ritz values that satisfy convergence
00918   ip(5) = 0; //ip(5) not referenced
00919   ip(6) = mode; // mode
00920   ip(7) = 0;
00921   ip(8) = 0;
00922   ip(9) = 0;
00923   ip(10) = 0;
00924   // ip(7) to ip(10) return values
00925 
00926   Array<octave_idx_type> iptr (dim_vector (14, 1));
00927   octave_idx_type *ipntr = iptr.fortran_vec ();
00928 
00929   octave_idx_type ido = 0;
00930   int iter = 0;
00931   octave_idx_type lwork = p * (p + 8);
00932 
00933   OCTAVE_LOCAL_BUFFER (double, v, n * p);
00934   OCTAVE_LOCAL_BUFFER (double, workl, lwork);
00935   OCTAVE_LOCAL_BUFFER (double, workd, 3 * n);
00936   double *presid = resid.fortran_vec ();
00937 
00938   do
00939     {
00940       F77_FUNC (dsaupd, DSAUPD)
00941         (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
00942          F77_CONST_CHAR_ARG2 ((typ.c_str()), 2),
00943          k, tol, presid, p, v, n, iparam,
00944          ipntr, workd, workl, lwork, info
00945          F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
00946 
00947       if (f77_exception_encountered)
00948         {
00949           (*current_liboctave_error_handler)
00950             ("eigs: unrecoverable exception encountered in dsaupd");
00951           return -1;
00952         }
00953 
00954       if (disp > 0 && !xisnan(workl[iptr(5)-1]))
00955         {
00956           if (iter++)
00957             {
00958               os << "Iteration " << iter - 1 <<
00959                 ": a few Ritz values of the " << p << "-by-" <<
00960                 p << " matrix\n";
00961               for (int i = 0 ; i < k; i++)
00962                 os << "    " << workl[iptr(5)+i-1] << "\n";
00963             }
00964 
00965           // This is a kludge, as ARPACK doesn't give its
00966           // iteration pointer. But as workl[iptr(5)-1] is
00967           // an output value updated at each iteration, setting
00968           // a value in this array to NaN and testing for it
00969           // is a way of obtaining the iteration counter.
00970           if (ido != 99)
00971             workl[iptr(5)-1] = octave_NaN;
00972         }
00973 
00974       if (ido == -1 || ido == 1 || ido == 2)
00975         {
00976           if (have_b)
00977             {
00978               Matrix mtmp (n,1);
00979               for (octave_idx_type i = 0; i < n; i++)
00980                 mtmp(i,0) = workd[i + iptr(0) - 1];
00981 
00982               mtmp = utsolve(bt, permB, m * ltsolve(b, permB, mtmp));
00983 
00984               for (octave_idx_type i = 0; i < n; i++)
00985                 workd[i+iptr(1)-1] = mtmp(i,0);
00986             }
00987           else if (!vector_product (m, workd + iptr(0) - 1,
00988                                     workd + iptr(1) - 1))
00989             break;
00990         }
00991       else
00992         {
00993           if (info < 0)
00994             {
00995               (*current_liboctave_error_handler)
00996                 ("eigs: error %d in dsaupd", info);
00997               return -1;
00998             }
00999           break;
01000         }
01001     }
01002   while (1);
01003 
01004   octave_idx_type info2;
01005 
01006   // We have a problem in that the size of the C++ bool
01007   // type relative to the fortran logical type. It appears
01008   // that fortran uses 4- or 8-bytes per logical and C++ 1-byte
01009   // per bool, though this might be system dependent. As
01010   // long as the HOWMNY arg is not "S", the logical array
01011   // is just workspace for ARPACK, so use int type to
01012   // avoid problems.
01013   Array<octave_idx_type> s (dim_vector (p, 1));
01014   octave_idx_type *sel = s.fortran_vec ();
01015 
01016   eig_vec.resize (n, k);
01017   double *z = eig_vec.fortran_vec ();
01018 
01019   eig_val.resize (k);
01020   double *d = eig_val.fortran_vec ();
01021 
01022   F77_FUNC (dseupd, DSEUPD)
01023     (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma,
01024      F77_CONST_CHAR_ARG2 (&bmat, 1), n,
01025      F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), k, tol, presid, p, v, n, iparam,
01026      ipntr, workd, workl, lwork, info2 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1)
01027      F77_CHAR_ARG_LEN(2));
01028 
01029   if (f77_exception_encountered)
01030     {
01031       (*current_liboctave_error_handler)
01032         ("eigs: unrecoverable exception encountered in dseupd");
01033       return -1;
01034     }
01035   else
01036     {
01037       if (info2 == 0)
01038         {
01039           octave_idx_type k2 = k / 2;
01040           if (typ != "SM" && typ != "BE")
01041             {
01042               for (octave_idx_type i = 0; i < k2; i++)
01043                 {
01044                   double dtmp = d[i];
01045                   d[i] = d[k - i - 1];
01046                   d[k - i - 1] = dtmp;
01047                 }
01048             }
01049 
01050           if (rvec)
01051             {
01052               if (typ != "SM" && typ != "BE")
01053                 {
01054                   OCTAVE_LOCAL_BUFFER (double, dtmp, n);
01055 
01056                   for (octave_idx_type i = 0; i < k2; i++)
01057                     {
01058                       octave_idx_type off1 = i * n;
01059                       octave_idx_type off2 = (k - i - 1) * n;
01060 
01061                       if (off1 == off2)
01062                         continue;
01063 
01064                       for (octave_idx_type j = 0; j < n; j++)
01065                         dtmp[j] = z[off1 + j];
01066 
01067                       for (octave_idx_type j = 0; j < n; j++)
01068                         z[off1 + j] = z[off2 + j];
01069 
01070                       for (octave_idx_type j = 0; j < n; j++)
01071                         z[off2 + j] = dtmp[j];
01072                     }
01073                 }
01074 
01075               if (note3)
01076                 eig_vec = ltsolve(b, permB, eig_vec);
01077             }
01078         }
01079       else
01080         {
01081           (*current_liboctave_error_handler)
01082             ("eigs: error %d in dseupd", info2);
01083           return -1;
01084         }
01085     }
01086 
01087   return ip(4);
01088 }
01089 
01090 template <class M>
01091 octave_idx_type
01092 EigsRealSymmetricMatrixShift (const M& m, double sigma,
01093                               octave_idx_type k, octave_idx_type p,
01094                               octave_idx_type &info, Matrix &eig_vec,
01095                               ColumnVector &eig_val, const M& _b,
01096                               ColumnVector &permB, ColumnVector &resid,
01097                               std::ostream& os, double tol, bool rvec,
01098                               bool cholB, int disp, int maxit)
01099 {
01100   M b(_b);
01101   octave_idx_type n = m.cols ();
01102   octave_idx_type mode = 3;
01103   bool have_b = ! b.is_empty();
01104   std::string typ = "LM";
01105 
01106   if (m.rows() != m.cols())
01107     {
01108       (*current_liboctave_error_handler) ("eigs: A must be square");
01109       return -1;
01110     }
01111   if (have_b && (m.rows() != b.rows() || m.rows() != b.cols()))
01112     {
01113       (*current_liboctave_error_handler)
01114         ("eigs: B must be square and the same size as A");
01115       return -1;
01116     }
01117 
01118   // FIXME: The "SM" type for mode 1 seems unstable though faster!!
01119   //if (! std::abs (sigma))
01120   //  return EigsRealSymmetricMatrix (m, "SM", k, p, info, eig_vec, eig_val,
01121   //                                _b, permB, resid, os, tol, rvec, cholB,
01122   //                                disp, maxit);
01123 
01124   if (resid.is_empty())
01125     {
01126       std::string rand_dist = octave_rand::distribution();
01127       octave_rand::distribution("uniform");
01128       resid = ColumnVector (octave_rand::vector(n));
01129       octave_rand::distribution(rand_dist);
01130     }
01131 
01132   if (n < 3)
01133     {
01134       (*current_liboctave_error_handler)
01135         ("eigs: n must be at least 3");
01136       return -1;
01137     }
01138 
01139   if (k <= 0 || k >= n - 1)
01140     {
01141       (*current_liboctave_error_handler)
01142         ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1-1).\n"
01143              "      Use 'eig(full(A))' instead");
01144       return -1;
01145     }
01146 
01147   if (p < 0)
01148     {
01149       p = k * 2;
01150 
01151       if (p < 20)
01152         p = 20;
01153 
01154       if (p > n - 1)
01155         p = n - 1 ;
01156     }
01157 
01158   if (p <= k || p >= n)
01159     {
01160       (*current_liboctave_error_handler)
01161         ("eigs: opts.p must be greater than k and less than n");
01162       return -1;
01163     }
01164 
01165   if (have_b && cholB && permB.length() != 0)
01166     {
01167       // Check the we really have a permutation vector
01168       if (permB.length() != n)
01169         {
01170           (*current_liboctave_error_handler) ("eigs: permB vector invalid");
01171           return -1;
01172         }
01173       else
01174         {
01175           Array<bool> checked (dim_vector (n, 1), false);
01176           for (octave_idx_type i = 0; i < n; i++)
01177             {
01178               octave_idx_type bidx =
01179                 static_cast<octave_idx_type> (permB(i));
01180               if (checked(bidx) || bidx < 0 ||
01181                   bidx >= n || D_NINT (bidx) != bidx)
01182                 {
01183                   (*current_liboctave_error_handler)
01184                     ("eigs: permB vector invalid");
01185                   return -1;
01186                 }
01187             }
01188         }
01189     }
01190 
01191   char bmat = 'I';
01192   if (have_b)
01193     bmat = 'G';
01194 
01195   Array<octave_idx_type> ip (dim_vector (11, 1));
01196   octave_idx_type *iparam = ip.fortran_vec ();
01197 
01198   ip(0) = 1; //ishift
01199   ip(1) = 0;   // ip(1) not referenced
01200   ip(2) = maxit; // mxiter, maximum number of iterations
01201   ip(3) = 1; // NB blocksize in recurrence
01202   ip(4) = 0; // nconv, number of Ritz values that satisfy convergence
01203   ip(5) = 0; //ip(5) not referenced
01204   ip(6) = mode; // mode
01205   ip(7) = 0;
01206   ip(8) = 0;
01207   ip(9) = 0;
01208   ip(10) = 0;
01209   // ip(7) to ip(10) return values
01210 
01211   Array<octave_idx_type> iptr (dim_vector (14, 1));
01212   octave_idx_type *ipntr = iptr.fortran_vec ();
01213 
01214   octave_idx_type ido = 0;
01215   int iter = 0;
01216   M L, U;
01217 
01218   OCTAVE_LOCAL_BUFFER (octave_idx_type, P, (have_b ? b.rows() : m.rows()));
01219   OCTAVE_LOCAL_BUFFER (octave_idx_type, Q, (have_b ? b.cols() : m.cols()));
01220 
01221   if (! LuAminusSigmaB(m, b, cholB, permB, sigma, L, U, P, Q))
01222     return -1;
01223 
01224   octave_idx_type lwork = p * (p + 8);
01225 
01226   OCTAVE_LOCAL_BUFFER (double, v, n * p);
01227   OCTAVE_LOCAL_BUFFER (double, workl, lwork);
01228   OCTAVE_LOCAL_BUFFER (double, workd, 3 * n);
01229   double *presid = resid.fortran_vec ();
01230 
01231   do
01232     {
01233       F77_FUNC (dsaupd, DSAUPD)
01234         (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
01235          F77_CONST_CHAR_ARG2 ((typ.c_str()), 2),
01236          k, tol, presid, p, v, n, iparam,
01237          ipntr, workd, workl, lwork, info
01238          F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
01239 
01240       if (f77_exception_encountered)
01241         {
01242           (*current_liboctave_error_handler)
01243             ("eigs: unrecoverable exception encountered in dsaupd");
01244           return -1;
01245         }
01246 
01247       if (disp > 0 && !xisnan(workl[iptr(5)-1]))
01248         {
01249           if (iter++)
01250             {
01251               os << "Iteration " << iter - 1 <<
01252                 ": a few Ritz values of the " << p << "-by-" <<
01253                 p << " matrix\n";
01254               for (int i = 0 ; i < k; i++)
01255                 os << "    " << workl[iptr(5)+i-1] << "\n";
01256             }
01257 
01258           // This is a kludge, as ARPACK doesn't give its
01259           // iteration pointer. But as workl[iptr(5)-1] is
01260           // an output value updated at each iteration, setting
01261           // a value in this array to NaN and testing for it
01262           // is a way of obtaining the iteration counter.
01263           if (ido != 99)
01264             workl[iptr(5)-1] = octave_NaN;
01265         }
01266 
01267       if (ido == -1 || ido == 1 || ido == 2)
01268         {
01269           if (have_b)
01270             {
01271               if (ido == -1)
01272                 {
01273                   OCTAVE_LOCAL_BUFFER (double, dtmp, n);
01274 
01275                   vector_product (m, workd+iptr(0)-1, dtmp);
01276 
01277                   Matrix tmp(n, 1);
01278 
01279                   for (octave_idx_type i = 0; i < n; i++)
01280                     tmp(i,0) = dtmp[P[i]];
01281 
01282                   lusolve (L, U, tmp);
01283 
01284                   double *ip2 = workd+iptr(1)-1;
01285                   for (octave_idx_type i = 0; i < n; i++)
01286                     ip2[Q[i]] = tmp(i,0);
01287                 }
01288               else if (ido == 2)
01289                 vector_product (b, workd+iptr(0)-1, workd+iptr(1)-1);
01290               else
01291                 {
01292                   double *ip2 = workd+iptr(2)-1;
01293                   Matrix tmp(n, 1);
01294 
01295                   for (octave_idx_type i = 0; i < n; i++)
01296                     tmp(i,0) = ip2[P[i]];
01297 
01298                   lusolve (L, U, tmp);
01299 
01300                   ip2 = workd+iptr(1)-1;
01301                   for (octave_idx_type i = 0; i < n; i++)
01302                     ip2[Q[i]] = tmp(i,0);
01303                 }
01304             }
01305           else
01306             {
01307               if (ido == 2)
01308                 {
01309                   for (octave_idx_type i = 0; i < n; i++)
01310                     workd[iptr(0) + i - 1] = workd[iptr(1) + i - 1];
01311                 }
01312               else
01313                 {
01314                   double *ip2 = workd+iptr(0)-1;
01315                   Matrix tmp(n, 1);
01316 
01317                   for (octave_idx_type i = 0; i < n; i++)
01318                     tmp(i,0) = ip2[P[i]];
01319 
01320                   lusolve (L, U, tmp);
01321 
01322                   ip2 = workd+iptr(1)-1;
01323                   for (octave_idx_type i = 0; i < n; i++)
01324                     ip2[Q[i]] = tmp(i,0);
01325                 }
01326             }
01327         }
01328       else
01329         {
01330           if (info < 0)
01331             {
01332               (*current_liboctave_error_handler)
01333                 ("eigs: error %d in dsaupd", info);
01334               return -1;
01335             }
01336           break;
01337         }
01338     }
01339   while (1);
01340 
01341   octave_idx_type info2;
01342 
01343   // We have a problem in that the size of the C++ bool
01344   // type relative to the fortran logical type. It appears
01345   // that fortran uses 4- or 8-bytes per logical and C++ 1-byte
01346   // per bool, though this might be system dependent. As
01347   // long as the HOWMNY arg is not "S", the logical array
01348   // is just workspace for ARPACK, so use int type to
01349   // avoid problems.
01350   Array<octave_idx_type> s (dim_vector (p, 1));
01351   octave_idx_type *sel = s.fortran_vec ();
01352 
01353   eig_vec.resize (n, k);
01354   double *z = eig_vec.fortran_vec ();
01355 
01356   eig_val.resize (k);
01357   double *d = eig_val.fortran_vec ();
01358 
01359   F77_FUNC (dseupd, DSEUPD)
01360     (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma,
01361      F77_CONST_CHAR_ARG2 (&bmat, 1), n,
01362      F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2),
01363      k, tol, presid, p, v, n, iparam, ipntr, workd, workl, lwork, info2
01364      F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
01365 
01366   if (f77_exception_encountered)
01367     {
01368       (*current_liboctave_error_handler)
01369         ("eigs: unrecoverable exception encountered in dseupd");
01370       return -1;
01371     }
01372   else
01373     {
01374       if (info2 == 0)
01375         {
01376           octave_idx_type k2 = k / 2;
01377           for (octave_idx_type i = 0; i < k2; i++)
01378             {
01379               double dtmp = d[i];
01380               d[i] = d[k - i - 1];
01381               d[k - i - 1] = dtmp;
01382             }
01383 
01384           if (rvec)
01385             {
01386               OCTAVE_LOCAL_BUFFER (double, dtmp, n);
01387 
01388               for (octave_idx_type i = 0; i < k2; i++)
01389                 {
01390                   octave_idx_type off1 = i * n;
01391                   octave_idx_type off2 = (k - i - 1) * n;
01392 
01393                   if (off1 == off2)
01394                     continue;
01395 
01396                   for (octave_idx_type j = 0; j < n; j++)
01397                     dtmp[j] = z[off1 + j];
01398 
01399                   for (octave_idx_type j = 0; j < n; j++)
01400                     z[off1 + j] = z[off2 + j];
01401 
01402                   for (octave_idx_type j = 0; j < n; j++)
01403                     z[off2 + j] = dtmp[j];
01404                 }
01405             }
01406         }
01407       else
01408         {
01409           (*current_liboctave_error_handler)
01410             ("eigs: error %d in dseupd", info2);
01411           return -1;
01412         }
01413     }
01414 
01415   return ip(4);
01416 }
01417 
01418 octave_idx_type
01419 EigsRealSymmetricFunc (EigsFunc fun, octave_idx_type n,
01420                        const std::string &_typ, double sigma,
01421                        octave_idx_type k, octave_idx_type p,
01422                        octave_idx_type &info, Matrix &eig_vec,
01423                        ColumnVector &eig_val, ColumnVector &resid,
01424                        std::ostream& os, double tol, bool rvec,
01425                        bool /* cholB */, int disp, int maxit)
01426 {
01427   std::string typ (_typ);
01428   bool have_sigma = (sigma ? true : false);
01429   char bmat = 'I';
01430   octave_idx_type mode = 1;
01431   int err = 0;
01432 
01433   if (resid.is_empty())
01434     {
01435       std::string rand_dist = octave_rand::distribution();
01436       octave_rand::distribution("uniform");
01437       resid = ColumnVector (octave_rand::vector(n));
01438       octave_rand::distribution(rand_dist);
01439     }
01440 
01441   if (n < 3)
01442     {
01443       (*current_liboctave_error_handler)
01444         ("eigs: n must be at least 3");
01445       return -1;
01446     }
01447 
01448   if (p < 0)
01449     {
01450       p = k * 2;
01451 
01452       if (p < 20)
01453         p = 20;
01454 
01455       if (p > n - 1)
01456         p = n - 1 ;
01457     }
01458 
01459   if (k <= 0 || k >= n - 1)
01460     {
01461       (*current_liboctave_error_handler)
01462         ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1).\n"
01463              "      Use 'eig(full(A))' instead");
01464       return -1;
01465     }
01466 
01467   if (p <= k || p >= n)
01468     {
01469       (*current_liboctave_error_handler)
01470         ("eigs: opts.p must be greater than k and less than n");
01471       return -1;
01472     }
01473 
01474   if (! have_sigma)
01475     {
01476       if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" &&
01477           typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" &&
01478           typ != "SI")
01479         (*current_liboctave_error_handler)
01480           ("eigs: unrecognized sigma value");
01481 
01482       if (typ == "LI" || typ == "SI" || typ == "LR" || typ == "SR")
01483         {
01484           (*current_liboctave_error_handler)
01485             ("eigs: invalid sigma value for real symmetric problem");
01486           return -1;
01487         }
01488 
01489       if (typ == "SM")
01490         {
01491           typ = "LM";
01492           sigma = 0.;
01493           mode = 3;
01494         }
01495     }
01496   else if (! std::abs (sigma))
01497     typ = "SM";
01498   else
01499     {
01500       typ = "LM";
01501       mode = 3;
01502     }
01503 
01504   Array<octave_idx_type> ip (dim_vector (11, 1));
01505   octave_idx_type *iparam = ip.fortran_vec ();
01506 
01507   ip(0) = 1; //ishift
01508   ip(1) = 0;   // ip(1) not referenced
01509   ip(2) = maxit; // mxiter, maximum number of iterations
01510   ip(3) = 1; // NB blocksize in recurrence
01511   ip(4) = 0; // nconv, number of Ritz values that satisfy convergence
01512   ip(5) = 0; //ip(5) not referenced
01513   ip(6) = mode; // mode
01514   ip(7) = 0;
01515   ip(8) = 0;
01516   ip(9) = 0;
01517   ip(10) = 0;
01518   // ip(7) to ip(10) return values
01519 
01520   Array<octave_idx_type> iptr (dim_vector (14, 1));
01521   octave_idx_type *ipntr = iptr.fortran_vec ();
01522 
01523   octave_idx_type ido = 0;
01524   int iter = 0;
01525   octave_idx_type lwork = p * (p + 8);
01526 
01527   OCTAVE_LOCAL_BUFFER (double, v, n * p);
01528   OCTAVE_LOCAL_BUFFER (double, workl, lwork);
01529   OCTAVE_LOCAL_BUFFER (double, workd, 3 * n);
01530   double *presid = resid.fortran_vec ();
01531 
01532   do
01533     {
01534       F77_FUNC (dsaupd, DSAUPD)
01535         (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
01536          F77_CONST_CHAR_ARG2 ((typ.c_str()), 2),
01537          k, tol, presid, p, v, n, iparam,
01538          ipntr, workd, workl, lwork, info
01539          F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
01540 
01541       if (f77_exception_encountered)
01542         {
01543           (*current_liboctave_error_handler)
01544             ("eigs: unrecoverable exception encountered in dsaupd");
01545           return -1;
01546         }
01547 
01548       if (disp > 0 && !xisnan(workl[iptr(5)-1]))
01549         {
01550           if (iter++)
01551             {
01552               os << "Iteration " << iter - 1 <<
01553                 ": a few Ritz values of the " << p << "-by-" <<
01554                 p << " matrix\n";
01555               for (int i = 0 ; i < k; i++)
01556                 os << "    " << workl[iptr(5)+i-1] << "\n";
01557             }
01558 
01559           // This is a kludge, as ARPACK doesn't give its
01560           // iteration pointer. But as workl[iptr(5)-1] is
01561           // an output value updated at each iteration, setting
01562           // a value in this array to NaN and testing for it
01563           // is a way of obtaining the iteration counter.
01564           if (ido != 99)
01565             workl[iptr(5)-1] = octave_NaN;
01566         }
01567 
01568 
01569       if (ido == -1 || ido == 1 || ido == 2)
01570         {
01571           double *ip2 = workd + iptr(0) - 1;
01572           ColumnVector x(n);
01573 
01574           for (octave_idx_type i = 0; i < n; i++)
01575             x(i) = *ip2++;
01576 
01577           ColumnVector y = fun (x, err);
01578 
01579           if (err)
01580             return false;
01581 
01582           ip2 = workd + iptr(1) - 1;
01583           for (octave_idx_type i = 0; i < n; i++)
01584             *ip2++ = y(i);
01585         }
01586       else
01587         {
01588           if (info < 0)
01589             {
01590               (*current_liboctave_error_handler)
01591                 ("eigs: error %d in dsaupd", info);
01592               return -1;
01593             }
01594           break;
01595         }
01596     }
01597   while (1);
01598 
01599   octave_idx_type info2;
01600 
01601   // We have a problem in that the size of the C++ bool
01602   // type relative to the fortran logical type. It appears
01603   // that fortran uses 4- or 8-bytes per logical and C++ 1-byte
01604   // per bool, though this might be system dependent. As
01605   // long as the HOWMNY arg is not "S", the logical array
01606   // is just workspace for ARPACK, so use int type to
01607   // avoid problems.
01608   Array<octave_idx_type> s (dim_vector (p, 1));
01609   octave_idx_type *sel = s.fortran_vec ();
01610 
01611   eig_vec.resize (n, k);
01612   double *z = eig_vec.fortran_vec ();
01613 
01614   eig_val.resize (k);
01615   double *d = eig_val.fortran_vec ();
01616 
01617   F77_FUNC (dseupd, DSEUPD)
01618     (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma,
01619      F77_CONST_CHAR_ARG2 (&bmat, 1), n,
01620      F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2),
01621      k, tol, presid, p, v, n, iparam, ipntr, workd, workl, lwork, info2
01622      F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
01623 
01624   if (f77_exception_encountered)
01625     {
01626       (*current_liboctave_error_handler)
01627         ("eigs: unrecoverable exception encountered in dseupd");
01628       return -1;
01629     }
01630   else
01631     {
01632       if (info2 == 0)
01633         {
01634           octave_idx_type k2 = k / 2;
01635           if (typ != "SM" && typ != "BE")
01636             {
01637               for (octave_idx_type i = 0; i < k2; i++)
01638                 {
01639                   double dtmp = d[i];
01640                   d[i] = d[k - i - 1];
01641                   d[k - i - 1] = dtmp;
01642                 }
01643             }
01644 
01645           if (rvec)
01646             {
01647               if (typ != "SM" && typ != "BE")
01648                 {
01649                   OCTAVE_LOCAL_BUFFER (double, dtmp, n);
01650 
01651                   for (octave_idx_type i = 0; i < k2; i++)
01652                     {
01653                       octave_idx_type off1 = i * n;
01654                       octave_idx_type off2 = (k - i - 1) * n;
01655 
01656                       if (off1 == off2)
01657                         continue;
01658 
01659                       for (octave_idx_type j = 0; j < n; j++)
01660                         dtmp[j] = z[off1 + j];
01661 
01662                       for (octave_idx_type j = 0; j < n; j++)
01663                         z[off1 + j] = z[off2 + j];
01664 
01665                       for (octave_idx_type j = 0; j < n; j++)
01666                         z[off2 + j] = dtmp[j];
01667                     }
01668                 }
01669             }
01670         }
01671       else
01672         {
01673           (*current_liboctave_error_handler)
01674             ("eigs: error %d in dseupd", info2);
01675           return -1;
01676         }
01677     }
01678 
01679   return ip(4);
01680 }
01681 
01682 template <class M>
01683 octave_idx_type
01684 EigsRealNonSymmetricMatrix (const M& m, const std::string typ,
01685                             octave_idx_type k, octave_idx_type p,
01686                             octave_idx_type &info, ComplexMatrix &eig_vec,
01687                             ComplexColumnVector &eig_val, const M& _b,
01688                             ColumnVector &permB, ColumnVector &resid,
01689                             std::ostream& os, double tol, bool rvec,
01690                             bool cholB, int disp, int maxit)
01691 {
01692   M b(_b);
01693   octave_idx_type n = m.cols ();
01694   octave_idx_type mode = 1;
01695   bool have_b = ! b.is_empty();
01696   bool note3 = false;
01697   char bmat = 'I';
01698   double sigmar = 0.;
01699   double sigmai = 0.;
01700   M bt;
01701 
01702   if (m.rows() != m.cols())
01703     {
01704       (*current_liboctave_error_handler) ("eigs: A must be square");
01705       return -1;
01706     }
01707   if (have_b && (m.rows() != b.rows() || m.rows() != b.cols()))
01708     {
01709       (*current_liboctave_error_handler)
01710         ("eigs: B must be square and the same size as A");
01711       return -1;
01712     }
01713 
01714   if (resid.is_empty())
01715     {
01716       std::string rand_dist = octave_rand::distribution();
01717       octave_rand::distribution("uniform");
01718       resid = ColumnVector (octave_rand::vector(n));
01719       octave_rand::distribution(rand_dist);
01720     }
01721 
01722   if (n < 3)
01723     {
01724       (*current_liboctave_error_handler)
01725         ("eigs: n must be at least 3");
01726       return -1;
01727     }
01728 
01729   if (p < 0)
01730     {
01731       p = k * 2 + 1;
01732 
01733       if (p < 20)
01734         p = 20;
01735 
01736       if (p > n - 1)
01737         p = n - 1 ;
01738     }
01739 
01740   if (k <= 0 || k >= n - 1)
01741     {
01742       (*current_liboctave_error_handler)
01743         ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1).\n"
01744          "      Use 'eig(full(A))' instead");
01745       return -1;
01746     }
01747 
01748   if (p <= k || p >= n)
01749     {
01750       (*current_liboctave_error_handler)
01751         ("eigs: opts.p must be greater than k and less than n");
01752       return -1;
01753     }
01754 
01755   if (have_b && cholB && permB.length() != 0)
01756     {
01757       // Check the we really have a permutation vector
01758       if (permB.length() != n)
01759         {
01760           (*current_liboctave_error_handler)
01761             ("eigs: permB vector invalid");
01762           return -1;
01763         }
01764       else
01765         {
01766           Array<bool> checked (dim_vector (n, 1), false);
01767           for (octave_idx_type i = 0; i < n; i++)
01768             {
01769               octave_idx_type bidx =
01770                 static_cast<octave_idx_type> (permB(i));
01771               if (checked(bidx) || bidx < 0 ||
01772                   bidx >= n || D_NINT (bidx) != bidx)
01773                 {
01774                   (*current_liboctave_error_handler)
01775                     ("eigs: permB vector invalid");
01776                   return -1;
01777                 }
01778             }
01779         }
01780     }
01781 
01782   if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" &&
01783       typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" &&
01784       typ != "SI")
01785     {
01786       (*current_liboctave_error_handler)
01787         ("eigs: unrecognized sigma value");
01788       return -1;
01789     }
01790 
01791   if (typ == "LA" || typ == "SA" || typ == "BE")
01792     {
01793       (*current_liboctave_error_handler)
01794         ("eigs: invalid sigma value for unsymmetric problem");
01795       return -1;
01796     }
01797 
01798   if (have_b)
01799     {
01800       // See Note 3 dsaupd
01801       note3 = true;
01802       if (cholB)
01803         {
01804           bt = b;
01805           b = b.transpose();
01806           if (permB.length() == 0)
01807             {
01808               permB = ColumnVector(n);
01809               for (octave_idx_type i = 0; i < n; i++)
01810                 permB(i) = i;
01811             }
01812         }
01813       else
01814         {
01815           if (! make_cholb(b, bt, permB))
01816             {
01817               (*current_liboctave_error_handler)
01818                 ("eigs: The matrix B is not positive definite");
01819               return -1;
01820             }
01821         }
01822     }
01823 
01824   Array<octave_idx_type> ip (dim_vector (11, 1));
01825   octave_idx_type *iparam = ip.fortran_vec ();
01826 
01827   ip(0) = 1; //ishift
01828   ip(1) = 0;   // ip(1) not referenced
01829   ip(2) = maxit; // mxiter, maximum number of iterations
01830   ip(3) = 1; // NB blocksize in recurrence
01831   ip(4) = 0; // nconv, number of Ritz values that satisfy convergence
01832   ip(5) = 0; //ip(5) not referenced
01833   ip(6) = mode; // mode
01834   ip(7) = 0;
01835   ip(8) = 0;
01836   ip(9) = 0;
01837   ip(10) = 0;
01838   // ip(7) to ip(10) return values
01839 
01840   Array<octave_idx_type> iptr (dim_vector (14, 1));
01841   octave_idx_type *ipntr = iptr.fortran_vec ();
01842 
01843   octave_idx_type ido = 0;
01844   int iter = 0;
01845   octave_idx_type lwork = 3 * p * (p + 2);
01846 
01847   OCTAVE_LOCAL_BUFFER (double, v, n * (p + 1));
01848   OCTAVE_LOCAL_BUFFER (double, workl, lwork + 1);
01849   OCTAVE_LOCAL_BUFFER (double, workd, 3 * n + 1);
01850   double *presid = resid.fortran_vec ();
01851 
01852   do
01853     {
01854       F77_FUNC (dnaupd, DNAUPD)
01855         (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
01856          F77_CONST_CHAR_ARG2 ((typ.c_str()), 2),
01857          k, tol, presid, p, v, n, iparam,
01858          ipntr, workd, workl, lwork, info
01859          F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
01860 
01861       if (f77_exception_encountered)
01862         {
01863           (*current_liboctave_error_handler)
01864             ("eigs: unrecoverable exception encountered in dnaupd");
01865           return -1;
01866         }
01867 
01868       if (disp > 0 && !xisnan(workl[iptr(5)-1]))
01869         {
01870           if (iter++)
01871             {
01872               os << "Iteration " << iter - 1 <<
01873                 ": a few Ritz values of the " << p << "-by-" <<
01874                 p << " matrix\n";
01875               for (int i = 0 ; i < k; i++)
01876                 os << "    " << workl[iptr(5)+i-1] << "\n";
01877             }
01878 
01879           // This is a kludge, as ARPACK doesn't give its
01880           // iteration pointer. But as workl[iptr(5)-1] is
01881           // an output value updated at each iteration, setting
01882           // a value in this array to NaN and testing for it
01883           // is a way of obtaining the iteration counter.
01884           if (ido != 99)
01885             workl[iptr(5)-1] = octave_NaN;
01886         }
01887 
01888       if (ido == -1 || ido == 1 || ido == 2)
01889         {
01890           if (have_b)
01891             {
01892               Matrix mtmp (n,1);
01893               for (octave_idx_type i = 0; i < n; i++)
01894                 mtmp(i,0) = workd[i + iptr(0) - 1];
01895 
01896               mtmp = utsolve(bt, permB, m * ltsolve(b, permB, mtmp));
01897 
01898               for (octave_idx_type i = 0; i < n; i++)
01899                 workd[i+iptr(1)-1] = mtmp(i,0);
01900             }
01901           else if (!vector_product (m, workd + iptr(0) - 1,
01902                                     workd + iptr(1) - 1))
01903             break;
01904         }
01905       else
01906         {
01907           if (info < 0)
01908             {
01909               (*current_liboctave_error_handler)
01910                 ("eigs: error %d in dnaupd", info);
01911               return -1;
01912             }
01913           break;
01914         }
01915     }
01916   while (1);
01917 
01918   octave_idx_type info2;
01919 
01920   // We have a problem in that the size of the C++ bool
01921   // type relative to the fortran logical type. It appears
01922   // that fortran uses 4- or 8-bytes per logical and C++ 1-byte
01923   // per bool, though this might be system dependent. As
01924   // long as the HOWMNY arg is not "S", the logical array
01925   // is just workspace for ARPACK, so use int type to
01926   // avoid problems.
01927   Array<octave_idx_type> s (dim_vector (p, 1));
01928   octave_idx_type *sel = s.fortran_vec ();
01929 
01930   // FIXME -- initialize eig_vec2 to zero; apparently dneupd can skip
01931   // the assignment to elements of Z that represent imaginary parts.
01932   // Found with valgrind and
01933   //
01934   //   A = [1,0,0,-1;0,1,0,0;0,0,1,0;0,0,2,1];
01935   //   [vecs, vals, f] = eigs (A, 1)
01936 
01937   Matrix eig_vec2 (n, k + 1, 0.0);
01938   double *z = eig_vec2.fortran_vec ();
01939 
01940   OCTAVE_LOCAL_BUFFER (double, dr, k + 1);
01941   OCTAVE_LOCAL_BUFFER (double, di, k + 1);
01942   OCTAVE_LOCAL_BUFFER (double, workev, 3 * p);
01943   for (octave_idx_type i = 0; i < k+1; i++)
01944     dr[i] = di[i] = 0.;
01945 
01946   F77_FUNC (dneupd, DNEUPD)
01947     (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, dr, di, z, n, sigmar,
01948      sigmai, workev,  F77_CONST_CHAR_ARG2 (&bmat, 1), n,
01949      F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), k, tol, presid, p, v, n, iparam,
01950      ipntr, workd, workl, lwork, info2 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1)
01951      F77_CHAR_ARG_LEN(2));
01952 
01953   if (f77_exception_encountered)
01954     {
01955       (*current_liboctave_error_handler)
01956         ("eigs: unrecoverable exception encountered in dneupd");
01957       return -1;
01958     }
01959   else
01960     {
01961       eig_val.resize (k+1);
01962       Complex *d = eig_val.fortran_vec ();
01963 
01964       if (info2 == 0)
01965         {
01966           octave_idx_type jj = 0;
01967           for (octave_idx_type i = 0; i < k+1; i++)
01968             {
01969               if (dr[i] == 0.0 && di[i] == 0.0 && jj == 0)
01970                 jj++;
01971               else
01972                 d [i-jj] = Complex (dr[i], di[i]);
01973             }
01974           if (jj == 0 && !rvec)
01975             for (octave_idx_type i = 0; i < k; i++)
01976               d[i] = d[i+1];
01977 
01978           octave_idx_type k2 = k / 2;
01979           for (octave_idx_type i = 0; i < k2; i++)
01980             {
01981               Complex dtmp = d[i];
01982               d[i] = d[k - i - 1];
01983               d[k - i - 1] = dtmp;
01984             }
01985           eig_val.resize(k);
01986 
01987           if (rvec)
01988             {
01989               OCTAVE_LOCAL_BUFFER (double, dtmp, n);
01990 
01991               for (octave_idx_type i = 0; i < k2; i++)
01992                 {
01993                   octave_idx_type off1 = i * n;
01994                   octave_idx_type off2 = (k - i - 1) * n;
01995 
01996                   if (off1 == off2)
01997                     continue;
01998 
01999                   for (octave_idx_type j = 0; j < n; j++)
02000                     dtmp[j] = z[off1 + j];
02001 
02002                   for (octave_idx_type j = 0; j < n; j++)
02003                     z[off1 + j] = z[off2 + j];
02004 
02005                   for (octave_idx_type j = 0; j < n; j++)
02006                     z[off2 + j] = dtmp[j];
02007                 }
02008 
02009               eig_vec.resize (n, k);
02010               octave_idx_type i = 0;
02011               while (i < k)
02012                 {
02013                   octave_idx_type off1 = i * n;
02014                   octave_idx_type off2 = (i+1) * n;
02015                   if (std::imag(eig_val(i)) == 0)
02016                     {
02017                       for (octave_idx_type j = 0; j < n; j++)
02018                         eig_vec(j,i) =
02019                           Complex(z[j+off1],0.);
02020                       i++;
02021                     }
02022                   else
02023                     {
02024                       for (octave_idx_type j = 0; j < n; j++)
02025                         {
02026                           eig_vec(j,i) =
02027                             Complex(z[j+off1],z[j+off2]);
02028                           if (i < k - 1)
02029                             eig_vec(j,i+1) =
02030                               Complex(z[j+off1],-z[j+off2]);
02031                         }
02032                       i+=2;
02033                     }
02034                 }
02035 
02036               if (note3)
02037                 eig_vec = ltsolve(M (b), permB, eig_vec);
02038             }
02039         }
02040       else
02041         {
02042           (*current_liboctave_error_handler)
02043             ("eigs: error %d in dneupd", info2);
02044           return -1;
02045         }
02046     }
02047 
02048   return ip(4);
02049 }
02050 
02051 template <class M>
02052 octave_idx_type
02053 EigsRealNonSymmetricMatrixShift (const M& m, double sigmar,
02054                                  octave_idx_type k, octave_idx_type p,
02055                                  octave_idx_type &info,
02056                                  ComplexMatrix &eig_vec,
02057                                  ComplexColumnVector &eig_val, const M& _b,
02058                                  ColumnVector &permB, ColumnVector &resid,
02059                                  std::ostream& os, double tol, bool rvec,
02060                                  bool cholB, int disp, int maxit)
02061 {
02062   M b(_b);
02063   octave_idx_type n = m.cols ();
02064   octave_idx_type mode = 3;
02065   bool have_b = ! b.is_empty();
02066   std::string typ = "LM";
02067   double sigmai = 0.;
02068 
02069   if (m.rows() != m.cols())
02070     {
02071       (*current_liboctave_error_handler) ("eigs: A must be square");
02072       return -1;
02073     }
02074   if (have_b && (m.rows() != b.rows() || m.rows() != b.cols()))
02075     {
02076       (*current_liboctave_error_handler)
02077         ("eigs: B must be square and the same size as A");
02078       return -1;
02079     }
02080 
02081   // FIXME: The "SM" type for mode 1 seems unstable though faster!!
02082   //if (! std::abs (sigmar))
02083   //  return EigsRealNonSymmetricMatrix (m, "SM", k, p, info, eig_vec, eig_val,
02084   //                                   _b, permB, resid, os, tol, rvec, cholB,
02085   //                                   disp, maxit);
02086 
02087   if (resid.is_empty())
02088     {
02089       std::string rand_dist = octave_rand::distribution();
02090       octave_rand::distribution("uniform");
02091       resid = ColumnVector (octave_rand::vector(n));
02092       octave_rand::distribution(rand_dist);
02093     }
02094 
02095   if (n < 3)
02096     {
02097       (*current_liboctave_error_handler)
02098         ("eigs: n must be at least 3");
02099       return -1;
02100     }
02101 
02102   if (p < 0)
02103     {
02104       p = k * 2 + 1;
02105 
02106       if (p < 20)
02107         p = 20;
02108 
02109       if (p > n - 1)
02110         p = n - 1 ;
02111     }
02112 
02113   if (k <= 0 || k >= n - 1)
02114     {
02115       (*current_liboctave_error_handler)
02116         ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1).\n"
02117              "      Use 'eig(full(A))' instead");
02118       return -1;
02119     }
02120 
02121   if (p <= k || p >= n)
02122     {
02123       (*current_liboctave_error_handler)
02124         ("eigs: opts.p must be greater than k and less than n");
02125       return -1;
02126     }
02127 
02128   if (have_b && cholB && permB.length() != 0)
02129     {
02130       // Check that we really have a permutation vector
02131       if (permB.length() != n)
02132         {
02133           (*current_liboctave_error_handler) ("eigs: permB vector invalid");
02134           return -1;
02135         }
02136       else
02137         {
02138           Array<bool> checked (dim_vector (n, 1), false);
02139           for (octave_idx_type i = 0; i < n; i++)
02140             {
02141               octave_idx_type bidx =
02142                 static_cast<octave_idx_type> (permB(i));
02143               if (checked(bidx) || bidx < 0 ||
02144                   bidx >= n || D_NINT (bidx) != bidx)
02145                 {
02146                   (*current_liboctave_error_handler)
02147                     ("eigs: permB vector invalid");
02148                   return -1;
02149                 }
02150             }
02151         }
02152     }
02153 
02154   char bmat = 'I';
02155   if (have_b)
02156     bmat = 'G';
02157 
02158   Array<octave_idx_type> ip (dim_vector (11, 1));
02159   octave_idx_type *iparam = ip.fortran_vec ();
02160 
02161   ip(0) = 1; //ishift
02162   ip(1) = 0;   // ip(1) not referenced
02163   ip(2) = maxit; // mxiter, maximum number of iterations
02164   ip(3) = 1; // NB blocksize in recurrence
02165   ip(4) = 0; // nconv, number of Ritz values that satisfy convergence
02166   ip(5) = 0; //ip(5) not referenced
02167   ip(6) = mode; // mode
02168   ip(7) = 0;
02169   ip(8) = 0;
02170   ip(9) = 0;
02171   ip(10) = 0;
02172   // ip(7) to ip(10) return values
02173 
02174   Array<octave_idx_type> iptr (dim_vector (14, 1));
02175   octave_idx_type *ipntr = iptr.fortran_vec ();
02176 
02177   octave_idx_type ido = 0;
02178   int iter = 0;
02179   M L, U;
02180 
02181   OCTAVE_LOCAL_BUFFER (octave_idx_type, P, (have_b ? b.rows() : m.rows()));
02182   OCTAVE_LOCAL_BUFFER (octave_idx_type, Q, (have_b ? b.cols() : m.cols()));
02183 
02184   if (! LuAminusSigmaB(m, b, cholB, permB, sigmar, L, U, P, Q))
02185     return -1;
02186 
02187   octave_idx_type lwork = 3 * p * (p + 2);
02188 
02189   OCTAVE_LOCAL_BUFFER (double, v, n * (p + 1));
02190   OCTAVE_LOCAL_BUFFER (double, workl, lwork + 1);
02191   OCTAVE_LOCAL_BUFFER (double, workd, 3 * n + 1);
02192   double *presid = resid.fortran_vec ();
02193 
02194   do
02195     {
02196       F77_FUNC (dnaupd, DNAUPD)
02197         (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
02198          F77_CONST_CHAR_ARG2 ((typ.c_str()), 2),
02199          k, tol, presid, p, v, n, iparam,
02200          ipntr, workd, workl, lwork, info
02201          F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
02202 
02203       if (f77_exception_encountered)
02204         {
02205           (*current_liboctave_error_handler)
02206             ("eigs: unrecoverable exception encountered in dsaupd");
02207           return -1;
02208         }
02209 
02210       if (disp > 0 && !xisnan(workl[iptr(5)-1]))
02211         {
02212           if (iter++)
02213             {
02214               os << "Iteration " << iter - 1 <<
02215                 ": a few Ritz values of the " << p << "-by-" <<
02216                 p << " matrix\n";
02217               for (int i = 0 ; i < k; i++)
02218                 os << "    " << workl[iptr(5)+i-1] << "\n";
02219             }
02220 
02221           // This is a kludge, as ARPACK doesn't give its
02222           // iteration pointer. But as workl[iptr(5)-1] is
02223           // an output value updated at each iteration, setting
02224           // a value in this array to NaN and testing for it
02225           // is a way of obtaining the iteration counter.
02226           if (ido != 99)
02227             workl[iptr(5)-1] = octave_NaN;
02228         }
02229 
02230       if (ido == -1 || ido == 1 || ido == 2)
02231         {
02232           if (have_b)
02233             {
02234               if (ido == -1)
02235                 {
02236                   OCTAVE_LOCAL_BUFFER (double, dtmp, n);
02237 
02238                   vector_product (m, workd+iptr(0)-1, dtmp);
02239 
02240                   Matrix tmp(n, 1);
02241 
02242                   for (octave_idx_type i = 0; i < n; i++)
02243                     tmp(i,0) = dtmp[P[i]];
02244 
02245                   lusolve (L, U, tmp);
02246 
02247                   double *ip2 = workd+iptr(1)-1;
02248                   for (octave_idx_type i = 0; i < n; i++)
02249                     ip2[Q[i]] = tmp(i,0);
02250                 }
02251               else if (ido == 2)
02252                 vector_product (b, workd+iptr(0)-1, workd+iptr(1)-1);
02253               else
02254                 {
02255                   double *ip2 = workd+iptr(2)-1;
02256                   Matrix tmp(n, 1);
02257 
02258                   for (octave_idx_type i = 0; i < n; i++)
02259                     tmp(i,0) = ip2[P[i]];
02260 
02261                   lusolve (L, U, tmp);
02262 
02263                   ip2 = workd+iptr(1)-1;
02264                   for (octave_idx_type i = 0; i < n; i++)
02265                     ip2[Q[i]] = tmp(i,0);
02266                 }
02267             }
02268           else
02269             {
02270               if (ido == 2)
02271                 {
02272                   for (octave_idx_type i = 0; i < n; i++)
02273                     workd[iptr(0) + i - 1] = workd[iptr(1) + i - 1];
02274                 }
02275               else
02276                 {
02277                   double *ip2 = workd+iptr(0)-1;
02278                   Matrix tmp(n, 1);
02279 
02280                   for (octave_idx_type i = 0; i < n; i++)
02281                     tmp(i,0) = ip2[P[i]];
02282 
02283                   lusolve (L, U, tmp);
02284 
02285                   ip2 = workd+iptr(1)-1;
02286                   for (octave_idx_type i = 0; i < n; i++)
02287                     ip2[Q[i]] = tmp(i,0);
02288                 }
02289             }
02290         }
02291       else
02292         {
02293           if (info < 0)
02294             {
02295               (*current_liboctave_error_handler)
02296                 ("eigs: error %d in dsaupd", info);
02297               return -1;
02298             }
02299           break;
02300         }
02301     }
02302   while (1);
02303 
02304   octave_idx_type info2;
02305 
02306   // We have a problem in that the size of the C++ bool
02307   // type relative to the fortran logical type. It appears
02308   // that fortran uses 4- or 8-bytes per logical and C++ 1-byte
02309   // per bool, though this might be system dependent. As
02310   // long as the HOWMNY arg is not "S", the logical array
02311   // is just workspace for ARPACK, so use int type to
02312   // avoid problems.
02313   Array<octave_idx_type> s (dim_vector (p, 1));
02314   octave_idx_type *sel = s.fortran_vec ();
02315 
02316   // FIXME -- initialize eig_vec2 to zero; apparently dneupd can skip
02317   // the assignment to elements of Z that represent imaginary parts.
02318   // Found with valgrind and
02319   //
02320   //   A = [1,0,0,-1;0,1,0,0;0,0,1,0;0,0,2,1];
02321   //   [vecs, vals, f] = eigs (A, 1)
02322 
02323   Matrix eig_vec2 (n, k + 1, 0.0);
02324   double *z = eig_vec2.fortran_vec ();
02325 
02326   OCTAVE_LOCAL_BUFFER (double, dr, k + 1);
02327   OCTAVE_LOCAL_BUFFER (double, di, k + 1);
02328   OCTAVE_LOCAL_BUFFER (double, workev, 3 * p);
02329   for (octave_idx_type i = 0; i < k+1; i++)
02330     dr[i] = di[i] = 0.;
02331 
02332   F77_FUNC (dneupd, DNEUPD)
02333     (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, dr, di, z, n, sigmar,
02334      sigmai, workev,  F77_CONST_CHAR_ARG2 (&bmat, 1), n,
02335      F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), k, tol, presid, p, v, n, iparam,
02336      ipntr, workd, workl, lwork, info2 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1)
02337      F77_CHAR_ARG_LEN(2));
02338 
02339   if (f77_exception_encountered)
02340     {
02341       (*current_liboctave_error_handler)
02342         ("eigs: unrecoverable exception encountered in dneupd");
02343       return -1;
02344     }
02345   else
02346     {
02347       eig_val.resize (k+1);
02348       Complex *d = eig_val.fortran_vec ();
02349 
02350       if (info2 == 0)
02351         {
02352           octave_idx_type jj = 0;
02353           for (octave_idx_type i = 0; i < k+1; i++)
02354             {
02355               if (dr[i] == 0.0 && di[i] == 0.0 && jj == 0)
02356                 jj++;
02357               else
02358                 d [i-jj] = Complex (dr[i], di[i]);
02359             }
02360           if (jj == 0 && !rvec)
02361             for (octave_idx_type i = 0; i < k; i++)
02362               d[i] = d[i+1];
02363 
02364           octave_idx_type k2 = k / 2;
02365           for (octave_idx_type i = 0; i < k2; i++)
02366             {
02367               Complex dtmp = d[i];
02368               d[i] = d[k - i - 1];
02369               d[k - i - 1] = dtmp;
02370             }
02371           eig_val.resize(k);
02372 
02373           if (rvec)
02374             {
02375               OCTAVE_LOCAL_BUFFER (double, dtmp, n);
02376 
02377               for (octave_idx_type i = 0; i < k2; i++)
02378                 {
02379                   octave_idx_type off1 = i * n;
02380                   octave_idx_type off2 = (k - i - 1) * n;
02381 
02382                   if (off1 == off2)
02383                     continue;
02384 
02385                   for (octave_idx_type j = 0; j < n; j++)
02386                     dtmp[j] = z[off1 + j];
02387 
02388                   for (octave_idx_type j = 0; j < n; j++)
02389                     z[off1 + j] = z[off2 + j];
02390 
02391                   for (octave_idx_type j = 0; j < n; j++)
02392                     z[off2 + j] = dtmp[j];
02393                 }
02394 
02395               eig_vec.resize (n, k);
02396               octave_idx_type i = 0;
02397               while (i < k)
02398                 {
02399                   octave_idx_type off1 = i * n;
02400                   octave_idx_type off2 = (i+1) * n;
02401                   if (std::imag(eig_val(i)) == 0)
02402                     {
02403                       for (octave_idx_type j = 0; j < n; j++)
02404                         eig_vec(j,i) =
02405                           Complex(z[j+off1],0.);
02406                       i++;
02407                     }
02408                   else
02409                     {
02410                       for (octave_idx_type j = 0; j < n; j++)
02411                         {
02412                           eig_vec(j,i) =
02413                             Complex(z[j+off1],z[j+off2]);
02414                           if (i < k - 1)
02415                             eig_vec(j,i+1) =
02416                               Complex(z[j+off1],-z[j+off2]);
02417                         }
02418                       i+=2;
02419                     }
02420                 }
02421             }
02422         }
02423       else
02424         {
02425           (*current_liboctave_error_handler)
02426             ("eigs: error %d in dneupd", info2);
02427           return -1;
02428         }
02429     }
02430 
02431   return ip(4);
02432 }
02433 
02434 octave_idx_type
02435 EigsRealNonSymmetricFunc (EigsFunc fun, octave_idx_type n,
02436                           const std::string &_typ, double sigmar,
02437                           octave_idx_type k, octave_idx_type p,
02438                           octave_idx_type &info, ComplexMatrix &eig_vec,
02439                           ComplexColumnVector &eig_val, ColumnVector &resid,
02440                           std::ostream& os, double tol, bool rvec,
02441                           bool /* cholB */, int disp, int maxit)
02442 {
02443   std::string typ (_typ);
02444   bool have_sigma = (sigmar ? true : false);
02445   char bmat = 'I';
02446   double sigmai = 0.;
02447   octave_idx_type mode = 1;
02448   int err = 0;
02449 
02450   if (resid.is_empty())
02451     {
02452       std::string rand_dist = octave_rand::distribution();
02453       octave_rand::distribution("uniform");
02454       resid = ColumnVector (octave_rand::vector(n));
02455       octave_rand::distribution(rand_dist);
02456     }
02457 
02458   if (n < 3)
02459     {
02460       (*current_liboctave_error_handler)
02461         ("eigs: n must be at least 3");
02462       return -1;
02463     }
02464 
02465   if (p < 0)
02466     {
02467       p = k * 2 + 1;
02468 
02469       if (p < 20)
02470         p = 20;
02471 
02472       if (p > n - 1)
02473         p = n - 1 ;
02474     }
02475 
02476   if (k <= 0 || k >= n - 1)
02477     {
02478       (*current_liboctave_error_handler)
02479         ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1).\n"
02480              "      Use 'eig(full(A))' instead");
02481       return -1;
02482     }
02483 
02484   if (p <= k || p >= n)
02485     {
02486       (*current_liboctave_error_handler)
02487         ("eigs: opts.p must be greater than k and less than n");
02488       return -1;
02489     }
02490 
02491 
02492   if (! have_sigma)
02493     {
02494       if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" &&
02495           typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" &&
02496           typ != "SI")
02497         (*current_liboctave_error_handler)
02498           ("eigs: unrecognized sigma value");
02499 
02500       if (typ == "LA" || typ == "SA" || typ == "BE")
02501         {
02502           (*current_liboctave_error_handler)
02503             ("eigs: invalid sigma value for unsymmetric problem");
02504           return -1;
02505         }
02506 
02507       if (typ == "SM")
02508         {
02509           typ = "LM";
02510           sigmar = 0.;
02511           mode = 3;
02512         }
02513     }
02514   else if (! std::abs (sigmar))
02515     typ = "SM";
02516   else
02517     {
02518       typ = "LM";
02519       mode = 3;
02520     }
02521 
02522   Array<octave_idx_type> ip (dim_vector (11, 1));
02523   octave_idx_type *iparam = ip.fortran_vec ();
02524 
02525   ip(0) = 1; //ishift
02526   ip(1) = 0;   // ip(1) not referenced
02527   ip(2) = maxit; // mxiter, maximum number of iterations
02528   ip(3) = 1; // NB blocksize in recurrence
02529   ip(4) = 0; // nconv, number of Ritz values that satisfy convergence
02530   ip(5) = 0; //ip(5) not referenced
02531   ip(6) = mode; // mode
02532   ip(7) = 0;
02533   ip(8) = 0;
02534   ip(9) = 0;
02535   ip(10) = 0;
02536   // ip(7) to ip(10) return values
02537 
02538   Array<octave_idx_type> iptr (dim_vector (14, 1));
02539   octave_idx_type *ipntr = iptr.fortran_vec ();
02540 
02541   octave_idx_type ido = 0;
02542   int iter = 0;
02543   octave_idx_type lwork = 3 * p * (p + 2);
02544 
02545   OCTAVE_LOCAL_BUFFER (double, v, n * (p + 1));
02546   OCTAVE_LOCAL_BUFFER (double, workl, lwork + 1);
02547   OCTAVE_LOCAL_BUFFER (double, workd, 3 * n + 1);
02548   double *presid = resid.fortran_vec ();
02549 
02550   do
02551     {
02552       F77_FUNC (dnaupd, DNAUPD)
02553         (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
02554          F77_CONST_CHAR_ARG2 ((typ.c_str()), 2),
02555          k, tol, presid, p, v, n, iparam,
02556          ipntr, workd, workl, lwork, info
02557          F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
02558 
02559       if (f77_exception_encountered)
02560         {
02561           (*current_liboctave_error_handler)
02562             ("eigs: unrecoverable exception encountered in dnaupd");
02563           return -1;
02564         }
02565 
02566       if (disp > 0 && !xisnan(workl[iptr(5)-1]))
02567         {
02568           if (iter++)
02569             {
02570               os << "Iteration " << iter - 1 <<
02571                 ": a few Ritz values of the " << p << "-by-" <<
02572                 p << " matrix\n";
02573               for (int i = 0 ; i < k; i++)
02574                 os << "    " << workl[iptr(5)+i-1] << "\n";
02575             }
02576 
02577           // This is a kludge, as ARPACK doesn't give its
02578           // iteration pointer. But as workl[iptr(5)-1] is
02579           // an output value updated at each iteration, setting
02580           // a value in this array to NaN and testing for it
02581           // is a way of obtaining the iteration counter.
02582           if (ido != 99)
02583             workl[iptr(5)-1] = octave_NaN;
02584         }
02585 
02586       if (ido == -1 || ido == 1 || ido == 2)
02587         {
02588           double *ip2 = workd + iptr(0) - 1;
02589           ColumnVector x(n);
02590 
02591           for (octave_idx_type i = 0; i < n; i++)
02592             x(i) = *ip2++;
02593 
02594           ColumnVector y = fun (x, err);
02595 
02596           if (err)
02597             return false;
02598 
02599           ip2 = workd + iptr(1) - 1;
02600           for (octave_idx_type i = 0; i < n; i++)
02601             *ip2++ = y(i);
02602         }
02603       else
02604         {
02605           if (info < 0)
02606             {
02607               (*current_liboctave_error_handler)
02608                 ("eigs: error %d in dsaupd", info);
02609               return -1;
02610             }
02611           break;
02612         }
02613     }
02614   while (1);
02615 
02616   octave_idx_type info2;
02617 
02618   // We have a problem in that the size of the C++ bool
02619   // type relative to the fortran logical type. It appears
02620   // that fortran uses 4- or 8-bytes per logical and C++ 1-byte
02621   // per bool, though this might be system dependent. As
02622   // long as the HOWMNY arg is not "S", the logical array
02623   // is just workspace for ARPACK, so use int type to
02624   // avoid problems.
02625   Array<octave_idx_type> s (dim_vector (p, 1));
02626   octave_idx_type *sel = s.fortran_vec ();
02627 
02628   // FIXME -- initialize eig_vec2 to zero; apparently dneupd can skip
02629   // the assignment to elements of Z that represent imaginary parts.
02630   // Found with valgrind and
02631   //
02632   //   A = [1,0,0,-1;0,1,0,0;0,0,1,0;0,0,2,1];
02633   //   [vecs, vals, f] = eigs (A, 1)
02634 
02635   Matrix eig_vec2 (n, k + 1, 0.0);
02636   double *z = eig_vec2.fortran_vec ();
02637 
02638   OCTAVE_LOCAL_BUFFER (double, dr, k + 1);
02639   OCTAVE_LOCAL_BUFFER (double, di, k + 1);
02640   OCTAVE_LOCAL_BUFFER (double, workev, 3 * p);
02641   for (octave_idx_type i = 0; i < k+1; i++)
02642     dr[i] = di[i] = 0.;
02643 
02644   F77_FUNC (dneupd, DNEUPD)
02645     (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, dr, di, z, n, sigmar,
02646      sigmai, workev,  F77_CONST_CHAR_ARG2 (&bmat, 1), n,
02647      F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), k, tol, presid, p, v, n, iparam,
02648      ipntr, workd, workl, lwork, info2 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1)
02649      F77_CHAR_ARG_LEN(2));
02650 
02651   if (f77_exception_encountered)
02652     {
02653       (*current_liboctave_error_handler)
02654         ("eigs: unrecoverable exception encountered in dneupd");
02655       return -1;
02656     }
02657   else
02658     {
02659       eig_val.resize (k+1);
02660       Complex *d = eig_val.fortran_vec ();
02661 
02662       if (info2 == 0)
02663         {
02664           octave_idx_type jj = 0;
02665           for (octave_idx_type i = 0; i < k+1; i++)
02666             {
02667               if (dr[i] == 0.0 && di[i] == 0.0 && jj == 0)
02668                 jj++;
02669               else
02670                 d [i-jj] = Complex (dr[i], di[i]);
02671             }
02672           if (jj == 0 && !rvec)
02673             for (octave_idx_type i = 0; i < k; i++)
02674               d[i] = d[i+1];
02675 
02676           octave_idx_type k2 = k / 2;
02677           for (octave_idx_type i = 0; i < k2; i++)
02678             {
02679               Complex dtmp = d[i];
02680               d[i] = d[k - i - 1];
02681               d[k - i - 1] = dtmp;
02682             }
02683           eig_val.resize(k);
02684 
02685           if (rvec)
02686             {
02687               OCTAVE_LOCAL_BUFFER (double, dtmp, n);
02688 
02689               for (octave_idx_type i = 0; i < k2; i++)
02690                 {
02691                   octave_idx_type off1 = i * n;
02692                   octave_idx_type off2 = (k - i - 1) * n;
02693 
02694                   if (off1 == off2)
02695                     continue;
02696 
02697                   for (octave_idx_type j = 0; j < n; j++)
02698                     dtmp[j] = z[off1 + j];
02699 
02700                   for (octave_idx_type j = 0; j < n; j++)
02701                     z[off1 + j] = z[off2 + j];
02702 
02703                   for (octave_idx_type j = 0; j < n; j++)
02704                     z[off2 + j] = dtmp[j];
02705                 }
02706 
02707               eig_vec.resize (n, k);
02708               octave_idx_type i = 0;
02709               while (i < k)
02710                 {
02711                   octave_idx_type off1 = i * n;
02712                   octave_idx_type off2 = (i+1) * n;
02713                   if (std::imag(eig_val(i)) == 0)
02714                     {
02715                       for (octave_idx_type j = 0; j < n; j++)
02716                         eig_vec(j,i) =
02717                           Complex(z[j+off1],0.);
02718                       i++;
02719                     }
02720                   else
02721                     {
02722                       for (octave_idx_type j = 0; j < n; j++)
02723                         {
02724                           eig_vec(j,i) =
02725                             Complex(z[j+off1],z[j+off2]);
02726                           if (i < k - 1)
02727                             eig_vec(j,i+1) =
02728                               Complex(z[j+off1],-z[j+off2]);
02729                         }
02730                       i+=2;
02731                     }
02732                 }
02733             }
02734         }
02735       else
02736         {
02737           (*current_liboctave_error_handler)
02738             ("eigs: error %d in dneupd", info2);
02739           return -1;
02740         }
02741     }
02742 
02743   return ip(4);
02744 }
02745 
02746 template <class M>
02747 octave_idx_type
02748 EigsComplexNonSymmetricMatrix (const M& m, const std::string typ,
02749                                octave_idx_type k, octave_idx_type p,
02750                                octave_idx_type &info, ComplexMatrix &eig_vec,
02751                                ComplexColumnVector &eig_val, const M& _b,
02752                                ColumnVector &permB,
02753                                ComplexColumnVector &cresid,
02754                                std::ostream& os, double tol, bool rvec,
02755                                bool cholB, int disp, int maxit)
02756 {
02757   M b(_b);
02758   octave_idx_type n = m.cols ();
02759   octave_idx_type mode = 1;
02760   bool have_b = ! b.is_empty();
02761   bool note3 = false;
02762   char bmat = 'I';
02763   Complex sigma = 0.;
02764   M bt;
02765 
02766   if (m.rows() != m.cols())
02767     {
02768       (*current_liboctave_error_handler) ("eigs: A must be square");
02769       return -1;
02770     }
02771   if (have_b && (m.rows() != b.rows() || m.rows() != b.cols()))
02772     {
02773       (*current_liboctave_error_handler)
02774         ("eigs: B must be square and the same size as A");
02775       return -1;
02776     }
02777 
02778   if (cresid.is_empty())
02779     {
02780       std::string rand_dist = octave_rand::distribution();
02781       octave_rand::distribution("uniform");
02782       Array<double> rr (octave_rand::vector(n));
02783       Array<double> ri (octave_rand::vector(n));
02784       cresid = ComplexColumnVector (n);
02785       for (octave_idx_type i = 0; i < n; i++)
02786         cresid(i) = Complex(rr(i),ri(i));
02787       octave_rand::distribution(rand_dist);
02788     }
02789 
02790   if (n < 3)
02791     {
02792       (*current_liboctave_error_handler)
02793         ("eigs: n must be at least 3");
02794       return -1;
02795     }
02796 
02797   if (p < 0)
02798     {
02799       p = k * 2 + 1;
02800 
02801       if (p < 20)
02802         p = 20;
02803 
02804       if (p > n - 1)
02805         p = n - 1 ;
02806     }
02807 
02808   if (k <= 0 || k >= n - 1)
02809     {
02810       (*current_liboctave_error_handler)
02811         ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1).\n"
02812          "      Use 'eig(full(A))' instead");
02813       return -1;
02814     }
02815 
02816   if (p <= k || p >= n)
02817     {
02818       (*current_liboctave_error_handler)
02819         ("eigs: opts.p must be greater than k and less than n");
02820       return -1;
02821     }
02822 
02823   if (have_b && cholB && permB.length() != 0)
02824     {
02825       // Check the we really have a permutation vector
02826       if (permB.length() != n)
02827         {
02828           (*current_liboctave_error_handler)
02829             ("eigs: permB vector invalid");
02830           return -1;
02831         }
02832       else
02833         {
02834           Array<bool> checked (dim_vector (n, 1), false);
02835           for (octave_idx_type i = 0; i < n; i++)
02836             {
02837               octave_idx_type bidx =
02838                 static_cast<octave_idx_type> (permB(i));
02839               if (checked(bidx) || bidx < 0 ||
02840                   bidx >= n || D_NINT (bidx) != bidx)
02841                 {
02842                   (*current_liboctave_error_handler)
02843                     ("eigs: permB vector invalid");
02844                   return -1;
02845                 }
02846             }
02847         }
02848     }
02849 
02850   if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" &&
02851       typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" &&
02852       typ != "SI")
02853     {
02854       (*current_liboctave_error_handler)
02855         ("eigs: unrecognized sigma value");
02856       return -1;
02857     }
02858 
02859   if (typ == "LA" || typ == "SA" || typ == "BE")
02860     {
02861       (*current_liboctave_error_handler)
02862         ("eigs: invalid sigma value for complex problem");
02863       return -1;
02864     }
02865 
02866   if (have_b)
02867     {
02868       // See Note 3 dsaupd
02869       note3 = true;
02870       if (cholB)
02871         {
02872           bt = b;
02873           b = b.hermitian();
02874           if (permB.length() == 0)
02875             {
02876               permB = ColumnVector(n);
02877               for (octave_idx_type i = 0; i < n; i++)
02878                 permB(i) = i;
02879             }
02880         }
02881       else
02882         {
02883           if (! make_cholb(b, bt, permB))
02884             {
02885               (*current_liboctave_error_handler)
02886                 ("eigs: The matrix B is not positive definite");
02887               return -1;
02888             }
02889         }
02890     }
02891 
02892   Array<octave_idx_type> ip (dim_vector (11, 1));
02893   octave_idx_type *iparam = ip.fortran_vec ();
02894 
02895   ip(0) = 1; //ishift
02896   ip(1) = 0;   // ip(1) not referenced
02897   ip(2) = maxit; // mxiter, maximum number of iterations
02898   ip(3) = 1; // NB blocksize in recurrence
02899   ip(4) = 0; // nconv, number of Ritz values that satisfy convergence
02900   ip(5) = 0; //ip(5) not referenced
02901   ip(6) = mode; // mode
02902   ip(7) = 0;
02903   ip(8) = 0;
02904   ip(9) = 0;
02905   ip(10) = 0;
02906   // ip(7) to ip(10) return values
02907 
02908   Array<octave_idx_type> iptr (dim_vector (14, 1));
02909   octave_idx_type *ipntr = iptr.fortran_vec ();
02910 
02911   octave_idx_type ido = 0;
02912   int iter = 0;
02913   octave_idx_type lwork = p * (3 * p + 5);
02914 
02915   OCTAVE_LOCAL_BUFFER (Complex, v, n * p);
02916   OCTAVE_LOCAL_BUFFER (Complex, workl, lwork);
02917   OCTAVE_LOCAL_BUFFER (Complex, workd, 3 * n);
02918   OCTAVE_LOCAL_BUFFER (double, rwork, p);
02919   Complex *presid = cresid.fortran_vec ();
02920 
02921   do
02922     {
02923       F77_FUNC (znaupd, ZNAUPD)
02924         (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
02925          F77_CONST_CHAR_ARG2 ((typ.c_str()), 2),
02926          k, tol, presid, p, v, n, iparam,
02927          ipntr, workd, workl, lwork, rwork, info
02928          F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
02929 
02930       if (f77_exception_encountered)
02931         {
02932           (*current_liboctave_error_handler)
02933             ("eigs: unrecoverable exception encountered in znaupd");
02934           return -1;
02935         }
02936 
02937       if (disp > 0 && !xisnan(workl[iptr(5)-1]))
02938         {
02939           if (iter++)
02940             {
02941               os << "Iteration " << iter - 1 <<
02942                 ": a few Ritz values of the " << p << "-by-" <<
02943                 p << " matrix\n";
02944               for (int i = 0 ; i < k; i++)
02945                 os << "    " << workl[iptr(5)+i-1] << "\n";
02946             }
02947 
02948           // This is a kludge, as ARPACK doesn't give its
02949           // iteration pointer. But as workl[iptr(5)-1] is
02950           // an output value updated at each iteration, setting
02951           // a value in this array to NaN and testing for it
02952           // is a way of obtaining the iteration counter.
02953           if (ido != 99)
02954             workl[iptr(5)-1] = octave_NaN;
02955         }
02956 
02957       if (ido == -1 || ido == 1 || ido == 2)
02958         {
02959           if (have_b)
02960             {
02961               ComplexMatrix mtmp (n,1);
02962               for (octave_idx_type i = 0; i < n; i++)
02963                 mtmp(i,0) = workd[i + iptr(0) - 1];
02964               mtmp = utsolve(bt, permB, m * ltsolve(b, permB, mtmp));
02965               for (octave_idx_type i = 0; i < n; i++)
02966                 workd[i+iptr(1)-1] = mtmp(i,0);
02967 
02968             }
02969           else if (!vector_product (m, workd + iptr(0) - 1,
02970                                     workd + iptr(1) - 1))
02971             break;
02972         }
02973       else
02974         {
02975           if (info < 0)
02976             {
02977               (*current_liboctave_error_handler)
02978                 ("eigs: error %d in znaupd", info);
02979               return -1;
02980             }
02981           break;
02982         }
02983     }
02984   while (1);
02985 
02986   octave_idx_type info2;
02987 
02988   // We have a problem in that the size of the C++ bool
02989   // type relative to the fortran logical type. It appears
02990   // that fortran uses 4- or 8-bytes per logical and C++ 1-byte
02991   // per bool, though this might be system dependent. As
02992   // long as the HOWMNY arg is not "S", the logical array
02993   // is just workspace for ARPACK, so use int type to
02994   // avoid problems.
02995   Array<octave_idx_type> s (dim_vector (p, 1));
02996   octave_idx_type *sel = s.fortran_vec ();
02997 
02998   eig_vec.resize (n, k);
02999   Complex *z = eig_vec.fortran_vec ();
03000 
03001   eig_val.resize (k+1);
03002   Complex *d = eig_val.fortran_vec ();
03003 
03004   OCTAVE_LOCAL_BUFFER (Complex, workev, 2 * p);
03005 
03006   F77_FUNC (zneupd, ZNEUPD)
03007     (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma, workev,
03008      F77_CONST_CHAR_ARG2 (&bmat, 1), n,
03009      F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2),
03010      k, tol, presid, p, v, n, iparam, ipntr, workd, workl, lwork, rwork, info2
03011      F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
03012 
03013   if (f77_exception_encountered)
03014     {
03015       (*current_liboctave_error_handler)
03016         ("eigs: unrecoverable exception encountered in zneupd");
03017       return -1;
03018     }
03019 
03020   if (info2 == 0)
03021     {
03022       octave_idx_type k2 = k / 2;
03023       for (octave_idx_type i = 0; i < k2; i++)
03024         {
03025           Complex ctmp = d[i];
03026           d[i] = d[k - i - 1];
03027           d[k - i - 1] = ctmp;
03028         }
03029       eig_val.resize(k);
03030 
03031       if (rvec)
03032         {
03033           OCTAVE_LOCAL_BUFFER (Complex, ctmp, n);
03034 
03035           for (octave_idx_type i = 0; i < k2; i++)
03036             {
03037               octave_idx_type off1 = i * n;
03038               octave_idx_type off2 = (k - i - 1) * n;
03039 
03040               if (off1 == off2)
03041                 continue;
03042 
03043               for (octave_idx_type j = 0; j < n; j++)
03044                 ctmp[j] = z[off1 + j];
03045 
03046               for (octave_idx_type j = 0; j < n; j++)
03047                 z[off1 + j] = z[off2 + j];
03048 
03049               for (octave_idx_type j = 0; j < n; j++)
03050                 z[off2 + j] = ctmp[j];
03051             }
03052 
03053           if (note3)
03054             eig_vec = ltsolve(b, permB, eig_vec);
03055         }
03056     }
03057   else
03058     {
03059       (*current_liboctave_error_handler)
03060         ("eigs: error %d in zneupd", info2);
03061       return -1;
03062     }
03063 
03064   return ip(4);
03065 }
03066 
03067 template <class M>
03068 octave_idx_type
03069 EigsComplexNonSymmetricMatrixShift (const M& m, Complex sigma,
03070                                     octave_idx_type k, octave_idx_type p,
03071                                     octave_idx_type &info,
03072                                     ComplexMatrix &eig_vec,
03073                                     ComplexColumnVector &eig_val, const M& _b,
03074                                     ColumnVector &permB,
03075                                     ComplexColumnVector &cresid,
03076                                     std::ostream& os, double tol, bool rvec,
03077                                     bool cholB, int disp, int maxit)
03078 {
03079   M b(_b);
03080   octave_idx_type n = m.cols ();
03081   octave_idx_type mode = 3;
03082   bool have_b = ! b.is_empty();
03083   std::string typ = "LM";
03084 
03085   if (m.rows() != m.cols())
03086     {
03087       (*current_liboctave_error_handler) ("eigs: A must be square");
03088       return -1;
03089     }
03090   if (have_b && (m.rows() != b.rows() || m.rows() != b.cols()))
03091     {
03092       (*current_liboctave_error_handler)
03093         ("eigs: B must be square and the same size as A");
03094       return -1;
03095     }
03096 
03097   // FIXME: The "SM" type for mode 1 seems unstable though faster!!
03098   //if (! std::abs (sigma))
03099   //  return EigsComplexNonSymmetricMatrix (m, "SM", k, p, info, eig_vec,
03100   //                                      eig_val, _b, permB, cresid, os, tol,
03101   //                                      rvec, cholB, disp, maxit);
03102 
03103   if (cresid.is_empty())
03104     {
03105       std::string rand_dist = octave_rand::distribution();
03106       octave_rand::distribution("uniform");
03107       Array<double> rr (octave_rand::vector(n));
03108       Array<double> ri (octave_rand::vector(n));
03109       cresid = ComplexColumnVector (n);
03110       for (octave_idx_type i = 0; i < n; i++)
03111         cresid(i) = Complex(rr(i),ri(i));
03112       octave_rand::distribution(rand_dist);
03113     }
03114 
03115   if (n < 3)
03116     {
03117       (*current_liboctave_error_handler)
03118         ("eigs: n must be at least 3");
03119       return -1;
03120     }
03121 
03122   if (p < 0)
03123     {
03124       p = k * 2 + 1;
03125 
03126       if (p < 20)
03127         p = 20;
03128 
03129       if (p > n - 1)
03130         p = n - 1 ;
03131     }
03132 
03133   if (k <= 0 || k >= n - 1)
03134     {
03135       (*current_liboctave_error_handler)
03136         ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1).\n"
03137              "      Use 'eig(full(A))' instead");
03138       return -1;
03139     }
03140 
03141   if (p <= k || p >= n)
03142     {
03143       (*current_liboctave_error_handler)
03144         ("eigs: opts.p must be greater than k and less than n");
03145       return -1;
03146     }
03147 
03148   if (have_b && cholB && permB.length() != 0)
03149     {
03150       // Check that we really have a permutation vector
03151       if (permB.length() != n)
03152         {
03153           (*current_liboctave_error_handler) ("eigs: permB vector invalid");
03154           return -1;
03155         }
03156       else
03157         {
03158           Array<bool> checked (dim_vector (n, 1), false);
03159           for (octave_idx_type i = 0; i < n; i++)
03160             {
03161               octave_idx_type bidx =
03162                 static_cast<octave_idx_type> (permB(i));
03163               if (checked(bidx) || bidx < 0 ||
03164                   bidx >= n || D_NINT (bidx) != bidx)
03165                 {
03166                   (*current_liboctave_error_handler)
03167                     ("eigs: permB vector invalid");
03168                   return -1;
03169                 }
03170             }
03171         }
03172     }
03173 
03174   char bmat = 'I';
03175   if (have_b)
03176     bmat = 'G';
03177 
03178   Array<octave_idx_type> ip (dim_vector (11, 1));
03179   octave_idx_type *iparam = ip.fortran_vec ();
03180 
03181   ip(0) = 1; //ishift
03182   ip(1) = 0;   // ip(1) not referenced
03183   ip(2) = maxit; // mxiter, maximum number of iterations
03184   ip(3) = 1; // NB blocksize in recurrence
03185   ip(4) = 0; // nconv, number of Ritz values that satisfy convergence
03186   ip(5) = 0; //ip(5) not referenced
03187   ip(6) = mode; // mode
03188   ip(7) = 0;
03189   ip(8) = 0;
03190   ip(9) = 0;
03191   ip(10) = 0;
03192   // ip(7) to ip(10) return values
03193 
03194   Array<octave_idx_type> iptr (dim_vector (14, 1));
03195   octave_idx_type *ipntr = iptr.fortran_vec ();
03196 
03197   octave_idx_type ido = 0;
03198   int iter = 0;
03199   M L, U;
03200 
03201   OCTAVE_LOCAL_BUFFER (octave_idx_type, P, (have_b ? b.rows() : m.rows()));
03202   OCTAVE_LOCAL_BUFFER (octave_idx_type, Q, (have_b ? b.cols() : m.cols()));
03203 
03204   if (! LuAminusSigmaB(m, b, cholB, permB, sigma, L, U, P, Q))
03205     return -1;
03206 
03207   octave_idx_type lwork = p * (3 * p + 5);
03208 
03209   OCTAVE_LOCAL_BUFFER (Complex, v, n * p);
03210   OCTAVE_LOCAL_BUFFER (Complex, workl, lwork);
03211   OCTAVE_LOCAL_BUFFER (Complex, workd, 3 * n);
03212   OCTAVE_LOCAL_BUFFER (double, rwork, p);
03213   Complex *presid = cresid.fortran_vec ();
03214 
03215   do
03216     {
03217       F77_FUNC (znaupd, ZNAUPD)
03218         (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
03219          F77_CONST_CHAR_ARG2 ((typ.c_str()), 2),
03220          k, tol, presid, p, v, n, iparam,
03221          ipntr, workd, workl, lwork, rwork, info
03222          F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
03223 
03224       if (f77_exception_encountered)
03225         {
03226           (*current_liboctave_error_handler)
03227             ("eigs: unrecoverable exception encountered in znaupd");
03228           return -1;
03229         }
03230 
03231       if (disp > 0 && !xisnan(workl[iptr(5)-1]))
03232         {
03233           if (iter++)
03234             {
03235               os << "Iteration " << iter - 1 <<
03236                 ": a few Ritz values of the " << p << "-by-" <<
03237                 p << " matrix\n";
03238               for (int i = 0 ; i < k; i++)
03239                 os << "    " << workl[iptr(5)+i-1] << "\n";
03240             }
03241 
03242           // This is a kludge, as ARPACK doesn't give its
03243           // iteration pointer. But as workl[iptr(5)-1] is
03244           // an output value updated at each iteration, setting
03245           // a value in this array to NaN and testing for it
03246           // is a way of obtaining the iteration counter.
03247           if (ido != 99)
03248             workl[iptr(5)-1] = octave_NaN;
03249         }
03250 
03251       if (ido == -1 || ido == 1 || ido == 2)
03252         {
03253           if (have_b)
03254             {
03255               if (ido == -1)
03256                 {
03257                   OCTAVE_LOCAL_BUFFER (Complex, ctmp, n);
03258 
03259                   vector_product (m, workd+iptr(0)-1, ctmp);
03260 
03261                   ComplexMatrix tmp(n, 1);
03262 
03263                   for (octave_idx_type i = 0; i < n; i++)
03264                     tmp(i,0) = ctmp[P[i]];
03265 
03266                   lusolve (L, U, tmp);
03267 
03268                   Complex *ip2 = workd+iptr(1)-1;
03269                   for (octave_idx_type i = 0; i < n; i++)
03270                     ip2[Q[i]] = tmp(i,0);
03271                 }
03272               else if (ido == 2)
03273                 vector_product (b, workd + iptr(0) - 1, workd + iptr(1) - 1);
03274               else
03275                 {
03276                   Complex *ip2 = workd+iptr(2)-1;
03277                   ComplexMatrix tmp(n, 1);
03278 
03279                   for (octave_idx_type i = 0; i < n; i++)
03280                     tmp(i,0) = ip2[P[i]];
03281 
03282                   lusolve (L, U, tmp);
03283 
03284                   ip2 = workd+iptr(1)-1;
03285                   for (octave_idx_type i = 0; i < n; i++)
03286                     ip2[Q[i]] = tmp(i,0);
03287                 }
03288             }
03289           else
03290             {
03291               if (ido == 2)
03292                 {
03293                   for (octave_idx_type i = 0; i < n; i++)
03294                     workd[iptr(0) + i - 1] =
03295                       workd[iptr(1) + i - 1];
03296                 }
03297               else
03298                 {
03299                   Complex *ip2 = workd+iptr(0)-1;
03300                   ComplexMatrix tmp(n, 1);
03301 
03302                   for (octave_idx_type i = 0; i < n; i++)
03303                     tmp(i,0) = ip2[P[i]];
03304 
03305                   lusolve (L, U, tmp);
03306 
03307                   ip2 = workd+iptr(1)-1;
03308                   for (octave_idx_type i = 0; i < n; i++)
03309                     ip2[Q[i]] = tmp(i,0);
03310                 }
03311             }
03312         }
03313       else
03314         {
03315           if (info < 0)
03316             {
03317               (*current_liboctave_error_handler)
03318                 ("eigs: error %d in dsaupd", info);
03319               return -1;
03320             }
03321           break;
03322         }
03323     }
03324   while (1);
03325 
03326   octave_idx_type info2;
03327 
03328   // We have a problem in that the size of the C++ bool
03329   // type relative to the fortran logical type. It appears
03330   // that fortran uses 4- or 8-bytes per logical and C++ 1-byte
03331   // per bool, though this might be system dependent. As
03332   // long as the HOWMNY arg is not "S", the logical array
03333   // is just workspace for ARPACK, so use int type to
03334   // avoid problems.
03335   Array<octave_idx_type> s (dim_vector (p, 1));
03336   octave_idx_type *sel = s.fortran_vec ();
03337 
03338   eig_vec.resize (n, k);
03339   Complex *z = eig_vec.fortran_vec ();
03340 
03341   eig_val.resize (k+1);
03342   Complex *d = eig_val.fortran_vec ();
03343 
03344   OCTAVE_LOCAL_BUFFER (Complex, workev, 2 * p);
03345 
03346   F77_FUNC (zneupd, ZNEUPD)
03347     (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma, workev,
03348      F77_CONST_CHAR_ARG2 (&bmat, 1), n,
03349      F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2),
03350      k, tol, presid, p, v, n, iparam, ipntr, workd, workl, lwork, rwork, info2
03351      F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
03352 
03353   if (f77_exception_encountered)
03354     {
03355       (*current_liboctave_error_handler)
03356         ("eigs: unrecoverable exception encountered in zneupd");
03357       return -1;
03358     }
03359 
03360   if (info2 == 0)
03361     {
03362       octave_idx_type k2 = k / 2;
03363       for (octave_idx_type i = 0; i < k2; i++)
03364         {
03365           Complex ctmp = d[i];
03366           d[i] = d[k - i - 1];
03367           d[k - i - 1] = ctmp;
03368         }
03369       eig_val.resize(k);
03370 
03371       if (rvec)
03372         {
03373           OCTAVE_LOCAL_BUFFER (Complex, ctmp, n);
03374 
03375           for (octave_idx_type i = 0; i < k2; i++)
03376             {
03377               octave_idx_type off1 = i * n;
03378               octave_idx_type off2 = (k - i - 1) * n;
03379 
03380               if (off1 == off2)
03381                 continue;
03382 
03383               for (octave_idx_type j = 0; j < n; j++)
03384                 ctmp[j] = z[off1 + j];
03385 
03386               for (octave_idx_type j = 0; j < n; j++)
03387                 z[off1 + j] = z[off2 + j];
03388 
03389               for (octave_idx_type j = 0; j < n; j++)
03390                 z[off2 + j] = ctmp[j];
03391             }
03392         }
03393     }
03394   else
03395     {
03396       (*current_liboctave_error_handler)
03397         ("eigs: error %d in zneupd", info2);
03398       return -1;
03399     }
03400 
03401   return ip(4);
03402 }
03403 
03404 octave_idx_type
03405 EigsComplexNonSymmetricFunc (EigsComplexFunc fun, octave_idx_type n,
03406                              const std::string &_typ, Complex sigma,
03407                              octave_idx_type k, octave_idx_type p,
03408                              octave_idx_type &info, ComplexMatrix &eig_vec,
03409                              ComplexColumnVector &eig_val,
03410                              ComplexColumnVector &cresid, std::ostream& os,
03411                              double tol, bool rvec, bool /* cholB */,
03412                              int disp, int maxit)
03413 {
03414   std::string typ (_typ);
03415   bool have_sigma = (std::abs(sigma) ? true : false);
03416   char bmat = 'I';
03417   octave_idx_type mode = 1;
03418   int err = 0;
03419 
03420   if (cresid.is_empty())
03421     {
03422       std::string rand_dist = octave_rand::distribution();
03423       octave_rand::distribution("uniform");
03424       Array<double> rr (octave_rand::vector(n));
03425       Array<double> ri (octave_rand::vector(n));
03426       cresid = ComplexColumnVector (n);
03427       for (octave_idx_type i = 0; i < n; i++)
03428         cresid(i) = Complex(rr(i),ri(i));
03429       octave_rand::distribution(rand_dist);
03430     }
03431 
03432   if (n < 3)
03433     {
03434       (*current_liboctave_error_handler)
03435         ("eigs: n must be at least 3");
03436       return -1;
03437     }
03438 
03439   if (p < 0)
03440     {
03441       p = k * 2 + 1;
03442 
03443       if (p < 20)
03444         p = 20;
03445 
03446       if (p > n - 1)
03447         p = n - 1 ;
03448     }
03449 
03450   if (k <= 0 || k >= n - 1)
03451     {
03452       (*current_liboctave_error_handler)
03453         ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1).\n"
03454              "      Use 'eig(full(A))' instead");
03455       return -1;
03456     }
03457 
03458   if (p <= k || p >= n)
03459     {
03460       (*current_liboctave_error_handler)
03461         ("eigs: opts.p must be greater than k and less than n");
03462       return -1;
03463     }
03464 
03465   if (! have_sigma)
03466     {
03467       if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" &&
03468           typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" &&
03469           typ != "SI")
03470         (*current_liboctave_error_handler)
03471           ("eigs: unrecognized sigma value");
03472 
03473       if (typ == "LA" || typ == "SA" || typ == "BE")
03474         {
03475           (*current_liboctave_error_handler)
03476             ("eigs: invalid sigma value for complex problem");
03477           return -1;
03478         }
03479 
03480       if (typ == "SM")
03481         {
03482           typ = "LM";
03483           sigma = 0.;
03484           mode = 3;
03485         }
03486     }
03487   else if (! std::abs (sigma))
03488     typ = "SM";
03489   else
03490     {
03491       typ = "LM";
03492       mode = 3;
03493     }
03494 
03495   Array<octave_idx_type> ip (dim_vector (11, 1));
03496   octave_idx_type *iparam = ip.fortran_vec ();
03497 
03498   ip(0) = 1; //ishift
03499   ip(1) = 0;   // ip(1) not referenced
03500   ip(2) = maxit; // mxiter, maximum number of iterations
03501   ip(3) = 1; // NB blocksize in recurrence
03502   ip(4) = 0; // nconv, number of Ritz values that satisfy convergence
03503   ip(5) = 0; //ip(5) not referenced
03504   ip(6) = mode; // mode
03505   ip(7) = 0;
03506   ip(8) = 0;
03507   ip(9) = 0;
03508   ip(10) = 0;
03509   // ip(7) to ip(10) return values
03510 
03511   Array<octave_idx_type> iptr (dim_vector (14, 1));
03512   octave_idx_type *ipntr = iptr.fortran_vec ();
03513 
03514   octave_idx_type ido = 0;
03515   int iter = 0;
03516   octave_idx_type lwork = p * (3 * p + 5);
03517 
03518   OCTAVE_LOCAL_BUFFER (Complex, v, n * p);
03519   OCTAVE_LOCAL_BUFFER (Complex, workl, lwork);
03520   OCTAVE_LOCAL_BUFFER (Complex, workd, 3 * n);
03521   OCTAVE_LOCAL_BUFFER (double, rwork, p);
03522   Complex *presid = cresid.fortran_vec ();
03523 
03524   do
03525     {
03526       F77_FUNC (znaupd, ZNAUPD)
03527         (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
03528          F77_CONST_CHAR_ARG2 ((typ.c_str()), 2),
03529          k, tol, presid, p, v, n, iparam,
03530          ipntr, workd, workl, lwork, rwork, info
03531          F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
03532 
03533       if (f77_exception_encountered)
03534         {
03535           (*current_liboctave_error_handler)
03536             ("eigs: unrecoverable exception encountered in znaupd");
03537           return -1;
03538         }
03539 
03540       if (disp > 0 && !xisnan(workl[iptr(5)-1]))
03541         {
03542           if (iter++)
03543             {
03544               os << "Iteration " << iter - 1 <<
03545                 ": a few Ritz values of the " << p << "-by-" <<
03546                 p << " matrix\n";
03547               for (int i = 0 ; i < k; i++)
03548                 os << "    " << workl[iptr(5)+i-1] << "\n";
03549             }
03550 
03551           // This is a kludge, as ARPACK doesn't give its
03552           // iteration pointer. But as workl[iptr(5)-1] is
03553           // an output value updated at each iteration, setting
03554           // a value in this array to NaN and testing for it
03555           // is a way of obtaining the iteration counter.
03556           if (ido != 99)
03557             workl[iptr(5)-1] = octave_NaN;
03558         }
03559 
03560       if (ido == -1 || ido == 1 || ido == 2)
03561         {
03562           Complex *ip2 = workd + iptr(0) - 1;
03563           ComplexColumnVector x(n);
03564 
03565           for (octave_idx_type i = 0; i < n; i++)
03566             x(i) = *ip2++;
03567 
03568           ComplexColumnVector y = fun (x, err);
03569 
03570           if (err)
03571             return false;
03572 
03573           ip2 = workd + iptr(1) - 1;
03574           for (octave_idx_type i = 0; i < n; i++)
03575             *ip2++ = y(i);
03576         }
03577       else
03578         {
03579           if (info < 0)
03580             {
03581               (*current_liboctave_error_handler)
03582                 ("eigs: error %d in dsaupd", info);
03583               return -1;
03584             }
03585           break;
03586         }
03587     }
03588   while (1);
03589 
03590   octave_idx_type info2;
03591 
03592   // We have a problem in that the size of the C++ bool
03593   // type relative to the fortran logical type. It appears
03594   // that fortran uses 4- or 8-bytes per logical and C++ 1-byte
03595   // per bool, though this might be system dependent. As
03596   // long as the HOWMNY arg is not "S", the logical array
03597   // is just workspace for ARPACK, so use int type to
03598   // avoid problems.
03599   Array<octave_idx_type> s (dim_vector (p, 1));
03600   octave_idx_type *sel = s.fortran_vec ();
03601 
03602   eig_vec.resize (n, k);
03603   Complex *z = eig_vec.fortran_vec ();
03604 
03605   eig_val.resize (k+1);
03606   Complex *d = eig_val.fortran_vec ();
03607 
03608   OCTAVE_LOCAL_BUFFER (Complex, workev, 2 * p);
03609 
03610   F77_FUNC (zneupd, ZNEUPD)
03611     (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma, workev,
03612      F77_CONST_CHAR_ARG2 (&bmat, 1), n,
03613      F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2),
03614      k, tol, presid, p, v, n, iparam, ipntr, workd, workl, lwork, rwork, info2
03615      F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
03616 
03617   if (f77_exception_encountered)
03618     {
03619       (*current_liboctave_error_handler)
03620         ("eigs: unrecoverable exception encountered in zneupd");
03621       return -1;
03622     }
03623 
03624   if (info2 == 0)
03625     {
03626       octave_idx_type k2 = k / 2;
03627       for (octave_idx_type i = 0; i < k2; i++)
03628         {
03629           Complex ctmp = d[i];
03630           d[i] = d[k - i - 1];
03631           d[k - i - 1] = ctmp;
03632         }
03633       eig_val.resize(k);
03634 
03635       if (rvec)
03636         {
03637           OCTAVE_LOCAL_BUFFER (Complex, ctmp, n);
03638 
03639           for (octave_idx_type i = 0; i < k2; i++)
03640             {
03641               octave_idx_type off1 = i * n;
03642               octave_idx_type off2 = (k - i - 1) * n;
03643 
03644               if (off1 == off2)
03645                 continue;
03646 
03647               for (octave_idx_type j = 0; j < n; j++)
03648                 ctmp[j] = z[off1 + j];
03649 
03650               for (octave_idx_type j = 0; j < n; j++)
03651                 z[off1 + j] = z[off2 + j];
03652 
03653               for (octave_idx_type j = 0; j < n; j++)
03654                 z[off2 + j] = ctmp[j];
03655             }
03656         }
03657     }
03658   else
03659     {
03660       (*current_liboctave_error_handler)
03661         ("eigs: error %d in zneupd", info2);
03662       return -1;
03663     }
03664 
03665   return ip(4);
03666 }
03667 
03668 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
03669 extern octave_idx_type
03670 EigsRealSymmetricMatrix (const Matrix& m, const std::string typ,
03671                          octave_idx_type k, octave_idx_type p,
03672                          octave_idx_type &info, Matrix &eig_vec,
03673                          ColumnVector &eig_val, const Matrix& b,
03674                          ColumnVector &permB, ColumnVector &resid,
03675                          std::ostream &os, double tol = DBL_EPSILON,
03676                          bool rvec = false, bool cholB = 0, int disp = 0,
03677                          int maxit = 300);
03678 
03679 extern octave_idx_type
03680 EigsRealSymmetricMatrix (const SparseMatrix& m, const std::string typ,
03681                          octave_idx_type k, octave_idx_type p,
03682                          octave_idx_type &info, Matrix &eig_vec,
03683                          ColumnVector &eig_val, const SparseMatrix& b,
03684                          ColumnVector &permB, ColumnVector &resid,
03685                          std::ostream& os, double tol = DBL_EPSILON,
03686                          bool rvec = false, bool cholB = 0, int disp = 0,
03687                          int maxit = 300);
03688 
03689 extern octave_idx_type
03690 EigsRealSymmetricMatrixShift (const Matrix& m, double sigma,
03691                               octave_idx_type k, octave_idx_type p,
03692                               octave_idx_type &info, Matrix &eig_vec,
03693                               ColumnVector &eig_val, const Matrix& b,
03694                               ColumnVector &permB, ColumnVector &resid,
03695                               std::ostream &os, double tol = DBL_EPSILON,
03696                               bool rvec = false, bool cholB = 0, int disp = 0,
03697                               int maxit = 300);
03698 
03699 extern octave_idx_type
03700 EigsRealSymmetricMatrixShift (const SparseMatrix& m, double sigma,
03701                               octave_idx_type k, octave_idx_type p,
03702                               octave_idx_type &info, Matrix &eig_vec,
03703                               ColumnVector &eig_val, const SparseMatrix& b,
03704                               ColumnVector &permB, ColumnVector &resid,
03705                               std::ostream &os, double tol = DBL_EPSILON,
03706                               bool rvec = false, bool cholB = 0, int disp = 0,
03707                               int maxit = 300);
03708 
03709 extern octave_idx_type
03710 EigsRealSymmetricFunc (EigsFunc fun, octave_idx_type n,
03711                        const std::string &typ, double sigma,
03712                        octave_idx_type k, octave_idx_type p,
03713                        octave_idx_type &info,
03714                        Matrix &eig_vec, ColumnVector &eig_val,
03715                        ColumnVector &resid, std::ostream &os,
03716                        double tol = DBL_EPSILON, bool rvec = false,
03717                        bool cholB = 0, int disp = 0, int maxit = 300);
03718 
03719 extern octave_idx_type
03720 EigsRealNonSymmetricMatrix (const Matrix& m, const std::string typ,
03721                             octave_idx_type k, octave_idx_type p,
03722                             octave_idx_type &info, ComplexMatrix &eig_vec,
03723                             ComplexColumnVector &eig_val, const Matrix& b,
03724                             ColumnVector &permB, ColumnVector &resid,
03725                             std::ostream &os, double tol = DBL_EPSILON,
03726                             bool rvec = false, bool cholB = 0, int disp = 0,
03727                             int maxit = 300);
03728 
03729 extern octave_idx_type
03730 EigsRealNonSymmetricMatrix (const SparseMatrix& m, const std::string typ,
03731                             octave_idx_type k, octave_idx_type p,
03732                             octave_idx_type &info, ComplexMatrix &eig_vec,
03733                             ComplexColumnVector &eig_val,
03734                             const SparseMatrix& b,
03735                             ColumnVector &permB, ColumnVector &resid,
03736                             std::ostream &os, double tol = DBL_EPSILON,
03737                             bool rvec = false, bool cholB = 0, int disp = 0,
03738                             int maxit = 300);
03739 
03740 extern octave_idx_type
03741 EigsRealNonSymmetricMatrixShift (const Matrix& m, double sigma,
03742                                  octave_idx_type k, octave_idx_type p,
03743                                  octave_idx_type &info,
03744                                  ComplexMatrix &eig_vec,
03745                                  ComplexColumnVector &eig_val, const Matrix& b,
03746                                  ColumnVector &permB, ColumnVector &resid,
03747                                  std::ostream &os, double tol = DBL_EPSILON,
03748                                  bool rvec = false, bool cholB = 0,
03749                                  int disp = 0, int maxit = 300);
03750 
03751 extern octave_idx_type
03752 EigsRealNonSymmetricMatrixShift (const SparseMatrix& m, double sigma,
03753                                  octave_idx_type k, octave_idx_type p,
03754                                  octave_idx_type &info,
03755                                  ComplexMatrix &eig_vec,
03756                                  ComplexColumnVector &eig_val,
03757                                  const SparseMatrix& b,
03758                                  ColumnVector &permB, ColumnVector &resid,
03759                                  std::ostream &os, double tol = DBL_EPSILON,
03760                                  bool rvec = false, bool cholB = 0,
03761                                  int disp = 0, int maxit = 300);
03762 
03763 extern octave_idx_type
03764 EigsRealNonSymmetricFunc (EigsFunc fun, octave_idx_type n,
03765                           const std::string &_typ, double sigma,
03766                           octave_idx_type k, octave_idx_type p,
03767                           octave_idx_type &info, ComplexMatrix &eig_vec,
03768                           ComplexColumnVector &eig_val,
03769                           ColumnVector &resid, std::ostream& os,
03770                           double tol = DBL_EPSILON, bool rvec = false,
03771                           bool cholB = 0, int disp = 0, int maxit = 300);
03772 
03773 extern octave_idx_type
03774 EigsComplexNonSymmetricMatrix (const ComplexMatrix& m, const std::string typ,
03775                                octave_idx_type k, octave_idx_type p,
03776                                octave_idx_type &info, ComplexMatrix &eig_vec,
03777                                ComplexColumnVector &eig_val,
03778                                const ComplexMatrix& b, ColumnVector &permB,
03779                                ComplexColumnVector &resid,
03780                                std::ostream &os, double tol = DBL_EPSILON,
03781                                bool rvec = false, bool cholB = 0, int disp = 0,
03782                                int maxit = 300);
03783 
03784 extern octave_idx_type
03785 EigsComplexNonSymmetricMatrix (const SparseComplexMatrix& m,
03786                                const std::string typ, octave_idx_type k,
03787                                octave_idx_type p, octave_idx_type &info,
03788                                ComplexMatrix &eig_vec,
03789                                ComplexColumnVector &eig_val,
03790                                const SparseComplexMatrix& b,
03791                                ColumnVector &permB,
03792                                ComplexColumnVector &resid,
03793                                std::ostream &os, double tol = DBL_EPSILON,
03794                                bool rvec = false, bool cholB = 0, int disp = 0,
03795                                int maxit = 300);
03796 
03797 extern octave_idx_type
03798 EigsComplexNonSymmetricMatrixShift (const ComplexMatrix& m, Complex sigma,
03799                                     octave_idx_type k, octave_idx_type p,
03800                                     octave_idx_type &info,
03801                                     ComplexMatrix &eig_vec,
03802                                     ComplexColumnVector &eig_val,
03803                                     const ComplexMatrix& b,
03804                                     ColumnVector &permB,
03805                                     ComplexColumnVector &resid,
03806                                     std::ostream &os, double tol = DBL_EPSILON,
03807                                     bool rvec = false, bool cholB = 0,
03808                                     int disp = 0, int maxit = 300);
03809 
03810 extern octave_idx_type
03811 EigsComplexNonSymmetricMatrixShift (const SparseComplexMatrix& m,
03812                                     Complex sigma,
03813                                     octave_idx_type k, octave_idx_type p,
03814                                     octave_idx_type &info,
03815                                     ComplexMatrix &eig_vec,
03816                                     ComplexColumnVector &eig_val,
03817                                     const SparseComplexMatrix& b,
03818                                     ColumnVector &permB,
03819                                     ComplexColumnVector &resid,
03820                                     std::ostream &os, double tol = DBL_EPSILON,
03821                                     bool rvec = false, bool cholB = 0,
03822                                     int disp = 0, int maxit = 300);
03823 
03824 extern octave_idx_type
03825 EigsComplexNonSymmetricFunc (EigsComplexFunc fun, octave_idx_type n,
03826                              const std::string &_typ, Complex sigma,
03827                              octave_idx_type k, octave_idx_type p,
03828                              octave_idx_type &info, ComplexMatrix &eig_vec,
03829                              ComplexColumnVector &eig_val,
03830                              ComplexColumnVector &resid, std::ostream& os,
03831                              double tol = DBL_EPSILON, bool rvec = false,
03832                              bool cholB = 0, int disp = 0, int maxit = 300);
03833 #endif
03834 
03835 #ifndef _MSC_VER
03836 template static octave_idx_type
03837 lusolve (const SparseMatrix&, const SparseMatrix&, Matrix&);
03838 
03839 template static octave_idx_type
03840 lusolve (const SparseComplexMatrix&, const SparseComplexMatrix&,
03841          ComplexMatrix&);
03842 
03843 template static octave_idx_type
03844 lusolve (const Matrix&, const Matrix&, Matrix&);
03845 
03846 template static octave_idx_type
03847 lusolve (const ComplexMatrix&, const ComplexMatrix&, ComplexMatrix&);
03848 
03849 template static ComplexMatrix
03850 ltsolve (const SparseComplexMatrix&, const ColumnVector&,
03851          const ComplexMatrix&);
03852 
03853 template static Matrix
03854 ltsolve (const SparseMatrix&, const ColumnVector&, const Matrix&);
03855 
03856 template static ComplexMatrix
03857 ltsolve (const ComplexMatrix&, const ColumnVector&, const ComplexMatrix&);
03858 
03859 template static Matrix
03860 ltsolve (const Matrix&, const ColumnVector&, const Matrix&);
03861 
03862 template static ComplexMatrix
03863 utsolve (const SparseComplexMatrix&, const ColumnVector&,
03864          const ComplexMatrix&);
03865 
03866 template static Matrix
03867 utsolve (const SparseMatrix&, const ColumnVector&, const Matrix&);
03868 
03869 template static ComplexMatrix
03870 utsolve (const ComplexMatrix&, const ColumnVector&, const ComplexMatrix&);
03871 
03872 template static Matrix
03873 utsolve (const Matrix&, const ColumnVector&, const Matrix&);
03874 #endif
03875 
03876 #endif