lapack.cpp

Go to the documentation of this file.

/* $Id: lapack.cpp,v 1.8.2.12 2019/05/28 11:13:02 garloff Exp $ */

//#define NO_NS

#define LAPACK_INLINE inline

#include "tbci/cplx.h"
#include "tbci/f_matrix.h"
#include "tbci/f_bandmatrix.h"
#include "tbci/lapack/lapack_lib.h"


NAMESPACE_TBCI


//********************************************************************************
//********************************************************************************
TVector<double> lu_solve_expert(const F_BandMatrix<double>& A, const Vector<double>& B, int equi)
{
  char fact='E';        // equlibrate
  if(!equi) fact='N';   // No equilibration
  char trans='N';       //No transpose
  integer n=A.columns();        // Dimension des Problems
  integer kl=A.numSub();        //echt genutzte ND
  integer ku=A.numSuper();      //echt genutzte ND
  integer nrhs=1;               // Anzahl rechter Seiten
  //ab=A
  integer ldab=A.ldab();
  double* afb=new double[n*(2*kl+ku+1)];
  integer ldafb=2*kl+ku+1;
  integer* ipiv=new integer[n];
  char equed;
  double* r=new double[n];
  double* c=new double[n];
  TVector<double> x(B.size());
  double rcond;
  double ferr;          //Vektor fuer mehrere rechte Seiten
  double berr;          //Vektor fuer mehrere rechte Seiten
  double* work=new double[3*n];
  integer* iwork=new integer[n];
  integer info = 0;

  dgbsvx_(&fact, &trans, &n, &kl, &ku, &nrhs, A.get_fortran_matrix(), &ldab,
          afb, &ldafb, ipiv, &equed,
          r, c, B.get_fortran_vector(), &n,
          x.get_fortran_vector(), &n,
          &rcond, &ferr, &berr, work, iwork, &info);

  delete[] afb;
  delete[] ipiv;
  delete[] r;
  delete[] c;
  delete[] work;
  delete[] iwork;

  if (info)
  {
    STD__ cerr << "Error inverting Matrix A! (" << info << ")" << STD__ endl;
    STD__ cerr << "\n Equalized= " << equed << STD__ endl;
    STD__ cerr << " 1/Cond= " << rcond << " epsilon= " << dlamch_ ("Epsilon")
        << " f_err= " << ferr
        << " b_err= " << berr << STD__ endl;
  }
  return x;
}


TVector<double> lu_solve(const F_BandMatrix<double>& A, const Vector<double>& B)
{
  integer n=A.columns();  // Dimension des Problems
  integer kl=A.numSub();  //echt genutzte ND
  integer ku=A.numSuper();//echt genutzte ND
  integer nrhs = 1;         // Anzahl rechter Seiten
  integer info = 0;
  TVector<double> invA(B);
  BVector<integer> ipiv(n);

  //Mit Kopie der Matrix, inclusive Platz fuer fill_in!
  {
    F_BandMatrix<double> temp(n,2*ku,kl);
    temp.clear();
    //Matrix kopieren
    for(int i=0; i<n; i++)
    {
      temp(i,i)=A(i,i);
      for(int j=1; j<=kl; j++) {if(i-j>=0) temp(i,i-j)=A(i,i-j);}
      for(int k=1; k<=ku; k++) {if(i+k<n)  temp(i,i+k)=A(i,i+k);}
    }
    integer ldab=temp.ldab();  // mindestens 2*kl+ku+1;
    dgbsv_(&n, &kl, &ku,
           &nrhs,  temp.get_fortran_matrix(), &ldab, ipiv.get_fortran_vector(),
           invA.get_fortran_vector(), &n, &info);
  }
  if (info)
        STD__ cerr << "Error inverting Matrix A! (" << info << ")" << STD__ endl;
  return invA;
}

