linear_algebra_double.h

/* Copyright (C) 2003-2005  Martin Förg
 
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * 
 * This library 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
 * Lesser General Public License for more details.
 * 
 * You should have received a copy of the GNU Lesser General Public
 * License along with this library; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
 *
 * Contact:
 *   martin.o.foerg@googlemail.com
 *
 */
 
#ifndef linear_algebrad_h
#define linear_algebrad_h
 
#include "square_matrix.h"
#include "vector.h"
#include "row_vector.h"
#include "diagonal_matrix.h"
#include "symmetric_matrix.h"
#include "band_matrix.h"
#include <complex>
 
//-------------------------------------
// Matrix operations
//-------------------------------------
 
namespace fmatvec {
 
  FMATVEC_EXPORT int svd(Matrix<General, Ref, Ref, double> &A,    Matrix<General, Ref, Ref, double> &S,  SquareMatrix<Ref, double> &U, SquareMatrix<Ref, double> &VT, int Rueckgabe);
 
  FMATVEC_EXPORT Vector<Ref, std::complex<double>> eigval(const SquareMatrix<Ref, double> &A);
 
  FMATVEC_EXPORT int eigvec(const SquareMatrix<Ref, double> &A, SquareMatrix<Ref, std::complex<double>> &V, Vector<Ref, std::complex<double>> &w);
 
  FMATVEC_EXPORT int eigvec(const Matrix<Symmetric, Ref, Ref, double> &A, const Matrix<Symmetric, Ref, Ref, double> &B, SquareMatrix<Ref, double> &eigenvectors, Vector<Ref, double> &eigenvalues);
 
  FMATVEC_EXPORT Vector<Ref, double> eigval(const Matrix<Symmetric, Ref, Ref, double> &A);
 
  FMATVEC_EXPORT Vector<Ref, double> eigvalSel(const Matrix<Symmetric, Ref, Ref, double> &A, int il, int iu, double abstol = 0);
 
  FMATVEC_EXPORT int Doolittle_LU_with_Pivoting_Solve(const double *A, double B[], const int pivot[], double x[], int n);
 
  FMATVEC_EXPORT Matrix<General, Ref, Ref, double> slvLU(const SquareMatrix<Ref, double> &A, const Matrix<General, Ref, Ref, double> &X);
  FMATVEC_EXPORT Matrix<General, Var, Var, double> slvLU(const SquareMatrix<Var, double> &A, const Matrix<General, Var, Var, double> &X, int & info);
 
  FMATVEC_EXPORT Vector<Ref, double> slvLU(const SquareMatrix<Ref, double> &A, const Vector<Ref, double> &x);
 
  FMATVEC_EXPORT Vector<Ref, double> slvLU(const SquareMatrix<Ref, double> &A, const Vector<Ref, double> &x, int & info);
  FMATVEC_EXPORT Vector<Var, double> slvLU(const SquareMatrix<Var, double> &A, const Vector<Var, double> &x, int & info);
 
  FMATVEC_EXPORT Matrix<General, Ref, Ref, double> slvLUFac(const SquareMatrix<Ref, double> &A, const Matrix<General, Ref, Ref, double> &X, const Vector<Ref, int> &ipiv);
  FMATVEC_EXPORT Matrix<General, Var, Var, double> slvLUFac(const SquareMatrix<Var, double> &A, const Matrix<General, Var, Var, double> &X, const Vector<Var, int> &ipiv);
 
  FMATVEC_EXPORT Vector<Ref, double> slvLUFac(const SquareMatrix<Ref, double> &A, const Vector<Ref, double> &x, const Vector<Ref, int> &ipiv);
  FMATVEC_EXPORT Vector<Var, double> slvLUFac(const SquareMatrix<Var, double> &A, const Vector<Var, double> &x, const Vector<Var, int> &ipiv);
 
  FMATVEC_EXPORT Matrix<General, Ref, Ref, double> slvLL(const Matrix<Symmetric, Ref, Ref, double> &A, const Matrix<General, Ref, Ref, double> &X);
 
  FMATVEC_EXPORT Vector<Ref, double> slvLL(const Matrix<Symmetric, Ref, Ref, double> &A, const Vector<Ref, double> &x);
 
  FMATVEC_EXPORT Matrix<General, Ref, Ref, double> slvQR(const SquareMatrix<Ref, double> &A, const Matrix<General, Ref, Ref, double> &X);
 
  FMATVEC_EXPORT Vector<Ref, double> slvQR(const SquareMatrix<Ref, double> &A, const Vector<Ref, double> &x);
 
  FMATVEC_EXPORT Vector<Ref, double> slvQRLQ(const Matrix<General, Ref, Ref, double> &A, const Vector<Ref, double> &b);
 
  FMATVEC_EXPORT Matrix<General, Ref, Ref, double> slvQRLQ(const Matrix<General, Ref, Ref, double> &A, const Matrix<General, Ref, Ref, double> &B);
 
  FMATVEC_EXPORT SquareMatrix<Ref, double> inv(const SquareMatrix<Ref, double> &A);
 
  FMATVEC_EXPORT Matrix<Symmetric, Ref, Ref, double> inv(const Matrix<Symmetric, Ref, Ref, double> &A);
 
  FMATVEC_EXPORT Matrix<Diagonal, Ref, Ref, double> inv(const Matrix<Diagonal, Ref, Ref, double> &A);
 
  FMATVEC_EXPORT Matrix<General, Ref, Ref, double> facLU(const Matrix<General, Ref, Ref, double> &A, Vector<Ref, int> &ipiv);
 
  FMATVEC_EXPORT SquareMatrix<Ref, double> facLU(const SquareMatrix<Ref, double> &A, Vector<Ref, int> &ipiv);
 
