dgematrix Class Reference

Real Double-precision General Dence Matrix Class. More...

#include <dgematrix.hpp>

Public Member Functions
	dgematrix ()
	dgematrix (const dgematrix &)
	dgematrix (const _dgematrix &)
	dgematrix (const CPPL_INT &, const CPPL_INT &)
	dgematrix (const char *)
	~dgematrix ()
_zgematrix	to_zgematrix () const
double &	operator() (const CPPL_INT &, const CPPL_INT &)
double	operator() (const CPPL_INT &, const CPPL_INT &) const
dgematrix &	set (const CPPL_INT &, const CPPL_INT &, const double &)
void	write (const char *) const
void	read (const char *)
void	clear ()
dgematrix &	zero ()
dgematrix &	identity ()
void	chsign ()
void	copy (const dgematrix &)
void	shallow_copy (const _dgematrix &)
dgematrix &	resize (const CPPL_INT &, const CPPL_INT &)
_drovector	row (const CPPL_INT &) const
_dcovector	col (const CPPL_INT &) const
CPPL_INT	dgesv (dgematrix &)
CPPL_INT	dgesv (dcovector &)
CPPL_INT	dgels (dgematrix &)
CPPL_INT	dgels (dcovector &)
CPPL_INT	dgels (dgematrix &, drovector &)
CPPL_INT	dgels (dcovector &, double &)
CPPL_INT	dgelss (dcovector &, dcovector &, CPPL_INT &, const double)
CPPL_INT	dgelss (dgematrix &, dcovector &, CPPL_INT &, const double)
CPPL_INT	dgelsd (dcovector &, dcovector &, CPPL_INT &, const double)
CPPL_INT	dgeev (std::vector< double > &, std::vector< double > &)
CPPL_INT	dgeev (zcovector &)
CPPL_INT	dgeev (std::vector< double > &, std::vector< double > &, std::vector< dcovector > &, std::vector< dcovector > &)
CPPL_INT	dgeev (std::vector< double > &, std::vector< double > &, std::vector< drovector > &, std::vector< drovector > &)
CPPL_INT	dggev (dgematrix &, std::vector< double > &, std::vector< double > &)
CPPL_INT	dggev (dgematrix &, std::vector< double > &, std::vector< double > &, std::vector< dcovector > &, std::vector< dcovector > &)
CPPL_INT	dggev (dgematrix &, std::vector< double > &, std::vector< double > &, std::vector< drovector > &, std::vector< drovector > &)
CPPL_INT	dgesvd (dgbmatrix &)
CPPL_INT	dgesvd (dcovector &, dgematrix &, dgematrix &)
CPPL_INT	dgglse (dgematrix &, dcovector &, dcovector &, dcovector &)
dgematrix &	operator= (const dgematrix &)
dgematrix &	operator= (const _dgematrix &)
dgematrix &	operator+= (const dgematrix &)
dgematrix &	operator+= (const _dgematrix &)
dgematrix &	operator+= (const dsymatrix &)
dgematrix &	operator+= (const _dsymatrix &)
dgematrix &	operator+= (const dgbmatrix &)
dgematrix &	operator+= (const _dgbmatrix &)
dgematrix &	operator+= (const dgsmatrix &)
dgematrix &	operator+= (const _dgsmatrix &)
dgematrix &	operator+= (const dssmatrix &)
dgematrix &	operator+= (const _dssmatrix &)
dgematrix &	operator-= (const dgematrix &)
dgematrix &	operator-= (const _dgematrix &)
dgematrix &	operator-= (const dsymatrix &)
dgematrix &	operator-= (const _dsymatrix &)
dgematrix &	operator-= (const dgbmatrix &)
dgematrix &	operator-= (const _dgbmatrix &)
dgematrix &	operator-= (const dgsmatrix &)
dgematrix &	operator-= (const _dgsmatrix &)
dgematrix &	operator-= (const dssmatrix &)
dgematrix &	operator-= (const _dssmatrix &)
dgematrix &	operator*= (const dgematrix &)
dgematrix &	operator*= (const _dgematrix &)
dgematrix &	operator*= (const dsymatrix &)
dgematrix &	operator*= (const _dsymatrix &)
dgematrix &	operator*= (const dgbmatrix &)
dgematrix &	operator*= (const _dgbmatrix &)
dgematrix &	operator*= (const dgsmatrix &)
dgematrix &	operator*= (const _dgsmatrix &)
dgematrix &	operator*= (const dssmatrix &)
dgematrix &	operator*= (const _dssmatrix &)
dgematrix &	operator*= (const double &)
dgematrix &	operator/= (const double &)

Public Attributes
CPPL_INT	m
	matrix row size More...
CPPL_INT	n
	matrix column size More...
double *	array
	1D array to store matrix data More...
double **	darray
	array of pointers of column head addresses More...

Friends
std::ostream &	operator<< (std::ostream &, const dgematrix &)
void	swap (dgematrix &, dgematrix &)
_dgematrix	_ (dgematrix &)
_dgematrix	t (const dgematrix &)
_dgematrix	i (const dgematrix &)
void	idamax (CPPL_INT &, CPPL_INT &, const dgematrix &)
double	damax (const dgematrix &)
const dgematrix &	operator+ (const dgematrix &)
_dgematrix	operator- (const dgematrix &)
_dgematrix	operator+ (const dgematrix &, const dgematrix &)
_dgematrix	operator+ (const dgematrix &, const _dgematrix &)
_dgematrix	operator+ (const dgematrix &, const dsymatrix &)
_dgematrix	operator+ (const dgematrix &, const _dsymatrix &)
_dgematrix	operator+ (const dgematrix &, const dgbmatrix &)
_dgematrix	operator+ (const dgematrix &, const _dgbmatrix &)
_dgematrix	operator+ (const dgematrix &, const dgsmatrix &)
_dgematrix	operator+ (const dgematrix &, const _dgsmatrix &)
_dgematrix	operator+ (const dgematrix &, const dssmatrix &)
_dgematrix	operator+ (const dgematrix &, const _dssmatrix &)
_dgematrix	operator- (const dgematrix &, const dgematrix &)
_dgematrix	operator- (const dgematrix &, const _dgematrix &)
_dgematrix	operator- (const dgematrix &, const dsymatrix &)
_dgematrix	operator- (const dgematrix &, const _dsymatrix &)
_dgematrix	operator- (const dgematrix &, const dgbmatrix &)
_dgematrix	operator- (const dgematrix &, const _dgbmatrix &)
_dgematrix	operator- (const dgematrix &, const dgsmatrix &)
_dgematrix	operator- (const dgematrix &, const _dgsmatrix &)
_dgematrix	operator- (const dgematrix &, const dssmatrix &)
_dgematrix	operator- (const dgematrix &, const _dssmatrix &)
_dcovector	operator* (const dgematrix &, const dcovector &)
_dcovector	operator* (const dgematrix &, const _dcovector &)
_dgematrix	operator* (const dgematrix &, const dgematrix &)
_dgematrix	operator* (const dgematrix &, const _dgematrix &)
_dgematrix	operator* (const dgematrix &, const dsymatrix &)
_dgematrix	operator* (const dgematrix &, const _dsymatrix &)
_dgematrix	operator* (const dgematrix &, const dgbmatrix &)
_dgematrix	operator* (const dgematrix &, const _dgbmatrix &)
_dgematrix	operator* (const dgematrix &, const dgsmatrix &)
_dgematrix	operator* (const dgematrix &, const _dgsmatrix &)
_dgematrix	operator* (const dgematrix &, const dssmatrix &)
_dgematrix	operator* (const dgematrix &, const _dssmatrix &)
_dgematrix	operator* (const dgematrix &, const double &)
_dgematrix	operator/ (const dgematrix &, const double &)
_drovector	operator% (const dgematrix &, const dgematrix &)
_dgematrix	operator* (const double &, const dgematrix &)
_dgematrix	hadamard (const dgematrix &, const dgematrix &)

Detailed Description

Real Double-precision General Dence Matrix Class.

Definition at line 3 of file dgematrix.hpp.

Constructor & Destructor Documentation

dgematrix::dgematrix ( )

inline

dgematrix constructor without arguments

Definition at line 3 of file dgematrix-constructor.hpp.

References array, darray, m, and n.

{CPPL_VERBOSE_REPORT;
  //////// initialize ////////
  m =0;
  n =0;
  array =NULL;
  darray =NULL;
}

dgematrix::dgematrix ( const dgematrix &  mat )

inline

dgematrix copy constructor

Definition at line 18 of file dgematrix-constructor.hpp.

References array, darray, i, m, and n.

{CPPL_VERBOSE_REPORT;
  //////// initialize ////////
  m =mat.m;
  n =mat.n;
  array =new double[m*n];
  darray =new double*[n];
  for(int i=0; i<n; i++){
    darray[i] =&array[i*m];
  }
  
  //////// copy ////////
  CPPL_INT mn =m*n;
  CPPL_INT inc =1;
  dcopy_(&mn, mat.array, &inc, array, &inc);
}

dgematrix::dgematrix ( const _dgematrix &  mat )

inline

dgematrix constructor to cast _dgematrix

Definition at line 37 of file dgematrix-constructor.hpp.

References array, _dgematrix::array, _dgematrix::darray, darray, m, _dgematrix::m, _dgematrix::n, n, and _dgematrix::nullify().

{CPPL_VERBOSE_REPORT;
  m =mat.m;
  n =mat.n;
  array =mat.array;
  darray =mat.darray;
  
  mat.nullify();
}

dgematrix::dgematrix ( const CPPL_INT &  _m, const CPPL_INT &  _n  )

inline

dgematrix constructor with size specification

Definition at line 53 of file dgematrix-constructor.hpp.

References array, darray, i, m, and n.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if( _m<0 || _n<0 ){
    ERROR_REPORT;
    std::cerr << "Matrix sizes must be positive integers. " << std::endl
              << "Your input was (" << _m << "," << _n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  //////// initialize ////////
  m =_m;
  n =_n;
  array =new double[m*n];
  darray  =new double*[n];
  for(int i=0; i<n; i++){
    darray[i] =&array[i*m];
  }
}

dgematrix::dgematrix ( const char *  filename )

inline

dgematrix constructor with filename

Definition at line 76 of file dgematrix-constructor.hpp.

References array, darray, and read().

{CPPL_VERBOSE_REPORT;
  array =NULL;
  darray =NULL;
  
  //// read ////
  read(filename);
}

dgematrix::~dgematrix ( )

inline

dgematrix destructor

Definition at line 91 of file dgematrix-constructor.hpp.

References array, and darray.

{CPPL_VERBOSE_REPORT;
  delete [] darray;
  delete [] array;
}

Member Function Documentation

_zgematrix dgematrix::to_zgematrix ( ) const

inline

cast to _zgematrix

Definition at line 3 of file dgematrix-cast.hpp.

References _, array, zgematrix::array, i, m, and n.

{CPPL_VERBOSE_REPORT;
  zgematrix newmat(m,n);
  for(CPPL_INT i=0; i<m*n; i++){
    newmat.array[i] =comple(array[i],0.0);
  }
  
  return _(newmat);
}

double & dgematrix::operator() ( const CPPL_INT &  i, const CPPL_INT &  j  )

inline

operator() for non-const object

Definition at line 3 of file dgematrix-io.hpp.

References darray, i, m, and n.

Referenced by identity(), operator*=(), operator+=(), operator-=(), read(), and write().

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if( i<0 || j<0 || m<=i || n<=j ){
    ERROR_REPORT;
    std::cerr << "The required component is out of the matrix size." << std::endl
              << "Your input is (" << i << "," << j << "), whereas the matrix size is " << m << "x" << n << "." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  //return array[i+m*j];
  return darray[j][i];
}

double dgematrix::operator() ( const CPPL_INT &  i, const CPPL_INT &  j  ) const

inline

operator() for const object

Definition at line 20 of file dgematrix-io.hpp.

References darray, i, m, and n.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if( i<0 || j<0 || m<=i || n<=j ){
    ERROR_REPORT;
    std::cerr << "The required component is out of the matrix size." << std::endl
              << "Your input is (" << i << "," << j << "), whereas the matrix size is " << m << "x" << n << "." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  //return array[i+m*j];
  return darray[j][i];
}