TVector<double> lu_solve(F_BandMatrix<double>& A, const Vector<double>& B, integer kl, integer ku)
{
  integer n=A.columns();        // Dimension des Problems
  //integer kl=A.numSub();      //echt genutzte ND
  //integer ku=A.numSuper();    //echt genutzte ND
  integer nrhs = 1;             // Anzahl rechter Seiten
  integer info = 0;
  TVector<double> invA(B);
  BVector<integer> ipiv(n);
  // Mit Ueberschreiben
  {
    if (ku > (long)A.numSuper()/2)
        STD__ cerr << "Matrix too small to store fill_in" << STD__ endl;
    integer ldab=A.ldab();  // mindestens 2*kl+ku+1;
    dgbsv_(&n, &kl, &ku,
           &nrhs,  A.get_fortran_matrix(), &ldab, ipiv.get_fortran_vector(),
           invA.get_fortran_vector(), &n, &info);
  }

  if (info)
        STD__ cerr << "Error inverting Matrix A! (" << info << ")" << STD__ endl;
  return invA;
}

TVector<double> lu_solve(const F_Matrix<double>& A, const Vector<double>& B, int overwriteA)
{
        integer n=A.columns(); // Dimension des Problems
        integer nrhs=1;
        integer info=0;
        TVector<double> invA(B);
        BVector<integer> ipiv(n);
        //not tested!!
        
        if (overwriteA)
        {
                dgesv_(&n, &nrhs, A.get_fortran_matrix(), &n, ipiv.get_fortran_vector(),
                       invA.get_fortran_vector(), &n, &info);
        }
        else
        {
                F_Matrix<double> temp(A);
                dgesv_(&n, &nrhs,  temp.get_fortran_matrix(), &n, ipiv.get_fortran_vector(),
                       invA.get_fortran_vector(), &n, &info);
        }

        if (info)
                STD__ cerr << "Error inverting Matrix A!" << STD__ endl;
        return invA;
}


F_TMatrix<double> lu_solve(const F_Matrix<double>& A, const F_Matrix<double>& B, int overwriteA)
{
        integer n=A.columns(); // Dimension des Problems
        integer info;
        F_TMatrix<double> invA(B);
        BVector<integer> ipiv(n);

        if (overwriteA)
        {
                dgesv_(&n, &n, A.get_fortran_matrix(), &n, ipiv.get_fortran_vector(),
                       invA.get_fortran_matrix(), &n, &info);
        }
        else
        {
                F_Matrix<double> temp(A);
                dgesv_(&n, &n,  temp.get_fortran_matrix(), &n, ipiv.get_fortran_vector(),
                       invA.get_fortran_matrix(), &n, &info);
        }
        
        if (info)
                STD__ cerr << "Error inverting Matrix A!" << STD__ endl;
        return invA;
}


F_TMatrix<double> inv(const F_Matrix<double>& A, int overwriteA)
{
        integer n=A.columns(); // Dimension des Problems
        integer info=0;
        // Rechte Seite ist die Einheitsmatrix!
        F_TMatrix<double> invA(0.0, n, n);
        for (int i=0; i<n; i++)
                invA.setval(1.0,i,i);
        BVector<integer> ipiv(n);
        F_Matrix<double> temp;
        if (!overwriteA)
        {
                temp.resize(n,n);
                temp=A;
        }

        if (overwriteA)
                dgesv_(&n, &n, A.get_fortran_matrix(), &n, ipiv.get_fortran_vector(),
                       invA.get_fortran_matrix(), &n, &info);
        else
                dgesv_(&n, &n,  temp.get_fortran_matrix(), &n, ipiv.get_fortran_vector(),
                       invA.get_fortran_matrix(), &n, &info);
        if (info)
                STD__ cerr << "Error inverting Matrix A!" << STD__ endl;
        return invA;
}

