lapac.h

//
// Created by milinda on 1/19/17.
//

#ifndef SFCSORTBENCH_LAPAC_H
#define SFCSORTBENCH_LAPAC_H

#include <cstring>
#include <iostream>
//#include "lapacke.h"

extern "C" void dgesv_( int* n, int* nrhs, double* a, int* lda, int* ipiv,double* b, int* ldb, int* info );
extern "C" void dsyev_( char* jobz, char* uplo, int* n, double* a, int* lda,double* w, double* work, int* lwork, int* info);

// LU decomoposition of a general matrix
extern "C" void dgetrf_(int* M, int *N, double* A, int* lda, int* IPIV, int* INFO);

// generate inverse of a matrix given its LU decomposition
extern "C" void dgetri_(int* N, double* A, int* lda, int* IPIV, double* WORK, int* lwork, int* INFO);


namespace lapack
{

inline void lapack_DGESV(int n , int nrhs, const double * A, int lda, double * B, double * X, int ldb, int info)
    {
        int * ipiv = new int [n] ;
        //memcpy(X,B,sizeof(double)*n*nrhs);

        for(unsigned int i=0;i<nrhs;i++)
            for(unsigned int j=0;j<n;j++)
                X[i*n+j]=B[j*nrhs+i];

        double * L=new double[n*n];

        for(unsigned int i=0;i<n;i++)
            for(unsigned int j=0;j<n;j++)
                L[i*n+j] =A [j*n+i];

        dgesv_(&n,&nrhs,L,&lda,ipiv,X,&lda,&info);

        if( info > 0 ) {
            printf( "The diagonal element of the triangular factor of A,\n" );
            printf( "U(%i,%i) is zero, so that A is singular;\n", info, info );
            printf( "the solution could not be computed.\n" );

        }

        memcpy(L,X,sizeof(double)*n*nrhs);

        for(unsigned int i=0;i<nrhs;i++)
            for(unsigned int j=0;j<n;j++)
                X[i*nrhs+j] = L[j*n+i];


        /*lapack_int * ipiv = (lapack_int *)malloc(n*sizeof(lapack_int)) ;
        memcpy(X,B,sizeof(double)*n*nrhs);

        double * L=new double[n*n];
        memcpy(L,A,sizeof(double)*n*n);

        info = LAPACKE_dgesv( LAPACK_ROW_MAJOR, n, nrhs, L, lda, ipiv, X,ldb);*/


        if (info <0) {printf(" lapack linear solve failed. \n"); }
        delete [] L;
        delete [] ipiv;

        return ;
    }

inline void lapack_DSYEV(int n, const double * A, int lda, double * wr,double * vs,int info )
{


    double * L=new double[n*n];

    for(unsigned int i=0;i<n;i++)
        for(unsigned int j=0;j<n;j++)
            L[i*n+j] =A [j*n+i];

    double wkopt;
    double* work;
    int lwork=-1;
    dsyev_((char*) "Vectors", (char*)"Upper", (int*)&n, L, (int*)&lda, wr, &wkopt, &lwork,(int*)&info );
    lwork = (int)wkopt;
    work = new double[lwork];
    /* Solve eigenproblem */
    dsyev_( (char*)"Vectors", (char*)"Upper", (int*)&n, L, (int*)&lda, wr, work, &lwork, (int*)&info );


    for(unsigned int i=0;i<n;i++)
        for(unsigned int j=0;j<n;j++)
            vs[i*n+j]=L[j*n+i];

   /* std::cout<<"lapack: "<<std::endl;

    for(unsigned int i=0;i<n;i++)
        std::cout<<" i: "<<i<<" eig : "<<wr[i]<<std::endl;

    for(unsigned int i=0;i<n;i++)
    {
        for(unsigned int j=0;j<n;j++)
        {
            std::cout<<vs[i*(n)+j]<<" ";
        }

        std::cout<<std::endl;
    }

    memcpy(vs,A,sizeof(double)*n*n);
    info=LAPACKE_dsyev(LAPACK_ROW_MAJOR,'V','U',n,vs,lda,wr);
    std::cout<<"lapacke: "<<std::endl;
    for(unsigned int i=0;i<n;i++)
        std::cout<<" i: "<<i<<" eig : "<<wr[i]<<std::endl;

    for(unsigned int i=0;i<n;i++)
    {
        for(unsigned int j=0;j<n;j++)
        {
            std::cout<<vs[i*(n)+j]<<" ";
        }

        std::cout<<std::endl;
    }*/

    delete [] work;
    delete [] L;
    if(info!=0) std::cout<<"lapack eigen solve failed. "<<std::endl;
    return;

}

inline void inverse(double* A, int N)
{
    int *IPIV = new int[N+1];
    int LWORK = N*N;
    double *WORK = new double[LWORK];
    int INFO;

    dgetrf_(&N,&N,A,&N,IPIV,&INFO);
    dgetri_(&N,A,&N,IPIV,WORK,&LWORK,&INFO);

    delete IPIV;
    delete WORK;
}





}// end of namespace


#endif //SFCSORTBENCH_LAPAC_H