dgematrix & dgematrix::set ( const CPPL_INT &  i, const CPPL_INT &  j, const double &  v  )

inline

set value for const object

Definition at line 41 of file dgematrix-io.hpp.

References darray, i, m, and n.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if( i<0 || j<0 || m<=i || n<=j ){
    ERROR_REPORT;
    std::cerr << "The required component is out of the matrix size." << std::endl
              << "Your input is (" << i << "," << j << "), whereas the matrix size is " << m << "x" << n << "." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  //array[i+m*j] =v;
  darray[j][i] =v;
  
  return *this;
}

void dgematrix::write ( const char *  filename ) const

inline

Definition at line 79 of file dgematrix-io.hpp.

References i, m, n, and operator()().

{CPPL_VERBOSE_REPORT;
  std::ofstream ofs(filename, std::ios::trunc);
  ofs.setf(std::cout.flags());
  ofs.precision(std::cout.precision());
  ofs.width(std::cout.width());
  ofs.fill(std::cout.fill());
  
  ofs << "#dgematrix" << " " << m << " " << n << std::endl;
  for(CPPL_INT i=0; i<m; i++){
    for(CPPL_INT j=0; j<n; j++ ){
      ofs << operator()(i,j) << " ";
    }
    ofs << std::endl;
  }
  
  ofs.close();
}

void dgematrix::read ( const char *  filename )

inline

Definition at line 99 of file dgematrix-io.hpp.

References i, m, n, operator()(), and resize().

Referenced by dgematrix().

{CPPL_VERBOSE_REPORT;
  std::ifstream s( filename );
  if(!s){
    ERROR_REPORT;
    std::cerr << "The file \"" << filename << "\" can not be opened." << std::endl;
    exit(1);
  }

  std::string id;
  s >> id;
  if( id != "dgematrix" && id != "#dgematrix" ){
    ERROR_REPORT;
    std::cerr << "The type name of the file \"" << filename << "\" is not dgematrix." << std::endl
              << "Its type name was " << id << " ." << std::endl;
    exit(1);
  }
  
  s >> m >> n;
  resize(m, n);
  for(CPPL_INT i=0; i<m; i++){
    for(CPPL_INT j=0; j<n; j++ ){
      s >> operator()(i,j);
    }
  }
  if(s.eof()){
    ERROR_REPORT;
    std::cerr << "There is something is wrong with the file \"" << filename << "\"." << std::endl
              << "Most likely, there is a lack of data components, or a linefeed code or space code is missing at the end of the last line." << std::endl;
    exit(1);
  }
  
  s >> id;
  if(!s.eof()){
    ERROR_REPORT;
    std::cerr << "There is something is wrong with the file \"" << filename << "\"." << std::endl
              << "Most likely, there are extra data components." << std::endl;
    exit(1);
  }

  
  s.close();
}

void dgematrix::clear ( )

inline

clear all the matrix data and set the sizes 0

Definition at line 3 of file dgematrix-misc.hpp.

References array, darray, m, and n.

Referenced by dgels(), and dgelss().

{CPPL_VERBOSE_REPORT;
  m =0;
  n =0;
  delete [] array;
  delete [] darray;
  array =NULL;
  darray =NULL;
}

dgematrix & dgematrix::zero ( )

inline

change the matrix into a zero matrix

Definition at line 15 of file dgematrix-misc.hpp.

References array, i, m, and n.

Referenced by operator*(), operator*=(), _dgsmatrix::to_dgematrix(), _dssmatrix::to_dgematrix(), dssmatrix::to_dgematrix(), dgrmatrix::to_dgematrix(), and dgsmatrix::to_dgematrix().

{CPPL_VERBOSE_REPORT;
  for(CPPL_INT i=0; i<m*n; i++){
    array[i] =0.0;
  }
  return *this;
}

dgematrix & dgematrix::identity ( )

inline

change the matrix into an identity matrix

Definition at line 25 of file dgematrix-misc.hpp.

References array, i, m, n, and operator()().

Referenced by i().

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(m!=n){
    ERROR_REPORT;
    std::cerr << "Only square matrix can be a identity matrix." << std::endl
              << "The matrix size was " << m << "x" << n << "." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  for(CPPL_INT i=0; i<m*n; i++){
    array[i] =0.0;
  }
  for(CPPL_INT i=0; i<m; i++){
    operator()(i,i) =1.0;
  }
  return *this;
}

void dgematrix::chsign ( )

inline

change sign(+/-) of the matrix

Definition at line 47 of file dgematrix-misc.hpp.

References array, i, m, and n.

{CPPL_VERBOSE_REPORT;
  for(CPPL_INT i=0; i<m*n; i++){
    array[i] =-array[i];
  }
}

void dgematrix::copy ( const dgematrix &  mat )

inline

make a deep copy of the matrix

Definition at line 56 of file dgematrix-misc.hpp.

References array, darray, i, m, and n.

Referenced by operator=().

{CPPL_VERBOSE_REPORT;
  m =mat.m;
  n =mat.n;
  delete [] array;
  array =new double[mat.m*mat.n];
  delete [] darray;
  darray =new double*[n];
  for(int i=0; i<n; i++){
    darray[i] =&array[i*m];
  }
  
  CPPL_INT mn =mat.m*mat.n;
  CPPL_INT inc =1;
  dcopy_(&mn, mat.array, &inc, array, &inc);
}

void dgematrix::shallow_copy ( const _dgematrix &  mat )

inline

make a shallow copy of the matrix
This function is not designed to be used in project codes.

Definition at line 76 of file dgematrix-misc.hpp.

References array, _dgematrix::array, _dgematrix::darray, darray, m, _dgematrix::m, _dgematrix::n, n, and _dgematrix::nullify().

Referenced by operator=().

{CPPL_VERBOSE_REPORT;
  m =mat.m;
  n =mat.n;
  delete [] array;
  array =mat.array;
  delete [] darray;
  darray =mat.darray;
  
  mat.nullify();
}

dgematrix & dgematrix::resize ( const CPPL_INT &  _m, const CPPL_INT &  _n  )

inline

resize the matrix

Definition at line 90 of file dgematrix-misc.hpp.

References array, darray, i, m, and n.

Referenced by dgesvd(), and read().

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if( _m<0 || _n<0 ){
    ERROR_REPORT;
    std::cerr << "Matrix sizes must be positive integers." << std::endl
              << "Your input was (" << _m << "," << _n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  m =_m;
  n =_n;
  delete [] array;
  array =new double[m*n];
  delete [] darray;
  darray =new double*[n];
  for(int i=0; i<n; i++){
    darray[i] =&array[i*m];
  }
  
  return *this;
}

_drovector dgematrix::row ( const CPPL_INT &  _m ) const

inline

get row of the matrix

Definition at line 120 of file dgematrix-misc.hpp.

References _, m, and n.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if( _m<0 || _m>m ){
    ERROR_REPORT;
    std::cerr << "Input row number must be between 0 and " << m << "." << std::endl
              << "Your input was " << _m << "." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  drovector v(n);
  for(CPPL_INT j=0; j<n; j++){
    v(j)=(*this)(_m,j);
  }
  return _(v);
}

_dcovector dgematrix::col ( const CPPL_INT &  _n ) const

inline

get column of the matrix

Definition at line 140 of file dgematrix-misc.hpp.

References _, i, m, and n.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if( _n<0 || _n>n ){
    ERROR_REPORT;
    std::cerr << "Input row number must be between 0 and " << n << "." << std::endl
              << "Your input was " << _n << "." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  dcovector v(m);
  for(CPPL_INT i=0; i<m; i++){
    v(i)=(*this)(i,_n);
  }
  return _(v);
}

CPPL_INT dgematrix::dgesv ( dgematrix &  mat )

inline

solve A*X=Y using dgesv
The argument is dgematrix Y. Y is overwritten and become the solution X. A is also overwritten and become P*l*U.

Definition at line 6 of file dgematrix-lapack.hpp.

References array, m, and n.

Referenced by i().

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(m!=n || m!=mat.m){
    ERROR_REPORT;
    std::cerr << "These two matrices cannot be solved." << std::endl
              << "Your input was (" << m << "x" << n << ") and (" << mat.m << "x" << mat.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG 
  
  CPPL_INT NRHS(mat.n), LDA(n), *IPIV(new CPPL_INT[n]), LDB(mat.m), INFO(1);
  dgesv_(&n, &NRHS, array, &LDA, IPIV, mat.array, &LDB, &INFO);
  delete [] IPIV;
  
  if(INFO!=0){
    WARNING_REPORT;
    std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
  }
  return INFO;
}

CPPL_INT dgematrix::dgesv ( dcovector &  vec )

inline

solve A*x=y using dgesv
The argument is dcovector y. y is overwritten and become the solution x. A is also overwritten and become P*l*U.

Definition at line 33 of file dgematrix-lapack.hpp.

References array, dcovector::array, dcovector::l, m, and n.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(m!=n || m!=vec.l){
    ERROR_REPORT;
    std::cerr << "These matrix and vector cannot be solved." << std::endl
              << "Your input was (" << m << "x" << n << ") and (" << vec.l << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG 
  
  CPPL_INT NRHS(1), LDA(n), *IPIV(new CPPL_INT[n]), LDB(vec.l), INFO(1);
  dgesv_(&n, &NRHS, array, &LDA, IPIV, vec.array, &LDB, &INFO);
  delete [] IPIV;
  
  if(INFO!=0){
    WARNING_REPORT;
    std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
  }
  return INFO;
}

CPPL_INT dgematrix::dgels ( dgematrix &  mat )

inline

solve overdetermined or underdetermined A*X=Y using dgels

Definition at line 62 of file dgematrix-lapack.hpp.

References array, clear(), i, m, n, and swap.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(m!=mat.m){
    ERROR_REPORT;
    std::cerr << "These two matrices cannot be solved." << std::endl
              << "Your input was (" << m << "x" << n << ") and (" << mat.m << "x" << mat.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG    
  
  if(m<n){ //underdetermined
    dgematrix tmp(n,mat.n);
    for(CPPL_INT i=0; i<mat.m; i++){
      for(CPPL_INT j=0; j<mat.n; j++){
        tmp(i,j) =mat(i,j);
      }
    }
    mat.clear();
    swap(mat,tmp);
  }
  
  char TRANS('n');
  CPPL_INT NRHS(mat.n), LDA(m), LDB(mat.m), LWORK(std::min(m,n)+std::max(std::min(m,n),NRHS)), INFO(1);
  double *WORK(new double[LWORK]);
  dgels_(&TRANS, &m, &n, &NRHS, array, &LDA, mat.array, &LDB, WORK, &LWORK, &INFO);
  delete [] WORK;
  
  if(m>n){ //overdetermined
    dgematrix tmp(n,mat.n);
    for(CPPL_INT i=0; i<tmp.m; i++){
      for(CPPL_INT j=0; j<tmp.n; j++){
        tmp(i,j) =mat(i,j);
      }
    }
    mat.clear();
    swap(mat,tmp);
  }
  
  if(INFO!=0){
    WARNING_REPORT;
    std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
  }
  return INFO;
}

CPPL_INT dgematrix::dgels ( dcovector &  vec )

inline

solve overdetermined or underdetermined A*x=y using dgels

Definition at line 111 of file dgematrix-lapack.hpp.

References array, dcovector::array, dcovector::clear(), i, dcovector::l, m, n, and swap.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(m!=vec.l){
    ERROR_REPORT;
    std::cerr << "These matrix and vector cannot be solved." << std::endl
              << "Your input was (" << m << "x" << n << ") and (" << vec.l << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG    
  
  if(m<n){ //underdetermined
    dcovector tmp(n);
    for(CPPL_INT i=0; i<vec.l; i++){
      tmp(i)=vec(i);
    }
    vec.clear();
    swap(vec,tmp);
  }
  
  char TRANS('n');
  CPPL_INT NRHS(1), LDA(m), LDB(vec.l), LWORK(std::min(m,n)+std::max(std::min(m,n),NRHS)), INFO(1);
  double *WORK(new double[LWORK]);
  dgels_(&TRANS, &m, &n, &NRHS, array, &LDA, vec.array, &LDB, WORK, &LWORK, &INFO);
  delete [] WORK;
  
  if(m>n){ //overdetermined
    dcovector tmp(n);
    for(CPPL_INT i=0; i<tmp.l; i++){
      tmp(i)=vec(i);
    }
    vec.clear();
    swap(vec,tmp);
  }
  
  if(INFO!=0){
    WARNING_REPORT;
    std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
  }
  return INFO;
}

