/* -*- c++ -*- *///// Copyright 1997, 1998, 1999 University of Notre Dame.// Authors: Andrew Lumsdaine, Jeremy G. Siek, Lie-Quan Lee//// This file is part of the Matrix Template Library//// You should have received a copy of the License Agreement for the// Matrix Template Library along with the software;  see the// file LICENSE.  If not, contact Office of Research, University of Notre// Dame, Notre Dame, IN  46556.//// Permission to modify the code and to distribute modified code is// granted, provided the text of this NOTICE is retained, a notice that// the code was modified is included with the above COPYRIGHT NOTICE and// with the COPYRIGHT NOTICE in the LICENSE file, and that the LICENSE// file is distributed with the modified code.//// LICENSOR MAKES NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED.// By way of example, but not limitation, Licensor MAKES NO// REPRESENTATIONS OR WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY// PARTICULAR PURPOSE OR THAT THE USE OF THE LICENSED SOFTWARE COMPONENTS// OR DOCUMENTATION WILL NOT INFRINGE ANY PATENTS, COPYRIGHTS, TRADEMARKS// OR OTHER RIGHTS.//#ifndef MTL_LAPACK_H#define MTL_LAPACK_H#include <iostream>#include "mtl/mtl_complex.h"#include "mtl/lapack_interface.h"#include "mtl/matrix.h"#define MTL_FCALL(x) x##_namespace mtl_lapack_dispatch {/* *  Dispatch routines used by the mtl2lapack functions *  (mtl2lapack functions are below) */  inline void gecon (const char & norm, const int & n,		     const double da[], const int & lda,		     const double & danorm, double & drcond,		     int & info)  {    double* dwork  = new double[4*n];    int*    iwork2 = new int[n];    MTL_FCALL(dgecon)(norm, n, da, lda, danorm, drcond, dwork, iwork2, info);    delete [] dwork;    delete [] iwork2;  }  inline void gecon (const char & norm, const int & n,		     const float sa[], const int & lda,		     const float & sanorm, float & srcond,		     int & info)  {    float* swork  = new float[4*n];    int*   iwork2 = new int[n];    MTL_FCALL(sgecon)(norm, n, sa, lda, sanorm, srcond, swork, iwork2, info);    delete [] swork;    delete [] iwork2;  }  inline void gecon(const char & norm, const int & n,		    const std::complex<double> za[], const int & lda,		    const double & danorm, double & drcond,		    int & info)  {    std::complex<double>* zwork  = new std::complex<double>[2*n];    double*          dwork2 = new double[2*n];    MTL_FCALL(zgecon)(norm, n, za, lda, danorm, drcond, zwork, dwork2, info);    delete [] zwork;    delete [] dwork2;  }  inline void gecon (const char & norm, const int & n,		     const std::complex<float> ca[], const int & lda,		     const float & canorm, float & crcond, int & info)  {    std::complex<float>* cwork  = new std::complex<float>[2*n];    float*          cwork2 = new float[2*n];    MTL_FCALL(cgecon)(norm, n, ca, lda, canorm, crcond, cwork, cwork2, info);    delete [] cwork;    delete [] cwork2;  }  inline void geev(const char& jobvl, const char& jobvr, const int&  n,		   double da[], const int& lda, std::complex<double> zw[],		   double dvl[], const int& ldvl, double dvr[],		   const int& ldvr, int& info)  {    int     ldwork = 4*n+1;    double* dwork  = new double[ldwork];    double* dwr    = new double[n];    double* dwi    = new double[n];    MTL_FCALL(dgeev)(jobvl, jobvr, n, da, lda, dwr, dwi, dvl, ldvl, dvr, ldvr, dwork,	   ldwork, info);    /*assemble dwr and dwi into zw */    for (int i = 0; i < n; i++)      zw[i] = std::complex<double>(dwr[i],dwi[i]);    delete [] dwork;    delete [] dwr;    delete [] dwi;  }  inline void geev(const char& jobvl, const char& jobvr, const int&  n,		   float sa[], const int& lda, std::complex<float> cw[],		   float svl[], const int& ldvl, float svr[],		   const int& ldvr, int& info)  {    int     ldwork = 4*n+1;    float*  swork  = new float[ldwork];    float*  swr    = new float[n];    float*  swi    = new float[n];    MTL_FCALL(sgeev)(jobvl, jobvr, n, sa, lda, swr, swi, svl, ldvl, svr, ldvr,	   swork, ldwork, info);    /* assemble dwr and dwi into zw */    for (int i = 0; i < n; i++)      cw[i] = std::complex<float>(swr[i],swi[i]);    delete [] swork;    delete [] swr;    delete [] swi;  }  inline void geev(const char& jobvl, const char& jobvr, const int&  n,		   std::complex<double> za[], const int& lda,		   std::complex<double> zw[], std::complex<double> zvl[],		   const int& ldvl, std::complex<double> zvr[],		   const int& ldvr, int& info)  {    int              ldwork  = 1+2*n;    std::complex<double>* zwork   = new std::complex<double>[ldwork];    double*          dwork2  = new double[2*n];    MTL_FCALL(zgeev)(jobvl, jobvr, n, za, lda, zw, zvl, ldvl, zvr, ldvr, zwork, ldwork,	   dwork2, info);    delete [] zwork;    delete [] dwork2;  }  inline void geev(const char& jobvl, const char& jobvr, const int&  n,		   std::complex<float> ca[], const int& lda,		   std::complex<float> cw[], std::complex<float> cvl[],		   const int& ldvl, std::complex<float> cvr[],		   const int& ldvr, int& info)  {    int              ldwork  = 1+2*n;    std::complex<float>*  cwork   = new std::complex<float>[ldwork];    float*           swork2  = new float[2*n];    MTL_FCALL(cgeev)(jobvl, jobvr, n, ca, lda, cw, cvl, ldvl, cvr, ldvr, cwork, ldwork,	   swork2, info);    delete [] cwork;    delete [] swork2;  }  inline void geqpf(const int & m, const int & n,		    double da[], const int & lda,		    int jpivot[], double dtau[],		    int & info)  {    double* dwork = new double[3*n];    MTL_FCALL(dgeqpf)(m, n, da, lda, jpivot, dtau, dwork, info);    delete [] dwork;  }  inline void geqpf(const int & m, const int & n,		    float sa[], const int & lda,		    int jpivot[], float stau[],		    int & info)  {    float* swork = new float[3*n];    MTL_FCALL(sgeqpf)(m, n, sa, lda, jpivot, stau, swork, info);    delete [] swork;  }  inline void geqpf(const int & m, const int & n,		    std::complex<double> za[], const int & lda,		    int jpivot[], std::complex<double> ztau[],		    int & info)  {    std::complex<double>* zwork = new std::complex<double>[n];    double * dwork2 = new double[2*n];    MTL_FCALL(zgeqpf)(m, n, za, lda, jpivot, ztau, zwork, dwork2, info);    delete [] zwork;    delete [] dwork2;  }  inline void geqpf(const int & m, const int & n,		    std::complex<float> ca[], const int & lda,		    int jpivot[], std::complex<float> ctau[],		    int & info)  {    std::complex<float>* cwork = new std::complex<float>[n];    float * swork2 = new float[2*n];    MTL_FCALL(cgeqpf)(m, n, ca, lda, jpivot, ctau, cwork, swork2, info);    delete [] swork2;  }  inline void geqrf(const int & m, const int & n,		    double da[], const int & lda,		    double dtau[], int & info)  {    int ldwork = n*n;    double* dwork   = new double[ldwork];    MTL_FCALL(dgeqrf)(m, n, da, lda, dtau, dwork, ldwork, info);    delete [] dwork;  }  inline void geqrf(const int & m, const int & n,		    float sa[], const int & lda,		    float stau[], int & info)  {    int ldwork = 4*n;    float* swork = new float[ldwork];    MTL_FCALL(sgeqrf)(m, n, sa, lda, stau, swork, ldwork, info);    delete [] swork;  }  inline void geqrf(const int & m, const int & n,		    std::complex<double> za[], const int & lda,		    std::complex<double> ztau[], int & info)  {    int ldwork = 4*n;    std::complex<double>* zwork = new std::complex<double>[ldwork];    MTL_FCALL(zgeqrf)(m, n, za, lda, ztau, zwork, ldwork, info);    delete [] zwork;  }  inline void geqrf(const int & m, const int & n,		    std::complex<float> ca[], const int & lda,		    std::complex<float> ctau[], int & info)  {    int ldwork = 4*n;    std::complex<float>* cwork = new std::complex<float>[ldwork];    MTL_FCALL(cgeqrf)(m, n, ca, lda, ctau, cwork, ldwork, info);    delete [] cwork;  }  inline void gesv(const int & n, const int & nrhs,		   double da[], const int & lda, int ipivot[],		   double db[], const int & ldb, int & info)  {    MTL_FCALL(dgesv)(n, nrhs, da, lda, ipivot, db, ldb, info);  }  inline void gesv(const int & n, const int & nrhs,		   float sa[], const int & lda,		   int ipivot[], float sb[],		   const int & ldb, int & info)  {    MTL_FCALL(sgesv)(n, nrhs, sa, lda, ipivot, sb, ldb, info);  }  inline void gesv(const int & n, const int & nrhs,		   std::complex<double> za[], const int & lda,		   int ipivot[], std::complex<double> zb[],		   const int & ldb, int & info)  {    MTL_FCALL(zgesv)(n, nrhs, za, lda, ipivot, zb, ldb, info);  }  inline void gesv(const int & n, const int & nrhs,		   std::complex<float> ca[], const int & lda,		   int ipivot[], std::complex<float> cb[],		   const int & ldb, int & info)  {    MTL_FCALL(cgesv)(n, nrhs, ca, lda, ipivot, cb, ldb, info);  }  inline void getrf (const int& m, const int& n,		     double da[], const int& lda,		     int ipivot[], int& info)  {    MTL_FCALL(dgetrf)(m, n, da, lda, ipivot, info);  }  inline void getrf (const int& m, const int& n,		     float sa[], const int& lda,		     int ipivot[], int& info)  {    MTL_FCALL(sgetrf)(m, n, sa, lda, ipivot, info);  }   inline void getrf (const int& m, const int& n,		     std::complex<double> za[], const int& lda,		     int ipivot[], int& info)  {    MTL_FCALL(zgetrf)(m, n, za, lda, ipivot, info);  }  inline void getrf (const int& m, const int& n,		     std::complex<float> ca[], const int& lda,		     int ipivot[], int& info)  {    MTL_FCALL(cgetrf)(m, n, ca, lda, ipivot, info);  }  inline void getrs (const char& transa, const int& n,		     const int& nrhs, const double da[],		     const int& lda, int ipivot[],		     double db[], const int& ldb,		     int& info)  {    MTL_FCALL(dgetrs)(transa, n, nrhs, da, lda, ipivot, db, ldb, info);  }  inline void getrs (const char& transa, const int& n,		     const int& nrhs, const float  sa[],		     const int& lda, int ipivot[],		     float sb[], const int& ldb,		     int& info)  {    MTL_FCALL(sgetrs)(transa, n, nrhs, sa, lda, ipivot, sb, ldb, info);  }  inline void getrs (const char& transa, const int& n,		     const int& nrhs, const std::complex<double> za[],		     const int& lda, int ipivot[],		     std::complex<double> zb[], const int& ldb,		     int& info)  {    MTL_FCALL(zgetrs)(transa, n, nrhs, za, lda, ipivot, zb, ldb, info);  }  inline void getrs (const char& transa, const int& n,		     const int& nrhs, const std::complex<float> ca[],		     const int& lda, int ipivot[],		     std::complex<float> cb[], const int& ldb,		     int& info)  {    MTL_FCALL(cgetrs)(transa, n, nrhs, ca, lda, ipivot, cb, ldb, info);  }  inline void geequ (const int& m, const int& n,		     const double da[], const int& lda,		     double r[], double c[], 		     double& rowcnd, double& colcnd,		     double& amax, int& info)  {    MTL_FCALL(dgeequ)(m, n, da, lda, r, c, rowcnd, colcnd, amax, info);  }  inline void geequ(const int& m, const int& n,		    const float da[], const int& lda,		    float r[], float c[], float& rowcnd,		    float& colcnd,		    float& amax, int& info)  {    MTL_FCALL(sgeequ)(m, n, da, lda, r,  c, rowcnd, colcnd, amax, info);  }  inline void geequ(const int& m, const int& n, 		    const std::complex<float> da[],		    const int& lda, std::complex<float> r[],		    std::complex<float> c[],		    float& rowcnd, float& colcnd,		    float& amax, int& info)  {    MTL_FCALL(cgeequ)(m, n, da, lda, r,  c, rowcnd, colcnd, amax, info);  }  inline void geequ(const int& m, const int& n,		    const std::complex<double> da[],		    const int& lda, std::complex<double> r[],		    std::complex<double> c[],		    double& rowcnd, double& colcnd,		    double& amax, int& info)  {    MTL_FCALL(zgeequ)(m, n, da, lda, r,  c, rowcnd, colcnd, amax, info);  }  inline void gelqf(const int& m, const int& n,		    double da[], const int& lda,		    double dtau[], int& info)