lapack_stdcplx.cpp

Go to the documentation of this file.

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

//#define NO_NS

#define LAPACK_INLINE inline

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


NAMESPACE_TBCI

//********************************************************************************
// Solution of linear eqution systems, partial pivoting
//********************************************************************************
TVector<double> lu_solve_expert(const F_BandMatrix<double>& A, const Vector<double>& B, int equi);
TVector<double> lu_solve(const F_BandMatrix<double>& A, const Vector<double>& B);
TVector<double> lu_solve(F_BandMatrix<double>& A, const Vector<double>& B, integer kl, integer ku);
TVector<double> lu_solve(const F_Matrix<double>& A, const Vector<double>& B, int overwriteA);
F_TMatrix<double> lu_solve(const F_Matrix<double>& A, const F_Matrix<double>& B, int overwriteA);

F_TMatrix<double> inv(const F_Matrix<double>& A, int overwriteA);


F_TMatrix<CPLX__ complex<double> > inv(const F_Matrix<CPLX__ complex<double> >& A, int overwriteA)
{
        integer info = 0;
        integer n=A.columns(); // Dimension des Problems
        BVector<integer> ipiv(n);
        F_Matrix<CPLX__ complex<double> > temp;
        if (!overwriteA)
        {
                temp.resize(n,n);
                temp=A;
        }
        // Rechte Seite ist die Einheitsmatrix!
        F_TMatrix<CPLX__ complex<double> > invA(0.0, n, n);
        for (int i=0; i<n; i++)
                invA.setval(CPLX__ complex<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__ complex<double> > lu_solve(const F_Matrix<CPLX__ complex<double> >& A, const F_Matrix<CPLX__ complex<double> >& B, int overwriteA)
{
        integer info = 0;
        integer n=A.columns(); // Dimension des Problems
        BVector<integer> ipiv(n);
        F_Matrix<CPLX__ complex<double> > temp;
        if (!overwriteA)
        {
                temp.resize(n,n);
                temp=A;
        }
        F_TMatrix<CPLX__ complex<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;
}


//********************************************************************************
// eigenvalues (and eigenvectors) \a A*EV=EW*EV of symmetric double Matrix \a A
//********************************************************************************
int eig(const F_Matrix<double>& A, Vector<double>& EW);

int eig(const F_Matrix<double>& A, Vector<double>& EW, F_Matrix<double>& EV);

//symmtric Band Matrix
int eig(const F_BandMatrix<double>& A, Vector<double>& EW, F_Matrix<double>& EV);

int eig(const F_BandMatrix<double>& A, Vector<double>& EW, F_Matrix<double>& EV,
        double ew_min, double ew_max);

//********************************************************************************
//********************************************************************************
int eig(const F_Matrix<CPLX__ complex<double> >& A, Vector<CPLX__ complex<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__ complex<double> >& A, Vector<CPLX__ complex<double> >& EW, F_Matrix<CPLX__ complex<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__ complex<double> >& A,
        const F_BandMatrix<CPLX__ complex<double> >& B,
        double EW_min, double EW_max, Vector<double>& Eigenwerte,
        F_Matrix<CPLX__ complex<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) = CPLX__ complex<double> (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__ complex<double> >& A, Vector<double>& EW, F_Matrix<CPLX__ complex<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__ complex<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