LapackDriver.cpp

Go to the documentation of this file.

//Copyright (c) 2004-2005, Baris Sumengen
//All rights reserved.
//
// CIMPL Matrix Performance Library
//
//Redistribution and use in source and binary
//forms, with or without modification, are
//permitted provided that the following
//conditions are met:
//
//    * No commercial use is allowed. 
//    This software can only be used
//    for non-commercial purposes. This 
//    distribution is mainly intended for
//    academic research and teaching.
//    * Redistributions of source code must
//    retain the above copyright notice, this
//    list of conditions and the following
//    disclaimer.
//    * Redistributions of binary form must
//    mention the above copyright notice, this
//    list of conditions and the following
//    disclaimer in a clearly visible part 
//    in associated product manual, 
//    readme, and web site of the redistributed 
//    software.
//    * Redistributions in binary form must
//    reproduce the above copyright notice,
//    this list of conditions and the
//    following disclaimer in the
//    documentation and/or other materials
//    provided with the distribution.
//    * The name of Baris Sumengen may not be
//    used to endorse or promote products
//    derived from this software without
//    specific prior written permission.
//
//THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT
//HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
//EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT
//NOT LIMITED TO, THE IMPLIED WARRANTIES OF
//MERCHANTABILITY AND FITNESS FOR A PARTICULAR
//PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
//CONTRIBUTORS BE LIABLE FOR ANY
//DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
//EXEMPLARY, OR CONSEQUENTIAL DAMAGES
//(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
//OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
//DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
//HOWEVER CAUSED AND ON ANY THEORY OF
//LIABILITY, WHETHER IN CONTRACT, STRICT
//LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
//OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
//OF THIS SOFTWARE, EVEN IF ADVISED OF THE
//POSSIBILITY OF SUCH DAMAGE.



#include "./Lapack.h"
#include "mkl.h"



//=========================================================
// Linear equations solver
//=========================================================

// Gesv
Matrix<float> Lapack::Gesv(Matrix<float>& A, Matrix<float>& B)
{
        if(A.Columns() != A.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        if(A.Columns() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }
        int info;
        Vector<int> ipiv(B.Rows());
        Matrix<float> LU = A.Clone();
        Matrix<float> X = B.Clone();
        int xrows = X.Rows();
        int xcols = X.Columns();
        int lurows = LU.Rows();
        SGESV(&xrows, &xcols, LU.Data(), &lurows, ipiv.Data(), X.Data(), &xrows, &info);
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of SGESV is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "U(%i,%i) is exactly zero. The LU factorization has been completed, but the factor U is exactly singular, so the solution could not be computed!", info, info);
                        Utility::RunTimeError(s);
                }
        }
        return X;
}

Matrix<double> Lapack::Gesv(Matrix<double>& A, Matrix<double>& B)
{
        if(A.Columns() != A.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        if(A.Columns() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }
        int info;
        Vector<int> ipiv(B.Rows());
        Matrix<double> LU = A.Clone();
        Matrix<double> X = B.Clone();
        int xrows = X.Rows();
        int xcols = X.Columns();
        int lurows = LU.Rows();
        DGESV(&xrows, &xcols, LU.Data(), &lurows, ipiv.Data(), X.Data(), &xrows, &info);
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of DGESV is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "U(%i,%i) is exactly zero. The LU factorization has been completed, but the factor U is exactly singular, so the solution could not be computed!", info, info);
                        Utility::RunTimeError(s);
                }
        }
        return X;
}

Matrix<ComplexFloat> Lapack::Gesv(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B)
{
        if(A.Columns() != A.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        if(A.Columns() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }
        int info;
        Vector<int> ipiv(B.Rows());
        Matrix<ComplexFloat> LU = A.Clone();
        Matrix<ComplexFloat> X = B.Clone();
        int xrows = X.Rows();
        int xcols = X.Columns();
        int lurows = LU.Rows();
        CGESV(&xrows, &xcols, reinterpret_cast<MKL_Complex8*>(LU.Data()), &lurows, ipiv.Data(), reinterpret_cast<MKL_Complex8*>(X.Data()), &xrows, &info);
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of CGESV is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "U(%i,%i) is exactly zero. The LU factorization has been completed, but the factor U is exactly singular, so the solution could not be computed!", info, info);
                        Utility::RunTimeError(s);
                }
        }
        return X;
}

Matrix<ComplexDouble> Lapack::Gesv(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B)
{
        if(A.Columns() != A.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        if(A.Columns() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }
        int info;
        Vector<int> ipiv(B.Rows());
        Matrix<ComplexDouble> LU = A.Clone();
        Matrix<ComplexDouble> X = B.Clone();
        int xrows = X.Rows();
        int xcols = X.Columns();
        int lurows = LU.Rows();
        ZGESV(&xrows, &xcols, reinterpret_cast<MKL_Complex16*>(LU.Data()), &lurows, ipiv.Data(), reinterpret_cast<MKL_Complex16*>(X.Data()), &xrows, &info);
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of ZGESV is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "U(%i,%i) is exactly zero. The LU factorization has been completed, but the factor U is exactly singular, so the solution could not be computed!", info, info);
                        Utility::RunTimeError(s);
                }
        }
        return X;
}


//Gesvx


//Posv
Matrix<float> Lapack::Posv(Matrix<float>& A, Matrix<float>& B)
{
        if(A.Columns() != A.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        if(A.Columns() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }
        int info;
        Matrix<float> LU = A.Clone();
        Matrix<float> X = B.Clone();
        int xrows = X.Rows();
        int xcols = X.Columns();
        int lurows = LU.Rows();
        SPOSV("U", &xrows, &xcols, LU.Data(), &lurows, X.Data(), &xrows, &info);
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of SPOSV is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "The leading minor of order %i (and hence the matrix A itself) is not positive definite, so the factorization could not be completed, and the solution has not been computed!", info);
                        Utility::RunTimeError(s);
                }
        }
        return X;
}

Matrix<double> Lapack::Posv(Matrix<double>& A, Matrix<double>& B)
{
        if(A.Columns() != A.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        if(A.Columns() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }
        int info;
        Matrix<double> LU = A.Clone();
        Matrix<double> X = B.Clone();
        int xrows = X.Rows();
        int xcols = X.Columns();
        int lurows = LU.Rows();
        DPOSV("U", &xrows, &xcols, LU.Data(), &lurows, X.Data(), &xrows, &info);
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of DPOSV is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "The leading minor of order %i (and hence the matrix A itself) is not positive definite, so the factorization could not be completed, and the solution has not been computed!", info);
                        Utility::RunTimeError(s);
                }
        }
        return X;
}

Matrix<ComplexFloat> Lapack::Posv(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B)
{
        if(A.Columns() != A.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        if(A.Columns() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }
        int info;
        Matrix<ComplexFloat> LU = A.Clone();
        Matrix<ComplexFloat> X = B.Clone();
        int xrows = X.Rows();
        int xcols = X.Columns();
        int lurows = LU.Rows();
        CPOSV("U", &xrows, &xcols, reinterpret_cast<MKL_Complex8*>(LU.Data()), &lurows, reinterpret_cast<MKL_Complex8*>(X.Data()), &xrows, &info);
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of CPOSV is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "The leading minor of order %i (and hence the matrix A itself) is not positive definite, so the factorization could not be completed, and the solution has not been computed!", info);
                        Utility::RunTimeError(s);
                }
        }
        return X;
}

Matrix<ComplexDouble> Lapack::Posv(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B)
{
        if(A.Columns() != A.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        if(A.Columns() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }
        int info;
        Matrix<ComplexDouble> LU = A.Clone();
        Matrix<ComplexDouble> X = B.Clone();
        int xrows = X.Rows();
        int xcols = X.Columns();
        int lurows = LU.Rows();
        ZPOSV("U", &xrows, &xcols, reinterpret_cast<MKL_Complex16*>(LU.Data()), &lurows, reinterpret_cast<MKL_Complex16*>(X.Data()), &xrows, &info);
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of ZPOSV is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "The leading minor of order %i (and hence the matrix A itself) is not positive definite, so the factorization could not be completed, and the solution has not been computed!", info);
                        Utility::RunTimeError(s);
                }
        }
        return X;
}

//Posvx


//Sysv
Matrix<float> Lapack::Sysv(Matrix<float>& A, Matrix<float>& B)
{
        if(A.Columns() != A.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        if(A.Columns() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }
        int info;
        Vector<int> ipiv(B.Rows());
        Matrix<float> LU = A.Clone();
        Matrix<float> X = B.Clone();
        int xrows = X.Rows();
        int xcols = X.Columns();
        int lurows = LU.Rows();
        int bs = 64*lurows;
        Vector<float> work(bs);
        SSYSV("U", &xrows, &xcols, LU.Data(), &lurows, ipiv.Data(), X.Data(), &xrows, work.Data(), &bs, &info);

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of SSYSV is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "d(%i,%i) is 0. The factorization has been completed, but D is exactly singular, so the solution could not be computed!", info, info);
                        Utility::RunTimeError(s);
                }
        }
        return X;
}

Matrix<double> Lapack::Sysv(Matrix<double>& A, Matrix<double>& B)
{
        if(A.Columns() != A.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        if(A.Columns() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }
        int info;
        Vector<int> ipiv(B.Rows());
        Matrix<double> LU = A.Clone();
        Matrix<double> X = B.Clone();
        int xrows = X.Rows();
        int xcols = X.Columns();
        int lurows = LU.Rows();
        int bs = 64*lurows;
        Vector<double> work(bs);
        DSYSV("U", &xrows, &xcols, LU.Data(), &lurows, ipiv.Data(), X.Data(), &xrows, work.Data(), &bs, &info);

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of DSYSV is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "d(%i,%i) is 0. The factorization has been completed, but D is exactly singular, so the solution could not be computed!", info, info);
                        Utility::RunTimeError(s);
                }
        }
        return X;
}

Matrix<ComplexFloat> Lapack::Sysv(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B)
{
        if(A.Columns() != A.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        if(A.Columns() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }
        int info;
        Vector<int> ipiv(B.Rows());
        Matrix<ComplexFloat> LU = A.Clone();
        Matrix<ComplexFloat> X = B.Clone();
        int xrows = X.Rows();
        int xcols = X.Columns();
        int lurows = LU.Rows();
        int bs = 64*lurows;
        Vector<ComplexFloat> work(bs);
        CSYSV("U", &xrows, &xcols, reinterpret_cast<MKL_Complex8*>(LU.Data()), &lurows, ipiv.Data(), reinterpret_cast<MKL_Complex8*>(X.Data()), &xrows, reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, &info);

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of CSYSV is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "d(%i,%i) is 0. The factorization has been completed, but D is exactly singular, so the solution could not be computed!", info, info);
                        Utility::RunTimeError(s);
                }
        }
        return X;
}

Matrix<ComplexDouble> Lapack::Sysv(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B)
{
        if(A.Columns() != A.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        if(A.Columns() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }
        int info;
        Vector<int> ipiv(B.Rows());
        Matrix<ComplexDouble> LU = A.Clone();
        Matrix<ComplexDouble> X = B.Clone();
        int xrows = X.Rows();
        int xcols = X.Columns();
        int lurows = LU.Rows();
        int bs = 64*lurows;
        Vector<ComplexDouble> work(bs);
        ZSYSV("U", &xrows, &xcols, reinterpret_cast<MKL_Complex16*>(LU.Data()), &lurows, ipiv.Data(), reinterpret_cast<MKL_Complex16*>(X.Data()), &xrows, reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, &info);

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of ZSYSV is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "d(%i,%i) is 0. The factorization has been completed, but D is exactly singular, so the solution could not be computed!", info, info);
                        Utility::RunTimeError(s);
                }
        }
        return X;
}

//Sysvx




//Hesv
Matrix<ComplexFloat> Lapack::Hesv(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B)
{
        if(A.Columns() != A.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        if(A.Columns() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }
        int info;
        Vector<int> ipiv(B.Rows());
        Matrix<ComplexFloat> LU = A.Clone();
        Matrix<ComplexFloat> X = B.Clone();
        int xrows = X.Rows();
        int xcols = X.Columns();
        int lurows = LU.Rows();
        int bs = 64*lurows;
        Vector<ComplexFloat> work(bs);
        CHESV("U", &xrows, &xcols, reinterpret_cast<MKL_Complex8*>(LU.Data()), &lurows, ipiv.Data(), reinterpret_cast<MKL_Complex8*>(X.Data()), &xrows, reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, &info);

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of CHESV is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "d(%i,%i) is 0. The factorization has been completed, but D is exactly singular, so the solution could not be computed!", info, info);
                        Utility::RunTimeError(s);
                }
        }
        return X;
}

Matrix<ComplexDouble> Lapack::Hesv(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B)
{
        if(A.Columns() != A.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        if(A.Columns() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }
        int info;
        Vector<int> ipiv(B.Rows());
        Matrix<ComplexDouble> LU = A.Clone();
        Matrix<ComplexDouble> X = B.Clone();
        int xrows = X.Rows();
        int xcols = X.Columns();
        int lurows = LU.Rows();
        int bs = 64*lurows;
        Vector<ComplexDouble> work(bs);
        ZHESV("U", &xrows, &xcols, reinterpret_cast<MKL_Complex16*>(LU.Data()), &lurows, ipiv.Data(), reinterpret_cast<MKL_Complex16*>(X.Data()), &xrows, reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, &info);

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of ZHESV is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "d(%i,%i) is 0. The factorization has been completed, but D is exactly singular, so the solution could not be computed!", info, info);
                        Utility::RunTimeError(s);
                }
        }
        return X;
}


//Hesvx















//=========================================================
// Linear Least Squares problems
//=========================================================



// Gels
Matrix<float> Lapack::Gels(Matrix<float>& A, Matrix<float>& B)
{
        if(A.Rows() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }

        int info;
        Matrix<float> F = A.Clone();
        Matrix<float> X;
        if(B.Rows() >= A.Columns())
        {
                X = B.Clone();
        }
        else
        {
                X = Matrix<float>(A.Columns(), B.Columns());
                X.ReadFromMatrix(B);
        }
        int m = F.Rows();
        int nrhs = X.Columns();
        int n = F.Columns();
        int ldb = X.Rows();

        int min = (m>n)?n:m;
        int max = (m>n)?m:n;
        max = (max>nrhs)?max:nrhs;
        int bs = min+64*max;
        Vector<float> work(bs);
        SGELS("N", &m, &n, &nrhs, F.Data(), &m, X.Data(), &ldb, work.Data(), &bs, &info);
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of SGELS is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        Utility::RunTimeError("Something went wrong when running SGELS!");
                }
        }
        if(B.Rows() > A.Columns())
        {
                X = X(0, A.Columns()-1, 0, nrhs-1);
        }
        return X;
}

Matrix<double> Lapack::Gels(Matrix<double>& A, Matrix<double>& B)
{
        if(A.Rows() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }

        int info;
        Matrix<double> F = A.Clone();
        Matrix<double> X;
        if(B.Rows() >= A.Columns())
        {
                X = B.Clone();
        }
        else
        {
                X = Matrix<double>(A.Columns(), B.Columns());
                X.ReadFromMatrix(B);
        }
        int m = F.Rows();
        int nrhs = X.Columns();
        int n = F.Columns();
        int ldb = X.Rows();

        int min = (m>n)?n:m;
        int max = (m>n)?m:n;
        max = (max>nrhs)?max:nrhs;
        int bs = min+64*max;
        Vector<double> work(bs);
        DGELS("N", &m, &n, &nrhs, F.Data(), &m, X.Data(), &ldb, work.Data(), &bs, &info);
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of DGELS is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        Utility::RunTimeError("Something went wrong when running DGELS!");
                }
        }
        if(B.Rows() > A.Columns())
        {
                X = X(0, A.Columns()-1, 0, nrhs-1);
        }
        return X;
}

Matrix<ComplexFloat> Lapack::Gels(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B)
{
        if(A.Rows() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }

        int info;
        Matrix<ComplexFloat> F = A.Clone();
        Matrix<ComplexFloat> X;
        if(B.Rows() >= A.Columns())
        {
                X = B.Clone();
        }
        else
        {
                X = Matrix<ComplexFloat>(A.Columns(), B.Columns());
                X.ReadFromMatrix(B);
        }
        int m = F.Rows();
        int nrhs = X.Columns();
        int n = F.Columns();
        int ldb = X.Rows();

        int min = (m>n)?n:m;
        int max = (m>n)?m:n;
        max = (max>nrhs)?max:nrhs;
        int bs = min+64*max;
        Vector<ComplexFloat> work(bs);
        CGELS("N", &m, &n, &nrhs, reinterpret_cast<MKL_Complex8*>(F.Data()), &m, reinterpret_cast<MKL_Complex8*>(X.Data()), &ldb, reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, &info);
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of CGELS is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        Utility::RunTimeError("Something went wrong when running CGELS!");
                }
        }
        if(B.Rows() > A.Columns())
        {
                X = X(0, A.Columns()-1, 0, nrhs-1);
        }
        return X;
}