CPPL_INT dgematrix::dgels ( dgematrix &  mat, drovector &  residual  )

inline

solve overdetermined or underdetermined A*X=Y using dgels with the sum of residual squares output
The residual is set as the columnwise sum of residual squares for overdetermined problems while it is always zero for underdetermined problems.

Definition at line 160 of file dgematrix-lapack.hpp.

References array, clear(), i, drovector::l, m, n, drovector::resize(), swap, and drovector::zero().

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(m!=mat.m){
    ERROR_REPORT;
    std::cerr << "These two matrices cannot be solved." << std::endl
              << "Your input was (" << m << "x" << n << ") and (" << mat.m << "x" << mat.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  residual.resize(mat.n);
  residual.zero();
  
  if(m<n){ //underdetermined
    dgematrix tmp(n,mat.n);
    for(CPPL_INT i=0; i<mat.m; i++){
      for(CPPL_INT j=0; j<mat.n; j++){
        tmp(i,j) =mat(i,j);
      }
    }
    mat.clear();
    swap(mat,tmp);
  }
  
  char TRANS('n');
  CPPL_INT NRHS(mat.n), LDA(m), LDB(mat.m), LWORK(std::min(m,n)+std::max(std::min(m,n),NRHS)), INFO(1);
  double *WORK(new double[LWORK]);
  dgels_(&TRANS, &m, &n, &NRHS, array, &LDA, mat.array, &LDB, WORK, &LWORK, &INFO);
  delete [] WORK;
  
  if(m>n){ //overdetermined
    for(CPPL_INT i=0; i<residual.l; i++){
      for(CPPL_INT j=0; j<m-n; j++){
        residual(i) += std::pow(mat(n+j,i), 2.0);
      }
    }
    
    dgematrix tmp(n,mat.n);
    for(CPPL_INT i=0; i<tmp.m; i++){
      for(CPPL_INT j=0; j<tmp.n; j++){
        tmp(i,j) =mat(i,j);
      }
    }
    mat.clear();
    swap(mat,tmp);
  }
  
  if(INFO!=0){
    WARNING_REPORT;
    std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
  }
  return INFO;
}

CPPL_INT dgematrix::dgels ( dcovector &  vec, double &  residual  )

inline

solve overdetermined or underdetermined A*x=y using dgels with the sum of residual squares output
The residual is set as the sum of residual squares for overdetermined problems while it is always zero for underdetermined problems.

Definition at line 222 of file dgematrix-lapack.hpp.

References array, dcovector::array, dcovector::clear(), i, dcovector::l, m, n, and swap.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(m!=vec.l){
    ERROR_REPORT;
    std::cerr << "These matrix and vector cannot be solved." << std::endl
              << "Your input was (" << m << "x" << n << ") and (" << vec.l << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG    
  
  residual=0.0;
  
  if(m<n){ //underdetermined
    dcovector tmp(n);
    for(CPPL_INT i=0; i<vec.l; i++){
      tmp(i)=vec(i);
    }
    vec.clear();
    swap(vec,tmp);
  }
  
  char TRANS('n');
  CPPL_INT NRHS(1), LDA(m), LDB(vec.l), LWORK(std::min(m,n)+std::max(std::min(m,n),NRHS)), INFO(1);
  double *WORK(new double[LWORK]);
  dgels_(&TRANS, &m, &n, &NRHS, array, &LDA, vec.array, &LDB, WORK, &LWORK, &INFO);
  delete [] WORK;
  
  if(m>n){ //overdetermined
    for(CPPL_INT i=0; i<m-n; i++){
      residual+=std::pow(vec(n+i),2.0);
    }
    
    dcovector tmp(n);
    for(CPPL_INT i=0; i<tmp.l; i++){
      tmp(i)=vec(i);
    }
    vec.clear();
    swap(vec,tmp);
  }
  
  if(INFO!=0){
    WARNING_REPORT;
    std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
  }
  return INFO;
}

CPPL_INT dgematrix::dgelss ( dcovector &  B, dcovector &  S, CPPL_INT &  RANK, const double  RCOND = -1.  )

inline

calculate the least-squares-least-norm solution for overdetermined or underdetermined A*x=y using dgelss

Definition at line 279 of file dgematrix-lapack.hpp.

References dcovector::array, dcovector::clear(), i(), dcovector::l, dcovector::resize(), and swap().

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(m!=B.l){
    ERROR_REPORT;
    std::cerr << "These matrix and vector cannot be solved." << std::endl
              << "Your input was (" << m << "x" << n << ") and (" << B.l << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG    
  
  if(m<n){ //underdetermined
    dcovector tmp(n);
    for(CPPL_INT i=0; i<B.l; i++){
      tmp(i)=B(i);
    }
    B.clear();
    swap(B,tmp);
  }
  
  S.resize(std::min(m,n));
  
  CPPL_INT NRHS(1), LDA(m), LDB(B.l), LWORK(3*std::min(m,n)+std::max(std::max(2*std::min(m,n),std::max(m,n)), NRHS)), INFO(1);
  double *WORK(new double[LWORK]);
  dgelss_(&m, &n, &NRHS, array, &LDA, B.array, &LDB, S.array, &RCOND, &RANK, WORK, &LWORK, &INFO);
  delete [] WORK;
  
  if(INFO!=0){
    WARNING_REPORT;
    std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
  }
  return INFO;
}

CPPL_INT dgematrix::dgelss ( dgematrix &  B, dcovector &  S, CPPL_INT &  RANK, const double  RCOND = -1.  )

inline

calculate the least-squares-least-norm solution for overdetermined or underdetermined A*x=y using dgelss

Definition at line 323 of file dgematrix-lapack.hpp.

References array, dcovector::array, clear(), i(), m, n, dcovector::resize(), and swap().

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(m!=B.m){
    ERROR_REPORT;
    std::cerr << "These matrix and vector cannot be solved." << std::endl
              << "Your input was (" << m << "x" << n << ") and (" << B.m << "x" << B.n  << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG    
  
  if(m<n){ //underdetermined
    dgematrix tmp(n,B.n);
    for(CPPL_INT i=0; i<B.m; i++){
      for(CPPL_INT j=0; j<B.n; j++){
        tmp(i,j)=B(i,j);
      }
    }
    B.clear();
    swap(B,tmp);
  }
  
  S.resize(std::min(m,n));
  
  CPPL_INT NRHS(B.n), LDA(m), LDB(B.m), LWORK(3*std::min(m,n)+std::max(std::max(2*std::min(m,n),std::max(m,n)), NRHS)), INFO(1);
  double *WORK(new double[LWORK]);
  dgelss_(&m, &n, &NRHS, array, &LDA, B.array, &LDB, S.array, &RCOND, &RANK, WORK, &LWORK, &INFO);
  delete [] WORK;

  if(INFO!=0){
    WARNING_REPORT;
    std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
  }
  return INFO;
}

CPPL_INT dgematrix::dgelsd ( dcovector &  B, dcovector &  S, CPPL_INT &  RANK, const double  RCOND = -1.  )

inline

calculate the least-squares-least-norm solution for overdetermined or underdetermined A*x=y using dgelsd

Definition at line 372 of file dgematrix-lapack.hpp.

References dcovector::array, dcovector::clear(), i(), dcovector::l, dcovector::resize(), and swap().

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(m!=B.l){
    ERROR_REPORT;
    std::cerr << "These matrix and vector cannot be solved." << std::endl
              << "Your input was (" << m << "x" << n << ") and (" << B.l << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG    
  
  //////// resize ////////
  if(m<n){ //underdetermined
    dcovector tmp(n);
    for(CPPL_INT i=0; i<B.l; i++){
      tmp(i)=B(i);
    }
    B.clear();
    swap(B,tmp);
  }
  
  S.resize(std::min(m,n));
  
  //////// prerun ////////
  CPPL_INT M =m;
  CPPL_INT N =n;
  CPPL_INT NRHS =1;
  CPPL_INT LDA =m;
  CPPL_INT LDB =B.l;
  std::vector<double> WORK(1);
  CPPL_INT LWORK =-1;
  CPPL_INT MINMN =std::min(M,N);
  CPPL_INT SMLSIZ =25;
  CPPL_INT NLVL =CPPL_INT( std::log(double(MINMN)/double(SMLSIZ+1))/std::log(2.) );
  CPPL_INT LIWORK =3*MINMN*NLVL+11*MINMN;
  std::vector<CPPL_INT> IWORK(LIWORK);
  CPPL_INT INFO =1;
  dgelsd_(&M, &N, &NRHS, array, &LDA, B.array, &LDB, S.array, &RCOND, &RANK, &WORK[0], &LWORK, &IWORK[0], &INFO);
  LWORK =CPPL_INT(WORK[0]);
  WORK.resize(LWORK);
  
  //////// run ////////
  INFO =1;
  RANK =1;
  dgelsd_(&M, &N, &NRHS, array, &LDA, B.array, &LDB, S.array, &RCOND, &RANK, &WORK[0], &LWORK, &IWORK[0], &INFO);
  
  //////// info ////////
  if(INFO!=0){
    WARNING_REPORT;
    std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
  }
  return INFO;
}

CPPL_INT dgematrix::dgeev ( std::vector< double > &  wr, std::vector< double > &  wi  )

inline

calculate eigenvalues
All of the arguments need not to be initialized. wr and wi are overwitten and become real and imaginary part of eigenvalues, respectively. This matrix is also overwritten.

Definition at line 442 of file dgematrix-lapack.hpp.

References array, m, and n.

Referenced by dgeev().

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(m!=n){
    ERROR_REPORT;
    std::cerr << "This matrix is not a square matrix." << std::endl
              << "This matrix is (" << m << "x" << n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  wr.resize(n);
  wi.resize(n);
  char JOBVL('n'), JOBVR('n');
  CPPL_INT LDA(n), LDVL(1), LDVR(1), LWORK(3*n), INFO(1);
  double *VL(NULL), *VR(NULL), *WORK(new double[LWORK]);
  dgeev_(&JOBVL, &JOBVR, &n, array, &LDA, &wr[0], &wi[0], VL, &LDVL, VR, &LDVR, WORK, &LWORK, &INFO);
  delete [] WORK;
  delete [] VL;
  delete [] VR;
  
  if(INFO!=0){
    WARNING_REPORT;
    std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
  }
  return INFO;
}

CPPL_INT dgematrix::dgeev ( zcovector &  w )

inline

calculate eigenvalues
All of the arguments need not to be initialized. w are overwitten and become eigenvalues, respectively. This matrix is also overwritten.

Definition at line 476 of file dgematrix-lapack.hpp.

References dgeev(), i, n, and zcovector::resize().

{CPPL_VERBOSE_REPORT;
  //// call dgeev ////
  std::vector<double> wr, wi;
  CPPL_INT INFO =dgeev(wr,wi);
  
  //// assign ////
  w.resize(n);
  for(CPPL_INT i=0; i<n; i++){
    w(i) =comple(wr[i],wi[i]);
  }
  
  return INFO;
}

CPPL_INT dgematrix::dgeev ( std::vector< double > &  wr, std::vector< double > &  wi, std::vector< dcovector > &  vrr, std::vector< dcovector > &  vri  )

inline

calculate right eigenvalues and right eigenvectors
All of the arguments need not to be initialized. wr, wi, vrr, vri are overwitten and become real and imaginary part of right eigenvalues and right eigenvectors, respectively. This matrix is also overwritten.

Definition at line 500 of file dgematrix-lapack.hpp.

