mtl2lapack.h

MTL C++ Numeric Library

/* -*- 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)