F_TMatrix<cplx<double> > inv(const F_Matrix<cplx<double> >& A, int overwriteA)
{
        integer n=A.columns(); // Dimension des Problems
        integer info=0;
        BVector<integer> ipiv(n);
        F_Matrix<cplx<double> > temp;
        if (!overwriteA)
        {
                temp.resize(n,n);
                temp=A;
        }
        // Rechte Seite ist die Einheitsmatrix!
        F_TMatrix<cplx<double> > invA(0.0, n, n);
        for (int i=0; i<n; i++)
                invA.setval(cplx<double>(1.0), i,i);

        if (overwriteA)
                zgesv_(&n, &n, (doublecomplex*) A.get_fortran_matrix(), &n, ipiv.get_fortran_vector(),
                (doublecomplex*) invA.get_fortran_matrix(), &n, &info);
        else
                zgesv_(&n, &n, (doublecomplex*) temp.get_fortran_matrix(), &n, ipiv.get_fortran_vector(),
                (doublecomplex*) invA.get_fortran_matrix(), &n, &info);
        if (info)
                STD__ cerr << "Error inverting Matrix A!" << STD__ endl;
        return invA;
}

F_TMatrix<cplx<double> > lu_solve(const F_Matrix<cplx<double> >& A, const F_Matrix<cplx<double> >& B, int overwriteA)
{
        integer n=A.columns(); // Dimension des Problems
        integer info=0;
        BVector<integer> ipiv(n);
        F_Matrix<cplx<double> > temp;
        if (!overwriteA)
        {
                temp.resize(n,n);
                temp=A;
        }
        F_TMatrix<cplx<double> > invA(B);

        if (overwriteA)
                zgesv_(&n, &n, (doublecomplex*) A.get_fortran_matrix(), &n, ipiv.get_fortran_vector(),
                (doublecomplex*) invA.get_fortran_matrix(), &n, &info);
        else
                zgesv_(&n, &n, (doublecomplex*) temp.get_fortran_matrix(), &n, ipiv.get_fortran_vector(),
                (doublecomplex*) invA.get_fortran_matrix(), &n, &info);
        if (info)
                STD__ cerr << "Error inverting Matrix A!" << STD__ endl;
        return invA;
}


//********************************************************************************
//********************************************************************************
int eig(const F_Matrix<double>& A, Vector<double>& EW)
{
        char jobz='N';         //Eigenwerte und Eigenvektoren
        char uplo='U';         //Obere Haelfte der Matrix lesen
        integer N=A.columns(); // Dimension des Problems
        integer lwork=3*N;
        BVector<double> work(lwork);
        integer info=0;
        F_Matrix<double> EV(A);  // Die Systemmatrix wird sonst ueberschrieben!!
        dsyev_(&jobz, &uplo, &N, EV.get_fortran_matrix(),
               &N, EW.get_fortran_vector(), work.get_fortran_vector(), &lwork,
               &info);
        return (int) info;
}


int eig(const F_Matrix<double>& A, Vector<double>& EW, F_Matrix<double>& EV)
{
        char jobz='V';         //Eigenwerte und Eigenvektoren
        char uplo='U';         //Obere Haelfte der Matrix lesen
        integer N=A.columns(); // Dimension des Problems
        integer lwork=3*N;
        BVector<double> work(lwork);
        integer info=0;
        EV = A;                 // Die Systemmatrix wird sonst ueberschrieben!!
        dsyev_(&jobz, &uplo, &N, EV.get_fortran_matrix(),
               &N, EW.get_fortran_vector(), work.get_fortran_vector(), &lwork,
               &info);
        return (int) info;
}