References array, and i().

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(m!=n){
    ERROR_REPORT;
    std::cerr << "This matrix is not a square matrix." << std::endl
              << "This matrix is (" << m << "x" << n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  wr.resize(n);
  wi.resize(n);
  vrr.resize(n);
  vri.resize(n);
  for(CPPL_INT i=0; i<n; i++){
    vrr[i].resize(n);
    vri[i].resize(n);
  }
  
  dgematrix VR(n,n);
  char JOBVL('n'), JOBVR('V');
  CPPL_INT LDA(n), LDVL(1), LDVR(n), LWORK(4*n), INFO(1);
  double *VL(NULL), *WORK(new double[LWORK]);
  dgeev_(&JOBVL, &JOBVR, &n, array, &LDA, &wr[0], &wi[0], VL, &LDVL, VR.array, &LDVR, WORK, &LWORK, &INFO);
  delete [] WORK;
  delete [] VL;
  
  //// forming ////
  for(CPPL_INT j=0; j<n; j++){
    if(fabs(wi[j])<DBL_MIN){
      for(CPPL_INT i=0; i<n; i++){
        vrr[j](i) = VR(i,j);
        vri[j](i) = 0.0;
      }
    }
    else{
      for(CPPL_INT i=0; i<n; i++){
        vrr[j](i)   = VR(i,j);
        vri[j](i)   = VR(i,j+1);
        vrr[j+1](i) = VR(i,j);
        vri[j+1](i) =-VR(i,j+1);
      }
      j++;
    }
  }
  
  if(INFO!=0){
    WARNING_REPORT;
    std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
  }
  return INFO;
}

CPPL_INT dgematrix::dgeev ( std::vector< double > &  wr, std::vector< double > &  wi, std::vector< drovector > &  vlr, std::vector< drovector > &  vli  )

inline

calculate left eigenvalues and left eigenvectors
All of the arguments need not to be initialized. wr, wi, vrr, vri are overwitten and become real and imaginary part of left eigenvalues and left eigenvectors, respectively. This matrix is also overwritten.

Definition at line 567 of file dgematrix-lapack.hpp.

References array, i, m, and n.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(m!=n){
    ERROR_REPORT;
    std::cerr << "This matrix is not a square matrix." << std::endl
              << "This matrix is (" << m << "x" << n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  wr.resize(n);
  wi.resize(n);
  vlr.resize(n);
  vli.resize(n);
  for(CPPL_INT i=0; i<n; i++){
    vlr[i].resize(n);
    vli[i].resize(n);
  }
  
  dgematrix VL(n,n);
  char JOBVL('V'), JOBVR('n');
  CPPL_INT LDA(n), LDVL(n), LDVR(1), LWORK(4*n), INFO(1);
  double *VR(NULL), *WORK(new double[LWORK]);
  dgeev_(&JOBVL, &JOBVR, &n, array, &LDA, &wr[0], &wi[0], VL.array, &LDVL, VR, &LDVR, WORK, &LWORK, &INFO);
  delete [] WORK;
  delete [] VR;

  //// forming ////
  for(CPPL_INT j=0; j<n; j++){
    if(fabs(wi[j])<DBL_MIN){
      for(CPPL_INT i=0; i<n; i++){
        vlr[j](i) = VL(i,j);
        vli[j](i) = 0.0;
      }
    }
    else{
      for(CPPL_INT i=0; i<n; i++){
        vlr[j](i)   = VL(i,j);
        vli[j](i)   =-VL(i,j+1);
        vlr[j+1](i) = VL(i,j);
        vli[j+1](i) = VL(i,j+1);
      }
      j++;
    }
  }
  

  if(INFO!=0){
    WARNING_REPORT;
    std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
  }
  return INFO;
}

CPPL_INT dgematrix::dggev ( dgematrix &  matB, std::vector< double > &  wr, std::vector< double > &  wi  )

inline

calculate generalized eigenvalues
All of the arguments don't need to be initialized. wr and wi are overwitten and become real and imaginary part of generalized eigenvalues, respectively. This matrix and matB are also overwritten.

Definition at line 637 of file dgematrix-lapack.hpp.

References array, i, m, and n.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(m!=n){
    ERROR_REPORT;
    std::cerr << "This matrix is not a square matrix." << std::endl
              << "This matrix is (" << m << "x" << n << ")." << std::endl;
    exit(1);
  }
  if(matB.m!=n || matB.n!=n){
    ERROR_REPORT;
    std::cerr << "The matrix B is not a square matrix having the same size as \"this\" matrix." << std::endl
              << "The B matrix is (" << matB.m << "x" << matB.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  wr.resize(n);
  wi.resize(n);
  char JOBVL('n'), JOBVR('n');
  CPPL_INT LDA(n), LDB(n), LDVL(1), LDVR(1), LWORK(8*n), INFO(1);
  double *BETA(new double[n]), *VL(NULL), *VR(NULL), *WORK(new double[LWORK]);
  dggev_(&JOBVL, &JOBVR, &n, array, &LDA, matB.array, &LDB, &wr[0], &wi[0], BETA, VL, &LDVL, VR, &LDVR, WORK, &LWORK, &INFO);
  delete [] WORK;
  delete [] VL;
  delete [] VR;
  
  //// reforming ////
  for(CPPL_INT i=0; i<n; i++){
    wr[i]/=BETA[i];
    wi[i]/=BETA[i];
  }
  delete [] BETA; 
  
  if(INFO!=0){
    WARNING_REPORT;
    std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
  }
  return INFO;
}

CPPL_INT dgematrix::dggev ( dgematrix &  matB, std::vector< double > &  wr, std::vector< double > &  wi, std::vector< dcovector > &  vrr, std::vector< dcovector > &  vri  )

inline

calculate generalized eigenvalues and generalized right eigenvectors
All of the arguments don't need to be initialized. wr, wi, vrr and vri are overwitten and become real and imaginary part of generalized eigenvalue and generalized right eigenvector, respectively. This matrix and matB are also overwritten.

Definition at line 687 of file dgematrix-lapack.hpp.

References array, i, m, and n.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(m!=n){
    ERROR_REPORT;
    std::cerr << "This matrix is not a square matrix." << std::endl
              << "This matrix is (" << m << "x" << n << ")." << std::endl;
    exit(1);
  }
  if(matB.m!=n || matB.n!=n){
    ERROR_REPORT;
    std::cerr << "The matrix B is not a square matrix having the same size as \"this\" matrix." << std::endl
              << "The B matrix is (" << matB.m << "x" << matB.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  wr.resize(n);
  wi.resize(n);
  vrr.resize(n);
  vri.resize(n);
  for(CPPL_INT i=0; i<n; i++){
    vrr[i].resize(n);
    vri[i].resize(n);
  }
  
  dgematrix VR(n,n);
  char JOBVL('n'), JOBVR('V');
  CPPL_INT LDA(n), LDB(n), LDVL(1), LDVR(n), LWORK(8*n), INFO(1);
  double *BETA(new double[n]), *VL(NULL), *WORK(new double[LWORK]);
  dggev_(&JOBVL, &JOBVR, &n, array, &LDA, matB.array, &LDB, &wr[0], &wi[0], BETA, VL, &LDVL, VR.array, &LDVR, WORK, &LWORK, &INFO);
  delete [] WORK;
  delete [] VL;
  
  //// reforming ////
  for(CPPL_INT i=0; i<n; i++){
    wr[i]/=BETA[i];
    wi[i]/=BETA[i];
  }
  delete [] BETA; 
  
  //// forming ////
  for(CPPL_INT j=0; j<n; j++){
    if(fabs(wi[j])<DBL_MIN){
      for(CPPL_INT i=0; i<n; i++){
        vrr[j](i) = VR(i,j);
        vri[j](i) = 0.0;
      }
    }
    else{
      for(CPPL_INT i=0; i<n; i++){
        vrr[j](i)   = VR(i,j);
        vri[j](i)   = VR(i,j+1);
        vrr[j+1](i) = VR(i,j);
        vri[j+1](i) =-VR(i,j+1);
      }
      j++;
    }
  }
  
  if(INFO!=0){
    WARNING_REPORT;
    std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
  }
  return INFO;
}

CPPL_INT dgematrix::dggev ( dgematrix &  matB, std::vector< double > &  wr, std::vector< double > &  wi, std::vector< drovector > &  vlr, std::vector< drovector > &  vli  )

inline

calculate generalized eigenvalues and generalized left eigenvectors
All of the arguments don't need to be initialized. wr, wi, vlr and vli are overwitten and become real and imaginary part of generalized eigenvalue and generalized left eigenvector, respectively. This matrix and matB are also overwritten.

Definition at line 765 of file dgematrix-lapack.hpp.

References array, i, m, and n.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(m!=n){
    ERROR_REPORT;
    std::cerr << "This matrix is not a square matrix." << std::endl
              << "This matrix is (" << m << "x" << n << ")." << std::endl;
    exit(1);
  }
  if(matB.m!=n || matB.n!=n){
    ERROR_REPORT;
    std::cerr << "The matrix B is not a square matrix having the same size as \"this\" matrix." << std::endl
              << "The B matrix is (" << matB.m << "x" << matB.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  wr.resize(n);
  wi.resize(n);
  vlr.resize(n);
  vli.resize(n);
  for(CPPL_INT i=0; i<n; i++){
    vlr[i].resize(n);
    vli[i].resize(n);
  }
  
  dgematrix VL(n,n);
  char JOBVL('V'), JOBVR('n');
  CPPL_INT LDA(n), LDB(n), LDVL(n), LDVR(1), LWORK(8*n), INFO(1);
  double *BETA(new double[n]), *VR(NULL), *WORK(new double[LWORK]);
  dggev_(&JOBVL, &JOBVR, &n, array, &LDA, matB.array, &LDB, &wr[0], &wi[0], BETA, VL.array, &LDVL, VR, &LDVR, WORK, &LWORK, &INFO);
  delete [] WORK;
  delete [] VR;
  
  //// reforming ////
  for(CPPL_INT i=0; i<n; i++){
    wr[i]/=BETA[i];
    wi[i]/=BETA[i];
  }
  delete [] BETA; 
  
  //// forming ////
  for(CPPL_INT j=0; j<n; j++){
    if(fabs(wi[j])<DBL_MIN){
      for(CPPL_INT i=0; i<n; i++){
        vlr[j](i) = VL(i,j);
        vli[j](i) = 0.0;
      }
    }
    else{
      for(CPPL_INT i=0; i<n; i++){
        vlr[j](i)   = VL(i,j);
        vli[j](i)   =-VL(i,j+1);
        vlr[j+1](i) = VL(i,j);
        vli[j+1](i) = VL(i,j+1);
      }
      j++;
    }
  }
  
  if(INFO!=0){
    WARNING_REPORT;
    std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
  }
  return INFO;
}

CPPL_INT dgematrix::dgesvd ( dgbmatrix &  S )

inline

compute the singular value decomposition (SVD)
The argument is dgbmatrix S. S doesn't need to be initialized. S is overwitten and become singular values. This matrix also overwritten.

Definition at line 847 of file dgematrix-lapack.hpp.

References dgbmatrix::array, and dgbmatrix::resize().

{CPPL_VERBOSE_REPORT;
  char JOBU ='N';
  char JOBVT ='N';
  // M
  // N
  // A
  CPPL_INT LDA =m;
  double* U =NULL;
  S.resize(m,n,0,0);
  CPPL_INT LDU =1;
  double* VT =NULL;
  CPPL_INT LDVT =1;
  CPPL_INT LWORK =std::max(3*std::min(m,n)+std::max(m,n),5*std::min(m,n));
  double *WORK =new double[LWORK];
  CPPL_INT INFO =1;
  
  dgesvd_(&JOBU, &JOBVT, &m, &n, array, &LDA, S.array, U, &LDU, VT, &LDVT, WORK, &LWORK, &INFO);
  delete [] WORK;
  
  if(INFO!=0){
    WARNING_REPORT;
    std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
  }
  return INFO;
}

CPPL_INT dgematrix::dgesvd ( dcovector &  S, dgematrix &  U, dgematrix &  VT  )

inline