Matrix<ComplexDouble> Lapack::Gels(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B)
{
        if(A.Rows() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }

        int info;
        Matrix<ComplexDouble> F = A.Clone();
        Matrix<ComplexDouble> X;
        if(B.Rows() >= A.Columns())
        {
                X = B.Clone();
        }
        else
        {
                X = Matrix<ComplexDouble>(A.Columns(), B.Columns());
                X.ReadFromMatrix(B);
        }
        int m = F.Rows();
        int nrhs = X.Columns();
        int n = F.Columns();
        int ldb = X.Rows();

        int min = (m>n)?n:m;
        int max = (m>n)?m:n;
        max = (max>nrhs)?max:nrhs;
        int bs = min+64*max;
        Vector<ComplexDouble> work(bs);
        ZGELS("N", &m, &n, &nrhs, reinterpret_cast<MKL_Complex16*>(F.Data()), &m, reinterpret_cast<MKL_Complex16*>(X.Data()), &ldb, reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, &info);
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of ZGELS is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        Utility::RunTimeError("Something went wrong when running ZGELS!");
                }
        }

        if(B.Rows() > A.Columns())
        {
                X = X(0, A.Columns()-1, 0, nrhs-1);
        }

        return X;
}


// Gelsy
Matrix<float> Lapack::Gelsy(Matrix<float>& A, Matrix<float>& B)
{
        if(A.Rows() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }

        int info;
        int rank;
        Matrix<float> F = A.Clone();
        Matrix<float> X;
        if(B.Rows() >= A.Columns())
        {
                X = B.Clone();
        }
        else
        {
                X = Matrix<float>(A.Columns(), B.Columns());
                X.ReadFromMatrix(B);
        }
        int m = F.Rows();
        int nrhs = X.Columns();
        int n = F.Columns();
        int ldb = X.Rows();

        int min = (m>n)?n:m;
        int max = (m>n)?m:n;
        int bs = min+2*n+64*(n+1);
        int temp = 2*min+64*nrhs;
        bs = (bs>temp)?bs:temp;
        Vector<float> work(bs);
        Vector<int> jpvt(n);
        float rcond = max*numeric_limits<float>::epsilon();

        SGELSY(&m, &n, &nrhs, F.Data(), &m, X.Data(), &ldb, jpvt.Data(), &rcond, &rank, work.Data(), &bs, &info);
        if(rank < min && info == 0)
        {
                        char s[500];
                        sprintf(s, "Matrix is rank deficient (rank is %i)!", rank);
                        Utility::Warning(s);
        }
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of SGELSY is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        Utility::RunTimeError("Something went wrong when running SGELSY!");
                }
        }
        if(B.Rows() > A.Columns())
        {
                X = X(0, A.Columns()-1, 0, nrhs-1);
        }
        return X;
}



Matrix<double> Lapack::Gelsy(Matrix<double>& A, Matrix<double>& B)
{
        if(A.Rows() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }

        int info;
        int rank;
        Matrix<double> F = A.Clone();
        Matrix<double> X;
        if(B.Rows() >= A.Columns())
        {
                X = B.Clone();
        }
        else
        {
                X = Matrix<double>(A.Columns(), B.Columns());
                X.ReadFromMatrix(B);
        }
        int m = F.Rows();
        int nrhs = X.Columns();
        int n = F.Columns();
        int ldb = X.Rows();

        int min = (m>n)?n:m;
        int max = (m>n)?m:n;
        int bs = min+2*n+64*(n+1);
        int temp = 2*min+64*nrhs;
        bs = (bs>temp)?bs:temp;
        Vector<double> work(bs);
        Vector<int> jpvt(n);
        double rcond = max*numeric_limits<double>::epsilon();

        DGELSY(&m, &n, &nrhs, F.Data(), &m, X.Data(), &ldb, jpvt.Data(), &rcond, &rank, work.Data(), &bs, &info);
        if(rank < min && info == 0)
        {
                        char s[500];
                        sprintf(s, "Matrix is rank deficient (rank is %i)!", rank);
                        Utility::Warning(s);
        }
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of DGELSY is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        Utility::RunTimeError("Something went wrong when running DGELSY!");
                }
        }
        if(B.Rows() > A.Columns())
        {
                X = X(0, A.Columns()-1, 0, nrhs-1);
        }
        return X;
}

Matrix<ComplexFloat> Lapack::Gelsy(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B)
{
        if(A.Rows() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }

        int info;
        int rank;
        Matrix<ComplexFloat> F = A.Clone();
        Matrix<ComplexFloat> X;
        if(B.Rows() >= A.Columns())
        {
                X = B.Clone();
        }
        else
        {
                X = Matrix<ComplexFloat>(A.Columns(), B.Columns());
                X.ReadFromMatrix(B);
        }
        int m = F.Rows();
        int nrhs = X.Columns();
        int n = F.Columns();
        int ldb = X.Rows();

        int min = (m>n)?n:m;
        int max = (m>n)?m:n;
        int bs = 2*min;
        int temp = 64*(n+1);
        bs = (bs>temp)?bs:temp;
        temp = min+64*min;
        bs = (bs>temp)?bs:temp;
        temp = min+64*nrhs;
        bs = (bs>temp)?bs:temp;
        bs += min;
        Vector<ComplexFloat> work(bs);
        Vector<float> rwork(2*n);
        Vector<int> jpvt(n);
        float rcond = max*numeric_limits<float>::epsilon();

        CGELSY(&m, &n, &nrhs, reinterpret_cast<MKL_Complex8*>(F.Data()), &m, reinterpret_cast<MKL_Complex8*>(X.Data()), &ldb, jpvt.Data(), &rcond, &rank, reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), &info);
        if(rank < min && info == 0)
        {
                        char s[500];
                        sprintf(s, "Matrix is rank deficient (rank is %i)!", rank);
                        Utility::Warning(s);
        }
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of CGELSY is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        Utility::RunTimeError("Something went wrong when running CGELSY!");
                }
        }
        if(B.Rows() > A.Columns())
        {
                X = X(0, A.Columns()-1, 0, nrhs-1);
        }
        return X;
}

Matrix<ComplexDouble> Lapack::Gelsy(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B)
{
        if(A.Rows() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }

        int info;
        int rank;
        Matrix<ComplexDouble> F = A.Clone();
        Matrix<ComplexDouble> X;
        if(B.Rows() >= A.Columns())
        {
                X = B.Clone();
        }
        else
        {
                X = Matrix<ComplexDouble>(A.Columns(), B.Columns());
                X.ReadFromMatrix(B);
        }
        int m = F.Rows();
        int nrhs = X.Columns();
        int n = F.Columns();
        int ldb = X.Rows();

        int min = (m>n)?n:m;
        int max = (m>n)?m:n;
        int bs = 2*min;
        int temp = 64*(n+1);
        bs = (bs>temp)?bs:temp;
        temp = min+64*min;
        bs = (bs>temp)?bs:temp;
        temp = min+64*nrhs;
        bs = (bs>temp)?bs:temp;
        bs += min;
        Vector<ComplexDouble> work(bs);
        Vector<double> rwork(2*n);
        Vector<int> jpvt(n);
        double rcond = max*numeric_limits<double>::epsilon();

        ZGELSY(&m, &n, &nrhs, reinterpret_cast<MKL_Complex16*>(F.Data()), &m, reinterpret_cast<MKL_Complex16*>(X.Data()), &ldb, jpvt.Data(), &rcond, &rank, reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), &info);
        if(rank < min && info == 0)
        {
                        char s[500];
                        sprintf(s, "Matrix is rank deficient (rank is %i)!", rank);
                        Utility::Warning(s);
        }
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of ZGELSY is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        Utility::RunTimeError("Something went wrong when running ZGELSY!");
                }
        }
        if(B.Rows() > A.Columns())
        {
                X = X(0, A.Columns()-1, 0, nrhs-1);
        }
        return X;
}


// Gelss
Matrix<float> Lapack::Gelss(Matrix<float>& A, Matrix<float>& B)
{
        if(A.Rows() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }

        int info;
        int rank;
        Matrix<float> F = A.Clone();
        Matrix<float> X;
        if(B.Rows() >= A.Columns())
        {
                X = B.Clone();
        }
        else
        {
                X = Matrix<float>(A.Columns(), B.Columns());
                X.ReadFromMatrix(B);
        }
        int m = F.Rows();
        int nrhs = X.Columns();
        int n = F.Columns();
        int ldb = X.Rows();

        int min = (m>n)?n:m;
        int max = (m>n)?m:n;
        int bs = max;
        int temp = 2*min;
        bs = (bs>temp)?bs:temp;
        bs = (bs>nrhs)?bs:nrhs;
        bs += 3*min;
        Vector<float> work(bs);
        Vector<float> s(min);
        float rcond = max*numeric_limits<float>::epsilon();

        SGELSS(&m, &n, &nrhs, F.Data(), &m, X.Data(), &ldb, s.Data(), &rcond, &rank, work.Data(), &bs, &info);
        if(rank < min && info == 0)
        {
                        char s[500];
                        sprintf(s, "Matrix is rank deficient (rank is %i)!", rank);
                        Utility::Warning(s);
        }
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of SGELSS is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "The algorithm for computing the SVD failed to converge; %i off-diagonal elements of an intermediate bidiagonal form did not converge to zero!", info);
                        Utility::RunTimeError(s);
                }
        }
        if(B.Rows() > A.Columns())
        {
                X = X(0, A.Columns()-1, 0, nrhs-1);
        }
        return X;
}

Matrix<double> Lapack::Gelss(Matrix<double>& A, Matrix<double>& B)
{
        if(A.Rows() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }

        int info;
        int rank;
        Matrix<double> F = A.Clone();
        Matrix<double> X;
        if(B.Rows() >= A.Columns())
        {
                X = B.Clone();
        }
        else
        {
                X = Matrix<double>(A.Columns(), B.Columns());
                X.ReadFromMatrix(B);
        }
        int m = F.Rows();
        int nrhs = X.Columns();
        int n = F.Columns();
        int ldb = X.Rows();

        int min = (m>n)?n:m;
        int max = (m>n)?m:n;
        int bs = max;
        int temp = 2*min;
        bs = (bs>temp)?bs:temp;
        bs = (bs>nrhs)?bs:nrhs;
        bs += 3*min;
        Vector<double> work(bs);
        Vector<double> s(min);
        double rcond = max*numeric_limits<double>::epsilon();

        DGELSS(&m, &n, &nrhs, F.Data(), &m, X.Data(), &ldb, s.Data(), &rcond, &rank, work.Data(), &bs, &info);
        if(rank < min && info == 0)
        {
                        char s[500];
                        sprintf(s, "Matrix is rank deficient (rank is %i)!", rank);
                        Utility::Warning(s);
        }
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of DGELSS is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "The algorithm for computing the SVD failed to converge; %i off-diagonal elements of an intermediate bidiagonal form did not converge to zero!", info);
                        Utility::RunTimeError(s);
                }
        }
        if(B.Rows() > A.Columns())
        {
                X = X(0, A.Columns()-1, 0, nrhs-1);
        }
        return X;
}

Matrix<ComplexFloat> Lapack::Gelss(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B)
{
        if(A.Rows() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }

        int info;
        int rank;
        Matrix<ComplexFloat> F = A.Clone();
        Matrix<ComplexFloat> X;
        if(B.Rows() >= A.Columns())
        {
                X = B.Clone();
        }
        else
        {
                X = Matrix<ComplexFloat>(A.Columns(), B.Columns());
                X.ReadFromMatrix(B);
        }
        int m = F.Rows();
        int nrhs = X.Columns();
        int n = F.Columns();
        int ldb = X.Rows();

        int min = (m>n)?n:m;
        int max = (m>n)?m:n;
        int bs = (max>nrhs)?max:nrhs;
        bs += 2*min;
        Vector<ComplexFloat> work(bs);
        Vector<float> s(min);
        Vector<float> rwork(5*min);
        float rcond = max*numeric_limits<float>::epsilon();

        CGELSS(&m, &n, &nrhs, reinterpret_cast<MKL_Complex8*>(F.Data()), &m, reinterpret_cast<MKL_Complex8*>(X.Data()), &ldb, s.Data(), &rcond, &rank, reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), &info);
        if(rank < min && info == 0)
        {
                        char s[500];
                        sprintf(s, "Matrix is rank deficient (rank is %i)!", rank);
                        Utility::Warning(s);
        }
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of CGELSS is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "The algorithm for computing the SVD failed to converge; %i off-diagonal elements of an intermediate bidiagonal form did not converge to zero!", info);
                        Utility::RunTimeError(s);
                }
        }
        if(B.Rows() > A.Columns())
        {
                X = X(0, A.Columns()-1, 0, nrhs-1);
        }
        return X;
}

Matrix<ComplexDouble> Lapack::Gelss(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B)
{
        if(A.Rows() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }

        int info;
        int rank;
        Matrix<ComplexDouble> F = A.Clone();
        Matrix<ComplexDouble> X;
        if(B.Rows() >= A.Columns())
        {
                X = B.Clone();
        }
        else
        {
                X = Matrix<ComplexDouble>(A.Columns(), B.Columns());
                X.ReadFromMatrix(B);
        }
        int m = F.Rows();
        int nrhs = X.Columns();
        int n = F.Columns();
        int ldb = X.Rows();

        int min = (m>n)?n:m;
        int max = (m>n)?m:n;
        int bs = (max>nrhs)?max:nrhs;
        bs += 2*min;
        Vector<ComplexDouble> work(bs);
        Vector<double> s(min);
        Vector<double> rwork(5*min);
        double rcond = max*numeric_limits<double>::epsilon();

        ZGELSS(&m, &n, &nrhs, reinterpret_cast<MKL_Complex16*>(F.Data()), &m, reinterpret_cast<MKL_Complex16*>(X.Data()), &ldb, s.Data(), &rcond, &rank, reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), &info);
        if(rank < min && info == 0)
        {
                        char s[500];
                        sprintf(s, "Matrix is rank deficient (rank is %i)!", rank);
                        Utility::Warning(s);
        }
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of ZGELSS is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "The algorithm for computing the SVD failed to converge; %i off-diagonal elements of an intermediate bidiagonal form did not converge to zero!", info);
                        Utility::RunTimeError(s);
                }
        }
        if(B.Rows() > A.Columns())
        {
                X = X(0, A.Columns()-1, 0, nrhs-1);
        }
        return X;
}

// Gelsd
Matrix<float> Lapack::Gelsd(Matrix<float>& A, Matrix<float>& B)
{
        if(A.Rows() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }

        int info;
        int rank;
        Matrix<float> F = A.Clone();
        Matrix<float> X;
        if(B.Rows() >= A.Columns())
        {
                X = B.Clone();
        }
        else
        {
                X = Matrix<float>(A.Columns(), B.Columns());
                X.ReadFromMatrix(B);
        }
        int m = F.Rows();
        int nrhs = X.Columns();
        int n = F.Columns();
        int ldb = X.Rows();

        int min = (m>n)?n:m;
        int max = (m>n)?m:n;
        int smlsiz = 25;
        int nlvl = (int)(log((double)min/(double)(smlsiz+1))/log(2.0))+1;
        nlvl = (nlvl>0)?nlvl:0;
        int bs = 0;
        if(m>=n)
        {
                bs = 12*n + 2*n*smlsiz + 8*n*nlvl + n*nrhs + (smlsiz+1)*(smlsiz+1);
        }
        else
        {
                bs = 12*m + 2*m*smlsiz + 8*m*nlvl + m*nrhs + (smlsiz+1)*(smlsiz+1);
        }
        Vector<float> work(bs);
        int iwork_size = 3*min*nlvl + 11*min;
        Vector<int> iwork(iwork_size);
        Vector<float> s(min);
        float rcond = max*numeric_limits<float>::epsilon();

        SGELSD(&m, &n, &nrhs, F.Data(), &m, X.Data(), &ldb, s.Data(), &rcond, &rank, work.Data(), &bs, iwork.Data(), &info);

        if(rank < min && info == 0)
        {
                        char s[500];
                        sprintf(s, "Matrix is rank deficient (rank is %i)!", rank);
                        Utility::Warning(s);
        }
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of SGELSD is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "The algorithm for computing the SVD failed to converge; %i off-diagonal elements of an intermediate bidiagonal form did not converge to zero!", info);
                        Utility::RunTimeError(s);
                }
        }
        if(B.Rows() > A.Columns())
        {
                X = X(0, A.Columns()-1, 0, nrhs-1);
        }
        return X;
}