//symmtric Band Matrix
int eig(const F_BandMatrix<double>& A, Vector<double>& EW, F_Matrix<double>& EV)
{
  char jobz='V';                // Eigenwerte und Eigenvektoren
  char uplo='U';                // Obere Haelfte der Matrix lesen
  integer N=A.columns();        // Dimension des Problems
  integer kd=A.numSuper();      // Anzahl der oberen Diags
  integer lwork=3*N-2;
  integer lda=1+A.numSuper()+A.numSub();
  integer ldz  =N;              //fuer Eigenvektoren
  BVector<double> work(lwork);
  integer info=0;

  /* int dsbev_(char *jobz, char *uplo, integer *n, integer *kd,
        doublereal *ab, integer *ldab, doublereal *w, doublereal *z, integer *
        ldz, doublereal *work, integer *info); */

  dsbev_(&jobz, &uplo, & N, &kd, A.get_fortran_matrix(), &lda,
         EW.get_fortran_vector(), EV.get_fortran_matrix(),
         &ldz, work.get_fortran_vector(), &info);

  return (int) info;
}


int eig(const F_BandMatrix<double>& A, Vector<double>& EW, F_Matrix<double>& EV,
        double ew_min, double ew_max)
{
  char jobz='V';                // Eigenwerte und Eigenvektoren
  char uplo='U';                // Obere Haelfte der Matrix lesen
  char range='V';               // Eigenwerte aus Intervall
  integer N=A.columns();        // Dimension des Problems
  integer kd=A.numSuper();      // Anzahl der oberen Diags
  integer lda=1+A.numSuper()+A.numSub();
  //integer lwork=7*N;
  integer ldz  =N;              //fuer Eigenvektoren
  BVector<double> work(7*N);
  BVector<integer> iwork(5*N);
  BVector<integer> ifail(N);
  integer ldq  =N;              //fuer Eigenvektoren
  BVector<double> q(N*ldq);
  integer info=0, dummy, m;
  double abstol=0;

  Vector<double>   EW_temp(N);
  F_Matrix<double> EV_temp(N,N);

/* int dsbevx_(char *jobz, char *range, char *uplo, integer *n,
        integer *kd, doublereal *ab, integer *ldab,
  doublereal *q, integer *ldq,
  doublereal *vl, doublereal *vu, integer *il, integer *iu,
        doublereal *abstol, integer *m, doublereal *w, doublereal *z, integer *ldz,
  doublereal *work, integer *iwork, integer *ifail, integer *info);
*/
  dsbevx_(&jobz, &range, &uplo, &N,
          &kd, A.get_fortran_matrix(), &lda,
          q.get_fortran_vector(), &ldq,
          &ew_min, &ew_max, &dummy, &dummy,
          &abstol, &m, EW_temp.get_fortran_vector(), EV_temp.get_fortran_matrix(), &ldz,
          work.get_fortran_vector(), iwork.get_fortran_vector(),
          ifail.get_fortran_vector(), &info);

  // Eigenwerte kopieren
  EV.resize(N,m);
  EW.resize(m);
  for(int i=0; i<m; i++)
  {
    EW(i)=EW_temp(i);
    EV.set_col(EV_temp.get_col(i), i);
  }

  return (int) info;
}



//********************************************************************************
//********************************************************************************
int eig(const F_Matrix<cplx<double> >& A, Vector<cplx<double> >& EW)
{
        char jobvl='N';         //linke Eigenwerte und Eigenvektoren
        char jobvr='N';         //rechte Eigenwerte und Eigenvektoren
        integer N=A.columns();  // Dimension des Problems
        doublecomplex* vl=0;      // Keine linken Eigenvektoren
        doublecomplex* vr=0;      // Keine linken Eigenvektoren
        integer lwork=2*N;
        doublecomplex* work =NEW(doublecomplex,lwork);
        doublereal*    rwork=NEW(doublereal,lwork);
        integer info=0;
        zgeev_(&jobvl, &jobvr, &N,
               (doublecomplex*) A.get_fortran_matrix(), &N,
               (doublecomplex*) EW.get_fortran_vector(), vl,
               &N, vr, &N, work, &lwork, rwork, &info);
        TBCIDELETE(doublecomplex, work, lwork);
        TBCIDELETE(doublereal, rwork, lwork);
        return (int) info;
}