compute the singular value decomposition (SVD)
The arguments are dcocector S, dgematrix U and VT. All of them need not to be initialized. S, U and VT are overwitten and become singular values, left singular vectors, and right singular vectors respectively. This matrix also overwritten.

Definition at line 885 of file dgematrix-lapack.hpp.

References array, dcovector::array, m, n, resize(), and dcovector::resize().

{CPPL_VERBOSE_REPORT;
  char JOBU('A'), JOBVT('A');
  CPPL_INT LDA(m), LDU(m), LDVT(n), LWORK(std::max(3*std::min(m,n)+std::max(m,n),5*std::min(m,n))), INFO(1);
  double *WORK(new double[LWORK]);
  S.resize(std::min(m,n));
  U.resize(LDU,m);
  VT.resize(LDVT,n);
  
  dgesvd_(&JOBU, &JOBVT, &m, &n, array, &LDA, S.array, U.array, &LDU, VT.array, &LDVT, WORK, &LWORK, &INFO);
  delete [] WORK;
  
  if(INFO!=0){
    WARNING_REPORT;
    std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
  }
  return INFO;
}

CPPL_INT dgematrix::dgglse ( dgematrix &  B, dcovector &  c, dcovector &  d, dcovector &  x  )

inline

solve the linear equality-constrained least squares (LSE) problem
Input matrix and vectors, B, c, and d, are overwitten. This matrix is also overwritten. The solution vector x is to be automatically resized.

Definition at line 915 of file dgematrix-lapack.hpp.

References array, dcovector::array, dcovector::clear(), dcovector::l, m, and dcovector::resize().

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(m!=c.l){
    ERROR_REPORT;
    std::cerr << "A.m and c.l should be the same." << std::endl
              << "Your input was A.m=" << m << " and c.l=" << c.l << std::endl;
    exit(1);
  }
  if(B.m!=d.l){
    ERROR_REPORT;
    std::cerr << "B.m and d.l should be the same." << std::endl
              << "Your input was B.m=" << B.m << " and d.l=" << d.l << std::endl;
    exit(1);
  }
  if( !(B.m<=n) || !(n<=m+B.m) ){
    ERROR_REPORT;
    std::cerr << "B.m<=A.n<=A.m+B.m should be satisfied." << std::endl
              << "Your input was B.m=" << B.m << ", A.n=" << n << ", and A.m+B.m=" << m+B.m << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  CPPL_INT lwork(-1), info(1);
  dcovector work(1);
  x.resize(n);
  
  //////// workspace query ////////
  CPPL_INT lda =std::max(CPPL_INT(1),m);
  CPPL_INT ldb =std::max(CPPL_INT(1),B.m);
  dgglse_(&m, &n, &B.m, array, &lda, B.array, &ldb, c.array, d.array, x.array, work.array, &lwork, &info);
  lwork =CPPL_INT(work(0));
  work.resize(lwork);
  info =1;
  
  //////// solve ////////
  dgglse_(&m, &n, &B.m, array, &lda, B.array, &ldb, c.array, d.array, x.array, work.array, &lwork, &info);
  work.clear();
  
  if(info!=0){
    WARNING_REPORT;
    std::cerr << "Serious trouble happend. info = " << info << "." << std::endl;
  }
  return info;
}

dgematrix & dgematrix::operator= ( const dgematrix &  mat )

inline

dgematrix=dgematrix operator

Definition at line 3 of file dgematrix-dgematrix.hpp.

References array, and copy().

{CPPL_VERBOSE_REPORT;
  if(array!=mat.array){ // if it is NOT self substitution
    copy(mat);
  }
  return *this;
}

dgematrix & dgematrix::operator= ( const _dgematrix &  mat )

inline

dgematrix=_dgematrix operator

Definition at line 3 of file dgematrix-_dgematrix.hpp.

References shallow_copy().

{CPPL_VERBOSE_REPORT;
  shallow_copy(mat);
  return *this;
}

dgematrix & dgematrix::operator+= ( const dgematrix &  mat )

inline

dgematrix+=dgematrix operator

Definition at line 17 of file dgematrix-dgematrix.hpp.

References array, i, m, and n.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(n!=mat.n || m!=mat.m){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a summation." << std::endl
              << "Your input was (" << m << "x" << n << ") += (" << mat.m << "x" << mat.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  const CPPL_INT mn =m*n;
  for(CPPL_INT i=0; i<mn; i++){
    array[i]+=mat.array[i];
  }
  return *this;
}

dgematrix & dgematrix::operator+= ( const _dgematrix &  mat )

inline

dgematrix+=_dgematrix operator

Definition at line 15 of file dgematrix-_dgematrix.hpp.

References array, _dgematrix::array, _dgematrix::destroy(), i, _dgematrix::m, m, n, and _dgematrix::n.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(n!=mat.n || m!=mat.m){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a summation." << std::endl
              << "Your input was (" << m << "x" << n << ") += (" << mat.m << "x" << mat.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  const CPPL_INT mn =m*n;
  for(CPPL_INT i=0; i<mn; i++){
    array[i]+=mat.array[i];
  }
  
  mat.destroy();
  return *this;
}

dgematrix & dgematrix::operator+= ( const dsymatrix &  mat )

inline

dgematrix+=dsymatrix operator

Definition at line 3 of file dgematrix-dsymatrix.hpp.

References i, m, dsymatrix::n, n, and operator()().

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(n!=mat.n || m!=mat.n){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a summation." << std::endl
              << "Your input was (" << m << "x" << n << ") += (" << mat.n << "x" << mat.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  for(CPPL_INT i=0; i<m; i++){
    for( CPPL_INT j=0; j<n; j++){
      operator() (i,j) += mat(i,j);
    }
  }
  
  return *this;
}

dgematrix & dgematrix::operator+= ( const _dsymatrix &  mat )

inline

dgematrix+=_dsymatrix operator

Definition at line 3 of file dgematrix-_dsymatrix.hpp.

References _dsymatrix::destroy(), i, m, n, _dsymatrix::n, and operator()().

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(n!=mat.n || m!=mat.n){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a summation." << std::endl
              << "Your input was (" << m << "x" << n << ") += (" << mat.n << "x" << mat.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  for(CPPL_INT i=0; i<m; i++){
    for(CPPL_INT j=0; j<n; j++){
      operator()(i,j) += mat(i,j);
    }
  }
  
  mat.destroy();
  return *this;
}

dgematrix & dgematrix::operator+= ( const dgbmatrix &  mat )

inline

dgematrix+=dgbmatrix operator

Definition at line 3 of file dgematrix-dgbmatrix.hpp.

References i, dgbmatrix::kl, dgbmatrix::ku, m, dgbmatrix::m, dgbmatrix::n, n, and operator()().

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(n!=mat.n || m!=mat.m){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a summation." << std::endl
              << "Your input was (" << m << "x" << n << ") += (" << mat.m << "x" << mat.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  for(CPPL_INT i=0; i<mat.m; i++){
    const CPPL_INT jmax =std::min(n,i+mat.ku+1);
    for(CPPL_INT j=std::max(CPPL_INT(0),i-mat.kl); j<jmax; j++){
      operator()(i,j) += mat(i,j);
    }
  }
  
  return *this;
}

dgematrix & dgematrix::operator+= ( const _dgbmatrix &  mat )

inline

dgematrix+=_dgbmatrix operator

Definition at line 3 of file dgematrix-_dgbmatrix.hpp.

References _dgbmatrix::destroy(), i, _dgbmatrix::kl, _dgbmatrix::ku, m, _dgbmatrix::m, _dgbmatrix::n, n, and operator()().

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(n!=mat.n || m!=mat.m){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a summation." << std::endl
              << "Your input was (" << m << "x" << n << ") += (" << mat.m << "x" << mat.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  for(CPPL_INT i=0; i<mat.m; i++){
    const CPPL_INT jmax =std::min(n,i+mat.ku+1);
    for(CPPL_INT j=std::max(CPPL_INT(0),i-mat.kl); j<jmax; j++){
      operator()(i,j) += mat(i,j);
    }
  }
  
  mat.destroy();
  return *this;
}

dgematrix& dgematrix::operator+= ( const dgsmatrix &  )

inline

dgematrix& dgematrix::operator+= ( const _dgsmatrix &  )

inline

dgematrix& dgematrix::operator+= ( const dssmatrix &  )

inline

dgematrix& dgematrix::operator+= ( const _dssmatrix &  )

inline

dgematrix & dgematrix::operator-= ( const dgematrix &  mat )

inline

dgematrix operator-=

Definition at line 37 of file dgematrix-dgematrix.hpp.

References array, i, m, and n.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(n!=mat.n || m!=mat.m){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a sutraction." << std::endl
              << "Your input was (" << m << "x" << n << ") -= (" << mat.m << "x" << mat.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  const CPPL_INT mn =m*n;
  for(CPPL_INT i=0; i<mn; i++){
    array[i]-=mat.array[i];
  }
  return *this;
}

dgematrix & dgematrix::operator-= ( const _dgematrix &  mat )

inline

dgematrix-=_dgematrix operator

Definition at line 37 of file dgematrix-_dgematrix.hpp.

References array, _dgematrix::array, _dgematrix::destroy(), i, _dgematrix::m, m, n, and _dgematrix::n.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(n!=mat.n || m!=mat.m){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a sutraction." << std::endl
              << "Your input was (" << m << "x" << n << ") -= (" << mat.m << "x" << mat.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  const CPPL_INT mn =m*n;
  for(CPPL_INT i=0; i<mn; i++){
    array[i]-=mat.array[i];
  }

  mat.destroy();
  return *this;
}

dgematrix & dgematrix::operator-= ( const dsymatrix &  mat )

inline

dgematrix-=dsymatrix operator

Definition at line 25 of file dgematrix-dsymatrix.hpp.

References i, m, dsymatrix::n, n, and operator()().

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(n!=mat.n || m!=mat.n){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a sutraction." << std::endl
              << "Your input was (" << m << "x" << n << ") -= (" << mat.n << "x" << mat.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  for(CPPL_INT i=0; i<m; i++){
    for(CPPL_INT j=0; j<n; j++){
      operator() (i,j) -= mat(i,j);
    }
  }
  
  return *this;
}

dgematrix & dgematrix::operator-= ( const _dsymatrix &  mat )

inline

dgematrix-=_dsymatrix operator

Definition at line 26 of file dgematrix-_dsymatrix.hpp.

References _dsymatrix::destroy(), i, m, n, _dsymatrix::n, and operator()().

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(n!=mat.n || m!=mat.n){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a sutraction." << std::endl
              << "Your input was (" << m << "x" << n << ") -= (" << mat.n << "x" << mat.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  for(CPPL_INT i=0; i<m;i++){
    for(CPPL_INT j=0; j<n; j++){
      operator()(i,j) -= mat(i,j);
    }
  }
  
  mat.destroy();
  return *this;
}

dgematrix & dgematrix::operator-= ( const dgbmatrix &  mat )

inline

dgematrix-=dgbmatrix operator

Definition at line 26 of file dgematrix-dgbmatrix.hpp.

References i, dgbmatrix::kl, dgbmatrix::ku, m, dgbmatrix::m, dgbmatrix::n, n, and operator()().

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(n!=mat.n || m!=mat.m){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a subtraction." << std::endl
              << "Your input was (" << m << "x" << n << ") -= (" << mat.m << "x" << mat.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  for(CPPL_INT i=0; i<mat.m; i++){
    const CPPL_INT jmax =std::min(n,i+mat.ku+1);
    for(CPPL_INT j=std::max(CPPL_INT(0),i-mat.kl); j<jmax; j++){
      operator()(i,j) -= mat(i,j);
    }
  }
  
  return *this;
}

dgematrix & dgematrix::operator-= ( const _dgbmatrix &  mat )

inline

dgematrix-=_dgbmatrix operator

Definition at line 27 of file dgematrix-_dgbmatrix.hpp.