Matrix<double> Lapack::Gelsd(Matrix<double>& A, Matrix<double>& B)
{
        if(A.Rows() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }

        int info;
        int rank;
        Matrix<double> F = A.Clone();
        Matrix<double> X;
        if(B.Rows() >= A.Columns())
        {
                X = B.Clone();
        }
        else
        {
                X = Matrix<double>(A.Columns(), B.Columns());
                X.ReadFromMatrix(B);
        }
        int m = F.Rows();
        int nrhs = X.Columns();
        int n = F.Columns();
        int ldb = X.Rows();

        int min = (m>n)?n:m;
        int max = (m>n)?m:n;
        int smlsiz = 25;
        int nlvl = (int)(log((double)min/(double)(smlsiz+1))/log(2.0))+1;
        nlvl = (nlvl>0)?nlvl:0;
        int bs = 0;
        if(m>=n)
        {
                bs = 12*n + 2*n*smlsiz + 8*n*nlvl + n*nrhs + (smlsiz+1)*(smlsiz+1);
        }
        else
        {
                bs = 12*m + 2*m*smlsiz + 8*m*nlvl + m*nrhs + (smlsiz+1)*(smlsiz+1);
        }
        Vector<double> work(bs);
        int iwork_size = 3*min*nlvl + 11*min;
        Vector<int> iwork(iwork_size);
        Vector<double> s(min);
        double rcond = max*numeric_limits<double>::epsilon();

        DGELSD(&m, &n, &nrhs, F.Data(), &m, X.Data(), &ldb, s.Data(), &rcond, &rank, work.Data(), &bs, iwork.Data(), &info);

        if(rank < min && info == 0)
        {
                        char s[500];
                        sprintf(s, "Matrix is rank deficient (rank is %i)!", rank);
                        Utility::Warning(s);
        }
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of DGELSD is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "The algorithm for computing the SVD failed to converge; %i off-diagonal elements of an intermediate bidiagonal form did not converge to zero!", info);
                        Utility::RunTimeError(s);
                }
        }
        if(B.Rows() > A.Columns())
        {
                X = X(0, A.Columns()-1, 0, nrhs-1);
        }
        return X;
}

Matrix<ComplexFloat> Lapack::Gelsd(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B)
{
        if(A.Rows() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }

        int info;
        int rank;
        Matrix<ComplexFloat> F = A.Clone();
        Matrix<ComplexFloat> X;
        if(B.Rows() >= A.Columns())
        {
                X = B.Clone();
        }
        else
        {
                X = Matrix<ComplexFloat>(A.Columns(), B.Columns());
                X.ReadFromMatrix(B);
        }
        int m = F.Rows();
        int nrhs = X.Columns();
        int n = F.Columns();
        int ldb = X.Rows();

        int min = (m>n)?n:m;
        int max = (m>n)?m:n;
        int smlsiz = 25;
        int nlvl = (int)(log((double)min/(double)(smlsiz+1))/log(2.0))+1;
        nlvl = (nlvl>0)?nlvl:0;
        int bs = -1;

        Vector<ComplexFloat> work(1);
        int lrwork = 0;
        if(m>=n)
        {
                lrwork = 10*n + 2*n*smlsiz + 8*n*nlvl + 3*smlsiz*nrhs + (smlsiz+1)*(smlsiz+1);
        }
        else
        {
                lrwork = 10*m + 2*m*smlsiz + 8*m*nlvl + 3*smlsiz*nrhs + (smlsiz+1)*(smlsiz+1);
        }
        Vector<float> rwork(lrwork);
        int iwork_size = 3*min*nlvl + 11*min;
        Vector<int> iwork(iwork_size);
        Vector<float> s(min);
        float rcond = max*numeric_limits<float>::epsilon();

        CGELSD(&m, &n, &nrhs, reinterpret_cast<MKL_Complex8*>(F.Data()), &m, reinterpret_cast<MKL_Complex8*>(X.Data()), &ldb, s.Data(), &rcond, &rank, reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), iwork.Data(), &info);
        bs = (int)real(work[0]);
        work = Vector<ComplexFloat>(bs);
        CGELSD(&m, &n, &nrhs, reinterpret_cast<MKL_Complex8*>(F.Data()), &m, reinterpret_cast<MKL_Complex8*>(X.Data()), &ldb, s.Data(), &rcond, &rank, reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), iwork.Data(), &info);

        if(rank < min && info == 0)
        {
                        char s[500];
                        sprintf(s, "Matrix is rank deficient (rank is %i)!", rank);
                        Utility::Warning(s);
        }
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of CGELSD is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "The algorithm for computing the SVD failed to converge; %i off-diagonal elements of an intermediate bidiagonal form did not converge to zero!", info);
                        Utility::RunTimeError(s);
                }
        }
        if(B.Rows() > A.Columns())
        {
                X = X(0, A.Columns()-1, 0, nrhs-1);
        }
        return X;
}

Matrix<ComplexDouble> Lapack::Gelsd(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B)
{
        if(A.Rows() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }

        int info;
        int rank;
        Matrix<ComplexDouble> F = A.Clone();
        Matrix<ComplexDouble> X;
        if(B.Rows() >= A.Columns())
        {
                X = B.Clone();
        }
        else
        {
                X = Matrix<ComplexDouble>(A.Columns(), B.Columns());
                X.ReadFromMatrix(B);
        }
        int m = F.Rows();
        int nrhs = X.Columns();
        int n = F.Columns();
        int ldb = X.Rows();

        int min = (m>n)?n:m;
        int max = (m>n)?m:n;
        int smlsiz = 25;
        int nlvl = (int)(log((double)min/(double)(smlsiz+1))/log(2.0))+1;
        nlvl = (nlvl>0)?nlvl:0;
        int bs = -1;

        Vector<ComplexDouble> work(1);
        int lrwork = 0;
        if(m>=n)
        {
                lrwork = 10*n + 2*n*smlsiz + 8*n*nlvl + 3*smlsiz*nrhs + (smlsiz+1)*(smlsiz+1);
        }
        else
        {
                lrwork = 10*m + 2*m*smlsiz + 8*m*nlvl + 3*smlsiz*nrhs + (smlsiz+1)*(smlsiz+1);
        }
        Vector<double> rwork(lrwork);
        int iwork_size = 3*min*nlvl + 11*min;
        Vector<int> iwork(iwork_size);
        Vector<double> s(min);
        double rcond = max*numeric_limits<double>::epsilon();
        
        ZGELSD(&m, &n, &nrhs, reinterpret_cast<MKL_Complex16*>(F.Data()), &m, reinterpret_cast<MKL_Complex16*>(X.Data()), &ldb, s.Data(), &rcond, &rank, reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), iwork.Data(), &info);
        bs = (int)real(work[0]);
        work = Vector<ComplexDouble>(bs);
        ZGELSD(&m, &n, &nrhs, reinterpret_cast<MKL_Complex16*>(F.Data()), &m, reinterpret_cast<MKL_Complex16*>(X.Data()), &ldb, s.Data(), &rcond, &rank, reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), iwork.Data(), &info);

        if(rank < min && info == 0)
        {
                        char s[500];
                        sprintf(s, "Matrix is rank deficient (rank is %i)!", rank);
                        Utility::Warning(s);
        }
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of ZGELSD is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "The algorithm for computing the SVD failed to converge; %i off-diagonal elements of an intermediate bidiagonal form did not converge to zero!", info);
                        Utility::RunTimeError(s);
                }
        }
        if(B.Rows() > A.Columns())
        {
                X = X(0, A.Columns()-1, 0, nrhs-1);
        }
        return X;
}






//=========================================================
// Generalized Linear Least Squares problems
//=========================================================


// Gglse
Vector<float> Lapack::Gglse(Matrix<float>& A, Matrix<float>& B, Vector<float>& c, Vector<float>& d)
{
        if(A.Columns() != B.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }
        if(A.Rows() != c.Length() || B.Rows() != d.Length())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix/Vector dimensions does not match!");
        }

        if(B.Rows() > A.Columns() || (B.Rows() + A.Rows()) < A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix/Vector dimensions have incompatibilities!");
        }

        int info;
        Matrix<float> A2 = A.Clone();
        Matrix<float> B2 = B.Clone();
        Vector<float> c2 = c.Clone();
        Vector<float> d2 = d.Clone();

        int m = A2.Rows();
        int n = A2.Columns();
        int p = B2.Rows();
        Vector<float> x(n);

        int min = (m>n)?n:m;
        int max = (m>n)?m:n;
        int nb = 64;
        int bs = p + min + max*nb;
        Vector<float> work(bs);

        SGGLSE(&m, &n, &p, A2.Data(), &m, B2.Data(), &p, c2.Data(), d2.Data(), x.Data(), work.Data(), &bs, &info);

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of SGGLSE is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "Problem experienced when running SGGLSE!");
                        Utility::RunTimeError(s);
                }
        }
        return x;
}



Vector<double> Lapack::Gglse(Matrix<double>& A, Matrix<double>& B, Vector<double>& c, Vector<double>& d)
{
        if(A.Columns() != B.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }
        if(A.Rows() != c.Length() || B.Rows() != d.Length())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix/Vector dimensions does not match!");
        }

        if(B.Rows() > A.Columns() || (B.Rows() + A.Rows()) < A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix/Vector dimensions have incompatibilities!");
        }

        int info;
        Matrix<double> A2 = A.Clone();
        Matrix<double> B2 = B.Clone();
        Vector<double> c2 = c.Clone();
        Vector<double> d2 = d.Clone();

        int m = A2.Rows();
        int n = A2.Columns();
        int p = B2.Rows();
        Vector<double> x(n);

        int min = (m>n)?n:m;
        int max = (m>n)?m:n;
        int nb = 64;
        int bs = p + min + max*nb;
        Vector<double> work(bs);

        DGGLSE(&m, &n, &p, A2.Data(), &m, B2.Data(), &p, c2.Data(), d2.Data(), x.Data(), work.Data(), &bs, &info);

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of DGGLSE is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "Problem experienced when running DGGLSE!");
                        Utility::RunTimeError(s);
                }
        }
        return x;
}


Vector<ComplexFloat> Lapack::Gglse(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B, Vector<ComplexFloat>& c, Vector<ComplexFloat>& d)
{
        if(A.Columns() != B.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }
        if(A.Rows() != c.Length() || B.Rows() != d.Length())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix/Vector dimensions does not match!");
        }

        if(B.Rows() > A.Columns() || (B.Rows() + A.Rows()) < A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix/Vector dimensions have incompatibilities!");
        }

        int info;
        Matrix<ComplexFloat> A2 = A.Clone();
        Matrix<ComplexFloat> B2 = B.Clone();
        Vector<ComplexFloat> c2 = c.Clone();
        Vector<ComplexFloat> d2 = d.Clone();

        int m = A2.Rows();
        int n = A2.Columns();
        int p = B2.Rows();
        Vector<ComplexFloat> x(n);

        int min = (m>n)?n:m;
        int max = (m>n)?m:n;
        int nb = 64;
        int bs = p + min + max*nb;
        Vector<ComplexFloat> work(bs);

        CGGLSE(&m, &n, &p, reinterpret_cast<MKL_Complex8*>(A2.Data()), &m, reinterpret_cast<MKL_Complex8*>(B2.Data()), &p, reinterpret_cast<MKL_Complex8*>(c2.Data()), reinterpret_cast<MKL_Complex8*>(d2.Data()), reinterpret_cast<MKL_Complex8*>(x.Data()), reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, &info);

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of CGGLSE is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "Problem experienced when running CGGLSE!");
                        Utility::RunTimeError(s);
                }
        }
        return x;
}


Vector<ComplexDouble> Lapack::Gglse(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B, Vector<ComplexDouble>& c, Vector<ComplexDouble>& d)
{
        if(A.Columns() != B.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }
        if(A.Rows() != c.Length() || B.Rows() != d.Length())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix/Vector dimensions does not match!");
        }

        if(B.Rows() > A.Columns() || (B.Rows() + A.Rows()) < A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix/Vector dimensions have incompatibilities!");
        }

        int info;
        Matrix<ComplexDouble> A2 = A.Clone();
        Matrix<ComplexDouble> B2 = B.Clone();
        Vector<ComplexDouble> c2 = c.Clone();
        Vector<ComplexDouble> d2 = d.Clone();

        int m = A2.Rows();
        int n = A2.Columns();
        int p = B2.Rows();
        Vector<ComplexDouble> x(n);

        int min = (m>n)?n:m;
        int max = (m>n)?m:n;
        int nb = 64;
        int bs = p + min + max*nb;
        Vector<ComplexDouble> work(bs);

        ZGGLSE(&m, &n, &p, reinterpret_cast<MKL_Complex16*>(A2.Data()), &m, reinterpret_cast<MKL_Complex16*>(B2.Data()), &p, reinterpret_cast<MKL_Complex16*>(c2.Data()), reinterpret_cast<MKL_Complex16*>(d2.Data()), reinterpret_cast<MKL_Complex16*>(x.Data()), reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, &info);

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of ZGGLSE is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "Problem experienced when running ZGGLSE!");
                        Utility::RunTimeError(s);
                }
        }
        return x;
}





// Ggglm
void Lapack::Ggglm(Matrix<float>& A, Matrix<float>& B, Vector<float>& d, Vector<float>& x, Vector<float>& y)
{
        if(A.Rows() != B.Rows() || B.Columns() < (A.Rows() - A.Columns()))
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }

        if(d.Length() != A.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix/Vector dimensions have incompatibilities!");
        }

        int info;
        Matrix<float> A2 = A.Clone();
        Matrix<float> B2 = B.Clone();
        Vector<float> d2 = d.Clone();

        int m = A2.Columns();
        int n = A2.Rows();
        int p = B2.Columns();
        x = Vector<float>(m);
        y = Vector<float>(p);

        int min = (n>p)?p:n;
        int max = (n>p)?n:p;
        int nb = 64;
        int bs = m + min + max*nb;
        Vector<float> work(bs);

        SGGGLM(&n, &m, &p, A2.Data(), &n, B2.Data(), &n, d2.Data(), x.Data(), y.Data(), work.Data(), &bs, &info);

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of SGGGLM is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "Problem experienced when running SGGGLM!");
                        Utility::RunTimeError(s);
                }
        }

}


void Lapack::Ggglm(Matrix<double>& A, Matrix<double>& B, Vector<double>& d, Vector<double>& x, Vector<double>& y)
{
        if(A.Rows() != B.Rows() || B.Columns() < (A.Rows() - A.Columns()))
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }

        if(d.Length() != A.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix/Vector dimensions have incompatibilities!");
        }

        int info;
        Matrix<double> A2 = A.Clone();
        Matrix<double> B2 = B.Clone();
        Vector<double> d2 = d.Clone();

        int m = A2.Columns();
        int n = A2.Rows();
        int p = B2.Columns();
        x = Vector<double>(m);
        y = Vector<double>(p);

        int min = (n>p)?p:n;
        int max = (n>p)?n:p;
        int nb = 64;
        int bs = m + min + max*nb;
        Vector<double> work(bs);

        DGGGLM(&n, &m, &p, A2.Data(), &n, B2.Data(), &n, d2.Data(), x.Data(), y.Data(), work.Data(), &bs, &info);

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of DGGGLM is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "Problem experienced when running DGGGLM!");
                        Utility::RunTimeError(s);
                }
        }
}


void Lapack::Ggglm(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B, Vector<ComplexFloat>& d, Vector<ComplexFloat>& x, Vector<ComplexFloat>& y)
{
        if(A.Rows() != B.Rows() || B.Columns() < (A.Rows() - A.Columns()))
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }

        if(d.Length() != A.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix/Vector dimensions have incompatibilities!");
        }

        int info;
        Matrix<ComplexFloat> A2 = A.Clone();
        Matrix<ComplexFloat> B2 = B.Clone();
        Vector<ComplexFloat> d2 = d.Clone();

        int m = A2.Columns();
        int n = A2.Rows();
        int p = B2.Columns();
        x = Vector<ComplexFloat>(m);
        y = Vector<ComplexFloat>(p);

        int min = (n>p)?p:n;
        int max = (n>p)?n:p;
        int nb = 64;
        int bs = m + min + max*nb;
        Vector<ComplexFloat> work(bs);

        CGGGLM(&n, &m, &p, reinterpret_cast<MKL_Complex8*>(A2.Data()), &n, reinterpret_cast<MKL_Complex8*>(B2.Data()), &n, reinterpret_cast<MKL_Complex8*>(d2.Data()), reinterpret_cast<MKL_Complex8*>(x.Data()), reinterpret_cast<MKL_Complex8*>(y.Data()), reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, &info);

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of CGGGLM is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "Problem experienced when running CGGGLM!");
                        Utility::RunTimeError(s);
                }
        }
}