int eig(const F_Matrix<cplx<double> >& A, Vector<cplx<double> >& EW, F_Matrix<cplx<double> >& EV)
{
        char jobvl='N';         //linke Eigenwerte und Eigenvektoren
        char jobvr='V';         //rechte Eigenwerte und Eigenvektoren
        integer N=A.columns();  // Dimension des Problems
        doublecomplex* vl=0;      // Keine linken Eigenvektoren
        integer lwork=2*N;
        doublecomplex* work =NEW(doublecomplex,lwork);
        doublereal*    rwork=NEW(doublereal,lwork);
        integer info=0;
        zgeev_(&jobvl, &jobvr, &N, (doublecomplex*) A.get_fortran_matrix(), &N,
               (doublecomplex*) EW.get_fortran_vector(), vl,
               &N, (doublecomplex*) EV.get_fortran_matrix(), &N, work, &lwork, rwork, &info);
        TBCIDELETE(doublecomplex, work, lwork);
        TBCIDELETE(doublereal, rwork, lwork);
        return (int) info;
}




//********************************************************************************
//********************************************************************************
int eig(const F_BandMatrix<cplx<double> >& A, const F_BandMatrix<cplx<double> >& B,
                                double EW_min, double EW_max,
                                Vector<double>& Eigenwerte, F_Matrix<cplx<double> >& Eigenvectoren)
{
        // Benutzt wird nur die obere Haelfte der Bandmatrix A!!!!
        long uplo=1;
        // Ordnung des EWP bestimmen
        long n=A.columns();
        // Anzahl der Oberen Diagonalen der Matrizen A und B bestimmen
        long ldab=A.numSuper()+1;
        long ldbb=B.numSuper()+1;
        // Allg. Variablen
        long Anzahl_Eigenwerte;
        long BerechneEV=1; // Eigenvektoren sollen berechnet werden

        // Laufvariablen
        int i,j, x,y;


        // *****************************************************************************
        // Arrays im Fortran-Stil deklarieren,
        // diese Arrays werden mit den Werten aus den Matrizen besetzt
        //cout << "ab.bb..." << flush;
        // *****************************************************************************
        doublecomplex *ab = new doublecomplex[ldab*n];
        doublecomplex *bb = new doublecomplex[ldbb*n];

        // Arrays kopieren
        for (i=0; i<n; i++)
        {
                // Array A nach ab kopieren
                for (j=i; j<i+ldab && j<n; j++)
                {
                        // x und y sind die Koordinaten in Lapack-Bandstruktur
                        x=ldab+i-j-1;
                        y=j;
                        ab[x+ldab*y].r = CPLX__ real(A(i,j));
                        ab[x+ldab*y].i = CPLX__ imag(A(i,j));
                }

                // Array B nach bb kopieren
                for (j=i; j<i+ldbb && j<n; j++)
                {
                        // x und y sind die Koordinaten in Lapack-Bandstruktur
                        x=ldbb+i-j-1;
                        y=j;
                        bb[x+ldbb*y].r = CPLX__ real(B(i,j));
                        bb[x+ldbb*y].i = CPLX__ imag(B(i,j));
                }
        }

        // *******************************************************************************
        // Hilfsvariablen fuer die Eigenwertroutine deklarieren

        double* ew = new double[n];  // Temporaerer Eigenwertspeicher
        long *iwork = new long[5*n];
        double *rwork = new double[7*n];
        doublecomplex *work = new doublecomplex[n];
        doublecomplex *q = new doublecomplex[n*n];

        // *******************************************************************************
        // Eigenwertroutine aufrufen:
        long info=own_ew_(&BerechneEV, &uplo,  &n, &Anzahl_Eigenwerte, ab, &ldab, bb, &ldbb,
                          &EW_min, &EW_max, ew, q, &n, work, rwork, iwork);
        if (!info==0)
        {
                STD__ cerr << " Error calculating eigenvalues: (status " << info << ")" << STD__ endl;
                return info;
        }
        delete [] ab ;
        delete [] bb ;
        // Wenn Eigenvektoren berechnen
        if (BerechneEV==1 && Anzahl_Eigenwerte>0)
        {
          // neue Hilfsvariablen initialisieren
          doublecomplex *EV = new doublecomplex[n*Anzahl_Eigenwerte];
          long *ifail = new long[Anzahl_Eigenwerte];
          // Eigenvektoren zu 0 setzen
          for (i=0;i<n*Anzahl_Eigenwerte;i++)
          {
            EV[i].r=0;
            EV[i].i=0;
          }
          info=own_ev_(&n, &Anzahl_Eigenwerte, ew, EV, &n, q, &n, work, rwork, iwork, ifail);
          if (!info==0)
          {
                   STD__ cerr << " Error calculating eigenvectors" << STD__ endl;
                   return info;
          }
          // *******************************************************************************
          // Eigenvektoren in die Matrix Eigenvektoren eintragen
          Eigenvectoren.resize(n,Anzahl_Eigenwerte);
          for (j=0;j<Anzahl_Eigenwerte;j++)
          {
                for (int i=0; i<n; i++)
                {
                Eigenvectoren(i,j).set(EV[i+n*j].r, EV[i+n*j].i);
                  }
          }
          delete[] ifail;
          delete[] EV;
        } //endif Eigenvektorberechnung

        // Eigenwerte in Eigenwertvektor eintragen
        Eigenwerte.resize(Anzahl_Eigenwerte);
        for(i=0; i<Anzahl_Eigenwerte; i++)      Eigenwerte(i) = ew[i];

        delete[] ew;
        delete[] q;
        delete[] iwork;
        delete[] rwork;
        delete[] work;
        return 0;
}