References _dgbmatrix::destroy(), i, _dgbmatrix::kl, _dgbmatrix::ku, m, _dgbmatrix::m, _dgbmatrix::n, n, and operator()().

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(n!=mat.n || m!=mat.m){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a subtraction." << std::endl
              << "Your input was (" << m << "x" << n << ") -= (" << mat.m << "x" << mat.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  for(CPPL_INT i=0; i<mat.m; i++){
    const CPPL_INT jmax =std::min(n,i+mat.ku+1);
    for(CPPL_INT j=std::max(CPPL_INT(0),i-mat.kl); j<jmax; j++){
      operator()(i,j) -= mat(i,j);
    }
  }
  
  mat.destroy();
  return *this;
}

dgematrix& dgematrix::operator-= ( const dgsmatrix &  )

inline

dgematrix& dgematrix::operator-= ( const _dgsmatrix &  )

inline

dgematrix& dgematrix::operator-= ( const dssmatrix &  )

inline

dgematrix& dgematrix::operator-= ( const _dssmatrix &  )

inline

dgematrix & dgematrix::operator*= ( const dgematrix &  mat )

inline

dgematrix operator*=

Definition at line 57 of file dgematrix-dgematrix.hpp.

References array, m, n, and swap.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(n!=mat.m){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a product." << std::endl
              << "Your input was (" << m << "x" << n << ") *= (" << mat.m << "x" << mat.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  dgematrix newmat( m, mat.n );
  char transa ='n';
  char transb ='n';
  double alpha =1.;
  double beta =0.;
  
  dgemm_( &transa, &transb, &m, &mat.n, &n, &alpha, array, &m, mat.array, &mat.m, &beta, newmat.array, &m );
  
  swap(*this,newmat);
  return *this;
}

dgematrix & dgematrix::operator*= ( const _dgematrix &  mat )

inline

dgematrix*=_dgematrix operator

Definition at line 59 of file dgematrix-_dgematrix.hpp.

References array, _dgematrix::array, _dgematrix::destroy(), _dgematrix::m, m, _dgematrix::n, n, and swap.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(n!=mat.m){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a product." << std::endl
              << "Your input was (" << m << "x" << n << ") *= (" << mat.m << "x" << mat.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  dgematrix newmat( m, mat.n );
  char transa ='n';
  char transb ='n';
  double alpha =1.;
  double beta =0.;
  
  dgemm_( &transa, &transb, &m, &mat.n, &n, &alpha, array, &m, mat.array, &mat.m, &beta, newmat.array, &m );
  
  swap(*this,newmat);
  mat.destroy();
  return *this;
}

dgematrix & dgematrix::operator*= ( const dsymatrix &  mat )

inline

dgematrix*=dsymatrix operator

Definition at line 47 of file dgematrix-dsymatrix.hpp.

References dsymatrix::array, array, m, dsymatrix::n, n, and swap.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(n!=mat.n){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a product." << std::endl
              << "Your input was (" << m << "x" << n << ") *= (" << mat.n << "x" << mat.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  dgematrix newmat( m, mat.n );
  char side ='R';
  char uplo ='l';
  double alpha =1.;
  double beta =0.;
  
  dsymm_( &side, &uplo, &mat.n, &n, &alpha, mat.array, &mat.n, array, &m, &beta, newmat.array, &newmat.m );
  
  swap(*this,newmat);
  return *this;
}

dgematrix & dgematrix::operator*= ( const _dsymatrix &  mat )

inline

dgematrix*=_dsymatrix operator

Definition at line 49 of file dgematrix-_dsymatrix.hpp.

References array, _dsymatrix::array, _dsymatrix::destroy(), m, n, _dsymatrix::n, and swap.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(n!=mat.n){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a product." << std::endl
              << "Your input was (" << m << "x" << n << ") *= (" << mat.n << "x" << mat.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  dgematrix newmat( m, mat.n );
  char side ='R';
  char uplo ='l';
  double alpha =1.;
  double beta =0.;
  
  dsymm_( &side, &uplo, &mat.n, &n, &alpha, mat.array, &mat.n, array, &m, &beta, newmat.array, &newmat.m );
  
  swap(*this,newmat);
  mat.destroy();
  return *this;
}

dgematrix & dgematrix::operator*= ( const dgbmatrix &  mat )

inline

dgematrix*=dgbmatrix operator

Definition at line 49 of file dgematrix-dgbmatrix.hpp.

References i, dgbmatrix::kl, dgbmatrix::ku, m, dgbmatrix::m, n, dgbmatrix::n, operator()(), swap, and zero().

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(n!=mat.m){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a product." << std::endl
              << "Your input was (" << m << "x" << n << ") *= (" << mat.m << "x" << mat.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  dgematrix newmat(m,mat.n);
  newmat.zero();
  
  for(CPPL_INT i=0; i<newmat.m; i++){
    for(CPPL_INT j=0; j<newmat.n; j++){
      const CPPL_INT kmax =std::min(mat.m,j+mat.kl+1);
      for(CPPL_INT k=std::max(CPPL_INT(0),j-mat.ku); k<kmax; k++){
        newmat(i,j) += operator()(i,k)*mat(k,j);
      }
    }
  }
  
  swap(*this,newmat);
  return *this;
}

dgematrix & dgematrix::operator*= ( const _dgbmatrix &  mat )

inline

dgematrix*=_dgbmatrix operator

Definition at line 50 of file dgematrix-_dgbmatrix.hpp.

References _dgbmatrix::destroy(), i, _dgbmatrix::kl, _dgbmatrix::ku, m, _dgbmatrix::m, n, _dgbmatrix::n, operator()(), swap, and zero().

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(n!=mat.m){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a product." << std::endl
              << "Your input was (" << m << "x" << n << ") *= (" << mat.m << "x" << mat.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  dgematrix newmat(m,mat.n);
  newmat.zero();
  
  for(CPPL_INT i=0; i<newmat.m; i++){
    for(CPPL_INT j=0; j<newmat.n; j++){
      const CPPL_INT kmax =std::min(mat.m,j+mat.kl+1);
      for(CPPL_INT k=std::max(CPPL_INT(0),j-mat.ku); k<kmax; k++){
        newmat(i,j) += operator()(i,k)*mat(k,j);
      }
    }
  }
  
  swap(*this, newmat);
  mat.destroy();
  return *this;
}

dgematrix& dgematrix::operator*= ( const dgsmatrix &  )

inline

dgematrix& dgematrix::operator*= ( const _dgsmatrix &  )

inline

dgematrix& dgematrix::operator*= ( const dssmatrix &  )

inline

dgematrix& dgematrix::operator*= ( const _dssmatrix &  )

inline

dgematrix & dgematrix::operator*= ( const double &  d )

inline

dgematrix*=double operator

Definition at line 3 of file dgematrix-double.hpp.

References array, m, and n.

{CPPL_VERBOSE_REPORT;
  CPPL_INT mn =m*n;
  CPPL_INT inc =1;
  dscal_(&mn, &d, array, &inc);
  return *this;
}

dgematrix & dgematrix::operator/= ( const double &  d )

inline

dgematrix/=double operator

Definition at line 13 of file dgematrix-double.hpp.

References array, m, and n.

{CPPL_VERBOSE_REPORT;
  CPPL_INT mn =m*n;
  double dinv =1./d;
  CPPL_INT inc =1;
  dscal_(&mn, &dinv, array, &inc);
  return *this;
}

Friends And Related Function Documentation

std::ostream& operator<< ( std::ostream &  s, const dgematrix &  mat  )

friend

Definition at line 63 of file dgematrix-io.hpp.

{CPPL_VERBOSE_REPORT;
  for(CPPL_INT i=0; i<mat.m; i++){
    for(CPPL_INT j=0; j<mat.n; j++){
      s << " " << mat(i,j);
    }
    s << std::endl;
  }
  return s;
}

void swap ( dgematrix &  A, dgematrix &  B  )

friend

swap two matrices

Definition at line 164 of file dgematrix-misc.hpp.

Referenced by dgels(), and operator*=().

{CPPL_VERBOSE_REPORT;
  CPPL_INT A_m =A.m, A_n =A.n;
  double* A_array =A.array;
  double** A_darray=A.darray;
  A.m=B.m; A.n=B.n; A.array=B.array; A.darray=B.darray;
  B.m=A_m; B.n=A_n; B.array=A_array; B.darray=A_darray;
}

_dgematrix _ ( dgematrix &  mat )

friend

convert user object to smart-temporary object

Definition at line 175 of file dgematrix-misc.hpp.

Referenced by col(), row(), and to_zgematrix().

{CPPL_VERBOSE_REPORT;
  _dgematrix newmat;
  
  //////// shallow copy ////////
  newmat.m =mat.m;
  newmat.n =mat.n;
  newmat.array =mat.array;
  newmat.darray =mat.darray;
  
  //////// nullify ////////
  mat.m =0;
  mat.n =0;
  mat.array =NULL;
  mat.darray =NULL;
  
  return newmat;
}

_dgematrix t ( const dgematrix &  mat )

friend

return transposed dgematrix

Definition at line 3 of file dgematrix-calc.hpp.

{CPPL_VERBOSE_REPORT;
  dgematrix newmat(mat.n,mat.m);
  
  for(CPPL_INT i=0; i<newmat.m; i++){
    for(CPPL_INT j=0; j<newmat.n; j++){
      newmat(i,j) =mat(j,i);
    }
  }
  
  return _(newmat);
}

_dgematrix i ( const dgematrix &  mat )

friend

return its inverse matrix

Definition at line 18 of file dgematrix-calc.hpp.

Referenced by chsign(), col(), copy(), dgeev(), dgels(), dgematrix(), dggev(), identity(), operator()(), operator*=(), operator+=(), operator-=(), read(), resize(), set(), to_zgematrix(), write(), and zero().

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(mat.m!=mat.n){
    ERROR_REPORT;
    std::cerr << "This matrix is not square and has no inverse matrix." << std::endl
              << "Your input was (" << mat.m << "x" << mat.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  dgematrix mat_cp(mat), mat_inv(mat.m,mat.n);
  mat_inv.identity();
  mat_cp.dgesv(mat_inv);
  
  return _(mat_inv);
}

void idamax ( CPPL_INT &  i, CPPL_INT &  j, const dgematrix &  mat  )

friend

search the index of element having the largest absolute value in 0-based numbering system

Definition at line 43 of file dgematrix-calc.hpp.

{CPPL_VERBOSE_REPORT;
  CPPL_INT mn =mat.m*mat.n;
  CPPL_INT inc =1;
  CPPL_INT index =idamax_(&mn, mat.array, &inc) -1;
  i =index%mat.m;
  j =index/mat.m;
}

double damax ( const dgematrix &  mat )

friend

return its largest absolute value

Definition at line 54 of file dgematrix-calc.hpp.

{CPPL_VERBOSE_REPORT;
  CPPL_INT mn =mat.m*mat.n;
  CPPL_INT inc =1;
  return mat.array[idamax_(&mn, mat.array, &inc) -1];
}

const dgematrix& operator+ ( const dgematrix &  mat )

friend

+dgematrix operator

Definition at line 3 of file dgematrix-unary.hpp.

{CPPL_VERBOSE_REPORT;
  return mat;
}

_dgematrix operator- ( const dgematrix &  mat )

friend

-dgematrix operator

Definition at line 10 of file dgematrix-unary.hpp.

{CPPL_VERBOSE_REPORT;
  dgematrix newmat(mat.m,mat.n);
  for(CPPL_INT i=0; i<newmat.m*newmat.n; i++){ newmat.array[i]=-mat.array[i]; }
  return _(newmat);
}

_dgematrix operator+ ( const dgematrix &  matA, const dgematrix &  matB  )

friend

dgematrix+dgematrix operator

Definition at line 86 of file dgematrix-dgematrix.hpp.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(matA.n!=matB.n || matA.m!=matB.m){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a summation." << std::endl
              << "Your input was (" << matA.m << "x" << matA.n << ") + (" << matB.m << "x" << matB.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG

  dgematrix newmat(matA.m,matA.n);
  
  const CPPL_INT mn =newmat.m*newmat.n;
  for(CPPL_INT i=0; i<mn; i++){
    newmat.array[i] =matA.array[i]+matB.array[i];
  }
  
  return _(newmat);
}

_dgematrix operator+ ( const dgematrix &  matA, const _dgematrix &  matB  )

friend

dgematrix+_dgematrix operator