void Lapack::Ggglm(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B, Vector<ComplexDouble>& d, Vector<ComplexDouble>& x, Vector<ComplexDouble>& y)
{
        if(A.Rows() != B.Rows() || B.Columns() < (A.Rows() - A.Columns()))
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }

        if(d.Length() != A.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix/Vector dimensions have incompatibilities!");
        }

        int info;
        Matrix<ComplexDouble> A2 = A.Clone();
        Matrix<ComplexDouble> B2 = B.Clone();
        Vector<ComplexDouble> d2 = d.Clone();

        int m = A2.Columns();
        int n = A2.Rows();
        int p = B2.Columns();
        x = Vector<ComplexDouble>(m);
        y = Vector<ComplexDouble>(p);

        int min = (n>p)?p:n;
        int max = (n>p)?n:p;
        int nb = 64;
        int bs = m + min + max*nb;
        Vector<ComplexDouble> work(bs);

        ZGGGLM(&n, &m, &p, reinterpret_cast<MKL_Complex16*>(A2.Data()), &n, reinterpret_cast<MKL_Complex16*>(B2.Data()), &n, reinterpret_cast<MKL_Complex16*>(d2.Data()), reinterpret_cast<MKL_Complex16*>(x.Data()), reinterpret_cast<MKL_Complex16*>(y.Data()), reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, &info);

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of ZGGGLM is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "Problem experienced when running ZGGGLM!");
                        Utility::RunTimeError(s);
                }
        }
}






//=========================================================
// Symmetric eigenproblem
//=========================================================



// Syevd or Heevd
Vector<float> Lapack::Syevd(Matrix<float>& A)
{
        if(A.Rows() != A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        
        int info;
        Matrix<float> F = A.Clone();
        int n = F.Rows();
        Vector<float> w(n);

        int bs = 2*n+1;
        Vector<float> work(bs);

        int iwork_size = 1;
        Vector<int> iwork(iwork_size);

        SSYEVD("N", "U", &n, F.Data(), &n, w.Data(), work.Data(), &bs, iwork.Data(), &iwork_size, &info);

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of SSYEVD is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "The algorithm failed to converge; %i elements of an intermediate tridiagonal form did not converge to zero!", info);
                        Utility::RunTimeError(s);
                }
        }
        return w;
}


Vector<float> Lapack::Syevd(Matrix<float>& A, Matrix<float>& EV)
{
        if(A.Rows() != A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        
        int info;
        EV = A.Clone();
        int n = EV.Rows();
        Vector<float> w(n);

        int k = 1;
        if(n>1)
        {
                k = (int)(log((double)(n-1))/log(2.0))+1;
        }
        int bs = 3*n*n+(5+2*k)*n+1;
        Vector<float> work(bs);

        int iwork_size = 5*n+3;
        Vector<int> iwork(iwork_size);

        SSYEVD("V", "U", &n, EV.Data(), &n, w.Data(), work.Data(), &bs, iwork.Data(), &iwork_size, &info);
        
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of SSYEVD is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "The algorithm failed to converge; %i elements of an intermediate tridiagonal form did not converge to zero!", info);
                        Utility::RunTimeError(s);
                }
        }
        return w;
}


Vector<double> Lapack::Syevd(Matrix<double>& A)
{
        if(A.Rows() != A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        
        int info;
        Matrix<double> F = A.Clone();
        int n = F.Rows();
        Vector<double> w(n);

        int bs = 2*n+1;
        Vector<double> work(bs);

        int iwork_size = 1;
        Vector<int> iwork(iwork_size);

        DSYEVD("N", "U", &n, F.Data(), &n, w.Data(), work.Data(), &bs, iwork.Data(), &iwork_size, &info);

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of DSYEVD is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "The algorithm failed to converge; %i elements of an intermediate tridiagonal form did not converge to zero!", info);
                        Utility::RunTimeError(s);
                }
        }
        return w;
}

Vector<double> Lapack::Syevd(Matrix<double>& A, Matrix<double>& EV)
{
        if(A.Rows() != A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        
        int info;
        EV = A.Clone();
        int n = EV.Rows();
        Vector<double> w(n);

        int k = 1;
        if(n>1)
        {
                k = (int)(log((double)(n-1))/log(2.0))+1;
        }
        int bs = 3*n*n+(5+2*k)*n+1;
        Vector<double> work(bs);

        int iwork_size = 5*n+3;
        Vector<int> iwork(iwork_size);

        DSYEVD("V", "U", &n, EV.Data(), &n, w.Data(), work.Data(), &bs, iwork.Data(), &iwork_size, &info);
        
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of DSYEVD is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "The algorithm failed to converge; %i elements of an intermediate tridiagonal form did not converge to zero!", info);
                        Utility::RunTimeError(s);
                }
        }
        return w;
}

Vector<float> Lapack::Heevd(Matrix<ComplexFloat>& A)
{
        if(A.Rows() != A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        
        int info;
        Matrix<ComplexFloat> F = A.Clone();
        int n = F.Rows();
        Vector<float> w(n);

        int bs = n+1;
        Vector<ComplexFloat> work(bs);

        int lrwork = n;
        Vector<float> rwork(lrwork);

        int iwork_size = 1;
        Vector<int> iwork(iwork_size);

        CHEEVD("N", "U", &n, reinterpret_cast<MKL_Complex8*>(F.Data()), &n, w.Data(), reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), &lrwork, iwork.Data(), &iwork_size, &info);

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of CHEEVD is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "The algorithm failed to converge; %i elements of an intermediate tridiagonal form did not converge to zero!", info);
                        Utility::RunTimeError(s);
                }
        }
        return w;
}

Vector<float> Lapack::Heevd(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& EV)
{
        if(A.Rows() != A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        
        int info;
        EV = A.Clone();
        int n = EV.Rows();
        Vector<float> w(n);

        int k = 1;
        if(n>1)
        {
                k = (int)(log((double)(n-1))/log(2.0))+1;
        }
        int bs = n*n+2*n;
        Vector<ComplexFloat> work(bs);

        int lrwork = 3*n*n+(4+2*k)*n+1;
        Vector<float> rwork(lrwork);

        int iwork_size = 5*n+3;
        Vector<int> iwork(iwork_size);

        CHEEVD("V", "U", &n, reinterpret_cast<MKL_Complex8*>(EV.Data()), &n, w.Data(), reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), &lrwork, iwork.Data(), &iwork_size, &info);
        
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of CHEEVD is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "The algorithm failed to converge; %i elements of an intermediate tridiagonal form did not converge to zero!", info);
                        Utility::RunTimeError(s);
                }
        }
        return w;
}

Vector<double> Lapack::Heevd(Matrix<ComplexDouble>& A)
{
        if(A.Rows() != A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        
        int info;
        Matrix<ComplexDouble> F = A.Clone();
        int n = F.Rows();
        Vector<double> w(n);

        int bs = n+1;
        Vector<ComplexDouble> work(bs);

        int lrwork = n;
        Vector<double> rwork(lrwork);

        int iwork_size = 1;
        Vector<int> iwork(iwork_size);

        ZHEEVD("N", "U", &n, reinterpret_cast<MKL_Complex16*>(F.Data()), &n, w.Data(), reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), &lrwork, iwork.Data(), &iwork_size, &info);

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of ZHEEVD is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "The algorithm failed to converge; %i elements of an intermediate tridiagonal form did not converge to zero!", info);
                        Utility::RunTimeError(s);
                }
        }
        return w;
}

Vector<double> Lapack::Heevd(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& EV)
{
        if(A.Rows() != A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        
        int info;
        EV = A.Clone();
        int n = EV.Rows();
        Vector<double> w(n);

        int k = 1;
        if(n>1)
        {
                k = (int)(log((double)(n-1))/log(2.0))+1;
        }
        int bs = n*n+2*n;
        Vector<ComplexDouble> work(bs);

        int lrwork = 3*n*n+(4+2*k)*n+1;
        Vector<double> rwork(lrwork);

        int iwork_size = 5*n+3;
        Vector<int> iwork(iwork_size);

        ZHEEVD("V", "U", &n, reinterpret_cast<MKL_Complex16*>(EV.Data()), &n, w.Data(), reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), &lrwork, iwork.Data(), &iwork_size, &info);
        
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of ZHEEVD is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "The algorithm failed to converge; %i elements of an intermediate tridiagonal form did not converge to zero!", info);
                        Utility::RunTimeError(s);
                }
        }
        return w;
}










// Syevr
Vector<float> Lapack::Syevr(Matrix<float>& A)
{
        if(A.Rows() != A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        
        int info;
        float abstol = SLAMCH("S");

        int m;
        float vl, vu;
        int il, iu;
        
        Matrix<float> F = A.Clone();
        int n = F.Rows();
        Vector<float> w(n);

        Vector<int> isuppz(2*n);

        int bs = (64+6)*n;
        Vector<float> work(bs);

        int iwork_size = 10*n;
        Vector<int> iwork(iwork_size);

        SSYEVR("N", "A", "U", &n, F.Data(), &n, &vl, &vu, &il, &iu, &abstol, &m,  w.Data(), 0, &n, isuppz.Data(), work.Data(), &bs, iwork.Data(), &iwork_size, &info);
        
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of SSYEVR is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "Internal error occured when running SSYEVR");
                        Utility::RunTimeError(s);
                }
        }
        return w;
}


Vector<float> Lapack::Syevr(Matrix<float>& A, int il, int iu)
{
        if(A.Rows() != A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }

        if(il > iu || il > A.Rows() || iu > A.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("The choice of parameters (il and iu) for selecting a subset of eigenvalues is not correct!");
        }
        
        int info;
        float abstol = SLAMCH("S");

        int m;
        float vl, vu;
        
        Matrix<float> F = A.Clone();
        int n = F.Rows();
        Vector<float> w(iu-il+1);

        Vector<int> isuppz(2*n);

        int bs = (64+6)*n;
        Vector<float> work(bs);

        int iwork_size = 10*n;
        Vector<int> iwork(iwork_size);

        SSYEVR("N", "I", "U", &n, F.Data(), &n, &vl, &vu, &il, &iu, &abstol, &m,  w.Data(), 0, &n, isuppz.Data(), work.Data(), &bs, iwork.Data(), &iwork_size, &info);
        
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of SSYEVR is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "Internal error occured when running SSYEVR");
                        Utility::RunTimeError(s);
                }
        }

        return w;
}

Vector<float> Lapack::Syevr(Matrix<float>& A, Matrix<float>& EV)
{
        if(A.Rows() != A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        
        int info;
        float abstol = SLAMCH("S");

        int m;
        float vl, vu;
        int il, iu;
        
        Matrix<float> F = A.Clone();
        int n = F.Rows();
        Vector<float> w(n);
        EV = Matrix<float>(n,n);

        Vector<int> isuppz(2*n);

        int bs = (64+6)*n;
        Vector<float> work(bs);

        int iwork_size = 10*n;
        Vector<int> iwork(iwork_size);

        SSYEVR("V", "A", "U", &n, F.Data(), &n, &vl, &vu, &il, &iu, &abstol, &m,  w.Data(), EV.Data(), &n, isuppz.Data(), work.Data(), &bs, iwork.Data(), &iwork_size, &info);
        
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of SSYEVR is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "Internal error occured when running SSYEVR");
                        Utility::RunTimeError(s);
                }
        }
        return w;
}

Vector<float> Lapack::Syevr(Matrix<float>& A, Matrix<float>& EV, int il, int iu)
{
        if(A.Rows() != A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }

        if(il > iu || il > A.Rows() || iu > A.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("The choice of parameters (il and iu) for selecting a subset of eigenvalues is not correct!");
        }
        
        int info;
        float abstol = SLAMCH("S");

        int m;
        float vl, vu;
        
        Matrix<float> F = A.Clone();
        int n = F.Rows();
        Vector<float> w(iu-il+1);
        EV = Matrix<float>(n,iu-il+1);

        Vector<int> isuppz(2*n);

        int bs = (64+6)*n;
        Vector<float> work(bs);

        int iwork_size = 10*n;
        Vector<int> iwork(iwork_size);

        SSYEVR("V", "I", "U", &n, F.Data(), &n, &vl, &vu, &il, &iu, &abstol, &m,  w.Data(), EV.Data(), &n, isuppz.Data(), work.Data(), &bs, iwork.Data(), &iwork_size, &info);
        
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of SSYEVR is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "Internal error occured when running SSYEVR");
                        Utility::RunTimeError(s);
                }
        }

        return w;
}






Vector<double> Lapack::Syevr(Matrix<double>& A)
{
        if(A.Rows() != A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        
        int info;
        double abstol = SLAMCH("S");

        int m;
        double vl, vu;
        int il, iu;
        
        Matrix<double> F = A.Clone();
        int n = F.Rows();
        Vector<double> w(n);

        Vector<int> isuppz(2*n);

        int bs = (64+6)*n;
        Vector<double> work(bs);

        int iwork_size = 10*n;
        Vector<int> iwork(iwork_size);

        DSYEVR("N", "A", "U", &n, F.Data(), &n, &vl, &vu, &il, &iu, &abstol, &m,  w.Data(), 0, &n, isuppz.Data(), work.Data(), &bs, iwork.Data(), &iwork_size, &info);
        
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of DSYEVR is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "Internal error occured when running DSYEVR");
                        Utility::RunTimeError(s);
                }
        }
        return w;
}


Vector<double> Lapack::Syevr(Matrix<double>& A, int il, int iu)
{
        if(A.Rows() != A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }

        if(il > iu || il > A.Rows() || iu > A.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("The choice of parameters (il and iu) for selecting a subset of eigenvalues is not correct!");
        }
        
        int info;
        double abstol = SLAMCH("S");

        int m;
        double vl, vu;
        
        Matrix<double> F = A.Clone();
        int n = F.Rows();
        Vector<double> w(iu-il+1);

        Vector<int> isuppz(2*n);

        int bs = (64+6)*n;
        Vector<double> work(bs);

        int iwork_size = 10*n;
        Vector<int> iwork(iwork_size);

        DSYEVR("N", "I", "U", &n, F.Data(), &n, &vl, &vu, &il, &iu, &abstol, &m,  w.Data(), 0, &n, isuppz.Data(), work.Data(), &bs, iwork.Data(), &iwork_size, &info);
        
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of DSYEVR is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "Internal error occured when running DSYEVR");
                        Utility::RunTimeError(s);
                }
        }

        return w;
}

Vector<double> Lapack::Syevr(Matrix<double>& A, Matrix<double>& EV)
{
        if(A.Rows() != A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        
        int info;
        double abstol = SLAMCH("S");

        int m;
        double vl, vu;
        int il, iu;
        
        Matrix<double> F = A.Clone();
        int n = F.Rows();
        Vector<double> w(n);
        EV = Matrix<double>(n,n);

        Vector<int> isuppz(2*n);

        int bs = (64+6)*n;
        Vector<double> work(bs);

        int iwork_size = 10*n;
        Vector<int> iwork(iwork_size);

        DSYEVR("V", "A", "U", &n, F.Data(), &n, &vl, &vu, &il, &iu, &abstol, &m,  w.Data(), EV.Data(), &n, isuppz.Data(), work.Data(), &bs, iwork.Data(), &iwork_size, &info);
        
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of DSYEVR is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "Internal error occured when running DSYEVR");
                        Utility::RunTimeError(s);
                }
        }
        return w;
}

Vector<double> Lapack::Syevr(Matrix<double>& A, Matrix<double>& EV, int il, int iu)
{
        if(A.Rows() != A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }

        if(il > iu || il > A.Rows() || iu > A.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("The choice of parameters (il and iu) for selecting a subset of eigenvalues is not correct!");
        }
        
        int info;
        double abstol = SLAMCH("S");

        int m;
        double vl, vu;
        
        Matrix<double> F = A.Clone();
        int n = F.Rows();
        Vector<double> w(iu-il+1);
        EV = Matrix<double>(n,iu-il+1);

        Vector<int> isuppz(2*n);

        int bs = (64+6)*n;
        Vector<double> work(bs);

        int iwork_size = 10*n;
        Vector<int> iwork(iwork_size);

        DSYEVR("V", "I", "U", &n, F.Data(), &n, &vl, &vu, &il, &iu, &abstol, &m,  w.Data(), EV.Data(), &n, isuppz.Data(), work.Data(), &bs, iwork.Data(), &iwork_size, &info);
        
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of DSYEVR is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "Internal error occured when running DSYEVR");
                        Utility::RunTimeError(s);
                }
        }

        return w;
}





