arpack/arpack.hpp

Go to the documentation of this file.

#ifndef _ARPACK_HPP
#define _ARPACK_HPP


//
// Copyright (c) 2009 Simen Kvaal
//
// This file is part of OpenFCI.
//
// OpenFCI is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// OpenFCI is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with OpenFCI. If not, see <http://www.gnu.org/licenses/>.
//



#include "saupp.h"
#include "seupp.h"
#include <cassert>
#include <iostream>

namespace arpack
{
  
  template<class matrix_t>
  class symmetric_eigensolver
  {
    private:
      typedef void (matrix_t::* operation_func_t)(double*, double*);

      bool debug;                 

      matrix_t *mat;              
      operation_func_t matrix_op; 
      
      ARint n;                    
      ARint nev;                  
      ARint ncv;                  
      double tol;                 

      ARint ido;                  
      char bmat;                  
      char which[3];              

      std::vector<double> resid;       
      std::vector<double> V;           
      ARint ldv;                  
      std::vector<ARint> iparam;       
      std::vector<ARint> ipntr;        
      std::vector<double> workd;       
      std::vector<double> workl;       
      ARint lworkl;               
      ARint info;                 

      ARint maxit;                

      ARint nconv;                
      std::vector<double> d;           
      std::vector<double> Z;           

    public:
      void set_matrix_op(matrix_t *the_mat, void (matrix_t::* the_matrix_op)(double*, double*))
      {
        mat = the_mat;
        matrix_op = the_matrix_op;
      }

      void set_parameters(ARint the_n, ARint the_nev,  double the_tol, ARint the_ncv = 0)
      {
        assert(the_n > 0);
        assert(the_nev > 0);
        assert(the_nev < the_n);
        
        n = the_n;
        nev = the_nev;
        tol = the_tol;
        if (the_ncv == 0)
        {
          ncv = min(2*nev+1,n);
        }
        else
          ncv = the_ncv;
        
      }
      

      symmetric_eigensolver() 
      { 
        maxit = 1000;
        debug = true;
      }

      symmetric_eigensolver(matrix_t *the_mat, void (matrix_t::* the_matrix_op)(double*, double*))  
      { 
        debug = true;
        maxit = 1000;
        set_matrix_op(the_mat, the_matrix_op);
      }

      void set_debug(bool d) { debug = d; }



      ARint solve()
      {
        

        bool finished;
        
        finished = false;

        ido = 0;
        bmat = 'I';
        strcpy(which,"SM");
        tol = 1e-8;
        resid.resize(n); // initial vector if info != 0.
        info = 0; // use random initial vector.
        V.resize(ncv*n+1);
        ldv = n;
        iparam.resize(12); // set them!
        iparam[1] = 1;
        iparam[3] = maxit;
        iparam[4] = 1;
        iparam[7] = 1;
        ipntr.resize(12);
        
        workd.resize(3*n + 1);
        lworkl = ncv*(ncv+8);
        workl.resize(lworkl + 1);
        
        
        int it = 1;

        while (finished == false)
        {
          it++; // increase iteration counter
          if (debug)
          {
            std::cerr << "Iteration number " << it << " commencing." << endl;
          }
          saupp(ido, bmat, n, which, nev, tol, &resid[0],
          ncv, &V[0], ldv, &iparam[0], &ipntr[0], &workd[0], &workl[0], lworkl, info);
          
          if (debug)
            std::cerr << "(ido, info) == " << ido << ", " << info << std::endl;
          
          if ((ido == -1) || (ido == 1))
          {
            (mat->*matrix_op)(&workd[ipntr[1]], &workd[ipntr[2]]);
          }
          else
            finished = true;
          
          
        }
        
        if (info != 0)
        {
          std::cerr << "Warning! Exit status of ARPACK was info == " << info << "." << std::endl;
        }
        

        nconv = iparam[5];
        bool rvec = true;
        double sigma = 0;

        if (debug)
        {
          std::cerr << "Found nconv == " << nconv << " out of nev == " << nev << " eigenvalues." << std::endl;
        }
        
        d.resize(n);
        Z.resize(nconv*n);
        
        seupp(rvec, 'A', &d[0], &Z[0],
          n, sigma, bmat, n,
          which, nev, tol, &resid[0],
          ncv, &V[0], ldv, &iparam[0],
          &ipntr[0], &workd[0], &workl[0],
          lworkl, info);
        
        
        return nconv;
        
      }
      

      void eigenvalues(std::vector<double>& evals)
      {
        evals.resize(nconv);
        for (size_t k = 0; k<(size_t)nconv; ++k)
          evals[k] = d[k];
      
      }

      double  get_eigenvalue(size_t k)
      {
        return d[k];
      }

      double get_eigenvector(size_t j, size_t k)
      {
        return Z[n*j + k];
      }

  };
  

} // namespace arpack


#endif

---