Definition at line 89 of file dgematrix-_dgematrix.hpp.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(matA.n!=matB.n || matA.m!=matB.m){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a summation." << std::endl
              << "Your input was (" << matA.m << "x" << matA.n << ") + (" << matB.m << "x" << matB.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  const CPPL_INT mn =matA.m*matA.n;
  for(CPPL_INT i=0; i<mn; i++){
    matB.array[i] +=matA.array[i];
  }
  
  return matB;
}

_dgematrix operator+ ( const dgematrix &  matA, const dsymatrix &  matB  )

friend

dgematrix+dsymatrix operator

Definition at line 76 of file dgematrix-dsymatrix.hpp.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(matA.n!=matB.n || matA.m!=matB.n){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a summation." << std::endl
              << "Your input was (" << matA.m << "x" << matA.n << ") + (" << matB.n << "x" << matB.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG

  dgematrix newmat(matA);
  for(CPPL_INT i=0; i<matA.m; i++){
    for(CPPL_INT j=0; j<matA.n; j++){
      newmat(i,j) += matB(i,j);
    }
  }
  
  return _(newmat);
}

_dgematrix operator+ ( const dgematrix &  matA, const _dsymatrix &  matB  )

friend

dgematrix+_dsymatrix operator

Definition at line 79 of file dgematrix-_dsymatrix.hpp.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(matA.n!=matB.n || matA.m!=matB.n){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a summation." << std::endl
              << "Your input was (" << matA.m << "x" << matA.n << ") + (" << matB.n << "x" << matB.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG

  dgematrix newmat(matA);
  for(CPPL_INT i=0; i<matA.m; i++){
    for(CPPL_INT j=0; j<matA.n; j++){
      newmat(i,j) += matB(i,j);
    }
  }
  
  matB.destroy();
  return _(newmat);
}

_dgematrix operator+ ( const dgematrix &  matA, const dgbmatrix &  matB  )

friend

dgematrix+dgbmatrix operator

Definition at line 81 of file dgematrix-dgbmatrix.hpp.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(matA.n!=matB.n || matA.m!=matB.m){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a summation." << std::endl
              << "Your input was (" << matA.m << "x" << matA.n << ") + (" << matB.m << "x" << matB.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  dgematrix newmat(matA);
  
  for(CPPL_INT i=0; i<matB.m; i++){
    const CPPL_INT jmax =std::min(matB.n,i+matB.ku+1);
    for(CPPL_INT j=std::max(CPPL_INT(0),i-matB.kl); j<jmax; j++){
      newmat(i,j) += matB(i,j);
    }
  }
  
  return _(newmat);
}

_dgematrix operator+ ( const dgematrix &  matA, const _dgbmatrix &  matB  )

friend

dgematrix+_dgbmatrix operator

Definition at line 84 of file dgematrix-_dgbmatrix.hpp.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(matA.n!=matB.n || matA.m!=matB.m){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a summation." << std::endl
              << "Your input was (" << matA.m << "x" << matA.n << ") + (" << matB.m << "x" << matB.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  dgematrix newmat(matA);
  
  for(CPPL_INT i=0; i<matB.m; i++){
    const CPPL_INT jmax =std::min(matB.n,i+matB.ku+1);
    for(CPPL_INT j=std::max(CPPL_INT(0),i-matB.kl); j<jmax; j++){
      newmat(i,j) += matB(i,j);
    }
  }
  
  matB.destroy();
  return _(newmat);
}

_dgematrix operator+ ( const dgematrix &  matA, const dgsmatrix &  matB  )

friend

dgematrix+dgsmatrix operator

Definition at line 3 of file dgematrix-dgsmatrix.hpp.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(matA.m!=matB.m || matA.n!=matB.n){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a summation." << std::endl
              << "Your input was (" << matA.m << "x" << matA.n << ") + (" << matB.m << "x" << matB.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  dgematrix newmat(matA);
  
  const std::vector<dcomponent>::const_iterator matB_data_end =matB.data.end();
  for(std::vector<dcomponent>::const_iterator it=matB.data.begin(); it!=matB_data_end; it++){
    newmat(it->i,it->j) += it->v;
  }
  
  return _(newmat);
}

_dgematrix operator+ ( const dgematrix &  matA, const _dgsmatrix &  matB  )

friend

dgematrix+_dgsmatrix operator

Definition at line 3 of file dgematrix-_dgsmatrix.hpp.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(matA.m!=matB.m || matA.n!=matB.n){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a summation." << std::endl
              << "Your input was (" << matA.m << "x" << matA.n << ") + (" << matB.m << "x" << matB.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  dgematrix newmat(matA);
  
  const std::vector<dcomponent>::const_iterator matB_data_end =matB.data.end();
  for(std::vector<dcomponent>::const_iterator it=matB.data.begin(); it!=matB_data_end; it++){
    newmat(it->i,it->j) += it->v;
  }
  
  matB.destroy();
  return _(newmat);
}

_dgematrix operator+ ( const dgematrix &  , const dssmatrix &    )

friend

_dgematrix operator+ ( const dgematrix &  , const _dssmatrix &    )

friend

_dgematrix operator- ( const dgematrix &  matA, const dgematrix &  matB  )

friend

dgematrix-dgematrix operator

Definition at line 109 of file dgematrix-dgematrix.hpp.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(matA.n!=matB.n || matA.m!=matB.m){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a subtraction." << std::endl
              << "Your input was (" << matA.m << "x" << matA.n << ") - (" << matB.m << "x" << matB.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG

  dgematrix newmat(matA.m,matA.n);
  
  const CPPL_INT mn =newmat.m*newmat.n;
  for(CPPL_INT i=0; i<mn; i++){
    newmat.array[i] =matA.array[i]-matB.array[i];
  }
  
  return _(newmat);
}

_dgematrix operator- ( const dgematrix &  matA, const _dgematrix &  matB  )

friend

dgematrix-_dgematrix operator

Definition at line 110 of file dgematrix-_dgematrix.hpp.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(matA.n!=matB.n || matA.m!=matB.m){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a subtraction." << std::endl
              << "Your input was (" << matA.m << "x" << matA.n << ") - (" << matB.m << "x" << matB.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG

  const CPPL_INT mn =matA.m*matA.n;
  for(CPPL_INT i=0; i<mn; i++){
    matB.array[i] =matA.array[i]-matB.array[i];
  }
  
  return matB;
}

_dgematrix operator- ( const dgematrix &  matA, const dsymatrix &  matB  )

friend

dgematrix-dsymatrix operator

Definition at line 99 of file dgematrix-dsymatrix.hpp.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(matA.n!=matB.n || matA.m!=matB.n){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a subtraction." << std::endl
              << "Your input was (" << matA.m << "x" << matA.n << ") - (" << matB.n << "x" << matB.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG

  dgematrix newmat(matA);
  for(CPPL_INT i=0; i<matA.m; i++){
    for(CPPL_INT j=0; j<matA.n; j++){
      newmat(i,j) -= matB(i,j);
    }
  }
  
  return _(newmat);
}

_dgematrix operator- ( const dgematrix &  matA, const _dsymatrix &  matB  )

friend

dgematrix-_dsymatrix operator

Definition at line 103 of file dgematrix-_dsymatrix.hpp.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(matA.n!=matB.n || matA.m!=matB.n){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a subtraction." << std::endl
              << "Your input was (" << matA.m << "x" << matA.n << ") - (" << matB.n << "x" << matB.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG

  dgematrix newmat(matA);
  for(CPPL_INT i=0; i<matA.m; i++){
    for(CPPL_INT j=0; j<matA.n; j++){
      newmat(i,j) -= matB(i,j);
    }
  }
  
  matB.destroy();
  return _(newmat);
}

_dgematrix operator- ( const dgematrix &  matA, const dgbmatrix &  matB  )

friend

dgematrix-dgbmatrix operator

Definition at line 106 of file dgematrix-dgbmatrix.hpp.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(matA.n!=matB.n || matA.m!=matB.m){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a summation." << std::endl
              << "Your input was (" << matA.m << "x" << matA.n << ") + (" << matB.m << "x" << matB.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  dgematrix newmat(matA);
  
  for(CPPL_INT i=0; i<matB.m; i++){
    const CPPL_INT jmax =std::min(matB.n,i+matB.ku+1);
    for(CPPL_INT j=std::max(CPPL_INT(0),i-matB.kl); j<jmax; j++){
      newmat(i,j) -= matB(i,j);
    }
  }
  
  return _(newmat);
}

_dgematrix operator- ( const dgematrix &  matA, const _dgbmatrix &  matB  )

friend

dgematrix-_dgbmatrix operator

Definition at line 110 of file dgematrix-_dgbmatrix.hpp.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(matA.n!=matB.n || matA.m!=matB.m){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a summation." << std::endl
              << "Your input was (" << matA.m << "x" << matA.n << ") + (" << matB.m << "x" << matB.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  dgematrix newmat(matA);
  
  for(CPPL_INT i=0; i<matB.m; i++){
    const CPPL_INT jmax =std::min(matB.n,i+matB.ku+1);
    for(CPPL_INT j=std::max(CPPL_INT(0),i-matB.kl); j<jmax; j++){
      newmat(i,j) -= matB(i,j);
    }
  }
  
  matB.destroy();
  return _(newmat);
}

_dgematrix operator- ( const dgematrix &  matA, const dgsmatrix &  matB  )

friend

dgematrix-dgsmatrix operator

Definition at line 26 of file dgematrix-dgsmatrix.hpp.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(matA.m!=matB.m || matA.n!=matB.n){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a subtraction." << std::endl
              << "Your input was (" << matA.m << "x" << matA.n << ") - (" << matB.m << "x" << matB.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  dgematrix newmat(matA);
  
  const std::vector<dcomponent>::const_iterator matB_data_end =matB.data.end();
  for(std::vector<dcomponent>::const_iterator it=matB.data.begin(); it!=matB_data_end; it++){
    newmat(it->i,it->j) -= it->v;
  }
  
  return _(newmat);
}

_dgematrix operator- ( const dgematrix &  matA, const _dgsmatrix &  matB  )

friend

dgematrix-_dgsmatrix operator

Definition at line 27 of file dgematrix-_dgsmatrix.hpp.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(matA.m!=matB.m || matA.n!=matB.n){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a subtraction." << std::endl
              << "Your input was (" << matA.m << "x" << matA.n << ") - (" << matB.m << "x" << matB.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG

  dgematrix newmat(matA);
  
  const std::vector<dcomponent>::const_iterator matB_data_end =matB.data.end();
  for(std::vector<dcomponent>::const_iterator it=matB.data.begin(); it!=matB_data_end; it++){
    newmat(it->i,it->j) -= it->v;
  }
  
  matB.destroy();
  return _(newmat);
}

_dgematrix operator- ( const dgematrix &  , const dssmatrix &    )

friend

_dgematrix operator- ( const dgematrix &  , const _dssmatrix &    )

friend

_dcovector operator* ( const dgematrix &  mat, const dcovector &  vec  )

friend

dgematrix*dcovector operator

Definition at line 3 of file dgematrix-dcovector.hpp.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(mat.n!=vec.l){
    ERROR_REPORT;
    std::cerr << "These matrix and vector can not make a product." << std::endl
              << "Your input was (" << mat.m << "x" << mat.n << ") * (" << vec.l << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  dcovector newvec(mat.m);
  char trans ='n';
  double alpha =1.;
  CPPL_INT inc =1;
  double beta =0.;
  
  dgemv_( &trans, &mat.m, &mat.n, &alpha, mat.array, &mat.m, vec.array, &inc, &beta, newvec.array, &inc );
  
  return _(newvec);
}

_dcovector operator* ( const dgematrix &  mat, const _dcovector &  vec  )

friend

dgematrix*_dcovector operator

Definition at line 3 of file dgematrix-_dcovector.hpp.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(mat.n!=vec.l){
    ERROR_REPORT;
    std::cerr << "These matrix and vector can not make a product." << std::endl
              << "Your input was (" << mat.m << "x" << mat.n << ") * (" << vec.l << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  dcovector newvec(mat.m);
  char trans ='n';
  double alpha =1.;
  CPPL_INT inc =1;
  double beta =0.;
  
  dgemv_( &trans, &mat.m, &mat.n, &alpha, mat.array, &mat.m, vec.array, &inc, &beta, newvec.array, &inc );
  
  vec.destroy();
  return _(newvec);
}