Vector<float> Lapack::Heevr(Matrix<ComplexFloat>& A)
{
        if(A.Rows() != A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        
        int info;
        float abstol = SLAMCH("S");

        int m;
        float vl, vu;
        int il, iu;
        
        Matrix<ComplexFloat> F = A.Clone();
        int n = F.Rows();
        Vector<float> w(n);

        Vector<int> isuppz(2*n);

        int bs = (64+1)*n;
        Vector<ComplexFloat> work(bs);

        int rwork_size = 24*n;
        Vector<float> rwork(rwork_size);

        int iwork_size = 10*n;
        Vector<int> iwork(iwork_size);

        CHEEVR("N", "A", "U", &n, reinterpret_cast<MKL_Complex8*>(F.Data()), &n, &vl, &vu, &il, &iu, &abstol, &m,  w.Data(), 0, &n, isuppz.Data(), reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), &rwork_size, iwork.Data(), &iwork_size, &info);
        
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of CHEEVR is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "Internal error occured when running CHEEVR");
                        Utility::RunTimeError(s);
                }
        }
        return w;
}


Vector<float> Lapack::Heevr(Matrix<ComplexFloat>& A, int il, int iu)
{
        if(A.Rows() != A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }

        if(il > iu || il > A.Rows() || iu > A.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("The choice of parameters (il and iu) for selecting a subset of eigenvalues is not correct!");
        }
        
        int info;
        float abstol = SLAMCH("S");

        int m;
        float vl, vu;
        
        Matrix<ComplexFloat> F = A.Clone();
        int n = F.Rows();
        Vector<float> w(iu-il+1);

        Vector<int> isuppz(2*n);

        int bs = (64+6)*n;
        Vector<ComplexFloat> work(bs);

        int rwork_size = 24*n;
        Vector<float> rwork(rwork_size);
        int iwork_size = 10*n;
        Vector<int> iwork(iwork_size);

        CHEEVR("N", "I", "U", &n, reinterpret_cast<MKL_Complex8*>(F.Data()), &n, &vl, &vu, &il, &iu, &abstol, &m,  w.Data(), 0, &n, isuppz.Data(), reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), &rwork_size, iwork.Data(), &iwork_size, &info);
        
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of CHEEVR is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "Internal error occured when running CHEEVR");
                        Utility::RunTimeError(s);
                }
        }

        return w;
}

Vector<float> Lapack::Heevr(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& EV)
{
        if(A.Rows() != A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        
        int info;
        float abstol = SLAMCH("S");

        int m;
        float vl, vu;
        int il, iu;
        
        Matrix<ComplexFloat> F = A.Clone();
        int n = F.Rows();
        Vector<float> w(n);
        EV = Matrix<ComplexFloat>(n,n);

        Vector<int> isuppz(2*n);

        int bs = (64+6)*n;
        Vector<ComplexFloat> work(bs);

        int rwork_size = 24*n;
        Vector<float> rwork(rwork_size);
        int iwork_size = 10*n;
        Vector<int> iwork(iwork_size);

        CHEEVR("V", "A", "U", &n, reinterpret_cast<MKL_Complex8*>(F.Data()), &n, &vl, &vu, &il, &iu, &abstol, &m,  w.Data(), reinterpret_cast<MKL_Complex8*>(EV.Data()), &n, isuppz.Data(), reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), &rwork_size, iwork.Data(), &iwork_size, &info);
        
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of CHEEVR is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "Internal error occured when running CHEEVR");
                        Utility::RunTimeError(s);
                }
        }
        return w;
}

Vector<float> Lapack::Heevr(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& EV, int il, int iu)
{
        if(A.Rows() != A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }

        if(il > iu || il > A.Rows() || iu > A.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("The choice of parameters (il and iu) for selecting a subset of eigenvalues is not correct!");
        }
        
        int info;
        float abstol = SLAMCH("S");

        int m;
        float vl, vu;
        
        Matrix<ComplexFloat> F = A.Clone();
        int n = F.Rows();
        Vector<float> w(iu-il+1);
        EV = Matrix<ComplexFloat>(n,iu-il+1);

        Vector<int> isuppz(2*n);

        int bs = (64+6)*n;
        Vector<ComplexFloat> work(bs);

        int rwork_size = 24*n;
        Vector<float> rwork(rwork_size);
        int iwork_size = 10*n;
        Vector<int> iwork(iwork_size);

        CHEEVR("V", "I", "U", &n, reinterpret_cast<MKL_Complex8*>(F.Data()), &n, &vl, &vu, &il, &iu, &abstol, &m,  w.Data(), reinterpret_cast<MKL_Complex8*>(EV.Data()), &n, isuppz.Data(), reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), &rwork_size, iwork.Data(), &iwork_size, &info);
        
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of CHEEVR is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "Internal error occured when running CHEEVR");
                        Utility::RunTimeError(s);
                }
        }

        return w;
}




Vector<double> Lapack::Heevr(Matrix<ComplexDouble>& A)
{
        if(A.Rows() != A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        
        int info;
        double abstol = SLAMCH("S");

        int m;
        double vl, vu;
        int il, iu;
        
        Matrix<ComplexDouble> F = A.Clone();
        int n = F.Rows();
        Vector<double> w(n);

        Vector<int> isuppz(2*n);

        int bs = (64+1)*n;
        Vector<ComplexDouble> work(bs);

        int rwork_size = 24*n;
        Vector<double> rwork(rwork_size);

        int iwork_size = 10*n;
        Vector<int> iwork(iwork_size);

        ZHEEVR("N", "A", "U", &n, reinterpret_cast<MKL_Complex16*>(F.Data()), &n, &vl, &vu, &il, &iu, &abstol, &m,  w.Data(), 0, &n, isuppz.Data(), reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), &rwork_size, iwork.Data(), &iwork_size, &info);
        
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of ZHEEVR is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "Internal error occured when running ZHEEVR");
                        Utility::RunTimeError(s);
                }
        }
        return w;
}


Vector<double> Lapack::Heevr(Matrix<ComplexDouble>& A, int il, int iu)
{
        if(A.Rows() != A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }

        if(il > iu || il > A.Rows() || iu > A.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("The choice of parameters (il and iu) for selecting a subset of eigenvalues is not correct!");
        }
        
        int info;
        double abstol = SLAMCH("S");

        int m;
        double vl, vu;
        
        Matrix<ComplexDouble> F = A.Clone();
        int n = F.Rows();
        Vector<double> w(iu-il+1);

        Vector<int> isuppz(2*n);

        int bs = (64+6)*n;
        Vector<ComplexDouble> work(bs);

        int rwork_size = 24*n;
        Vector<double> rwork(rwork_size);
        int iwork_size = 10*n;
        Vector<int> iwork(iwork_size);

        ZHEEVR("N", "I", "U", &n, reinterpret_cast<MKL_Complex16*>(F.Data()), &n, &vl, &vu, &il, &iu, &abstol, &m,  w.Data(), 0, &n, isuppz.Data(), reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), &rwork_size, iwork.Data(), &iwork_size, &info);
        
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of ZHEEVR is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "Internal error occured when running ZHEEVR");
                        Utility::RunTimeError(s);
                }
        }

        return w;
}

Vector<double> Lapack::Heevr(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& EV)
{
        if(A.Rows() != A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        
        int info;
        double abstol = SLAMCH("S");

        int m;
        double vl, vu;
        int il, iu;
        
        Matrix<ComplexDouble> F = A.Clone();
        int n = F.Rows();
        Vector<double> w(n);
        EV = Matrix<ComplexDouble>(n,n);

        Vector<int> isuppz(2*n);

        int bs = (64+6)*n;
        Vector<ComplexDouble> work(bs);

        int rwork_size = 24*n;
        Vector<double> rwork(rwork_size);
        int iwork_size = 10*n;
        Vector<int> iwork(iwork_size);

        ZHEEVR("V", "A", "U", &n, reinterpret_cast<MKL_Complex16*>(F.Data()), &n, &vl, &vu, &il, &iu, &abstol, &m,  w.Data(), reinterpret_cast<MKL_Complex16*>(EV.Data()), &n, isuppz.Data(), reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), &rwork_size, iwork.Data(), &iwork_size, &info);
        
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of ZHEEVR is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "Internal error occured when running ZHEEVR");
                        Utility::RunTimeError(s);
                }
        }
        return w;
}

Vector<double> Lapack::Heevr(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& EV, int il, int iu)
{
        if(A.Rows() != A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }

        if(il > iu || il > A.Rows() || iu > A.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("The choice of parameters (il and iu) for selecting a subset of eigenvalues is not correct!");
        }
        
        int info;
        double abstol = SLAMCH("S");

        int m;
        double vl, vu;
        
        Matrix<ComplexDouble> F = A.Clone();
        int n = F.Rows();
        Vector<double> w(iu-il+1);
        EV = Matrix<ComplexDouble>(n,iu-il+1);

        Vector<int> isuppz(2*n);

        int bs = (64+6)*n;
        Vector<ComplexDouble> work(bs);

        int rwork_size = 24*n;
        Vector<double> rwork(rwork_size);
        int iwork_size = 10*n;
        Vector<int> iwork(iwork_size);

        ZHEEVR("V", "I", "U", &n, reinterpret_cast<MKL_Complex16*>(F.Data()), &n, &vl, &vu, &il, &iu, &abstol, &m,  w.Data(), reinterpret_cast<MKL_Complex16*>(EV.Data()), &n, isuppz.Data(), reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), &rwork_size, iwork.Data(), &iwork_size, &info);
        
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of ZHEEVR is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "Internal error occured when running ZHEEVR");
                        Utility::RunTimeError(s);
                }
        }

        return w;
}











//=========================================================
// Non-Symmetric eigenvalue problem
//=========================================================


// Gees
Vector<ComplexFloat> Lapack::Gees(Matrix<float>& A, Matrix<float>& T)
{
        if(A.Rows() != A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        
        int info;
        T = A.Clone();
        int n = T.Rows();
        Vector<float> wr(n);
        Vector<float> wi(n);

        int bs = 3*n;
        Vector<float> work(bs);
        int sdim;
        
        SGEES("N", "N", 0, &n, T.Data(), &n, &sdim, wr.Data(), wi.Data(), 0, &n, work.Data(), &bs, 0, &info);

        Vector<ComplexFloat> w(n);
        for(int i=0;i<n;i++)
        {
                w[i] = ComplexFloat(wr[i],wi[i]);
        }

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of SGEES is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "The QR algorithm failed to compute all the eigenvalues!");
                        Utility::RunTimeError(s);
                }
        }
        return w;
}

Vector<ComplexFloat> Lapack::Gees(Matrix<float>& A, Matrix<float>& T, Matrix<float>& Z)
{
        if(A.Rows() != A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        
        int info;
        T = A.Clone();
        int n = T.Rows();
        Vector<float> wr(n);
        Vector<float> wi(n);

        Z = Matrix<float>(n,n);

        int bs = 3*n;
        Vector<float> work(bs);
        int sdim;

        SGEES("V", "N", 0, &n, T.Data(), &n, &sdim, wr.Data(), wi.Data(), Z.Data(), &n, work.Data(), &bs, 0, &info);

        Vector<ComplexFloat> w(n);
        for(int i=0;i<n;i++)
        {
                w[i] = ComplexFloat(wr[i],wi[i]);
        }

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of SGEES is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "The QR algorithm failed to compute all the eigenvalues!");
                        Utility::RunTimeError(s);
                }
        }
        return w;
}

Vector<ComplexDouble> Lapack::Gees(Matrix<double>& A, Matrix<double>& T)
{
        if(A.Rows() != A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        
        int info;
        T = A.Clone();
        int n = T.Rows();
        Vector<double> wr(n);
        Vector<double> wi(n);

        int bs = 3*n;
        Vector<double> work(bs);
        int sdim;
        
        DGEES("N", "N", 0, &n, T.Data(), &n, &sdim, wr.Data(), wi.Data(), 0, &n, work.Data(), &bs, 0, &info);

        Vector<ComplexDouble> w(n);
        for(int i=0;i<n;i++)
        {
                w[i] = ComplexDouble(wr[i],wi[i]);
        }

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of DGEES is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "The QR algorithm failed to compute all the eigenvalues!");
                        Utility::RunTimeError(s);
                }
        }
        return w;
}

Vector<ComplexDouble> Lapack::Gees(Matrix<double>& A, Matrix<double>& T, Matrix<double>& Z)
{
        if(A.Rows() != A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        
        int info;
        T = A.Clone();
        int n = T.Rows();
        Vector<double> wr(n);
        Vector<double> wi(n);

        Z = Matrix<double>(n,n);

        int bs = 3*n;
        Vector<double> work(bs);
        int sdim;

        DGEES("V", "N", 0, &n, T.Data(), &n, &sdim, wr.Data(), wi.Data(), Z.Data(), &n, work.Data(), &bs, 0, &info);

        Vector<ComplexDouble> w(n);
        for(int i=0;i<n;i++)
        {
                w[i] = ComplexDouble(wr[i],wi[i]);
        }

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of DGEES is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "The QR algorithm failed to compute all the eigenvalues!");
                        Utility::RunTimeError(s);
                }
        }
        return w;
}



Vector<ComplexFloat> Lapack::Gees(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& T)
{
        if(A.Rows() != A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        
        int info;
        T = A.Clone();
        int n = T.Rows();
        Vector<ComplexFloat> w(n);

        int bs = 2*n;
        Vector<ComplexFloat> work(bs);
        Vector<float> rwork(n);
        int sdim;

        CGEES("N", "N", 0, &n, reinterpret_cast<MKL_Complex8*>(T.Data()), &n, &sdim, reinterpret_cast<MKL_Complex8*>(w.Data()), 0, &n, reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), 0, &info);

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of CGEES is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "The QR algorithm failed to compute all the eigenvalues!");
                        Utility::RunTimeError(s);
                }
        }
        return w;
}


Vector<ComplexFloat> Lapack::Gees(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& T, Matrix<ComplexFloat>& Z)
{
        if(A.Rows() != A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        
        int info;
        T = A.Clone();
        int n = T.Rows();
        Vector<ComplexFloat> w(n);
        Z = Matrix<ComplexFloat>(n,n);

        int bs = 2*n;
        Vector<ComplexFloat> work(bs);
        Vector<float> rwork(n);
        int sdim;

        CGEES("V", "N", 0, &n, reinterpret_cast<MKL_Complex8*>(T.Data()), &n, &sdim, reinterpret_cast<MKL_Complex8*>(w.Data()), reinterpret_cast<MKL_Complex8*>(Z.Data()), &n, reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), 0, &info);

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of CGEES is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "The QR algorithm failed to compute all the eigenvalues!");
                        Utility::RunTimeError(s);
                }
        }
        return w;
}



Vector<ComplexDouble> Lapack::Gees(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& T)
{
        if(A.Rows() != A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        
        int info;
        T = A.Clone();
        int n = T.Rows();
        Vector<ComplexDouble> w(n);

        int bs = 2*n;
        Vector<ComplexDouble> work(bs);
        Vector<double> rwork(n);
        int sdim;

        ZGEES("N", "N", 0, &n, reinterpret_cast<MKL_Complex16*>(T.Data()), &n, &sdim, reinterpret_cast<MKL_Complex16*>(w.Data()), 0, &n, reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), 0, &info);

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of ZGEES is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "The QR algorithm failed to compute all the eigenvalues!");
                        Utility::RunTimeError(s);
                }
        }
        return w;
}


Vector<ComplexDouble> Lapack::Gees(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& T, Matrix<ComplexDouble>& Z)
{
        if(A.Rows() != A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        
        int info;
        T = A.Clone();
        int n = T.Rows();
        Vector<ComplexDouble> w(n);
        Z = Matrix<ComplexDouble>(n,n);

        int bs = 2*n;
        Vector<ComplexDouble> work(bs);
        Vector<double> rwork(n);
        int sdim;

        ZGEES("V", "N", 0, &n, reinterpret_cast<MKL_Complex16*>(T.Data()), &n, &sdim, reinterpret_cast<MKL_Complex16*>(w.Data()), reinterpret_cast<MKL_Complex16*>(Z.Data()), &n, reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), 0, &info);

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of ZGEES is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "The QR algorithm failed to compute all the eigenvalues!");
                        Utility::RunTimeError(s);
                }
        }
        return w;
}






// Geev