  FMATVEC_EXPORT SquareMatrix<Var, double> facLU(const SquareMatrix<Var, double> &A, Vector<Var, int> &ipiv);
 
  FMATVEC_EXPORT int facLU(double *A, int pivot[], int n);
  
  FMATVEC_EXPORT Matrix<Symmetric, Ref, Ref, double> facLL(const Matrix<Symmetric, Ref, Ref, double> &A);
 
  FMATVEC_EXPORT double nrm1(const Vector<Ref, double> &x);
 
  FMATVEC_EXPORT double nrm2(const Vector<Ref, double> &x);
 
  FMATVEC_EXPORT double nrmInf(const Vector<Ref, double> &x);
  
  FMATVEC_EXPORT double nrm1(const Matrix<General, Ref, Ref, double> &A);
 
  FMATVEC_EXPORT double nrm2(const Matrix<General, Ref, Ref, double> &A);
 
  FMATVEC_EXPORT double nrmInf(const Matrix<General, Ref, Ref, double> &A);
 
  FMATVEC_EXPORT double nrmInf(const Matrix<Symmetric, Ref, Ref, double> &A);
 
  FMATVEC_EXPORT double nrmFro(const Matrix<General, Ref, Ref, double> &A);
 
  FMATVEC_EXPORT double rho(const SquareMatrix<Ref, double> &A);
 
  FMATVEC_EXPORT double rho(const Matrix<Symmetric, Ref, Ref, double> &A);
 
  FMATVEC_EXPORT Vector<Ref, double> slvLLFac(const Matrix<Symmetric, Ref, Ref, double> &A, const Vector<Ref, double> &x);
 
  FMATVEC_EXPORT Matrix<General, Ref, Ref, double> slvLLFac(const Matrix<Symmetric, Ref, Ref, double> &A, const Matrix<General, Ref, Ref, double> &X);
 
  FMATVEC_EXPORT Matrix<General, Ref, Ref, double> slvLS(const Matrix<General, Ref, Ref, double> &A, const Matrix<General, Ref, Ref, double> &B, double rcond = -1);
  
  FMATVEC_EXPORT Vector<Ref, double> slvLS(const Matrix<General, Ref, Ref, double> &A, const Vector<Ref, double> &b, double rcond = -1);
 
//Matrix<General,Ref,Ref,double> slvLU(CBLAS_SIDE side, CBLAS_UPLO uplo, CBLAS_DIAG unit, const SquareMatrix<Ref,double> &A, const Matrix<General,Ref,Ref,double> &X, const Vector<Ref,int> &ipiv );
 
//Matrix<General,Ref,Ref,double> swap(const Matrix<General,Ref,Ref,double> &A, const Vector<Ref,int> &ipiv );
 
  template <int size>
  Vector<Fixed<size>, double> slvLU(const SquareMatrix<Fixed<size>, double> &A, const Vector<Fixed<size>, double> &b, int & info) {
 
    if (size == 0)
      return b;
 
    int ipiv[size];
 
    SquareMatrix<Fixed<size>, double> LU = A;
 
    info = facLU(LU(), ipiv, size);
 
    Vector<Fixed<size>, double> x;
    Vector<Fixed<size>, double> y = b;
 
    info = Doolittle_LU_with_Pivoting_Solve(LU(), y(), ipiv, x(), size);
 
    return x;
 
  }
 
  template <int size>
  SquareMatrix<Fixed<size>, double> facLU(const SquareMatrix<Fixed<size>, double> &A, Vector<Fixed<size>, int> &ipiv) {
    SquareMatrix<Fixed<size>, double> LU = A;
    facLU(LU(), ipiv(), size);
    return LU;
  }
 
  template <int size>
  Vector<Fixed<size>, double> slvLUFac(const SquareMatrix<Fixed<size>, double> &ALU, const Vector<Fixed<size>, double> &b, const Vector<Fixed<size>, int> &ipiv, int &info) {
    Vector<Fixed<size>, double> x;
    Vector<Fixed<size>, double> y = b;
    info = Doolittle_LU_with_Pivoting_Solve(ALU(), y(), ipiv(), x(), size);
    return x;
  }
 
  template <int size, int nrrhs>
  Matrix<General, Fixed<size>, Fixed<nrrhs>, double> slvLUFac(const SquareMatrix<Fixed<size>, double> &ALU, const Matrix<General, Fixed<size>, Fixed<nrrhs>, double> &B, const Vector<Fixed<size>, int> &ipiv, int &info) {
    Matrix<General, Fixed<size>, Fixed<nrrhs>, double> X;
    Vector<Fixed<size>, double> x;
    info = 0;
    for(int c=0; c<nrrhs; c++) {
      Vector<Fixed<size>, double> y = B.col(c);
      int localInfo = Doolittle_LU_with_Pivoting_Solve(ALU(), y(), ipiv(), x(), size);
      if(localInfo!=0)
        info = localInfo;
      X.set(c, x);
    }
    return X;
  }
 
  template <int size>
  Matrix<General, Fixed<size>, Var, double> slvLUFac(const SquareMatrix<Fixed<size>, double> &ALU, const Matrix<General, Fixed<size>, Var, double> &B, const Vector<Fixed<size>, int> &ipiv, int &info) {
    Matrix<General, Fixed<size>, Var, double> X(size,B.cols());
    Vector<Fixed<size>, double> x;
    info = 0;
    for(int c=0; c<B.cols(); c++) {
      Vector<Fixed<size>, double> y = B.col(c);
      int localInfo = Doolittle_LU_with_Pivoting_Solve(ALU(), y(), ipiv(), x(), size);
      if(localInfo!=0)
        info = localInfo;
      X.set(c, x);
    }
    return X;
  }
 
}
 
#endif