_dgematrix operator* ( const dgematrix &  matA, const dgematrix &  matB  )

friend

dgematrix*dgematrix operator

Definition at line 132 of file dgematrix-dgematrix.hpp.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(matA.n!=matB.m){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a product." << std::endl
              << "Your input was (" << matA.m << "x" << matA.n << ") * (" << matB.m << "x" << matB.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  dgematrix newmat( matA.m, matB.n );
  char transa ='n';
  char transb ='n';
  double alpha =1.;
  double beta =0.;
  
  dgemm_( &transa, &transb, &matA.m, &matB.n, &matA.n, &alpha, matA.array, &matA.m, matB.array, &matB.m, &beta, newmat.array, &matA.m );
  
  return _(newmat);
}

_dgematrix operator* ( const dgematrix &  matA, const _dgematrix &  matB  )

friend

dgematrix*_dgematrix operator

Definition at line 131 of file dgematrix-_dgematrix.hpp.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(matA.n!=matB.m){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a product." << std::endl
              << "Your input was (" << matA.m << "x" << matA.n << ") * (" << matB.m << "x" << matB.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  dgematrix newmat( matA.m, matB.n );
  char transa ='n';
  char transb ='n';
  double alpha =1.;
  double beta =0.;
  
  dgemm_( &transa, &transb, &matA.m, &matB.n, &matA.n, &alpha, matA.array, &matA.m, matB.array, &matB.m, &beta, newmat.array, &matA.m );
  
  matB.destroy();
  return _(newmat);
}

_dgematrix operator* ( const dgematrix &  matA, const dsymatrix &  matB  )

friend

dgematrix*dsymatrix operator

Definition at line 122 of file dgematrix-dsymatrix.hpp.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(matA.n!=matB.n){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a product." << std::endl
              << "Your input was (" << matA.m << "x" << matA.n << ") * (" << matB.n << "x" << matB.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  dgematrix newmat( matA.m, matB.n );
  char side ='R';
  char uplo ='l';
  double alpha =1.;
  double beta =0.;
  
  dsymm_( &side, &uplo, &newmat.m, &newmat.n, &alpha, matB.array, &matB.n, matA.array, &matA.m, &beta, newmat.array, &newmat.m );
  
  return _(newmat);
}

_dgematrix operator* ( const dgematrix &  matA, const _dsymatrix &  matB  )

friend

dgematrix*_dsymatrix operator

Definition at line 127 of file dgematrix-_dsymatrix.hpp.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(matA.n!=matB.n){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a product." << std::endl
              << "Your input was (" << matA.m << "x" << matA.n << ") * (" << matB.n << "x" << matB.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  dgematrix newmat( matA.m, matB.n );
  char side ='R';
  char uplo ='l';
  double alpha =1.;
  double beta =0.;
  
  dsymm_( &side, &uplo, &newmat.m, &newmat.n, &alpha, matB.array, &matB.n, matA.array, &matA.m, &beta, newmat.array, &newmat.m );
  
  matB.destroy();
  return _(newmat);
}

_dgematrix operator* ( const dgematrix &  matA, const dgbmatrix &  matB  )

friend

dgematrix*dgbmatrix operator

Definition at line 131 of file dgematrix-dgbmatrix.hpp.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(matA.n!=matB.m){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a product." << std::endl
              << "Your input was (" << matA.m << "x" << matA.n << ") * (" << matB.m << "x" << matB.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  dgematrix newmat( matA.m, matB.n );
  newmat.zero();
  
  for(CPPL_INT i=0; i<newmat.m; i++){
    for(CPPL_INT j=0; j<newmat.n; j++){
      const CPPL_INT kmax =std::min(matB.m,j+matB.kl+1);
      for(CPPL_INT k=std::max(CPPL_INT(0),j-matB.ku); k<kmax; k++){
        newmat(i,j) += matA(i,k)*matB(k,j);
      }
    }
  }
  
  return _(newmat);
}

_dgematrix operator* ( const dgematrix &  matA, const _dgbmatrix &  matB  )

friend

dgematrix*_dgbmatrix operator

Definition at line 136 of file dgematrix-_dgbmatrix.hpp.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(matA.n!=matB.m){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a product." << std::endl
              << "Your input was (" << matA.m << "x" << matA.n << ") * (" << matB.m << "x" << matB.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  dgematrix newmat( matA.m, matB.n );
  newmat.zero();
  
  for(CPPL_INT i=0; i<newmat.m; i++){
    for(CPPL_INT j=0; j<newmat.n; j++){
      const CPPL_INT kmax =std::min(matB.m,j+matB.kl+1);
      for(CPPL_INT k=std::max(CPPL_INT(0),j-matB.ku); k<kmax; k++){
        newmat(i,j) += matA(i,k)*matB(k,j);
      }
    }
  }
  
  matB.destroy();
  return _(newmat);
}

_dgematrix operator* ( const dgematrix &  matA, const dgsmatrix &  matB  )

friend

dgematrix*dgsmatrix operator

Definition at line 49 of file dgematrix-dgsmatrix.hpp.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(matA.n!=matB.m){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a product." << std::endl
              << "Your input was (" << matA.m << "x" << matA.n << ") * (" << matB.m << "x" << matB.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  dgematrix newmat(matA.m, matB.n);
  newmat.zero();
  
  const std::vector<dcomponent>::const_iterator matB_data_end =matB.data.end();
  for(std::vector<dcomponent>::const_iterator it=matB.data.begin(); it!=matB_data_end; it++){
    for(CPPL_INT i=0; i<matA.m; i++){
      newmat(i,it->j) += matA(i,it->i)*it->v;
    }
  }
  
  return _(newmat);
}

_dgematrix operator* ( const dgematrix &  matA, const _dgsmatrix &  matB  )

friend

dgematrix*_dgsmatrix operator

Definition at line 51 of file dgematrix-_dgsmatrix.hpp.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(matA.n!=matB.m){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a product." << std::endl
              << "Your input was (" << matA.m << "x" << matA.n << ") * (" << matB.m << "x" << matB.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  dgematrix newmat(matA.m, matB.n);
  newmat.zero();
  
  const std::vector<dcomponent>::const_iterator matB_data_end =matB.data.end();
  for(std::vector<dcomponent>::const_iterator it=matB.data.begin(); it!=matB_data_end; it++){
    for(CPPL_INT i=0; i<matA.m; i++){
      newmat(i,it->j) += matA(i,it->i)*it->v;
    }
  }
  
  matB.destroy();
  return _(newmat);
}

_dgematrix operator* ( const dgematrix &  , const dssmatrix &    )

friend

_dgematrix operator* ( const dgematrix &  , const _dssmatrix &    )

friend

_dgematrix operator* ( const dgematrix &  mat, const double &  d  )

friend

dgematrix*double operator

Definition at line 28 of file dgematrix-double.hpp.

{CPPL_VERBOSE_REPORT;
  dgematrix newmat(mat.m, mat.n);
  
  const CPPL_INT mn =mat.m*mat.n;
  for(CPPL_INT i=0; i<mn; i++){
    newmat.array[i] =mat.array[i]*d;
  }
  
  return _(newmat);
}

_dgematrix operator/ ( const dgematrix &  mat, const double &  d  )

friend

dgematrix/double operator

Definition at line 42 of file dgematrix-double.hpp.

{CPPL_VERBOSE_REPORT;
  dgematrix newmat(mat.m, mat.n);
  
  const CPPL_INT mn =mat.m*mat.n;
  for(CPPL_INT i=0; i<mn; i++){
    newmat.array[i] =mat.array[i]/d;
  }
  
  return _(newmat);
}

_drovector operator% ( const dgematrix &  matA, const dgematrix &  matB  )

friend

dgematrixdgematrix operator

Definition at line 156 of file dgematrix-dgematrix.hpp.

{CPPL_VERBOSE_REPORT;
#ifdef  CPPL_DEBUG
  if(matA.m!=matB.m || matA.n!=matB.n){
    ERROR_REPORT;
    std::cerr << "These two matrises can not make a product." << std::endl
              << "Your input was (" << matA.m << "x" << matA.n << ") % (" << matB.m << "x" << matB.n << ")." << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  drovector newvec( matA.n );
  
  newvec.zero();
  for(CPPL_INT j=0; j<matA.n; j++){
    for(CPPL_INT i=0; i<matA.m; i++){
      newvec(j) +=matA(i,j)*matB(i,j);
    }
  }
  
  return _(newvec);
}

_dgematrix operator* ( const double &  d, const dgematrix &  mat  )

friend

double*dgematrix operator

Definition at line 3 of file double-dgematrix.hpp.

{CPPL_VERBOSE_REPORT;
  dgematrix newmat(mat.m, mat.n);
  
  const CPPL_INT size =mat.m*mat.n;
  for(CPPL_INT i=0; i<size; i++){
    newmat.array[i] =d*mat.array[i];
  }
  
  return _(newmat);
}

_dgematrix hadamard ( const dgematrix &  , const dgematrix &    )

friend

Member Data Documentation

CPPL_INT dgematrix::m

matrix row size

Definition at line 9 of file dgematrix.hpp.

Referenced by _(), chsign(), clear(), col(), copy(), damax(), dgbmatrix::dgbsv(), dgeev(), dgels(), dgelss(), dgematrix(), dgematrix_small< m, n >::dgematrix_small(), dgesv(), dgesvd(), dggev(), dgglse(), dsymatrix::dsysv(), hadamerd(), i(), idamax(), identity(), operator%(), operator()(), operator*(), operator*=(), operator+(), operator+=(), operator-(), operator-=(), operator/(), operator/=(), operator<<(), read(), resize(), row(), set(), shallow_copy(), swap(), t(), to_zgematrix(), write(), and zero().

CPPL_INT dgematrix::n

matrix column size

Definition at line 10 of file dgematrix.hpp.

Referenced by _(), chsign(), clear(), col(), copy(), damax(), dgbmatrix::dgbsv(), dgeev(), dgels(), dgelss(), dgematrix(), dgematrix_small< m, n >::dgematrix_small(), dgesv(), dgesvd(), dggev(), dsymatrix::dsysv(), hadamerd(), i(), idamax(), identity(), operator%(), operator()(), operator*(), operator*=(), operator+(), operator+=(), operator-(), operator-=(), operator/(), operator/=(), operator<<(), read(), resize(), row(), set(), shallow_copy(), swap(), t(), to_zgematrix(), write(), and zero().

double* dgematrix::array

1D array to store matrix data

Definition at line 11 of file dgematrix.hpp.

Referenced by _(), chsign(), clear(), copy(), damax(), dgbmatrix::dgbsv(), dgeev(), dgels(), dgelss(), dgematrix(), dgematrix_small< m, n >::dgematrix_small(), dgesv(), dgesvd(), dggev(), dgglse(), dsymatrix::dsysv(), idamax(), identity(), operator*(), operator*=(), operator+(), operator+=(), operator-(), operator-=(), operator/(), operator/=(), operator=(), resize(), shallow_copy(), swap(), to_zgematrix(), zero(), and ~dgematrix().

double** dgematrix::darray

array of pointers of column head addresses

Definition at line 12 of file dgematrix.hpp.

Referenced by _(), clear(), copy(), dgematrix(), operator()(), resize(), set(), shallow_copy(), swap(), and ~dgematrix().

---

The documentation for this class was generated from the following files:

• dgematrix.hpp
• dgematrix-_dgbmatrix.hpp
• dgematrix-_dgematrix.hpp
• dgematrix-_dsymatrix.hpp
• dgematrix-cast.hpp
• dgematrix-constructor.hpp
• dgematrix-dgbmatrix.hpp
• dgematrix-dgematrix.hpp
• dgematrix-double.hpp
• dgematrix-dsymatrix.hpp
• dgematrix-io.hpp
• dgematrix-lapack.hpp
• dgematrix-misc.hpp