Vector<ComplexFloat> Lapack::Geev(Matrix<float>& A)
{
        if(A.Rows() != A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        
        int info;
        Matrix<float> F = A.Clone();
        int n = F.Rows();
        Vector<float> wr(n);
        Vector<float> wi(n);

        int bs = 4*n;
        Vector<float> work(bs);

        SGEEV("N", "N", &n, F.Data(), &n, wr.Data(), wi.Data(), 0, &n, 0, &n, work.Data(), &bs, &info);

        Vector<ComplexFloat> w(n);
        for(int i=0;i<n;i++)
        {
                w[i] = ComplexFloat(wr[i],wi[i]);
        }

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of SGEEV is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "The QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed!");
                        Utility::RunTimeError(s);
                }
        }
        return w;
}

Vector<ComplexFloat> Lapack::Geev(Matrix<float>& A, Matrix<ComplexFloat>& Vl, Matrix<ComplexFloat>& Vr)
{
        if(A.Rows() != A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        
        int info;
        Matrix<float> F = A.Clone();
        int n = F.Rows();
        Vector<float> wr(n);
        Vector<float> wi(n);

        Vl = Matrix<ComplexFloat>(n,n);
        Vr = Matrix<ComplexFloat>(n,n);
        Matrix<float> vl(n,n);
        Matrix<float> vr(n,n);

        int bs = 4*n;
        Vector<float> work(bs);

        SGEEV("V", "V", &n, F.Data(), &n, wr.Data(), wi.Data(), vl.Data(), &n, vr.Data(), &n, work.Data(), &bs, &info);

        Vector<ComplexFloat> w(n);
        bool conj = false;
        for(int i=0;i<n;i++)
        {
                w[i] = ComplexFloat(wr[i],wi[i]);
                if(fabs(wi[i]) < numeric_limits<float>::epsilon())
                {
                        for(int j=0;j<n;j++)
                        {
                                Vl.ElemNC(j,i) = ComplexFloat(vl.ElemNC(j,i),0);
                                Vr.ElemNC(j,i) = ComplexFloat(vr.ElemNC(j,i),0);
                        }
                }
                else
                {
                        for(int j=0;j<n;j++)
                        {
                                if(conj)
                                {
                                        Vl.ElemNC(j,i) = ComplexFloat(vl.ElemNC(j,i-1),-vl.ElemNC(j,i));
                                        Vr.ElemNC(j,i) = ComplexFloat(vr.ElemNC(j,i-1),-vr.ElemNC(j,i));
                                }
                                else
                                {
                                        Vl.ElemNC(j,i) = ComplexFloat(vl.ElemNC(j,i),vl.ElemNC(j,i+1));
                                        Vr.ElemNC(j,i) = ComplexFloat(vr.ElemNC(j,i),vr.ElemNC(j,i+1));
                                }
                        }
                        conj = !conj;
                }
        }

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of SGEEV is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "The QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed!");
                        Utility::RunTimeError(s);
                }
        }
        return w;
}

Vector<ComplexDouble> Lapack::Geev(Matrix<double>& A)
{
        if(A.Rows() != A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        
        int info;
        Matrix<double> F = A.Clone();
        int n = F.Rows();
        Vector<double> wr(n);
        Vector<double> wi(n);

        int bs = 4*n;
        Vector<double> work(bs);

        DGEEV("N", "N", &n, F.Data(), &n, wr.Data(), wi.Data(), 0, &n, 0, &n, work.Data(), &bs, &info);

        Vector<ComplexDouble> w(n);
        for(int i=0;i<n;i++)
        {
                w[i] = ComplexDouble(wr[i],wi[i]);
        }

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of DGEEV is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "The QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed!");
                        Utility::RunTimeError(s);
                }
        }
        return w;
}

Vector<ComplexDouble> Lapack::Geev(Matrix<double>& A, Matrix<ComplexDouble>& Vl, Matrix<ComplexDouble>& Vr)
{
        if(A.Rows() != A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        
        int info;
        Matrix<double> F = A.Clone();
        int n = F.Rows();
        Vector<double> wr(n);
        Vector<double> wi(n);

        Vl = Matrix<ComplexDouble>(n,n);
        Vr = Matrix<ComplexDouble>(n,n);
        Matrix<double> vl(n,n);
        Matrix<double> vr(n,n);

        int bs = 4*n;
        Vector<double> work(bs);

        DGEEV("V", "V", &n, F.Data(), &n, wr.Data(), wi.Data(), vl.Data(), &n, vr.Data(), &n, work.Data(), &bs, &info);

        Vector<ComplexDouble> w(n);
        bool conj = false;
        for(int i=0;i<n;i++)
        {
                w[i] = ComplexDouble(wr[i],wi[i]);
                if(fabs(wi[i]) < numeric_limits<double>::epsilon())
                {
                        for(int j=0;j<n;j++)
                        {
                                Vl.ElemNC(j,i) = ComplexDouble(vl.ElemNC(j,i),0);
                                Vr.ElemNC(j,i) = ComplexDouble(vr.ElemNC(j,i),0);
                        }
                }
                else
                {
                        for(int j=0;j<n;j++)
                        {
                                if(conj)
                                {
                                        Vl.ElemNC(j,i) = ComplexDouble(vl.ElemNC(j,i-1),-vl.ElemNC(j,i));
                                        Vr.ElemNC(j,i) = ComplexDouble(vr.ElemNC(j,i-1),-vr.ElemNC(j,i));
                                }
                                else
                                {
                                        Vl.ElemNC(j,i) = ComplexDouble(vl.ElemNC(j,i),vl.ElemNC(j,i+1));
                                        Vr.ElemNC(j,i) = ComplexDouble(vr.ElemNC(j,i),vr.ElemNC(j,i+1));
                                }
                        }
                        conj = !conj;
                }
        }

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of DGEEV is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "The QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed!");
                        Utility::RunTimeError(s);
                }
        }
        return w;
}

Vector<ComplexFloat> Lapack::Geev(Matrix<ComplexFloat>& A)
{
        if(A.Rows() != A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        
        int info;
        Matrix<ComplexFloat> F = A.Clone();
        int n = F.Rows();
        Vector<ComplexFloat> w(n);

        int bs = 4*n;
        Vector<ComplexFloat> work(bs);
        
        Vector<float> rwork(2*n);

        CGEEV("N", "N", &n, reinterpret_cast<MKL_Complex8*>(F.Data()), &n, reinterpret_cast<MKL_Complex8*>(w.Data()), 0, &n, 0, &n, reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), &info);

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of CGEEV is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "The QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed!");
                        Utility::RunTimeError(s);
                }
        }
        return w;
}

Vector<ComplexFloat> Lapack::Geev(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& Vl, Matrix<ComplexFloat>& Vr)
{
        if(A.Rows() != A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        
        int info;
        Matrix<ComplexFloat> F = A.Clone();
        int n = F.Rows();
        Vector<ComplexFloat> w(n);

        Vl = Matrix<ComplexFloat>(n,n);
        Vr = Matrix<ComplexFloat>(n,n);

        int bs = 4*n;
        Vector<ComplexFloat> work(bs);

        Vector<float> rwork(2*n);

        CGEEV("V", "V", &n, reinterpret_cast<MKL_Complex8*>(F.Data()), &n, reinterpret_cast<MKL_Complex8*>(w.Data()), reinterpret_cast<MKL_Complex8*>(Vl.Data()), &n, reinterpret_cast<MKL_Complex8*>(Vr.Data()), &n, reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), &info);

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of CGEEV is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "The QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed!");
                        Utility::RunTimeError(s);
                }
        }
        return w;
}

Vector<ComplexDouble> Lapack::Geev(Matrix<ComplexDouble>& A)
{
        if(A.Rows() != A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        
        int info;
        Matrix<ComplexDouble> F = A.Clone();
        int n = F.Rows();
        Vector<ComplexDouble> w(n);

        int bs = 4*n;
        Vector<ComplexDouble> work(bs);
        
        Vector<double> rwork(2*n);

        ZGEEV("N", "N", &n, reinterpret_cast<MKL_Complex16*>(F.Data()), &n, reinterpret_cast<MKL_Complex16*>(w.Data()), 0, &n, 0, &n, reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), &info);

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of ZGEEV is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "The QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed!");
                        Utility::RunTimeError(s);
                }
        }
        return w;
}

Vector<ComplexDouble> Lapack::Geev(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& Vl, Matrix<ComplexDouble>& Vr)
{
        if(A.Rows() != A.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        
        int info;
        Matrix<ComplexDouble> F = A.Clone();
        int n = F.Rows();
        Vector<ComplexDouble> w(n);

        Vl = Matrix<ComplexDouble>(n,n);
        Vr = Matrix<ComplexDouble>(n,n);

        int bs = 4*n;
        Vector<ComplexDouble> work(bs);

        Vector<double> rwork(2*n);

        ZGEEV("V", "V", &n, reinterpret_cast<MKL_Complex16*>(F.Data()), &n, reinterpret_cast<MKL_Complex16*>(w.Data()), reinterpret_cast<MKL_Complex16*>(Vl.Data()), &n, reinterpret_cast<MKL_Complex16*>(Vr.Data()), &n, reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), &info);

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of ZGEEV is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "The QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed!");
                        Utility::RunTimeError(s);
                }
        }
        return w;
}






//=========================================================
// Singular Value Decomposition (SVD)
//=========================================================




// Gesdd
Vector<float> Lapack::Gesdd(Matrix<float>& A)
{
        int info;
        Matrix<float> F = A.Clone();
        int m = F.Rows();
        int n = F.Columns();
        int min = (m>n)?n:m;
        int max = (m>n)?m:n;

        Vector<float> s(min);

        int temp = (max>6*min) ? max : 6*min;
        int bs = 3*min + temp;
        Vector<float> work(bs);
        Vector<int> iwork(8*min);
        SGESDD("N", &m, &n, F.Data(), &m, s.Data(), 0, &m, 0, &n, work.Data(), &bs, iwork.Data(), &info);
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of SGESDD is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "_bdsdc did not converge, updating process failed!");
                        Utility::RunTimeError(s);
                }
        }
        return s;
}

Vector<float> Lapack::Gesdd(Matrix<float>& A, Matrix<float>& U, Matrix<float>& V)
{
        int info;
        Matrix<float> F = A.Clone();
        int m = F.Rows();
        int n = F.Columns();
        int min = (m>n)?n:m;
        int max = (m>n)?m:n;
        U = Matrix<float>(m,m);
        V = Matrix<float>(n,n);

        Vector<float> s(min);

        int temp = (max > 4*min*min+4*min) ? max : 4*min*min+4*min;
        int bs = 3*min*min + temp;
        Vector<float> work(bs);
        Vector<int> iwork(8*min);
        SGESDD("A", &m, &n, F.Data(), &m, s.Data(), U.Data(), &m, V.Data(), &n, work.Data(), &bs, iwork.Data(), &info);
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of SGESDD is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "_bdsdc did not converge, updating process failed!");
                        Utility::RunTimeError(s);
                }
        }
        
        V.Transpose();
                
        return s;
}

Vector<double> Lapack::Gesdd(Matrix<double>& A)
{
        int info;
        Matrix<double> F = A.Clone();
        int m = F.Rows();
        int n = F.Columns();
        int min = (m>n)?n:m;
        int max = (m>n)?m:n;

        Vector<double> s(min);

        int temp = (max>6*min) ? max : 6*min;
        int bs = 3*min + temp;
        Vector<double> work(bs);
        Vector<int> iwork(8*min);
        DGESDD("N", &m, &n, F.Data(), &m, s.Data(), 0, &m, 0, &n, work.Data(), &bs, iwork.Data(), &info);
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of DGESDD is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "_bdsdc did not converge, updating process failed!");
                        Utility::RunTimeError(s);
                }
        }
        return s;
}

Vector<double> Lapack::Gesdd(Matrix<double>& A, Matrix<double>& U, Matrix<double>& V)
{
        int info;
        Matrix<double> F = A.Clone();
        int m = F.Rows();
        int n = F.Columns();
        int min = (m>n)?n:m;
        int max = (m>n)?m:n;
        U = Matrix<double>(m,m);
        V = Matrix<double>(n,n);

        Vector<double> s(min);

        int temp = (max > 4*min*min+4*min) ? max : 4*min*min+4*min;
        int bs = 3*min*min + temp;
        Vector<double> work(bs);
        Vector<int> iwork(8*min);
        DGESDD("A", &m, &n, F.Data(), &m, s.Data(), U.Data(), &m, V.Data(), &n, work.Data(), &bs, iwork.Data(), &info);
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of DGESDD is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "_bdsdc did not converge, updating process failed!");
                        Utility::RunTimeError(s);
                }
        }
        
        V.Transpose();
                
        return s;
}

Vector<float> Lapack::Gesdd(Matrix<ComplexFloat>& A)
{
        int info;
        Matrix<ComplexFloat> F = A.Clone();
        int m = F.Rows();
        int n = F.Columns();
        int min = (m>n)?n:m;
        int max = (m>n)?m:n;

        Vector<float> s(min);

        int bs = 2*min + max;
        Vector<ComplexFloat> work(bs);
        Vector<int> iwork(8*min);
        Vector<float> rwork(5*min);

        CGESDD("N", &m, &n, reinterpret_cast<MKL_Complex8*>(F.Data()), &m, s.Data(), 0, &m, 0, &n, reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), iwork.Data(), &info);
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of CGESDD is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "_bdsdc did not converge, updating process failed!");
                        Utility::RunTimeError(s);
                }
        }
        return s;
}

Vector<float> Lapack::Gesdd(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& U, Matrix<ComplexFloat>& V)
{
        int info;
        Matrix<ComplexFloat> F = A.Clone();
        int m = F.Rows();
        int n = F.Columns();
        int min = (m>n)?n:m;
        int max = (m>n)?m:n;
        U = Matrix<ComplexFloat>(m,m);
        V = Matrix<ComplexFloat>(n,n);

        Vector<float> s(min);

        int bs = min*min + max + 2*min;
        Vector<ComplexFloat> work(bs);
        Vector<int> iwork(8*min);
        Vector<float> rwork(5*min*min + 7*min);
        CGESDD("A", &m, &n, reinterpret_cast<MKL_Complex8*>(F.Data()), &m, s.Data(), reinterpret_cast<MKL_Complex8*>(U.Data()), &m, reinterpret_cast<MKL_Complex8*>(V.Data()), &n, reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), iwork.Data(), &info);
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of CGESDD is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "_bdsdc did not converge, updating process failed!");
                        Utility::RunTimeError(s);
                }
        }
        
        V.Transpose();
        for(int i=0; i<V.Length(); i++)
        {
                V.ElemNC(i) = conj(V.ElemNC(i));
        }
                
        return s;
}

Vector<double> Lapack::Gesdd(Matrix<ComplexDouble>& A)
{
        int info;
        Matrix<ComplexDouble> F = A.Clone();
        int m = F.Rows();
        int n = F.Columns();
        int min = (m>n)?n:m;
        int max = (m>n)?m:n;

        Vector<double> s(min);

        int bs = 2*min + max;
        Vector<ComplexDouble> work(bs);
        Vector<int> iwork(8*min);
        Vector<double> rwork(5*min);

        ZGESDD("N", &m, &n, reinterpret_cast<MKL_Complex16*>(F.Data()), &m, s.Data(), 0, &m, 0, &n, reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), iwork.Data(), &info);
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of ZGESDD is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "_bdsdc did not converge, updating process failed!");
                        Utility::RunTimeError(s);
                }
        }
        return s;
}

Vector<double> Lapack::Gesdd(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& U, Matrix<ComplexDouble>& V)
{
        int info;
        Matrix<ComplexDouble> F = A.Clone();
        int m = F.Rows();
        int n = F.Columns();
        int min = (m>n)?n:m;
        int max = (m>n)?m:n;
        U = Matrix<ComplexDouble>(m,m);
        V = Matrix<ComplexDouble>(n,n);

        Vector<double> s(min);

        int bs = min*min + max + 2*min;
        Vector<ComplexDouble> work(bs);
        Vector<int> iwork(8*min);
        Vector<double> rwork(5*min*min + 7*min);
        ZGESDD("A", &m, &n, reinterpret_cast<MKL_Complex16*>(F.Data()), &m, s.Data(), reinterpret_cast<MKL_Complex16*>(U.Data()), &m, reinterpret_cast<MKL_Complex16*>(V.Data()), &n, reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), iwork.Data(), &info);
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of ZGESDD is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "_bdsdc did not converge, updating process failed!");
                        Utility::RunTimeError(s);
                }
        }
        
        V.Transpose();
        for(int i=0; i<V.Length(); i++)
        {
                V.ElemNC(i) = conj(V.ElemNC(i));
        }

        return s;
}