// ***************************************************************************
// ***************************************************************************
int eig(const F_Matrix<cplx<double> >& A, Vector<double>& EW, F_Matrix<cplx<double> >& EVec)
{
        char jobz='V';                  // Eigenvektoren berechnen
        char uplo='U';                  // obere Haelfte der Matrix wird gespeichert
        integer n=A.rows();             // Ordnung des EWP bestimmen

        integer lwork =2*n-1;
        doublecomplex *work = new doublecomplex[lwork];
        double *rwork       = new double[3*n-2];
        integer info=0;

        EVec=A;
        zheev_(&jobz, &uplo, &n, (doublecomplex*) EVec.get_fortran_matrix(), &n,
               EW.get_fortran_vector(), work, &lwork, rwork, &info);

        if (!info==0)
        {
                STD__ cerr << " Error calculating eigenvalues: (status " << info << ")" << STD__ endl;
                return info;
        }

        // Alle alten Arrays loeschen
        delete[] work;
        delete[] rwork;
        return 0;
}


int eig(F_Matrix<cplx<double> > A, Vector<double>& EW)
{
        char jobz='N';                  // Eigenvektoren berechnen
        char uplo='U';                  // obere Haelfte der Matrix wird gespeichert
        integer n=A.rows();             // Ordnung des EWP bestimmen

        integer lwork =2*n-1;
        doublecomplex *work = new doublecomplex[lwork];
        double *rwork       = new double[3*n-2];
        integer info=0;

        zheev_(&jobz, &uplo, &n, (doublecomplex*) A.get_fortran_matrix(), &n,
               EW.get_fortran_vector(), work, &lwork, rwork, &info);

        if (!info==0)
        {
                STD__ cerr << " Error calculating eigenvalues: (status " << info << ")" << STD__ endl;
                return info;
        }

        // Alle alten Arrays loeschen
        delete[] work;
        delete[] rwork;
        return 0;
}

template class tbci_memalloc<doublecomplex>;
template class tbci_memalloc_cache<doublecomplex>;

NAMESPACE_END