//=========================================================
// Generalized Symmetric Definite eigenproblem
//=========================================================


        // Sygv
        
        
        
        
        // Sygvx
        
        
        
        
        
        // Sygvd,Hegvd
Vector<float> Lapack::Sygvd(Matrix<float>& A, Matrix<float>& B, int ptype)
{
        if(ptype < 1 || ptype > 3 )
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Incorrext problem type. ptype can be 1, 2, or 3!");
        }
        
        if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        if(A.Rows() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }
        
        int info;
        Matrix<float> A2 = A.Clone();
        Matrix<float> B2 = B.Clone();
        int n = A2.Rows();
        Vector<float> w(n);

        int bs = 2*n+1;
        Vector<float> work(bs);

        int iwork_size = 1;
        Vector<int> iwork(iwork_size);

        SSYGVD(&ptype, "N", "U", &n, A2.Data(), &n, B2.Data(), &n, w.Data(), work.Data(), &bs, iwork.Data(), &iwork_size, &info);
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of SSYGVD is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        cerr << "spotrf/dpotrf and ssyev/dsyev returned an error code: ";
                        if(info <= n)
                        {
                                cerr << "ssyev/dsyev failed to converge, and " << info <<  " off-diagonal elements of an intermediate tridiagonal did not converge to zero!";
                        }
                        else if(info > n && info <=  2*n)
                        {
                                cerr << "The leading minor of order " << info-n << " of B is not positive-definite. The factorization of B could not be completed and no eigenvalues or eigenvectors were computed!";
                        }
                        Utility::RunTimeError("");
                }
        }
        return w;
}

Vector<float> Lapack::Sygvd(Matrix<float>& A, Matrix<float>& B, Matrix<float>& EV, int ptype)
{
        if(ptype < 1 || ptype > 3 )
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Incorrext problem type. ptype can be 1, 2, or 3!");
        }
        
        if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        if(A.Rows() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }
        
        int info;
        EV = A.Clone();
        Matrix<float> B2 = B.Clone();
        int n = B2.Rows();
        Vector<float> w(n);

        int bs = 2*n*n + 6*n +1;
        Vector<float> work(bs);

        int iwork_size = 5*n + 3;
        Vector<int> iwork(iwork_size);

        SSYGVD(&ptype, "V", "U", &n, EV.Data(), &n, B2.Data(), &n, w.Data(), work.Data(), &bs, iwork.Data(), &iwork_size, &info);
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of SSYGVD is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        cerr << "spotrf/dpotrf and ssyev/dsyev returned an error code: ";
                        if(info <= n)
                        {
                                cerr << "ssyev/dsyev failed to converge, and " << info <<  " off-diagonal elements of an intermediate tridiagonal did not converge to zero!";
                        }
                        else if(info > n && info <=  2*n)
                        {
                                cerr << "The leading minor of order " << info-n << " of B is not positive-definite. The factorization of B could not be completed and no eigenvalues or eigenvectors were computed!";
                        }
                        Utility::RunTimeError("");
                }
        }
        return w;
}

Vector<double> Lapack::Sygvd(Matrix<double>& A, Matrix<double>& B, int ptype)
{
        if(ptype < 1 || ptype > 3 )
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Incorrext problem type. ptype can be 1, 2, or 3!");
        }
        
        if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        if(A.Rows() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }
        
        int info;
        Matrix<double> A2 = A.Clone();
        Matrix<double> B2 = B.Clone();
        int n = A2.Rows();
        Vector<double> w(n);

        int bs = 2*n+1;
        Vector<double> work(bs);

        int iwork_size = 1;
        Vector<int> iwork(iwork_size);

        DSYGVD(&ptype, "N", "U", &n, A2.Data(), &n, B2.Data(), &n, w.Data(), work.Data(), &bs, iwork.Data(), &iwork_size, &info);
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of DSYGVD is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        cerr << "spotrf/dpotrf and ssyev/dsyev returned an error code: ";
                        if(info <= n)
                        {
                                cerr << "ssyev/dsyev failed to converge, and " << info <<  " off-diagonal elements of an intermediate tridiagonal did not converge to zero!";
                        }
                        else if(info > n && info <=  2*n)
                        {
                                cerr << "The leading minor of order " << info-n << " of B is not positive-definite. The factorization of B could not be completed and no eigenvalues or eigenvectors were computed!";
                        }
                        Utility::RunTimeError("");
                }
        }
        return w;
}

Vector<double> Lapack::Sygvd(Matrix<double>& A, Matrix<double>& B, Matrix<double>& EV, int ptype)
{
        if(ptype < 1 || ptype > 3 )
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Incorrext problem type. ptype can be 1, 2, or 3!");
        }
        
        if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        if(A.Rows() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }
        
        int info;
        EV = A.Clone();
        Matrix<double> B2 = B.Clone();
        int n = B2.Rows();
        Vector<double> w(n);

        int bs = 2*n*n + 6*n +1;
        Vector<double> work(bs);

        int iwork_size = 5*n + 3;
        Vector<int> iwork(iwork_size);

        DSYGVD(&ptype, "V", "U", &n, EV.Data(), &n, B2.Data(), &n, w.Data(), work.Data(), &bs, iwork.Data(), &iwork_size, &info);
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of DSYGVD is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        cerr << "spotrf/dpotrf and ssyev/dsyev returned an error code: ";
                        if(info <= n)
                        {
                                cerr << "ssyev/dsyev failed to converge, and " << info <<  " off-diagonal elements of an intermediate tridiagonal did not converge to zero!";
                        }
                        else if(info > n && info <=  2*n)
                        {
                                cerr << "The leading minor of order " << info-n << " of B is not positive-definite. The factorization of B could not be completed and no eigenvalues or eigenvectors were computed!";
                        }
                        Utility::RunTimeError("");
                }
        }
        return w;
}


// Hegvd
Vector<float> Lapack::Hegvd(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B, int ptype)
{
        if(ptype < 1 || ptype > 3 )
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Incorrext problem type. ptype can be 1, 2, or 3!");
        }
        
        if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        if(A.Rows() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }
        
        int info;
        Matrix<ComplexFloat> A2 = A.Clone();
        Matrix<ComplexFloat> B2 = B.Clone();
        int n = A2.Rows();
        Vector<float> w(n);

        int bs = n+1;
        Vector<ComplexFloat> work(bs);

        int rwork_size = n;
        Vector<float> rwork(rwork_size);

        int iwork_size = 1;
        Vector<int> iwork(iwork_size);

        CHEGVD(&ptype, "N", "U", &n, reinterpret_cast<MKL_Complex8*>(A2.Data()), &n, reinterpret_cast<MKL_Complex8*>(B2.Data()), &n, w.Data(), reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), &rwork_size, iwork.Data(), &iwork_size, &info);
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of CHEGVD is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        cerr << "cpotrf/zpotrf and cheev/zheev returned an error code: ";
                        if(info <= n)
                        {
                                cerr << "cheev/zheev failed to converge, and " << info <<  " off-diagonal elements of an intermediate tridiagonal did not converge to zero!";
                        }
                        else if(info > n && info <=  2*n)
                        {
                                cerr << "The leading minor of order " << info-n << " of B is not positive-definite. The factorization of B could not be completed and no eigenvalues or eigenvectors were computed!";
                        }
                        Utility::RunTimeError("");
                }
        }
        return w;
}

Vector<float> Lapack::Hegvd(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B, Matrix<ComplexFloat>& EV, int ptype)
{
        if(ptype < 1 || ptype > 3 )
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Incorrext problem type. ptype can be 1, 2, or 3!");
        }
        
        if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        if(A.Rows() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }
        
        int info;
        EV = A.Clone();
        Matrix<ComplexFloat> B2 = B.Clone();
        int n = B2.Rows();
        Vector<float> w(n);

        int bs = n*n+2*n;
        Vector<ComplexFloat> work(bs);

        int rwork_size = 2*n*n+5*n+1;
        Vector<float> rwork(rwork_size);

        int iwork_size = 5*n+3;
        Vector<int> iwork(iwork_size);

        CHEGVD(&ptype, "V", "U", &n, reinterpret_cast<MKL_Complex8*>(EV.Data()), &n, reinterpret_cast<MKL_Complex8*>(B2.Data()), &n, w.Data(), reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), &rwork_size, iwork.Data(), &iwork_size, &info);
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of CHEGVD is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        cerr << "cpotrf/zpotrf and cheev/zheev returned an error code: ";
                        if(info <= n)
                        {
                                cerr << "cheev/zheev failed to converge, and " << info <<  " off-diagonal elements of an intermediate tridiagonal did not converge to zero!";
                        }
                        else if(info > n && info <=  2*n)
                        {
                                cerr << "The leading minor of order " << info-n << " of B is not positive-definite. The factorization of B could not be completed and no eigenvalues or eigenvectors were computed!";
                        }
                        Utility::RunTimeError("");
                }
        }
        return w;
}


Vector<double> Lapack::Hegvd(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B, int ptype)
{
        if(ptype < 1 || ptype > 3 )
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Incorrext problem type. ptype can be 1, 2, or 3!");
        }
        
        if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        if(A.Rows() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }
        
        int info;
        Matrix<ComplexDouble> A2 = A.Clone();
        Matrix<ComplexDouble> B2 = B.Clone();
        int n = A2.Rows();
        Vector<double> w(n);

        int bs = n+1;
        Vector<ComplexDouble> work(bs);

        int rwork_size = n;
        Vector<double> rwork(rwork_size);

        int iwork_size = 1;
        Vector<int> iwork(iwork_size);

        ZHEGVD(&ptype, "N", "U", &n, reinterpret_cast<MKL_Complex16*>(A2.Data()), &n, reinterpret_cast<MKL_Complex16*>(B2.Data()), &n, w.Data(), reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), &rwork_size, iwork.Data(), &iwork_size, &info);
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of ZHEGVD is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        cerr << "cpotrf/zpotrf and cheev/zheev returned an error code: ";
                        if(info <= n)
                        {
                                cerr << "cheev/zheev failed to converge, and " << info <<  " off-diagonal elements of an intermediate tridiagonal did not converge to zero!";
                        }
                        else if(info > n && info <=  2*n)
                        {
                                cerr << "The leading minor of order " << info-n << " of B is not positive-definite. The factorization of B could not be completed and no eigenvalues or eigenvectors were computed!";
                        }
                        Utility::RunTimeError("");
                }
        }
        return w;
}

Vector<double> Lapack::Hegvd(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B, Matrix<ComplexDouble>& EV, int ptype)
{
        if(ptype < 1 || ptype > 3 )
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Incorrext problem type. ptype can be 1, 2, or 3!");
        }
        
        if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        if(A.Rows() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }
        
        int info;
        EV = A.Clone();
        Matrix<ComplexDouble> B2 = B.Clone();
        int n = B2.Rows();
        Vector<double> w(n);

        int bs = n*n+2*n;
        Vector<ComplexDouble> work(bs);

        int rwork_size = 2*n*n+5*n+1;
        Vector<double> rwork(rwork_size);

        int iwork_size = 5*n+3;
        Vector<int> iwork(iwork_size);

        ZHEGVD(&ptype, "V", "U", &n, reinterpret_cast<MKL_Complex16*>(EV.Data()), &n, reinterpret_cast<MKL_Complex16*>(B2.Data()), &n, w.Data(), reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), &rwork_size, iwork.Data(), &iwork_size, &info);
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of ZHEGVD is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        cerr << "cpotrf/zpotrf and cheev/zheev returned an error code: ";
                        if(info <= n)
                        {
                                cerr << "cheev/zheev failed to converge, and " << info <<  " off-diagonal elements of an intermediate tridiagonal did not converge to zero!";
                        }
                        else if(info > n && info <=  2*n)
                        {
                                cerr << "The leading minor of order " << info-n << " of B is not positive-definite. The factorization of B could not be completed and no eigenvalues or eigenvectors were computed!";
                        }
                        Utility::RunTimeError("");
                }
        }
        return w;
}





//=========================================================
// Generalized Non-Symmetric eigenvalue problem
//=========================================================


// Gges
void Lapack::Gges(Matrix<float>& A, Matrix<float>& B, Vector<ComplexFloat>& alpha, Vector<float>& beta, Matrix<float>& S, Matrix<float>& T)
{
        if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        if(A.Rows() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }

        
        int info;
        S = A.Clone();
        T = B.Clone();
        int n = A.Rows();
        Vector<float> wr(n);
        Vector<float> wi(n);
        beta = Vector<float>(n);

        int bs = 39*n+23;
        Vector<float> work(bs);

        int sdim;

        SGGES("N", "N", "N", 0, &n, S.Data(), &n, T.Data(), &n, &sdim, wr.Data(), wi.Data(), beta.Data(), 0, &n, 0, &n, work.Data(), &bs, 0, &info);

        alpha = Vector<ComplexFloat>(n);
        for(int i=0;i<n;i++)
        {
                alpha[i] = ComplexFloat(wr[i],wi[i]);
        }

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of SGGES is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "Lapack error!");
                        Utility::RunTimeError(s);
                }
        }
}

void Lapack::Gges(Matrix<float>& A, Matrix<float>& B, Vector<ComplexFloat>& alpha, Vector<float>& beta, Matrix<float>& S, Matrix<float>& T, Matrix<float>& Vl, Matrix<float>& Vr)
{
        if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        if(A.Rows() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }

        
        int info;
        S = A.Clone();
        T = B.Clone();
        int n = A.Rows();
        Vector<float> wr(n);
        Vector<float> wi(n);
        beta = Vector<float>(n);

        Vl = Matrix<float>(n,n);
        Vr = Matrix<float>(n,n);

        int bs = 39*n+23;
        Vector<float> work(bs);

        int sdim;

        SGGES("V", "V", "N", 0, &n, S.Data(), &n, T.Data(), &n, &sdim, wr.Data(), wi.Data(), beta.Data(), Vl.Data(), &n, Vr.Data(), &n, work.Data(), &bs, 0, &info);

        alpha = Vector<ComplexFloat>(n);
        for(int i=0;i<n;i++)
        {
                alpha[i] = ComplexFloat(wr[i],wi[i]);
        }

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of SGGES is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "Lapack error!");
                        Utility::RunTimeError(s);
                }
        }
}

void Lapack::Gges(Matrix<double>& A, Matrix<double>& B, Vector<ComplexDouble>& alpha, Vector<double>& beta, Matrix<double>& S, Matrix<double>& T)
{
        if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        if(A.Rows() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }

        
        int info;
        S = A.Clone();
        T = B.Clone();
        int n = A.Rows();
        Vector<double> wr(n);
        Vector<double> wi(n);
        beta = Vector<double>(n);

        int bs = 39*n+23;
        Vector<double> work(bs);

        int sdim;

        DGGES("N", "N", "N", 0, &n, S.Data(), &n, T.Data(), &n, &sdim, wr.Data(), wi.Data(), beta.Data(), 0, &n, 0, &n, work.Data(), &bs, 0, &info);

        alpha = Vector<ComplexDouble>(n);
        for(int i=0;i<n;i++)
        {
                alpha[i] = ComplexDouble(wr[i],wi[i]);
        }

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of DGGES is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "Lapack error!");
                        Utility::RunTimeError(s);
                }
        }
}


void Lapack::Gges(Matrix<double>& A, Matrix<double>& B, Vector<ComplexDouble>& alpha, Vector<double>& beta, Matrix<double>& S, Matrix<double>& T, Matrix<double>& Vl, Matrix<double>& Vr)
{
        if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        if(A.Rows() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }

        
        int info;
        S = A.Clone();
        T = B.Clone();
        int n = A.Rows();
        Vector<double> wr(n);
        Vector<double> wi(n);
        beta = Vector<double>(n);

        Vl = Matrix<double>(n,n);
        Vr = Matrix<double>(n,n);

        int bs = 39*n+23;
        Vector<double> work(bs);

        int sdim;

        DGGES("V", "V", "N", 0, &n, S.Data(), &n, T.Data(), &n, &sdim, wr.Data(), wi.Data(), beta.Data(), Vl.Data(), &n, Vr.Data(), &n, work.Data(), &bs, 0, &info);

        alpha = Vector<ComplexDouble>(n);
        for(int i=0;i<n;i++)
        {
                alpha[i] = ComplexDouble(wr[i],wi[i]);
        }

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of DGGES is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "Lapack error!");
                        Utility::RunTimeError(s);
                }
        }
}



void Lapack::Gges(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B, Vector<ComplexFloat>& alpha, Vector<ComplexFloat>& beta, Matrix<ComplexFloat>& S, Matrix<ComplexFloat>& T)
{
        if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        if(A.Rows() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }

        
        int info;
        S = A.Clone();
        T = B.Clone();
        int n = A.Rows();
        alpha = Vector<ComplexFloat>(n);
        beta = Vector<ComplexFloat>(n);

        int bs = 33*n;
        Vector<ComplexFloat> work(bs);

        Vector<float> rwork(8*n);
        int sdim;

        CGGES("N", "N", "N", 0, &n, reinterpret_cast<MKL_Complex8*>(S.Data()), &n, reinterpret_cast<MKL_Complex8*>(T.Data()), &n, &sdim, reinterpret_cast<MKL_Complex8*>(alpha.Data()), reinterpret_cast<MKL_Complex8*>(beta.Data()), 0, &n, 0, &n, reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), 0, &info);


        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of CGGES is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "Lapack error!");
                        Utility::RunTimeError(s);
                }
        }
}

void Lapack::Gges(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B, Vector<ComplexFloat>& alpha, Vector<ComplexFloat>& beta, Matrix<ComplexFloat>& S, Matrix<ComplexFloat>& T, Matrix<ComplexFloat>& Vl, Matrix<ComplexFloat>& Vr)
{
        if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        if(A.Rows() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }

        
        int info;
        S = A.Clone();
        T = B.Clone();
        int n = A.Rows();
        alpha = Vector<ComplexFloat>(n);
        beta = Vector<ComplexFloat>(n);

        Vl = Matrix<ComplexFloat>(n,n);
        Vr = Matrix<ComplexFloat>(n,n);

        int bs = 33*n;
        Vector<ComplexFloat> work(bs);

        Vector<float> rwork(8*n);
        int sdim;

        CGGES("V", "V", "N", 0, &n, reinterpret_cast<MKL_Complex8*>(S.Data()), &n, reinterpret_cast<MKL_Complex8*>(T.Data()), &n, &sdim, reinterpret_cast<MKL_Complex8*>(alpha.Data()), reinterpret_cast<MKL_Complex8*>(beta.Data()), reinterpret_cast<MKL_Complex8*>(Vl.Data()), &n, reinterpret_cast<MKL_Complex8*>(Vr.Data()), &n, reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), 0, &info);

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of CGGES is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "Lapack error!");
                        Utility::RunTimeError(s);
                }
        }
}

void Lapack::Gges(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B, Vector<ComplexDouble>& alpha, Vector<ComplexDouble>& beta, Matrix<ComplexDouble>& S, Matrix<ComplexDouble>& T)
{
        if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        if(A.Rows() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }

        
        int info;
        S = A.Clone();
        T = B.Clone();
        int n = A.Rows();
        alpha = Vector<ComplexDouble>(n);
        beta = Vector<ComplexDouble>(n);

        int bs = 33*n;
        Vector<ComplexDouble> work(bs);

        Vector<double> rwork(8*n);
        int sdim;

        ZGGES("N", "N", "N", 0, &n, reinterpret_cast<MKL_Complex16*>(S.Data()), &n, reinterpret_cast<MKL_Complex16*>(T.Data()), &n, &sdim, reinterpret_cast<MKL_Complex16*>(alpha.Data()), reinterpret_cast<MKL_Complex16*>(beta.Data()), 0, &n, 0, &n, reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), 0, &info);


        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of ZGGES is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "Lapack error!");
                        Utility::RunTimeError(s);
                }
        }
}

void Lapack::Gges(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B, Vector<ComplexDouble>& alpha, Vector<ComplexDouble>& beta, Matrix<ComplexDouble>& S, Matrix<ComplexDouble>& T, Matrix<ComplexDouble>& Vl, Matrix<ComplexDouble>& Vr)
{
        if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        if(A.Rows() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }

        
        int info;
        S = A.Clone();
        T = B.Clone();
        int n = A.Rows();
        alpha = Vector<ComplexDouble>(n);
        beta = Vector<ComplexDouble>(n);

        Vl = Matrix<ComplexDouble>(n,n);
        Vr = Matrix<ComplexDouble>(n,n);

        int bs = 33*n;
        Vector<ComplexDouble> work(bs);

        Vector<double> rwork(8*n);
        int sdim;

        ZGGES("V", "V", "N", 0, &n, reinterpret_cast<MKL_Complex16*>(S.Data()), &n, reinterpret_cast<MKL_Complex16*>(T.Data()), &n, &sdim, reinterpret_cast<MKL_Complex16*>(alpha.Data()), reinterpret_cast<MKL_Complex16*>(beta.Data()), reinterpret_cast<MKL_Complex16*>(Vl.Data()), &n, reinterpret_cast<MKL_Complex16*>(Vr.Data()), &n, reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), 0, &info);


        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of CGGES is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "Lapack error!");
                        Utility::RunTimeError(s);
                }
        }
}







// Ggev

void Lapack::Ggev(Matrix<float>& A, Matrix<float>& B, Vector<ComplexFloat>& alpha, Vector<float>& beta)
{
        if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        if(A.Rows() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }

        
        int info;
        Matrix<float> A2 = A.Clone();
        Matrix<float> B2 = B.Clone();
        int n = A2.Rows();
        Vector<float> wr(n);
        Vector<float> wi(n);
        beta = Vector<float>(n);

        int bs = 8*n+16;
        Vector<float> work(bs);

        SGGEV("N", "N", &n, A2.Data(), &n, B2.Data(), &n, wr.Data(), wi.Data(), beta.Data(), 0, &n, 0, &n, work.Data(), &bs, &info);
        
        alpha = Vector<ComplexFloat>(n);
        for(int i=0;i<n;i++)
        {
                alpha[i] = ComplexFloat(wr[i],wi[i]);
        }

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of SGGEV is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "Lapack error!");
                        Utility::RunTimeError(s);
                }
        }
}

void Lapack::Ggev(Matrix<float>& A, Matrix<float>& B, Vector<ComplexFloat>& alpha, Vector<float>& beta, Matrix<ComplexFloat>& Vl, Matrix<ComplexFloat>& Vr)
{
        if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        if(A.Rows() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }
        
        int info;
        Matrix<float> A2 = A.Clone();
        Matrix<float> B2 = B.Clone();
        int n = A2.Rows();
        Vector<float> wr(n);
        Vector<float> wi(n);
        beta = Vector<float>(n);

        Vl = Matrix<ComplexFloat>(n,n);
        Vr = Matrix<ComplexFloat>(n,n);
        Matrix<float> vl(n,n);
        Matrix<float> vr(n,n);

        int bs = 8*n+16;
        Vector<float> work(bs);

        SGGEV("V", "V", &n, A2.Data(), &n, B2.Data(), &n, wr.Data(), wi.Data(), beta.Data(), vl.Data(), &n, vr.Data(), &n, work.Data(), &bs, &info);

        alpha = Vector<ComplexFloat>(n);
        bool conj = false;
        for(int i=0;i<n;i++)
        {
                alpha[i] = ComplexFloat(wr[i],wi[i]);
                if(fabs(wi[i]) < numeric_limits<float>::epsilon())
                {
                        for(int j=0;j<n;j++)
                        {
                                Vl.ElemNC(j,i) = ComplexFloat(vl.ElemNC(j,i),0);
                                Vr.ElemNC(j,i) = ComplexFloat(vr.ElemNC(j,i),0);
                        }
                }
                else
                {
                        for(int j=0;j<n;j++)
                        {
                                if(conj)
                                {
                                        Vl.ElemNC(j,i) = ComplexFloat(vl.ElemNC(j,i-1),-vl.ElemNC(j,i));
                                        Vr.ElemNC(j,i) = ComplexFloat(vr.ElemNC(j,i-1),-vr.ElemNC(j,i));
                                }
                                else
                                {
                                        Vl.ElemNC(j,i) = ComplexFloat(vl.ElemNC(j,i),vl.ElemNC(j,i+1));
                                        Vr.ElemNC(j,i) = ComplexFloat(vr.ElemNC(j,i),vr.ElemNC(j,i+1));
                                }
                        }
                        conj = !conj;
                }
        }

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of SGGEV is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "Lapack error!");
                        Utility::RunTimeError(s);
                }
        }
}

void Lapack::Ggev(Matrix<double>& A, Matrix<double>& B, Vector<ComplexDouble>& alpha, Vector<double>& beta)
{
        if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        if(A.Rows() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }

        
        int info;
        Matrix<double> A2 = A.Clone();
        Matrix<double> B2 = B.Clone();
        int n = A2.Rows();
        Vector<double> wr(n);
        Vector<double> wi(n);
        beta = Vector<double>(n);

        int bs = 8*n+16;
        Vector<double> work(bs);

        DGGEV("N", "N", &n, A2.Data(), &n, B2.Data(), &n, wr.Data(), wi.Data(), beta.Data(), 0, &n, 0, &n, work.Data(), &bs, &info);
        alpha = Vector<ComplexDouble>(n);
        for(int i=0;i<n;i++)
        {
                alpha[i] = ComplexDouble(wr[i],wi[i]);
        }

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of DGGEV is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "Lapack error!");
                        Utility::RunTimeError(s);
                }
        }
}

void Lapack::Ggev(Matrix<double>& A, Matrix<double>& B, Vector<ComplexDouble>& alpha, Vector<double>& beta, Matrix<ComplexDouble>& Vl, Matrix<ComplexDouble>& Vr)
{
        if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        if(A.Rows() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }
        
        int info;
        Matrix<double> A2 = A.Clone();
        Matrix<double> B2 = B.Clone();
        int n = A2.Rows();
        Vector<double> wr(n);
        Vector<double> wi(n);
        beta = Vector<double>(n);

        Vl = Matrix<ComplexDouble>(n,n);
        Vr = Matrix<ComplexDouble>(n,n);
        Matrix<double> vl(n,n);
        Matrix<double> vr(n,n);

        int bs = 8*n+16;
        Vector<double> work(bs);

        DGGEV("V", "V", &n, A2.Data(), &n, B2.Data(), &n, wr.Data(), wi.Data(), beta.Data(), vl.Data(), &n, vr.Data(), &n, work.Data(), &bs, &info);

        alpha = Vector<ComplexDouble>(n);
        bool conj = false;
        for(int i=0;i<n;i++)
        {
                alpha[i] = ComplexDouble(wr[i],wi[i]);
                if(fabs(wi[i]) < numeric_limits<double>::epsilon())
                {
                        for(int j=0;j<n;j++)
                        {
                                Vl.ElemNC(j,i) = ComplexDouble(vl.ElemNC(j,i),0);
                                Vr.ElemNC(j,i) = ComplexDouble(vr.ElemNC(j,i),0);
                        }
                }
                else
                {
                        for(int j=0;j<n;j++)
                        {
                                if(conj)
                                {
                                        Vl.ElemNC(j,i) = ComplexDouble(vl.ElemNC(j,i-1),-vl.ElemNC(j,i));
                                        Vr.ElemNC(j,i) = ComplexDouble(vr.ElemNC(j,i-1),-vr.ElemNC(j,i));
                                }
                                else
                                {
                                        Vl.ElemNC(j,i) = ComplexDouble(vl.ElemNC(j,i),vl.ElemNC(j,i+1));
                                        Vr.ElemNC(j,i) = ComplexDouble(vr.ElemNC(j,i),vr.ElemNC(j,i+1));
                                }
                        }
                        conj = !conj;
                }
        }

        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of DGGEV is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "Lapack error!");
                        Utility::RunTimeError(s);
                }
        }
}

void Lapack::Ggev(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B, Vector<ComplexFloat>& alpha, Vector<ComplexFloat>& beta)
{
        if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        if(A.Rows() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }

        
        int info;
        Matrix<ComplexFloat> A2 = A.Clone();
        Matrix<ComplexFloat> B2 = B.Clone();
        int n = A2.Rows();
        alpha = Vector<ComplexFloat>(n);
        beta = Vector<ComplexFloat>(n);

        int bs = 2*n;
        Vector<ComplexFloat> work(bs);
        int rwork_size = 8*n;
        Vector<float> rwork(rwork_size);

        CGGEV("N", "N", &n, reinterpret_cast<MKL_Complex8*>(A2.Data()), &n, reinterpret_cast<MKL_Complex8*>(B2.Data()), &n, reinterpret_cast<MKL_Complex8*>(alpha.Data()), reinterpret_cast<MKL_Complex8*>(beta.Data()), 0, &n, 0, &n, reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), &info);
        
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of CGGEV is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "Lapack error!");
                        Utility::RunTimeError(s);
                }
        }
}

void Lapack::Ggev(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B, Vector<ComplexFloat>& alpha, Vector<ComplexFloat>& beta, Matrix<ComplexFloat>& Vl, Matrix<ComplexFloat>& Vr)
{
        if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        if(A.Rows() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }

        
        int info;
        Matrix<ComplexFloat> A2 = A.Clone();
        Matrix<ComplexFloat> B2 = B.Clone();
        int n = A2.Rows();
        alpha = Vector<ComplexFloat>(n);
        beta = Vector<ComplexFloat>(n);

        Vl = Matrix<ComplexFloat>(n,n);
        Vr = Matrix<ComplexFloat>(n,n);

        int bs = 2*n;
        Vector<ComplexFloat> work(bs);
        int rwork_size = 8*n;
        Vector<float> rwork(rwork_size);

        CGGEV("V", "V", &n, reinterpret_cast<MKL_Complex8*>(A2.Data()), &n, reinterpret_cast<MKL_Complex8*>(B2.Data()), &n, reinterpret_cast<MKL_Complex8*>(alpha.Data()), reinterpret_cast<MKL_Complex8*>(beta.Data()), reinterpret_cast<MKL_Complex8*>(Vl.Data()), &n, reinterpret_cast<MKL_Complex8*>(Vr.Data()), &n, reinterpret_cast<MKL_Complex8*>(work.Data()), &bs, rwork.Data(), &info);
        
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of CGGEV is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "Lapack error!");
                        Utility::RunTimeError(s);
                }
        }
}

void Lapack::Ggev(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B, Vector<ComplexDouble>& alpha, Vector<ComplexDouble>& beta)
{
        if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        if(A.Rows() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }

        
        int info;
        Matrix<ComplexDouble> A2 = A.Clone();
        Matrix<ComplexDouble> B2 = B.Clone();
        int n = A2.Rows();
        alpha = Vector<ComplexDouble>(n);
        beta = Vector<ComplexDouble>(n);

        int bs = 2*n;
        Vector<ComplexDouble> work(bs);
        int rwork_size = 8*n;
        Vector<double> rwork(rwork_size);

        ZGGEV("N", "N", &n, reinterpret_cast<MKL_Complex16*>(A2.Data()), &n, reinterpret_cast<MKL_Complex16*>(B2.Data()), &n, reinterpret_cast<MKL_Complex16*>(alpha.Data()), reinterpret_cast<MKL_Complex16*>(beta.Data()), 0, &n, 0, &n, reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), &info);
        
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of ZGGEV is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "Lapack error!");
                        Utility::RunTimeError(s);
                }
        }
}

void Lapack::Ggev(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B, Vector<ComplexDouble>& alpha, Vector<ComplexDouble>& beta, Matrix<ComplexDouble>& Vl, Matrix<ComplexDouble>& Vr)
{
        if(A.Rows() != A.Columns() || B.Rows() != B.Columns())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix is not square!");
        }
        if(A.Rows() != B.Rows())
        {
                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                Utility::RunTimeError("Matrix dimensions does not match!");
        }

        
        int info;
        Matrix<ComplexDouble> A2 = A.Clone();
        Matrix<ComplexDouble> B2 = B.Clone();
        int n = A2.Rows();
        alpha = Vector<ComplexDouble>(n);
        beta = Vector<ComplexDouble>(n);

        Vl = Matrix<ComplexDouble>(n,n);
        Vr = Matrix<ComplexDouble>(n,n);

        int bs = 2*n;
        Vector<ComplexDouble> work(bs);
        int rwork_size = 8*n;
        Vector<double> rwork(rwork_size);

        ZGGEV("V", "V", &n, reinterpret_cast<MKL_Complex16*>(A2.Data()), &n, reinterpret_cast<MKL_Complex16*>(B2.Data()), &n, reinterpret_cast<MKL_Complex16*>(alpha.Data()), reinterpret_cast<MKL_Complex16*>(beta.Data()), reinterpret_cast<MKL_Complex16*>(Vl.Data()), &n, reinterpret_cast<MKL_Complex16*>(Vr.Data()), &n, reinterpret_cast<MKL_Complex16*>(work.Data()), &bs, rwork.Data(), &info);
        
        if(info != 0)
        {
                if(info < 0)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "%i'th argument of ZGGEV is an illegal value!", -info);
                        Utility::RunTimeError(s);
                }
                else if(info > 0)  
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        char s[500];
                        sprintf(s, "Lapack error!");
                        Utility::RunTimeError(s);
                }
        }
}

---