Complex version of $function
", @_);
}
#pp_bless('PDL::Complex'); pp_addpm({At=>'Top'},<<'EOD'); use strict; use PDL::Complex; use PDL::LinearAlgebra::Real;
{ package PDL; my $warningFlag; BEGIN{ $warningFlag = $^W; $^W = 0; } use overload ( 'x' => sub {UNIVERSAL::isa($_[1],'PDL::Complex') ? PDL::cmmult(PDL::Complex::r2C($_[0]), $_[1]): PDL::mmult($_[0], $_[1]); }); BEGIN{ $^W = $warningFlag ; } } { package PDL::Complex; my $warningFlag; BEGIN{ $warningFlag = $^W; $^W = 0; } use overload ( 'x' => sub {UNIVERSAL::isa($_[1],'PDL::Complex') ? PDL::cmmult($_[0], $_[1]) : PDL::cmmult($_[0], PDL::Complex::r2C($_[1])); }, ); BEGIN{ $^W = $warningFlag ; } }
PDL::LinearAlgebra::Complex - PDL interface to the lapack linear algebra programming library (complex number)
use PDL::Complex use PDL::LinearAlgebra::Complex; $a = r2C random (100,100); $s = r2C zeroes(100); $u = r2C zeroes(100,100); $v = r2C zeroes(100,100); $info = 0; $job = 0; cgesdd($a, $job, $info, $s , $u, $v);
This module provides an interface to parts of the lapack library (complex numbers). These routines accept either float or double piddles.
EOD
pp_defc("gesvd", HandleBad => 0, RedoDimsCode => '$SIZE(r) = $PDL(A)->ndims > 2 ? min($PDL(A)->dims[1], $PDL(A)->dims[2]) : 1;', Pars => '[io,phys]A(2,m,n); int jobu(); int jobvt(); [o,phys]s(r); [o,phys]U(2,p,q); [o,phys]VT(2,s,t); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code '
integer lwork;
char trau, travt;
types(F) %{
extern int cgesvd_(char *jobu, char *jobvt, integer *m, integer *n, float *a,
integer *lda, float *s, float *u, int *ldu,
float *vt, integer *ldvt, float *work, integer *lwork, float *rwork,
integer *info);
float *rwork;
float tmp_work[2];
%}
types(D) %{
extern int zgesvd_(char *jobz,char *jobvt, integer *m, integer *n,
double *a, integer *lda, double *s, double *u, int *ldu,
double *vt, integer *ldvt, double *work, integer *lwork, double *rwork,
integer *info);
double *rwork;
double tmp_work[2];
%}
lwork = ($PRIV(__m_size) < $PRIV(__n_size)) ? 5*$PRIV(__m_size) : 5*$PRIV(__n_size);
types(F) %{
rwork = (float *)malloc(lwork * sizeof(float));
%}
types(D) %{
rwork = (double *)malloc(lwork * sizeof(double));
%}
lwork = -1;
switch ($jobu())
{
case 1: trau = \'A\';
break;
case 2: trau = \'S\';
break;
case 3: trau = \'O\';
break;
default: trau = \'N\';
}
switch ($jobvt())
{
case 1: travt = \'A\';
break;
case 2: travt = \'S\';
break;
case 3: travt = \'O\';
break;
default: travt = \'N\';
}
$TFD(cgesvd_,zgesvd_)(
&trau,
&travt,
&$PRIV(__m_size),
&$PRIV(__n_size),
$P(A),
&$PRIV(__m_size),
$P(s),
$P(U),
&$PRIV(__p_size),
$P(VT),
&$PRIV(__s_size),
&tmp_work[0],
&lwork,
rwork,
$P(info));
lwork = (integer )tmp_work[0];
{
types(F) %{
float *work = (float *)malloc(2*lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2*lwork * sizeof(double));
%}
$TFD(cgesvd_,zgesvd_)(
&trau,
&travt,
&$PRIV(__m_size),
&$PRIV(__n_size),
$P(A),
&$PRIV(__m_size),
$P(s),
$P(U),
&$PRIV(__p_size),
$P(VT),
&$PRIV(__s_size),
work,
&lwork,
rwork,
$P(info));
free(work);
}
free(rwork);
',
Doc=>'
Complex version of gesvd.
The SVD is written
A = U * SIGMA * ConjugateTranspose(V)
');
pp_defc("gesdd", HandleBad => 0, RedoDimsCode => '$SIZE(r) = $PDL(A)->ndims > 2 ? min($PDL(A)->dims[1], $PDL(A)->dims[2]) : 1;', Pars => '[io,phys]A(2,m,n); int job(); [o,phys]s(r); [o,phys]U(2,p,q); [o,phys]VT(2,s,t); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code ' integer lwork; integer *iwork; char tra; types(F) %{
extern int cgesdd_(char *jobz, integer *m, integer *n, float *
a, integer *lda, float *s, float *u, int *ldu,
float *vt, integer *ldvt, float *work, integer *lwork,
float *rwork, integer *iwork, integer *info);
float *rwork;
float tmp_work[2];
%}
types(D) %{
extern int zgesdd_(char *jobz, integer *m, integer *n, double *
a, integer *lda, double *s, double *u, int *ldu,
double *vt, integer *ldvt, double *work, integer *lwork,
double *rwork, integer *iwork, integer *info);
double *rwork;
double tmp_work[2];
%}
lwork = ($PRIV(__m_size) < $PRIV(__n_size)) ? $PRIV(__m_size) : $PRIV(__n_size);
iwork = (integer *)malloc(lwork * 8 * sizeof(integer));
types(F) %{
switch ($job())
{
case 1: tra = \'A\';
rwork = (float *)malloc( (5 * lwork *lwork + 5 * lwork) * sizeof(float));
break;
case 2: tra = \'S\';
rwork = (float *)malloc( (5 * lwork *lwork + 5 * lwork) * sizeof(float));
break;
case 3: tra = \'O\';
rwork = (float *)malloc( (5 * lwork *lwork + 5 * lwork) * sizeof(float));
break;
default: tra = \'N\';
rwork = (float *)malloc( 7 * lwork * sizeof(float));
break;
}
%}
types(D) %{
switch ($job())
{
case 1: tra = \'A\';
rwork = (double *)malloc( (5 * lwork *lwork + 5 * lwork) * sizeof(double));
break;
case 2: tra = \'S\';
rwork = (double *)malloc( (5 * lwork *lwork + 5 * lwork) * sizeof(double));
break;
case 3: tra = \'O\';
rwork = (double *)malloc( (5 * lwork *lwork + 5 * lwork) * sizeof(double));
break;
default: tra = \'N\';
rwork = (double *)malloc( 7 * lwork * sizeof(double));
break;
}
%}
lwork = -1;
$TFD(cgesdd_,zgesdd_)(
&tra,
&$PRIV(__m_size),
&$PRIV(__n_size),
$P(A),
&$PRIV(__m_size),
$P(s),
$P(U),
&$PRIV(__p_size),
$P(VT),
&$PRIV(__s_size),
&tmp_work[0],
&lwork,
rwork,
iwork,
$P(info));
lwork = (integer )tmp_work[0];
{
types(F) %{
float *work = (float *)malloc(2 * lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2 * lwork * sizeof(double));
%}
$TFD(cgesdd_,zgesdd_)(
&tra,
&$PRIV(__m_size),
&$PRIV(__n_size),
$P(A),
&$PRIV(__m_size),
$P(s),
$P(U),
&$PRIV(__p_size),
$P(VT),
&$PRIV(__s_size),
work,
&lwork,
rwork,
iwork,
$P(info));
free(work);
}
free(iwork);
free(rwork);
',
Doc=>'
Complex version of gesdd.
The SVD is written
A = U * SIGMA * ConjugateTranspose(V)
');
pp_defc("ggsvd", HandleBad => 0, Pars => '[io,phys]A(2,m,n); int jobu(); int jobv(); int jobq(); [io,phys]B(2,p,n); int [o,phys]k(); int [o,phys]l();[o,phys]alpha(n);[o,phys]beta(n); [o,phys]U(2,q,r); [o,phys]V(2,s,t); [o,phys]Q(2,u,v); int [o,phys]iwork(n); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code ' char pjobu = \'N\'; char pjobv = \'N\'; char pjobq = \'N\';
types(F) %{
extern int cggsvd_(char *jobu, char *jobv, char *jobq, integer *m,
integer *n, integer *p, integer *k, integer *l, float *a,
integer *lda, float *b, integer *ldb, float *alpha,
float *beta, float *u, integer *ldu, float *v, integer
*ldv, float *q, integer *ldq, float *work, float *rwork, integer *iwork,
integer *info);
float *work, *rwork;
%}
types(D) %{
extern int zggsvd_(char *jobu, char *jobv, char *jobq, integer *m,
integer *n, integer *p, integer *k, integer *l, double *a,
integer *lda, double *b, integer *ldb, double *alpha,
double *beta, double *u, integer *ldu, double *v, integer
*ldv, double *q, integer *ldq, double *work, double *rwork, integer *iwork,
integer *info);
double *work, *rwork;
%}
integer lwork = ($SIZE (m) < $SIZE (n)) ? $SIZE (n): $SIZE (m);
if ($SIZE (p) > lwork)
lwork = $SIZE (p);
types(F) %{
work = (float *)malloc(2*(3*lwork + $SIZE (n))* sizeof(float));
rwork = (float *)malloc(2 * $SIZE (n) * sizeof(float));
%}
types(D) %{
work = (double *)malloc(2*(3*lwork + $SIZE (n)) * sizeof(double));
rwork = (double *)malloc(2 * $SIZE (n) * sizeof(double));
%}
if ($jobu())
pjobu = \'U\';
if ($jobv())
pjobv = \'V\';
if ($jobq())
pjobq = \'Q\';
$TFD(cggsvd_,zggsvd_)(
&pjobu,
&pjobv,
&pjobq,
&$PRIV(__m_size),
&$PRIV(__n_size),
&$PRIV(__p_size),
$P(k),
$P(l),
$P(A),
&$PRIV(__m_size),
$P(B),
&$PRIV(__p_size),
$P(alpha),
$P(beta),
$P(U),
&$PRIV(__q_size),
$P(V),
&$PRIV(__s_size),
$P(Q),
&$PRIV(__u_size),
work,
rwork,
$P(iwork),
$P(info));
free(rwork);
free(work);
');
pp_defc("geev", HandleBad => 0, Pars => '[phys]A(2,n,n); int jobvl(); int jobvr(); [o,phys]w(2,n); [o,phys]vl(2,m,m); [o,phys]vr(2,p,p); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code '
char jvl = \'N\';
char jvr = \'N\';
types(F) %{
extern int cgeev_(char *jobvl, char *jobvr, integer *n, float *a,
integer *lda, float *w, float *vl, integer *ldvl, float *vr,
integer *ldvr, float *work, integer *lwork, float *rwork, integer *info);
float tmp_work[2], *rwork;
%}
types(D) %{
extern int zgeev_(char *jobvl, char *jobvr, integer *n, double *
a, integer *lda, double *w, double *vl,
integer *ldvl, double *vr, integer *ldvr, double *work,
integer *lwork, double *rwork, integer *info);
double tmp_work[2], *rwork;
%}
integer lwork = -1;
if ($jobvl())
jvl = \'V\';
if ($jobvr())
jvr = \'V\';
types(F) %{
rwork = (float *)malloc( 2 * $PRIV(__n_size) * sizeof(float));
%}
types(D) %{
rwork = (double *)malloc(2 * $PRIV(__n_size) * sizeof(double));
%}
$TFD(cgeev_,zgeev_)(
&jvl,
&jvr,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(w),
$P(vl),
&$PRIV(__m_size),
$P(vr),
&$PRIV(__p_size),
&tmp_work[0],
&lwork,
rwork,
$P(info));
lwork = (integer )tmp_work[0];
{
types(F) %{
float *work = (float *)malloc(2*lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2*lwork * sizeof(double));
%}
$TFD(cgeev_,zgeev_)(
&jvl,
&jvr,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(w),
$P(vl),
&$PRIV(__m_size),
$P(vr),
&$PRIV(__p_size),
work,
&lwork,
rwork,
$P(info));
free(work);
}
free(rwork);
');
pp_defc("geevx", HandleBad => 0, Pars => '[io,phys]A(2,n,n); int jobvl(); int jobvr(); int balance(); int sense(); [o,phys]w(2,n); [o,phys]vl(2,m,m); [o,phys]vr(2,p,p); int [o,phys]ilo(); int [o,phys]ihi(); [o,phys]scale(n); [o,phys]abnrm(); [o,phys]rconde(q); [o,phys]rcondv(r); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code '
char jvl = \'N\';
char jvr = \'N\';
char balanc, sens;
integer lwork = -1;
types(F) %{
extern int cgeevx_(char *balanc, char *jobvl, char *jobvr, char *
sense, integer *n, float *a, integer *lda, float *w,
float *vl, integer *ldvl, float *vr,
integer *ldvr, integer *ilo, integer *ihi, float *scale,
float *abnrm, float *rconde, float *rcondv, float
*work, integer *lwork, float *rwork, integer *info);
float tmp_work[2], *rwork;
%}
types(D) %{
extern int zgeevx_(char *balanc, char *jobvl, char *jobvr, char *
sense, integer *n, double *a, integer *lda, double *w,
double *vl, integer *ldvl, double *vr,
integer *ldvr, integer *ilo, integer *ihi, double *scale,
double *abnrm, double *rconde, double *rcondv, double
*work, integer *lwork, double *rwork, integer *info);
double tmp_work[2], *rwork;
%}
if ($jobvl())
jvl = \'V\';
if ($jobvr())
jvr = \'V\';
switch ($balance())
{
case 1: balanc = \'P\';
break;
case 2: balanc = \'S\';
break;
case 3: balanc = \'B\';
break;
default: balanc = \'N\';
}
switch ($sense())
{
case 1: sens = \'E\';
break;
case 2: sens = \'V\';
break;
case 3: sens = \'B\';
break;
default: sens = \'N\';
}
types(F) %{
rwork = (float *)malloc( 2 * $PRIV(__n_size) * sizeof(float));
%}
types(D) %{
rwork = (double *)malloc(2 * $PRIV(__n_size) * sizeof(double));
%}
$TFD(cgeevx_,zgeevx_)(
&balanc,
&jvl,
&jvr,
&sens,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(w),
$P(vl),
&$PRIV(__m_size),
$P(vr),
&$PRIV(__p_size),
$P(ilo),
$P(ihi),
$P(scale),
$P(abnrm),
$P(rconde),
$P(rcondv),
&tmp_work[0],
&lwork,
rwork,
$P(info));
lwork = (integer )tmp_work[0];
{
types(F) %{
float *work = (float *)malloc(2 * lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2 * lwork * sizeof(double));
%}
$TFD(cgeevx_,zgeevx_)(
&balanc,
&jvl,
&jvr,
&sens,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(w),
$P(vl),
&$PRIV(__m_size),
$P(vr),
&$PRIV(__p_size),
$P(ilo),
$P(ihi),
$P(scale),
$P(abnrm),
$P(rconde),
$P(rcondv),
work,
&lwork,
rwork,
$P(info));
free(work);
}
free(rwork);
');
pp_defc("ggev", HandleBad => 0, Pars => '[phys]A(2,n,n); int jobvl();int jobvr();[phys]B(2,n,n);[o,phys]alpha(2,n);[o,phys]beta(2,n);[o,phys]VL(2,m,m);[o,phys]VR(2,p,p);int [o,phys]info()', GenericTypes => [F,D], Code => generate_code '
integer lwork = -1;
char pjobvl = \'N\', pjobvr = \'N\';
types(F) %{
extern int cggev_(char *jobvl, char *jobvr, integer *n, float *
a, integer *lda, float *b, integer *ldb, float *alpha,
float *beta, float *vl, integer *ldvl,
float *vr, integer *ldvr, float *work, integer *lwork,
float *rwork, integer *info);
float tmp_work[2], *rwork;
%}
types(D) %{
extern int zggev_(char *jobvl, char *jobvr, integer *n, double *
a, integer *lda, double *b, integer *ldb, double *alpha,
double *beta, double *vl, integer *ldvl,
double *vr, integer *ldvr, double *work, integer *lwork,
double *rwork, integer *info);
double tmp_work[2], *rwork;
%}
if ($jobvl())
pjobvl = \'V\';
if ($jobvr())
pjobvr = \'V\';
types(F) %{
rwork = (float *)malloc(8 * $SIZE(n) * sizeof(float));
%}
types(D) %{
rwork = (double *)malloc(8 * $SIZE(n) * sizeof(double));
%}
$TFD(cggev_,zggev_)(
&pjobvl,
&pjobvr,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(B),
&$PRIV(__n_size),
$P(alpha),
$P(beta),
$P(VL),
&$PRIV(__m_size),
$P(VR),
&$PRIV(__p_size),
&tmp_work[0],
&lwork,
rwork,
$P(info));
lwork = (integer )tmp_work[0];
{
types(F) %{
float *work = (float *)malloc(2 * lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2 * lwork * sizeof(double));
%}
$TFD(cggev_,zggev_)(
&pjobvl,
&pjobvr,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(B),
&$PRIV(__n_size),
$P(alpha),
$P(beta),
$P(VL),
&$PRIV(__m_size),
$P(VR),
&$PRIV(__p_size),
work,
&lwork,
rwork,
$P(info));
free(work);
}
free(rwork);
');
pp_defc("ggevx", HandleBad => 0, Pars => '[io,phys]A(2,n,n);int balanc();int jobvl();int jobvr();int sense();[io,phys]B(2,n,n);[o,phys]alpha(2,n);[o,phys]beta(2,n);[o,phys]VL(2,m,m);[o,phys]VR(2,p,p);int [o,phys]ilo();int [o,phys]ihi();[o,phys]lscale(n);[o,phys]rscale(n);[o,phys]abnrm();[o,phys]bbnrm();[o,phys]rconde(r);[o,phys]rcondv(s);int [o,phys]info()', GenericTypes => [F,D], Code => generate_code '
integer lwork = -1, *iwork, *bwork;
char pjobvl = \'N\', pjobvr = \'N\';
char pbalanc, psens;
types(F) %{
int cggevx_(char *balanc, char *jobvl, char *jobvr, char *
sense, integer *n, float *a, integer *lda, float *b,
integer *ldb, float *alpha, float *
beta, float *vl, integer *ldvl, float *vr, integer *ldvr,
integer *ilo, integer *ihi, float *lscale, float *rscale,
float *abnrm, float *bbnrm, float *rconde, float *
rcondv, float *work, integer *lwork, float *rwork, integer *iwork, logical *
bwork, integer *info);
float tmp_work[2], *rwork;
rwork = (float *)malloc(6 * $SIZE(n) * sizeof(float));
%}
types(D) %{
extern int zggevx_(char *balanc, char *jobvl, char *jobvr, char *
sense, integer *n, double *a, integer *lda, double *b,
integer *ldb, double *alpha, double *
beta, double *vl, integer *ldvl, double *vr, integer *ldvr,
integer *ilo, integer *ihi, double *lscale, double *rscale,
double *abnrm, double *bbnrm, double *rconde, double *
rcondv, double *work, integer *lwork, double *rwork, integer *iwork, logical *
bwork, integer *info);
double tmp_work[2], *rwork;
rwork = (double *)malloc(6 * $SIZE(n) * sizeof(double));
%}
if ($jobvl())
pjobvl = \'V\';
if ($jobvr())
pjobvr = \'V\';
switch ($balanc())
{
case 1: pbalanc = \'P\';
break;
case 2: pbalanc = \'S\';
break;
case 3: pbalanc = \'B\';
break;
default: pbalanc = \'N\';
}
switch ($sense())
{
case 1: psens = \'E\';
bwork = (integer *)malloc($SIZE(n) * sizeof(integer));
break;
case 2: psens = \'V\';
iwork = (integer *)malloc(($SIZE(n) + 2) * sizeof(integer));
bwork = (integer *)malloc($SIZE(n) * sizeof(integer));
break;
case 3: psens = \'B\';
iwork = (integer *)malloc(($SIZE(n) + 2) * sizeof(integer));
bwork = (integer *)malloc($SIZE(n) * sizeof(integer));
break;
default: psens = \'N\';
iwork = (integer *)malloc(($SIZE(n) + 2) * sizeof(integer));
}
$TFD(cggevx_,zggevx_)(
&pbalanc,
&pjobvl,
&pjobvr,
&psens,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(B),
&$PRIV(__n_size),
$P(alpha),
$P(beta),
$P(VL),
&$PRIV(__m_size),
$P(VR),
&$PRIV(__p_size),
$P(ilo),
$P(ihi),
$P(lscale),
$P(rscale),
$P(abnrm),
$P(bbnrm),
$P(rconde),
$P(rcondv),
&tmp_work[0],
&lwork,
rwork,
iwork,
bwork,
$P(info));
lwork = (integer )tmp_work[0];
{
types(F) %{
float *work = (float *)malloc(2 * lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2 * lwork * sizeof(double));
%}
$TFD(cggevx_,zggevx_)(
&pbalanc,
&pjobvl,
&pjobvr,
&psens,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(B),
&$PRIV(__n_size),
$P(alpha),
$P(beta),
$P(VL),
&$PRIV(__m_size),
$P(VR),
&$PRIV(__p_size),
$P(ilo),
$P(ihi),
$P(lscale),
$P(rscale),
$P(abnrm),
$P(bbnrm),
$P(rconde),
$P(rcondv),
work,
&lwork,
rwork,
iwork,
bwork,
$P(info));
free(work);
}
free(rwork);
if ($sense())
free(bwork);
if ($sense() != 1)
free(iwork);
');
pp_addhdr(' static SV* fselect_func; PDL_Long fselect_wrapper(float *p) { dSP ;
int count;
long ret;
SV *pdl1;
HV *bless_stash;
pdl *pdl;
PDL_Long odims[1];
PDL_Long dims[] = {2,1};
pdl = PDL->pdlnew();
PDL->setdims (pdl, dims, 2);
pdl->datatype = PDL_F;
pdl->data = p;
pdl->state |= PDL_DONTTOUCHDATA | PDL_ALLOCATED;
bless_stash = gv_stashpv("PDL::Complex", 0);
ENTER ; SAVETMPS ; PUSHMARK(sp) ;
pdl1 = sv_newmortal();
PDL->SetSV_PDL(pdl1, pdl);
pdl1 = sv_bless(pdl1, bless_stash); /* bless in PDL::Complex */
XPUSHs(pdl1);
PUTBACK ;
count = perl_call_sv(fselect_func, G_SCALAR);
SPAGAIN;
if (count !=1)
croak("Error calling perl function\n");
// For pdl_free
odims[0] = 0;
PDL->setdims (pdl, odims, 0);
pdl->state &= ~ (PDL_ALLOCATED |PDL_DONTTOUCHDATA);
pdl->data=NULL;
ret = (long ) POPl ;
PUTBACK ; FREETMPS ; LEAVE ;
return ret;
}
static SV* dselect_func; PDL_Long dselect_wrapper(double *p) { dSP ;
int count;
long ret;
SV *pdl1;
HV *bless_stash;
pdl *pdl;
PDL_Long odims[1];
PDL_Long dims[] = {2,1};
pdl = PDL->pdlnew();
PDL->setdims (pdl, dims, 2);
pdl->datatype = PDL_D;
pdl->data = p;
pdl->state |= PDL_DONTTOUCHDATA | PDL_ALLOCATED;
bless_stash = gv_stashpv("PDL::Complex", 0);
ENTER ; SAVETMPS ; PUSHMARK(sp) ;
pdl1 = sv_newmortal();
PDL->SetSV_PDL(pdl1, pdl);
pdl1 = sv_bless(pdl1, bless_stash); /* bless in PDL::Complex */
XPUSHs(pdl1);
PUTBACK ;
count = perl_call_sv(dselect_func, G_SCALAR);
SPAGAIN;
if (count !=1)
croak("Error calling perl function\n");
// For pdl_free
odims[0] = 0;
PDL->setdims (pdl, odims, 0);
pdl->state &= ~ (PDL_ALLOCATED |PDL_DONTTOUCHDATA);
pdl->data=NULL;
ret = (long ) POPl ;
PUTBACK ; FREETMPS ; LEAVE ;
return ret;
}
');
pp_defc("gees", HandleBad => 0, Pars => '[io,phys]A(2,n,n); int jobvs(); int sort(); [o,phys]w(2,n); [o,phys]vs(2,p,p); int [o,phys]sdim(); int [o,phys]info()', OtherPars => "SV* select_func" , GenericTypes => [F,D], Code => generate_code '
char jvs = \'N\';
char psort = \'N\';
integer *bwork;
integer lwork = -1;
types(F) %{
extern int cgees_(char *jobvs, char *sort, L_fp select, integer *n,
float *a, integer *lda, integer *sdim, float *w,
float *vs, integer *ldvs, float *work,
integer *lwork, float *rwork, integer *bwork, integer *info);
float tmp_work[2];
float *rwork, *work;
rwork = (float *) malloc ($PRIV(__n_size) * sizeof (float));
fselect_func = $PRIV(select_func);
%}
types(D) %{
extern int zgees_(char *jobvs, char *sort, L_fp select, integer *n,
double *a, integer *lda, integer *sdim, double *w,
double *vs, integer *ldvs, double *work,
integer *lwork, double *rwork, integer *bwork, integer *info);
double *rwork, *work;
double tmp_work[2];
dselect_func = $PRIV(select_func);
rwork = (double *) malloc ($PRIV(__n_size) * sizeof (double));
%}
if ($jobvs())
jvs = \'V\';
if ($sort()){
psort = \'S\';
bwork = (integer * ) malloc ($PRIV(__n_size) * sizeof (integer));
}
types(F) %{
cgees_(
&jvs,
&psort,
fselect_wrapper,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(sdim),
$P(w),
$P(vs),
&$PRIV(__p_size),
&tmp_work[0],
&lwork,
rwork,
bwork,
$P(info));
%}
types(D) %{
zgees_(
&jvs,
&psort,
dselect_wrapper,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(sdim),
$P(w),
$P(vs),
&$PRIV(__p_size),
&tmp_work[0],
&lwork,
rwork,
bwork,
$P(info));
%}
lwork = (integer )tmp_work[0];
types(F) %{
work = (float *) malloc(2 * lwork * sizeof(float));
cgees_(
&jvs,
&psort,
fselect_wrapper,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(sdim),
$P(w),
$P(vs),
&$PRIV(__p_size),
work,
&lwork,
rwork,
bwork,
$P(info));
free(work);
%}
types(D) %{
work = (double *) malloc(2*lwork * sizeof(double));
zgees_(
&jvs,
&psort,
dselect_wrapper,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(sdim),
$P(w),
$P(vs),
&$PRIV(__p_size),
work,
&lwork,
rwork,
bwork,
$P(info));
free(work);
%}
if ($sort())
free(bwork);
free(rwork);
',
Doc=>'
Complex version of gees
select_func:
If sort = 1, select_func is used to select eigenvalues to sort
to the top left of the Schur form.
If sort = 0, select_func is not referenced.
An complex eigenvalue w is selected if
select_func(PDL::Complex(w)) is true;
Note that a selected complex eigenvalue may no longer
satisfy select_func(PDL::Complex(w)) = 1 after ordering, since
ordering may change the value of complex eigenvalues
(especially if the eigenvalue is ill-conditioned); in this
case info is set to N+2.
');
pp_defc("geesx", HandleBad => 0, Pars => '[io,phys]A(2,n,n); int jobvs(); int sort(); int sense(); [o,phys]w(2,n);[o,phys]vs(2,p,p); int [o,phys]sdim(); [o,phys]rconde();[o,phys]rcondv(); int [o,phys]info()', OtherPars => "SV* select_func" , GenericTypes => [F,D], Code => generate_code ' char jvs = \'N\'; char psort = \'N\'; integer *bwork; integer lwork = 0; char sens;
types(F) %{
extern int cgeesx_(char *jobvs, char *sort, L_fp select, char * sense,
integer *n, float *a, integer *lda, integer *sdim, float *w,
float *vs, integer *ldvs, float *rconde, float *rcondv,
float *work, integer *lwork, float *rwork,
integer *bwork, integer *info);
float *work, *rwork;
rwork = (float *) malloc ($PRIV(__n_size) * sizeof (float));
fselect_func = $PRIV(select_func);
%}
types(D) %{
extern int zgeesx_(char *jobvs, char *sort, L_fp select, char * sense,
integer *n, double *a, integer *lda, integer *sdim, double *w,
double *vs, integer *ldvs, double *rconde, double *rcondv,
double *work, integer *lwork, double *rwork,
integer *bwork, integer *info);
double *work, *rwork;
dselect_func = $PRIV(select_func);
rwork = (double *) malloc ($PRIV(__n_size) * sizeof (double));
%}
if ($jobvs())
jvs = \'V\';
if ($sort()){
psort = \'S\';
bwork = (integer * ) malloc ($PRIV(__n_size) * sizeof (integer));
}
switch ($sense())
{
case 1: sens = \'E\';
lwork = (integer ) ($PRIV(__n_size) * ($PRIV(__n_size)/2));
break;
case 2: sens = \'V\';
lwork = (integer ) ($PRIV(__n_size) * ($PRIV(__n_size)/2));
break;
case 3: sens = \'B\';
lwork = (integer ) ($PRIV(__n_size) * ($PRIV(__n_size)/2));
break;
default: sens = \'N\';
lwork = (integer ) ($PRIV(__n_size) * 2);
}
types(D) %{
work = (double * )malloc(2*lwork * sizeof (double));
%}
types(F) %{
work = (float * )malloc(2*lwork * sizeof (float));
%}
types(F) %{
cgeesx_(
&jvs,
&psort,
fselect_wrapper,
&sens,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(sdim),
$P(w),
$P(vs),
&$PRIV(__p_size),
$P(rconde),
$P(rcondv),
work,
&lwork,
rwork,
bwork,
$P(info));
%}
types(D) %{
zgeesx_(
&jvs,
&psort,
dselect_wrapper,
&sens,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(sdim),
$P(w),
$P(vs),
&$PRIV(__p_size),
$P(rconde),
$P(rcondv),
work,
&lwork,
rwork,
bwork,
$P(info));
%}
free(work);
if ($sort())
free(bwork);
free(rwork);
',
Doc => '
Complex version of geesx
select_func:
If sort = 1, select_func is used to select eigenvalues to sort
to the top left of the Schur form.
If sort = 0, select_func is not referenced.
An complex eigenvalue w is selected if
select_func(PDL::Complex(w)) is true;
Note that a selected complex eigenvalue may no longer
satisfy select_func(PDL::Complex(w)) = 1 after ordering, since
ordering may change the value of complex eigenvalues
(especially if the eigenvalue is ill-conditioned); in this
case info is set to N+2.
');
pp_addhdr(' static SV* fgselect_func; PDL_Long fgselect_wrapper(float *p, float *q) { dSP ;
int count;
long ret;
SV *svpdl1, *svpdl2;
HV *bless_stash;
pdl *pdl1, *pdl2;
PDL_Long odims[1];
PDL_Long dims[] = {2,1};
pdl1 = PDL->pdlnew();
PDL->setdims (pdl1, dims, 2);
pdl1->datatype = PDL_F;
pdl1->data = p;
pdl1->state |= PDL_DONTTOUCHDATA | PDL_ALLOCATED;
pdl2 = PDL->pdlnew();
PDL->setdims (pdl2, dims, 2);
pdl2->datatype = PDL_F;
pdl2->data = q;
pdl2->state |= PDL_DONTTOUCHDATA | PDL_ALLOCATED;
bless_stash = gv_stashpv("PDL::Complex", 0);
ENTER ; SAVETMPS ; PUSHMARK(sp) ;
svpdl1 = sv_newmortal();
PDL->SetSV_PDL(svpdl1, pdl1);
svpdl1 = sv_bless(svpdl1, bless_stash); /* bless in PDL::Complex */
svpdl2 = sv_newmortal();
PDL->SetSV_PDL(svpdl2, pdl2);
svpdl2 = sv_bless(svpdl2, bless_stash); /* bless in PDL::Complex */
XPUSHs(svpdl1);
XPUSHs(svpdl2);
PUTBACK ;
count = perl_call_sv(fgselect_func, G_SCALAR);
SPAGAIN;
if (count !=1)
croak("Error calling perl function\n");
// For pdl_free
odims[0] = 0;
PDL->setdims (pdl1, odims, 0);
pdl1->state &= ~ (PDL_ALLOCATED |PDL_DONTTOUCHDATA);
pdl1->data=NULL;
PDL->setdims (pdl2, odims, 0);
pdl1->state &= ~ (PDL_ALLOCATED |PDL_DONTTOUCHDATA);
pdl1->data=NULL;
ret = (long ) POPl ;
PUTBACK ; FREETMPS ; LEAVE ;
return ret;
}
static SV* dgselect_func; PDL_Long dgselect_wrapper(double *p, double *q) { dSP ;
int count;
long ret;
SV *svpdl1, *svpdl2;
HV *bless_stash;
pdl *pdl1, *pdl2;
PDL_Long odims[1];
PDL_Long dims[] = {2,1};
pdl1 = PDL->pdlnew();
PDL->setdims (pdl1, dims, 2);
pdl1->datatype = PDL_D;
pdl1->data = p;
pdl1->state |= PDL_DONTTOUCHDATA | PDL_ALLOCATED;
pdl2 = PDL->pdlnew();
PDL->setdims (pdl2, dims, 2);
pdl2->datatype = PDL_D;
pdl2->data = q;
pdl2->state |= PDL_DONTTOUCHDATA | PDL_ALLOCATED;
bless_stash = gv_stashpv("PDL::Complex", 0);
ENTER ; SAVETMPS ; PUSHMARK(sp) ;
svpdl1 = sv_newmortal();
PDL->SetSV_PDL(svpdl1, pdl1);
svpdl1 = sv_bless(svpdl1, bless_stash); /* bless in PDL::Complex */
svpdl2 = sv_newmortal();
PDL->SetSV_PDL(svpdl2, pdl2);
svpdl2 = sv_bless(svpdl2, bless_stash); /* bless in PDL::Complex */
XPUSHs(svpdl1);
XPUSHs(svpdl2);
PUTBACK ;
count = perl_call_sv(dgselect_func, G_SCALAR);
SPAGAIN;
if (count !=1)
croak("Error calling perl function\n");
// For pdl_free
odims[0] = 0;
PDL->setdims (pdl1, odims, 0);
pdl1->state &= ~ (PDL_ALLOCATED |PDL_DONTTOUCHDATA);
pdl1->data=NULL;
PDL->setdims (pdl2, odims, 0);
pdl2->state &= ~ (PDL_ALLOCATED |PDL_DONTTOUCHDATA);
pdl2->data=NULL;
ret = (long ) POPl ;
PUTBACK ; FREETMPS ; LEAVE ;
return ret;
}
');
pp_defc("gges", HandleBad => 0, Pars => '[io,phys]A(2,n,n); int jobvsl();int jobvsr();int sort();[io,phys]B(2,n,n);[o,phys]alpha(2,n);[o,phys]beta(2,n);[o,phys]VSL(2,m,m);[o,phys]VSR(2,p,p);int [o,phys]sdim();int [o,phys]info()', OtherPars => "SV* select_func" , GenericTypes => [F,D], Code => generate_code '
integer lwork = -1;
char pjobvsl = \'N\', pjobvsr = \'N\', psort = \'N\';
integer *bwork;
types(F) %{
extern int cgges_(char *jobvsl, char *jobvsr, char *sort, L_fp
delctg, integer *n, float *a, integer *lda, float *b,
integer *ldb, integer *sdim, float *alpha,
float *beta, float *vsl, integer *ldvsl, float *vsr,
integer *ldvsr, float *work, integer *lwork, float *rwork,
logical *bwork, integer *info);
float tmp_work[2], *rwork;
fgselect_func = $PRIV(select_func);
rwork = (float *)malloc(8 * $SIZE(n) * sizeof(float));
%}
types(D) %{
extern int zgges_(char *jobvsl, char *jobvsr, char *sort, L_fp
delctg, integer *n, double *a, integer *lda, double *b,
integer *ldb, integer *sdim, double *alpha,
double *beta, double *vsl, integer *ldvsl, double *vsr,
integer *ldvsr, double *work, integer *lwork, double *rwork,
logical *bwork, integer *info);
double tmp_work[2], *rwork;
dgselect_func = $PRIV(select_func);
rwork = (double *)malloc(8 * $SIZE(n) * sizeof(double));
%}
if ($jobvsl())
pjobvsl = \'V\';
if ($jobvsr())
pjobvsr = \'V\';
if ($sort()){
psort = \'S\';
bwork = (integer *)malloc($PRIV(__n_size) * sizeof(integer));
}
types(F) %{
cgges_(
&pjobvsl,
&pjobvsr,
&psort,
fgselect_wrapper,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(B),
&$PRIV(__n_size),
$P(sdim),
$P(alpha),
$P(beta),
$P(VSL),
&$PRIV(__m_size),
$P(VSR),
&$PRIV(__p_size),
&tmp_work[0],
&lwork,
rwork,
bwork,
$P(info));
%}
types(D) %{
zgges_(
&pjobvsl,
&pjobvsr,
&psort,
dgselect_wrapper,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(B),
&$PRIV(__n_size),
$P(sdim),
$P(alpha),
$P(beta),
$P(VSL),
&$PRIV(__m_size),
$P(VSR),
&$PRIV(__p_size),
&tmp_work[0],
&lwork,
rwork,
bwork,
$P(info));
%}
lwork = (integer )tmp_work[0];
{
types(F) %{
float *work = (float *)malloc(2 * lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2 * lwork * sizeof(double));
%}
types(F) %{
cgges_(
&pjobvsl,
&pjobvsr,
&psort,
fgselect_wrapper,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(B),
&$PRIV(__n_size),
$P(sdim),
$P(alpha),
$P(beta),
$P(VSL),
&$PRIV(__m_size),
$P(VSR),
&$PRIV(__p_size),
work,
&lwork,
rwork,
bwork,
$P(info));
%}
types(D) %{
zgges_(
&pjobvsl,
&pjobvsr,
&psort,
dgselect_wrapper,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(B),
&$PRIV(__n_size),
$P(sdim),
$P(alpha),
$P(beta),
$P(VSL),
&$PRIV(__m_size),
$P(VSR),
&$PRIV(__p_size),
work,
&lwork,
rwork,
bwork,
$P(info));
%}
free(work);
}
if ($sort())
free (bwork);
free(rwork);
', Doc=>'
Complex version of ggees
select_func:
If sort = 1, select_func is used to select eigenvalues to sort
to the top left of the Schur form.
If sort = 0, select_func is not referenced.
An eigenvalue w = w/beta is selected if
select_func(PDL::Complex(w), PDL::Complex(beta)) is true;
Note that a selected complex eigenvalue may no longer
satisfy select_func(PDL::Complex(w),PDL::Complex(beta)) = 1 after ordering, since
ordering may change the value of complex eigenvalues
(especially if the eigenvalue is ill-conditioned); in this
case info is set to N+2.
');
pp_defc("ggesx", HandleBad => 0, Pars => '[io,phys]A(2,n,n); int jobvsl();int jobvsr();int sort();int sense();[io,phys]B(2,n,n);[o,phys]alpha(2,n);[o,phys]beta(2,n);[o,phys]VSL(2,m,m);[o,phys]VSR(2,p,p);int [o,phys]sdim();[o,phys]rconde(q);[o,phys]rcondv(r);int [o,phys]info()', OtherPars => "SV* select_func" , GenericTypes => [F,D], Code => generate_code '
integer lwork, maxwrk;
integer liwork = 1;
integer minwrk = 1;
static integer c__0 = 0;
static integer c__1 = 1;
static integer c_n1 = -1;
char pjobvsl = \'N\';
char pjobvsr = \'N\';
char psort = \'N\';
char psens = \'N\';
integer *bwork;
integer *iwork;
extern integer ilaenv_(integer *ispec, char *name__, char *opts, integer *n1,
integer *n2, integer *n3, integer *n4, ftnlen name_len, ftnlen
opts_len);
types(F) %{
extern int cggesx_(char *jobvsl, char *jobvsr, char *sort, L_fp
delctg, char *sense, integer *n, float *a, integer *lda, float *b,
integer *ldb, integer *sdim, float *alpha,
float *beta, float *vsl, integer *ldvsl, float *vsr,
integer *ldvsr, float *rconde, float *rcondv, float *work,
integer *lwork, float *rwork, integer *iwork, integer *liwork, logical *bwork,
integer *info);
float *rwork = (float *)malloc(8 * $SIZE(n) * sizeof(float));
fgselect_func = $PRIV(select_func);
%}
types(D) %{
extern int zggesx_(char *jobvsl, char *jobvsr, char *sort, L_fp
delctg, char *sense, integer *n, double *a, integer *lda, double *b,
integer *ldb, integer *sdim, double *alpha,
double *beta, double *vsl, integer *ldvsl, double *vsr,
integer *ldvsr, double *rconde, double *rcondv, double *work,
integer *lwork, double *rwork, integer *iwork, integer *liwork, logical *bwork,
integer *info);
double *rwork = (double *)malloc(8 * $SIZE(n) * sizeof(double));
dgselect_func = $PRIV(select_func);
%}
if ($jobvsr())
pjobvsr = \'V\';
if ($sort()){
psort = \'S\';
bwork = (integer *)malloc($PRIV(__n_size) * sizeof(integer));
}
switch ($sense())
{
case 1: psens = \'E\';
break;
case 2: psens = \'V\';
break;
case 3: psens = \'B\';
break;
default: psens = \'N\';
}
// if (!$sense())
// liwork = 1;
// else
// {
liwork = $SIZE(n) + 2;
iwork = (integer *)malloc(liwork * sizeof(integer));
// }
// Code modified from Lapack
// TODO other shur form above
// The actual updated release (clapack 09/20/2000) do not allow
// the workspace query. See release notes of Lapack
// for this feature.
minwrk = $SIZE(n) << 1;
maxwrk = $SIZE(n) + $SIZE(n) * ilaenv_(&c__1, "ZGEQRF", " ", &$PRIV(__n_size), &c__1,
&$PRIV(__n_size), &c__0, (ftnlen)6, (ftnlen)1);
if ($jobvsl())
{
integer i__1 = maxwrk;
integer i__2 = $SIZE(n) + $SIZE(n) * ilaenv_(&c__1, "ZUNGQR"
, " ", &$PRIV(__n_size), &c__1, &$PRIV(__n_size), &c_n1, (ftnlen)6, (ftnlen)1);
maxwrk = max(i__1,i__2);
pjobvsl = \'V\';
}
lwork = max(maxwrk,minwrk);
{
types(F) %{
float *work = (float *)malloc( 2 * lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2 * lwork * sizeof(double));
%}
types(F) %{
cggesx_(
&pjobvsl,
&pjobvsr,
&psort,
fgselect_wrapper,
&psens,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(B),
&$PRIV(__n_size),
$P(sdim),
$P(alpha),
$P(beta),
$P(VSL),
&$PRIV(__m_size),
$P(VSR),
&$PRIV(__p_size),
$P(rconde),
$P(rcondv),
work,
&lwork,
rwork,
iwork,
&liwork,
bwork,
$P(info));
%}
types(D) %{
zggesx_(
&pjobvsl,
&pjobvsr,
&psort,
dgselect_wrapper,
&psens,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(B),
&$PRIV(__n_size),
$P(sdim),
$P(alpha),
$P(beta),
$P(VSL),
&$PRIV(__m_size),
$P(VSR),
&$PRIV(__p_size),
$P(rconde),
$P(rcondv),
work,
&lwork,
rwork,
iwork,
&liwork,
bwork,
$P(info));
%}
free(work);
}
if ($sort())
free (bwork);
// if ($sense())
free(iwork);
free(rwork);
',
Doc=>'
Complex version of ggeesx
select_func:
If sort = 1, select_func is used to select eigenvalues to sort
to the top left of the Schur form.
If sort = 0, select_func is not referenced.
An eigenvalue w = w/beta is selected if
select_func(PDL::Complex(w), PDL::Complex(beta)) is true;
Note that a selected complex eigenvalue may no longer
satisfy select_func(PDL::Complex(w),PDL::Complex(beta)) = 1 after ordering, since
ordering may change the value of complex eigenvalues
(especially if the eigenvalue is ill-conditioned); in this
case info is set to N+3.
');
pp_defc("heev", HandleBad => 0, Pars => '[io,phys]A(2,n,n); int jobz(); int uplo(); [o,phys]w(n); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code ' char jz = \'N\'; char puplo = \'U\'; integer lwork = -1; types(F) %{ extern int cheev_(char *jobz, char *uplo, integer *n, float *a, integer *lda, float *w, float *work, integer *lwork, float *rwork, integer *info); float *rwork; float tmp_work[2]; rwork = (float *) malloc ((3*$PRIV(__n_size)-2) * sizeof(float)); %} types(D) %{ extern int zheev_(char *jobz, char *uplo, integer *n, double *a, integer *lda, double *w, double *work, integer *lwork, double *rwork, integer *info); double *rwork; double tmp_work[2]; rwork = (double *) malloc ((3*$PRIV(__n_size)-2) * sizeof(double)); %}
if ($jobz())
jz = \'V\';
if ($uplo())
puplo = \'L\';
$TFD(cheev_,zheev_)(
&jz,
&puplo,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(w),
&tmp_work[0],
&lwork,
rwork,
$P(info));
lwork = (integer )tmp_work[0];
{
types(F) %{
float *work = (float *)malloc(2*lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2*lwork * sizeof(double));
%}
$TFD(cheev_,zheev_)(
&jz,
&puplo,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(w),
work,
&lwork,
rwork,
$P(info));
free(work);
}
free(rwork);
',
Doc=>'
Complex version of syev for Hermitian matrix
');
pp_defc("heevd", HandleBad => 0, Pars => '[io,phys]A(2,n,n); int jobz(); int uplo(); [o,phys]w(n); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code '
char jz = \'N\';
char puplo = \'U\';
integer lwork = -1;
integer lrwork, liwork;
integer tmpi_work;
integer *iwork;
types(F) %{
extern int cheevd_(char *jobz, char *uplo, integer *n, float *a,
integer *lda, float *w, float *work, integer *lwork, float *rwork, integer *lrwork,
integer *iwork, integer *liwork, integer *info);
float tmp_work[2];
float tmpr_work;
%}
types(D) %{
extern int zheevd_(char *jobz, char *uplo, integer *n, double *a,
integer *lda, double *w, double *work, integer *lwork, double *rwork, integer *lrwork,
integer *iwork, integer *liwork, integer *info);
double tmp_work[2];
double tmpr_work;
%}
if ($jobz())
jz = \'V\';
if ($uplo())
puplo = \'L\';
$TFD(cheevd_,zheevd_)(
&jz,
&puplo,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(w),
&tmp_work[0],
&lwork,
&tmpr_work,
&lwork,
&tmpi_work,
&lwork,
$P(info));
lwork = (integer )tmp_work[0];
lrwork = (integer )tmpr_work;
liwork = (integer )tmpi_work;
iwork = (integer *)malloc(liwork * sizeof(integer));
{
types(F) %{
float *work = (float *)malloc(2*lwork * sizeof(float));
float *rwork = (float *)malloc(lrwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2*lwork * sizeof(double));
double *rwork = (double *)malloc(lrwork * sizeof(double));
%}
$TFD(cheevd_,zheevd_)(
&jz,
&puplo,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(w),
work,
&lwork,
rwork,
&lrwork,
iwork,
&liwork,
$P(info));
free(rwork);
free(work);
}
free(iwork);
',
Doc=>'
Complex version of syevd for Hermitian matrix
');
pp_defc("heevx", HandleBad => 0, Pars => '[phys]A(2,n,n); int jobz(); int range(); int uplo(); [phys]vl(); [phys]vu(); int [phys]il(); int [phys]iu(); [phys]abstol(); int [o,phys]m(); [o,phys]w(n); [o,phys]z(2,p,q);int [o,phys]ifail(r); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code ' char jz = \'N\'; char puplo = \'U\'; char prange = \'A\'; integer lwork = -1; integer *iwork;
types(F) %{
extern int cheevx_(char *jobz, char *range, char *uplo, integer *n,
float *a, integer *lda, float *vl, float *vu, integer *
il, integer *iu, float *abstol, integer *m, float *w,
float *z__, integer *ldz, float *work, integer *lwork,
float *rwork, integer *iwork, integer *ifail, integer *info);
float *rwork;
float tmp_work[2];
rwork = (float *)malloc(7 * $SIZE(n) * sizeof(float));
%}
types(D) %{
extern int zheevx_(char *jobz, char *range, char *uplo, integer *n,
double *a, integer *lda, double *vl, double *vu, integer *
il, integer *iu, double *abstol, integer *m, double *w,
double *z__, integer *ldz, double *work, integer *lwork,
double *rwork, integer *iwork, integer *ifail, integer *info);
double *rwork;
double tmp_work[2];
rwork = (double *)malloc(7 * $SIZE(n) * sizeof(double));
%}
if ($jobz())
jz = \'V\';
if ($uplo())
puplo = \'L\';
switch ($range())
{
case 1: prange = \'V\';
break;
case 2: prange = \'I\';
break;
default: prange = \'A\';
}
iwork = (integer *)malloc(5 * $SIZE (n) * sizeof(integer));
$TFD(cheevx_,zheevx_)(
&jz,
&prange,
&puplo,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(vl),
$P(vu),
$P(il),
$P(iu),
$P(abstol),
$P(m),
$P(w),
$P(z),
&$PRIV(__p_size),
&tmp_work[0],
&lwork,
rwork,
iwork,
$P(ifail),
$P(info));
lwork = (integer )tmp_work[0];
{
types(F) %{
float *work = (float *)malloc(2* lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2 * lwork * sizeof(double));
%}
$TFD(cheevx_,zheevx_)(
&jz,
&prange,
&puplo,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(vl),
$P(vu),
$P(il),
$P(iu),
$P(abstol),
$P(m),
$P(w),
$P(z),
&$PRIV(__p_size),
work,
&lwork,
rwork,
iwork,
$P(ifail),
$P(info));
free(work);
}
free(iwork);
free(rwork);
',
Doc=>'
Complex version of syevx for Hermitian matrix
');
pp_defc("heevr", HandleBad => 0, Pars => '[phys]A(2,n,n); int jobz(); int range(); int uplo(); [phys]vl(); [phys]vu(); int [phys]il(); int [phys]iu(); [phys]abstol(); int [o,phys]m(); [o,phys]w(n); [o,phys]z(2,p,q);int [o,phys]isuppz(r); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code ' char jz = \'N\'; char puplo = \'U\'; char prange = \'A\'; integer lwork = -1; integer liwork,lrwork; integer tmpi_work;
types(F) %{
extern int cheevr_(char *jobz, char *range, char *uplo, integer *n,
float *a, integer *lda, float *vl, float *vu, integer *
il, integer *iu, float *abstol, integer *m, float *w,
float *z__, integer *ldz, integer *isuppz, float *work, integer *lwork, float *rwork, integer *lrwork,
integer *iwork, integer *liwork, integer *info);
float tmp_work[2];
float tmpr_work;
%}
types(D) %{
extern int zheevr_(char *jobz, char *range, char *uplo, integer *n,
double *a, integer *lda, double *vl, double *vu, integer *
il, integer *iu, double *abstol, integer *m, double *w,
double *z__, integer *ldz, integer *isuppz, double *work, integer *lwork, double *rwork, integer *lrwork,
integer *iwork, integer *liwork, integer *info);
double tmp_work[2];
double tmpr_work;
%}
if ($jobz())
jz = \'V\';
if ($uplo())
puplo = \'L\';
switch ($range())
{
case 1: prange = \'V\';
break;
case 2: prange = \'I\';
break;
default: prange = \'A\';
}
$TFD(cheevr_,zheevr_)(
&jz,
&prange,
&puplo,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(vl),
$P(vu),
$P(il),
$P(iu),
$P(abstol),
$P(m),
$P(w),
$P(z),
&$PRIV(__p_size),
$P(isuppz),
&tmp_work[0],
&lwork,
&tmpr_work,
&lwork,
&tmpi_work,
&lwork,
$P(info));
lwork = (integer )tmp_work[0];
lrwork = (integer )tmpr_work;
liwork = (integer )tmpi_work;
{
types(F) %{
float *work = (float *)malloc(2* lwork * sizeof(float));
float *rwork = (float *)malloc(lrwork * sizeof(float));
integer *iwork = (integer *)malloc(liwork * sizeof(integer));
%}
types(D) %{
double *work = (double *)malloc(2 * lwork * sizeof(double));
double *rwork = (double *)malloc(lrwork * sizeof(double));
integer *iwork = (integer *)malloc(liwork * sizeof(integer));
%}
$TFD(cheevr_,zheevr_)(
&jz,
&prange,
&puplo,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(vl),
$P(vu),
$P(il),
$P(iu),
$P(abstol),
$P(m),
$P(w),
$P(z),
&$PRIV(__p_size),
$P(isuppz),
work,
&lwork,
rwork,
&lrwork,
iwork,
&liwork,
$P(info));
free(work);
free(iwork);
free(rwork);
}
', Doc=>'
Complex version of syevr for Hermitian matrix
');
pp_defc("hegv", HandleBad => 0, Pars => '[io,phys]A(2,n,n);int [phys]itype();int jobz(); int uplo();[io,phys]B(2,n,n);[o,phys]w(n); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code '
char jz = \'N\';
char puplo = \'U\';
integer lwork = -1;
types(F) %{
extern int chegv_(integer *itype, char *jobz, char *uplo, integer *
n, float *a, integer *lda, float *b, integer *ldb,
float *w, float *work, integer *lwork, float *rwork, integer *info);
float tmp_work[2], *rwork;
rwork = (float *) malloc( (3 * $SIZE(n) - 2 ) * sizeof(float));
%}
types(D) %{
extern int zhegv_(integer *itype, char *jobz, char *uplo, integer *
n, double *a, integer *lda, double *b, integer *ldb,
double *w, double *work, integer *lwork, double *rwork, integer *info);
double tmp_work[2], *rwork;
rwork = (double *) malloc( (3 * $SIZE(n) - 2 ) * sizeof(double));
%}
if ($jobz())
jz = \'V\';
if ($uplo())
puplo = \'L\';
$TFD(chegv_,zhegv_)(
$P(itype),
&jz,
&puplo,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(B),
&$PRIV(__n_size),
$P(w),
&tmp_work[0],
&lwork,
rwork,
$P(info));
lwork = (integer )tmp_work[0];
{
types(F) %{
float *work = (float *)malloc(2 * lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2 * lwork * sizeof(double));
%}
$TFD(chegv_,zhegv_)(
$P(itype),
&jz,
&puplo,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(B),
&$PRIV(__n_size),
$P(w),
work,
&lwork,
rwork,
$P(info));
free(work);
}
free(rwork);
',
Doc=>'
Complex version of sygv for Hermitian matrix
');
pp_defc("hegvd", HandleBad => 0, Pars => '[io,phys]A(2,n,n);int [phys]itype();int jobz(); int uplo();[io,phys]B(2,n,n);[o,phys]w(n); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code '
char jz = \'N\';
char puplo = \'U\';
integer lwork = -1;
integer liwork = -1;
integer lrwork = -1;
integer *iwork;
integer tmp_iwork;
types(F) %{
extern int chegvd_(integer *itype, char *jobz, char *uplo, integer *
n, float *a, integer *lda, float *b, integer *ldb,
float *w, float *work, integer *lwork, float *rwork, integer *lrwork,
integer *iwork, integer *liwork, integer *info);
float tmp_work[2], tmp_rwork;
%}
types(D) %{
extern int zhegvd_(integer *itype, char *jobz, char *uplo, integer *
n, double *a, integer *lda, double *b, integer *ldb,
double *w, double *work, integer *lwork, double *rwork, integer *lrwork,
integer *iwork, integer *liwork, integer *info);
double tmp_work[2], tmp_rwork;
%}
if ($jobz())
jz = \'V\';
if ($uplo())
puplo = \'L\';
$TFD(chegvd_,zhegvd_)(
$P(itype),
&jz,
&puplo,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(B),
&$PRIV(__n_size),
$P(w),
&tmp_work[0],
&lwork,
&tmp_rwork,
&lrwork,
&tmp_iwork,
&liwork,
$P(info));
lwork = (integer )tmp_work[0];
lrwork = (integer )tmp_rwork;
liwork = (integer )tmp_iwork;
iwork = (integer *)malloc(liwork * sizeof(integer));
{
types(F) %{
float *work = (float *)malloc(2 * lwork * sizeof(float));
float *rwork = (float *)malloc(lrwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2 * lwork * sizeof(double));
double *rwork = (double *)malloc(lrwork * sizeof(double));
%}
$TFD(chegvd_,zhegvd_)(
$P(itype),
&jz,
&puplo,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(B),
&$PRIV(__n_size),
$P(w),
work,
&lwork,
rwork,
&lrwork,
iwork,
&liwork,
$P(info));
free(work);
free(rwork);
}
free(iwork);
',
Doc=>'
Complex version of sygvd for Hermitian matrix
');
pp_defc("hegvx", HandleBad => 0, Pars => '[io,phys]A(2,n,n);int [phys]itype();int jobz();int range(); int uplo();[io,phys]B(2,n,n);[phys]vl();[phys]vu();int [phys]il();int [phys]iu();[phys]abstol();int [o,phys]m();[o,phys]w(n); [o,phys]Z(2,p,q);int [o,phys]ifail(r);int [o,phys]info()', GenericTypes => [F,D], Code => generate_code ' char jz = \'N\'; char puplo = \'U\'; char prange; integer lwork = -1; integer *iwork;
types(F) %{
extern int chegvx_(integer *itype, char *jobz, char *range, char *
uplo, integer *n, float *a, integer *lda, float *b, integer
*ldb, float *vl, float *vu, integer *il, integer *iu,
float *abstol, integer *m, float *w, float *z__,
integer *ldz, float *work, integer *lwork, float *rwork,
integer *iwork, integer *ifail, integer *info);
float tmp_work[2], *rwork;
rwork = (float *)malloc(7 * $SIZE(n) * sizeof(float));
%}
types(D) %{
extern int zhegvx_(integer *itype, char *jobz, char *range, char *
uplo, integer *n, double *a, integer *lda, double *b, integer
*ldb, double *vl, double *vu, integer *il, integer *iu,
double *abstol, integer *m, double *w, double *z__,
integer *ldz, double *work, integer *lwork, double *rwork,
integer *iwork, integer *ifail, integer *info);
double tmp_work[2], *rwork;
rwork = (double *)malloc(7 * $SIZE(n) * sizeof(double));
%}
if ($jobz())
jz = \'V\';
if ($uplo())
puplo = \'L\';
switch ($range())
{
case 1: prange = \'V\';
break;
case 2: prange = \'I\';
break;
default: prange = \'A\';
}
iwork = (integer *)malloc((5 * $SIZE(n)) * sizeof(integer));
$TFD(chegvx_,zhegvx_)(
$P(itype),
&jz,
&prange,
&puplo,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(B),
&$PRIV(__n_size),
$P(vl),
$P(vu),
$P(il),
$P(iu),
$P(abstol),
$P(m),
$P(w),
$P(Z),
&$PRIV(__p_size),
&tmp_work[0],
&lwork,
rwork,
iwork,
$P(ifail),
$P(info));
lwork = (integer )tmp_work[0];
{
types(F) %{
float *work = (float *)malloc( 2 * lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc( 2 * lwork * sizeof(double));
%}
$TFD(chegvx_,zhegvx_)(
$P(itype),
&jz,
&prange,
&puplo,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(B),
&$PRIV(__n_size),
$P(vl),
$P(vu),
$P(il),
$P(iu),
$P(abstol),
$P(m),
$P(w),
$P(Z),
&$PRIV(__p_size),
work,
&lwork,
rwork,
iwork,
$P(ifail),
$P(info));
free(work);
}
free(iwork);
free(rwork);
',
Doc=>'
Complex version of sygvx for Hermitian matrix
');
pp_defc("gesv", HandleBad => 0, Pars => '[io,phys]A(2,n,n); [io,phys]B(2,n,m); int [o,phys]ipiv(n); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code '
types(F) %{
extern int cgesv_(integer *n, integer *nrhs, float *a, integer
*lda, integer *ipiv, float *b, integer *ldb, integer *info);
%}
types(D) %{
extern int zgesv_(integer *n, integer *nrhs, double *a, integer
*lda, integer *ipiv, double *b, integer *ldb, integer *info);
%}
$TFD(cgesv_,zgesv_)(
&$PRIV(__n_size),
&$PRIV(__m_size),
$P(A),
&$PRIV(__n_size),
$P(ipiv),
$P(B),
&$PRIV(__n_size),
$P(info));
');
pp_defc("gesvx",
HandleBad => 0,
Pars => '[io,phys]A(2,n,n); int trans(); int fact(); [io,phys]B(2,n,m); [io,phys]af(2,n,n); int [io,phys]ipiv(n); int [io]equed(); [io,phys]r(n); [io,phys]c(n); [o,phys]X(2,n,m); [o,phys]rcond(); [o,phys]ferr(m); [o,phys]berr(m); [o,phys]rpvgrw(); int [o,phys]info()',
GenericTypes => [F,D],
Code => generate_code '
char ptrans, pfact, pequed;
types(F) %{
extern int cgesvx_(char *fact, char *trans, integer *n, integer *
nrhs, float *a, integer *lda, float *af, integer *ldaf,
integer *ipiv, char *equed, float *r__, float *c__,
float *b, integer *ldb, float *x, integer *ldx, float *
rcond, float *ferr, float *berr, float *work, float *
rwork, integer *info);
float *work = (float *) malloc(4 * $PRIV(__n_size) * sizeof(float));
float *rwork = (float *) malloc(4 * $PRIV(__n_size) * sizeof(float));
%}
types(D) %{
extern int zgesvx_(char *fact, char *trans, integer *n, integer *
nrhs, double *a, integer *lda, double *af, integer *ldaf,
integer *ipiv, char *equed, double *r__, double *c__,
double *b, integer *ldb, double *x, integer *ldx, double *
rcond, double *ferr, double *berr, double *work, double *
rwork, integer *info);
double *work = (double *) malloc(4 * $PRIV(__n_size) * sizeof(double));
double *rwork = (double *) malloc(4 * $PRIV(__n_size) * sizeof(double));
%}
switch ($trans())
{
case 1: ptrans = \'T\';
break;
case 2: ptrans = \'C\';
break;
default: ptrans = \'N\';
}
switch ($fact())
{
case 1: pfact = \'N\';
break;
case 2: pfact = \'E\';
break;
default: pfact = \'F\';
}
switch ($equed())
{
case 1: pequed = \'R\';
break;
case 2: pequed = \'C\';
break;
case 3: pequed = \'B\';
break;
default: pequed = \'N\';
}
$TFD(cgesvx_,zgesvx_)(
&pfact,
&ptrans,
&$PRIV(__n_size),
&$PRIV(__m_size),
$P(A),
&$PRIV(__n_size),
$P(af),
&$PRIV(__n_size),
$P(ipiv),
&pequed,
$P(r),
$P(c),
$P(B),
&$PRIV(__n_size),
$P(X),
&$PRIV(__n_size),
$P(rcond),
$P(ferr),
$P(berr),
work,
rwork,
$P(info));
free(work);
free(rwork);
switch (pequed)
{
case \'R\': $equed() = 1;
break;
case \'C\': $equed() = 2;
break;
case \'B\': $equed() = 3;
break;
default: $equed()= 0;
}
$rpvgrw() = rwork[0];
',
Doc => '
Complex version of gesvx.
trans: Specifies the form of the system of equations:
= 0: A * X = B (No transpose)
= 1: A\' * X = B (Transpose)
= 2: A**H * X = B (Conjugate transpose)
');
pp_defc("sysv", HandleBad => 0, Pars => '[io,phys]A(2,n,n); int uplo(); [io,phys]B(2,n,m); int [o,phys]ipiv(n); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code '
char puplo = \'U\';
integer lwork = -1;
types(F) %{
extern int csysv_(char *uplo, integer *n, integer *nrhs, float
*a, integer *lda, integer *ipiv, float *b, integer *ldb,
float *work, integer *lwork, integer *info);
float tmp_work[2];
%}
types(D) %{
extern int zsysv_(char *uplo, integer *n, integer *nrhs, double
*a, integer *lda, integer *ipiv, double *b, integer *ldb,
double *work, integer *lwork, integer *info);
double tmp_work[2];
%}
if ($uplo())
puplo = \'L\';
$TFD(csysv_,zsysv_)(
&puplo,
&$PRIV(__n_size),
&$PRIV(__m_size),
$P(A),
&$PRIV(__n_size),
$P(ipiv),
$P(B),
&$PRIV(__n_size),
&tmp_work[0],
&lwork,
$P(info));
lwork = (integer )tmp_work[0];
{
types(F) %{
float *work = (float *)malloc(2*lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2*lwork * sizeof(double));
%}
$TFD(csysv_,zsysv_)(
&puplo,
&$PRIV(__n_size),
&$PRIV(__m_size),
$P(A),
&$PRIV(__n_size),
$P(ipiv),
$P(B),
&$PRIV(__n_size),
work,
&lwork,
$P(info));
}
');
pp_defc("sysvx", HandleBad => 0, Pars => '[phys]A(2,n,n); int uplo(); int fact(); [phys]B(2,n,m); [io,phys]af(2,n,n); int [io,phys]ipiv(n); [o,phys]X(2,n,m); [o,phys]rcond(); [o,phys]ferr(m); [o,phys]berr(m); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code '
char pfact = \'N\';
char puplo = \'U\';
integer lwork = -1;
types(F) %{
extern int csysvx_(char *fact, char *uplo, integer *n, integer *
nrhs, float *a, integer *lda, float *af, integer *ldaf,
integer *ipiv, float *b, integer *ldb, float *x, integer *
ldx, float *rcond, float *ferr, float *berr,
float *work, integer *lwork, float *rwork, integer *info);
float *rwork = (float * )malloc ($PRIV(__n_size)* sizeof (float));
float tmp_work[2];
%}
types(D) %{
extern int zsysvx_(char *fact, char *uplo, integer *n, integer *
nrhs, double *a, integer *lda, double *af, integer *ldaf,
integer *ipiv, double *b, integer *ldb, double *x, integer *
ldx, double *rcond, double *ferr, double *berr,
double *work, integer *lwork, double *rwork, integer *info);
double *rwork = (double * )malloc ($PRIV(__n_size)* sizeof (double));
double tmp_work[2];
%}
if($fact())
pfact = \'F\';
if ($uplo())
puplo = \'L\';
$TFD(csysvx_,zsysvx_)(
&pfact,
&puplo,
&$PRIV(__n_size),
&$PRIV(__m_size),
$P(A),
&$PRIV(__n_size),
$P(af),
&$PRIV(__n_size),
$P(ipiv),
$P(B),
&$PRIV(__n_size),
$P(X),
&$PRIV(__n_size),
$P(rcond),
$P(ferr),
$P(berr),
&tmp_work[0],
&lwork,
rwork,
$P(info));
lwork = (integer )tmp_work[0];
{
types(F) %{
float *work = (float *)malloc(2*lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2*lwork * sizeof(double));
%}
$TFD(csysvx_,zsysvx_)(
&pfact,
&puplo,
&$PRIV(__n_size),
&$PRIV(__m_size),
$P(A),
&$PRIV(__n_size),
$P(af),
&$PRIV(__n_size),
$P(ipiv),
$P(B),
&$PRIV(__n_size),
$P(X),
&$PRIV(__n_size),
$P(rcond),
$P(ferr),
$P(berr),
work,
&lwork,
rwork,
$P(info));
free(work);
}
free(rwork);
');
pp_defc("hesv", HandleBad => 0, Pars => '[io,phys]A(2,n,n); int uplo(); [io,phys]B(2,n,m); int [o,phys]ipiv(n); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code '
char puplo = \'U\';
integer lwork = -1;
types(F) %{
extern int chesv_(char *uplo, integer *n, integer *nrhs, float
*a, integer *lda, integer *ipiv, float *b, integer *ldb,
float *work, integer *lwork, integer *info);
float tmp_work[2];
%}
types(D) %{
extern int zhesv_(char *uplo, integer *n, integer *nrhs, double
*a, integer *lda, integer *ipiv, double *b, integer *ldb,
double *work, integer *lwork, integer *info);
double tmp_work[2];
%}
if ($uplo())
puplo = \'L\';
$TFD(chesv_,zhesv_)(
&puplo,
&$PRIV(__n_size),
&$PRIV(__m_size),
$P(A),
&$PRIV(__n_size),
$P(ipiv),
$P(B),
&$PRIV(__n_size),
&tmp_work[0],
&lwork,
$P(info));
lwork = (integer )tmp_work[0];
{
types(F) %{
float *work = (float *)malloc(2*lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2*lwork * sizeof(double));
%}
$TFD(chesv_,zhesv_)(
&puplo,
&$PRIV(__n_size),
&$PRIV(__m_size),
$P(A),
&$PRIV(__n_size),
$P(ipiv),
$P(B),
&$PRIV(__n_size),
work,
&lwork,
$P(info));
}
',
Doc=>'
Complex version of sysv for Hermitian matrix
');
pp_defc("hesvx", HandleBad => 0, Pars => '[phys]A(2,n,n); int uplo(); int fact(); [phys]B(2,n,m); [io,phys]af(2,n,n); int [io,phys]ipiv(n); [o,phys]X(2,n,m); [o,phys]rcond(); [o,phys]ferr(m); [o,phys]berr(m); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code '
char pfact = \'N\';
char puplo = \'U\';
integer lwork = -1;
types(F) %{
extern int chesvx_(char *fact, char *uplo, integer *n, integer *
nrhs, float *a, integer *lda, float *af, integer *ldaf,
integer *ipiv, float *b, integer *ldb, float *x, integer *
ldx, float *rcond, float *ferr, float *berr,
float *work, integer *lwork, float *rwork, integer *info);
float *rwork = (float * )malloc ($PRIV(__n_size)* sizeof (float));
float tmp_work[2];
%}
types(D) %{
extern int zhesvx_(char *fact, char *uplo, integer *n, integer *
nrhs, double *a, integer *lda, double *af, integer *ldaf,
integer *ipiv, double *b, integer *ldb, double *x, integer *
ldx, double *rcond, double *ferr, double *berr,
double *work, integer *lwork, double *rwork, integer *info);
double *rwork = (double * )malloc ($PRIV(__n_size)* sizeof (double));
double tmp_work[2];
%}
if($fact())
pfact = \'F\';
if ($uplo())
puplo = \'L\';
$TFD(chesvx_,zhesvx_)(
&pfact,
&puplo,
&$PRIV(__n_size),
&$PRIV(__m_size),
$P(A),
&$PRIV(__n_size),
$P(af),
&$PRIV(__n_size),
$P(ipiv),
$P(B),
&$PRIV(__n_size),
$P(X),
&$PRIV(__n_size),
$P(rcond),
$P(ferr),
$P(berr),
&tmp_work[0],
&lwork,
rwork,
$P(info));
lwork = (integer )tmp_work[0];
{
types(F) %{
float *work = (float *)malloc(2*lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2*lwork * sizeof(double));
%}
$TFD(chesvx_,zhesvx_)(
&pfact,
&puplo,
&$PRIV(__n_size),
&$PRIV(__m_size),
$P(A),
&$PRIV(__n_size),
$P(af),
&$PRIV(__n_size),
$P(ipiv),
$P(B),
&$PRIV(__n_size),
$P(X),
&$PRIV(__n_size),
$P(rcond),
$P(ferr),
$P(berr),
work,
&lwork,
rwork,
$P(info));
free(work);
}
free(rwork);
', Doc=>'
Complex version of sysvx for Hermitian matrix
');
pp_defc("posv", HandleBad => 0, Pars => '[io,phys]A(2,n,n); int uplo(); [io,phys]B(2,n,m); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code '
char puplo = \'U\';
types(F) %{
extern int cposv_(char *uplo, integer *n, integer *nrhs, float
*a, integer *lda, float *b, integer *ldb, integer *info);
%}
types(D) %{
extern int zposv_(char *uplo, integer *n, integer *nrhs, double
*a, integer *lda, double *b, integer *ldb, integer *info);
%}
if ($uplo())
puplo = \'L\';
$TFD(cposv_,zposv_)(
&puplo,
&$PRIV(__n_size),
&$PRIV(__m_size),
$P(A),
&$PRIV(__n_size),
$P(B),
&$PRIV(__n_size),
$P(info));
',
Doc=>'
Complex version of posv for Hermitian positive definite matrix
'); pp_defc("posvx", HandleBad => 0, Pars => '[io,phys]A(2,n,n); int uplo(); int fact(); [io,phys]B(2,n,m); [io,phys]af(2,n,n); int [io]equed(); [io,phys]s(n); [o,phys]X(2,n,m); [o,phys]rcond(); [o,phys]ferr(m); [o,phys]berr(m); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code '
char pfact;
char pequed = \'N\';
char puplo = \'U\';
types(F) %{
extern int cposvx_(char *fact, char *uplo, integer *n, integer *
nrhs, float *a, integer *lda, float *af, integer *ldaf,
char *equed, float *s, float *b, integer *ldb, float *
x, integer *ldx, float *rcond, float *ferr, float *
berr, float *work, float *rwork, integer *info);
float *work, *rwork;
%}
types(D) %{
extern int zposvx_(char *fact, char *uplo, integer *n, integer *
nrhs, double *a, integer *lda, double *af, integer *ldaf,
char *equed, double *s, double *b, integer *ldb, double *
x, integer *ldx, double *rcond, double *ferr, double *
berr, double *work, double *rwork, integer *info);
double *work, *rwork;
%}
switch ($fact())
{
case 1: pfact = \'N\';
break;
case 2: pfact = \'E\';
break;
default: pfact = \'F\';
}
if ($equed())
pequed = \'Y\';
if ($uplo())
puplo = \'L\';
types(F) %{
work = (float *) malloc(4 * $PRIV(__n_size) * sizeof(float));
rwork = (float *) malloc(2 * $PRIV(__n_size) * sizeof(float));
%}
types(D) %{
work = (double *) malloc(4 * $PRIV(__n_size) * sizeof(double));
rwork = (double *) malloc(2 * $PRIV(__n_size) * sizeof(double));
%}
$TFD(cposvx_,zposvx_)(
&pfact,
&puplo,
&$PRIV(__n_size),
&$PRIV(__m_size),
$P(A),
&$PRIV(__n_size),
$P(af),
&$PRIV(__n_size),
&pequed,
$P(s),
$P(B),
&$PRIV(__n_size),
$P(X),
&$PRIV(__n_size),
$P(rcond),
$P(ferr),
$P(berr),
work,
rwork,
$P(info));
free(work);
free(rwork);
switch (pequed)
{
case \'Y\': $equed() = 1;
break;
default: $equed()= 0;
}
', Doc => '
Complex version of posvx for Hermitian positive definite matrix
');
pp_defc("gels", HandleBad => 0, Pars => '[io,phys]A(2,m,n); int trans(); [io,phys]B(2,p,q);int [o,phys]info()', GenericTypes => [F,D], Code => generate_code '
char ptrans = \'N\';
integer lwork = -1;
types(F) %{
extern int cgels_(char *trans, integer *m, integer *n, integer *
nrhs, float *a, integer *lda, float *b, integer *ldb,
float *work, integer *lwork, integer *info);
float tmp_work[2];
%}
types(D) %{
extern int zgels_(char *trans, integer *m, integer *n, integer *
nrhs, double *a, integer *lda, double *b, integer *ldb,
double *work, integer *lwork, integer *info);
double tmp_work[2];
%}
if($trans())
ptrans = \'C\';
$TFD(cgels_,zgels_)(
&ptrans,
&$PRIV(__m_size),
&$PRIV(__n_size),
&$PRIV(__q_size),
$P(A),
&$PRIV(__m_size),
$P(B),
&$PRIV(__p_size),
&tmp_work[0],
&lwork,
$P(info));
lwork = (integer )tmp_work[0];
{
types(F) %{
float *work = (float *)malloc(2*lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2*lwork * sizeof(double));
%}
$TFD(cgels_,zgels_)(
&ptrans,
&$PRIV(__m_size),
&$PRIV(__n_size),
&$PRIV(__q_size),
$P(A),
&$PRIV(__m_size),
$P(B),
&$PRIV(__p_size),
work,
&lwork,
$P(info));
free(work);
}
', Doc=>'
Solves overdetermined or underdetermined complex linear systems involving an M-by-N matrix A, or its conjugate-transpose. Complex version of gels.
trans: = 0: the linear system involves A;
= 1: the linear system involves A**H.
');
pp_defc("gelsy", HandleBad => 0, Pars => '[io,phys]A(2,m,n); [io,phys]B(2,p,q); [phys]rcond(); int [io,phys]jpvt(n); int [o,phys]rank();int [o,phys]info()', GenericTypes => [F,D], Code => generate_code '
integer lwork = -1;
types(F) %{
extern int cgelsy_(integer *m, integer *n, integer *nrhs,
float *a, integer *lda, float *b, integer *ldb, integer *
jpvt, float *rcond, integer *rank, float *work, integer *
lwork, float *rwork, integer *info);
float tmp_work[2];
float *rwork;
%}
types(D) %{
extern int zgelsy_(integer *m, integer *n, integer *nrhs,
double *a, integer *lda, double *b, integer *ldb, integer *
jpvt, double *rcond, integer *rank, double *work, integer *
lwork, double *rwork, integer *info);
double tmp_work[2];
double *rwork;
%}
types(F) %{
rwork = (float *)malloc( 2 * $PRIV(__m_size) * sizeof(float));
%}
types(D) %{
rwork = (double *)malloc(2 * $PRIV(__m_size) * sizeof(double));
%}
$TFD(cgelsy_,zgelsy_)(
&$PRIV(__m_size),
&$PRIV(__n_size),
&$PRIV(__q_size),
$P(A),
&$PRIV(__m_size),
$P(B),
&$PRIV(__p_size),
$P(jpvt),
$P(rcond),
$P(rank),
&tmp_work[0],
&lwork,
rwork,
$P(info));
lwork = (integer )tmp_work[0];
{
types(F) %{
float *work = (float *)malloc(2 * lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2 * lwork * sizeof(double));
%}
$TFD(cgelsy_,zgelsy_)(
&$PRIV(__m_size),
&$PRIV(__n_size),
&$PRIV(__q_size),
$P(A),
&$PRIV(__m_size),
$P(B),
&$PRIV(__p_size),
$P(jpvt),
$P(rcond),
$P(rank),
work,
&lwork,
rwork,
$P(info));
free(work);
}
free(rwork);
');
pp_defc("gelss", HandleBad => 0, Pars => '[io,phys]A(2,m,n); [io,phys]B(2,p,q); [phys]rcond(); [o,phys]s(r); int [o,phys]rank();int [o,phys]info()', GenericTypes => [F,D], Code => generate_code '
integer lwork = -1;
integer lrwork;
types(F) %{
extern int cgelss_(integer *m, integer *n, integer *nrhs,
float *a, integer *lda, float *b, integer *ldb, float *s,
float *rcond, integer *rank, float *work, integer *
lwork, float *rwork, integer *info);
float tmp_work[2];
float *rwork;
%}
types(D) %{
extern int zgelss_(integer *m, integer *n, integer *nrhs,
double *a, integer *lda, double *b, integer *ldb,
double *s,double *rcond, integer *rank, double *work, integer *
lwork, double *rwork, integer *info);
double tmp_work[2];
double *rwork;
%}
lrwork = min($PRIV(__m_size), $PRIV(__n_size));
types(F) %{
rwork = (float *)malloc(5 * lrwork * sizeof(float));
%}
types(D) %{
rwork = (double *)malloc(5 * lrwork * sizeof(double));
%}
$TFD(cgelss_,zgelss_)(
&$PRIV(__m_size),
&$PRIV(__n_size),
&$PRIV(__q_size),
$P(A),
&$PRIV(__m_size),
$P(B),
&$PRIV(__p_size),
$P(s),
$P(rcond),
$P(rank),
&tmp_work[0],
&lwork,
rwork,
$P(info));
lwork = (integer )tmp_work[0];
{
types(F) %{
float *work = (float *)malloc(2 * lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2 * lwork * sizeof(double));
%}
$TFD(cgelss_,zgelss_)(
&$PRIV(__m_size),
&$PRIV(__n_size),
&$PRIV(__q_size),
$P(A),
&$PRIV(__m_size),
$P(B),
&$PRIV(__p_size),
$P(s),
$P(rcond),
$P(rank),
work,
&lwork,
rwork,
$P(info));
free(work);
}
free(rwork);
');
pp_defc("gelsd", HandleBad => 0, Pars => '[io,phys]A(2,m,n); [io,phys]B(2,p,q); [phys]rcond(); [o,phys]s(r); int [o,phys]rank();int [o,phys]info()', GenericTypes => [F,D], Code => generate_code '
integer lwork = -1;
integer smlsiz, size_i, nlvl, *iwork;
integer minmn = min( $SIZE(m), $SIZE(n) );
extern integer ilaenv_(integer *ispec, char *name__, char *opts, integer *n1,
integer *n2, integer *n3, integer *n4, ftnlen name_len, ftnlen
opts_len);
types(F) %{
extern int cgelsd_(integer *m, integer *n, integer *nrhs,
float *a, integer *lda, float *b, integer *ldb, float *s,
float *rcond, integer *rank, float *work, integer *
lwork, float *rwork, integer *iwork, integer *info);
float *rwork;
float tmp_work[2];
%}
types(D) %{
extern int zgelsd_(integer *m, integer *n, integer *nrhs,
double *a, integer *lda, double *b, integer *ldb,
double *s,double *rcond, integer *rank, double *work, integer *
lwork, double *rwork, integer *iwork,integer *info);
double *rwork;
double tmp_work[2];
%}
minmn = max(1,minmn);
types(F) %{
smlsiz = ilaenv_(&c_nine, "CGELSD", " ", &c_zero, &c_zero, &c_zero, &c_zero, (ftnlen)6, (ftnlen)1);
size_i = (integer) (log((float) minmn / (float) (smlsiz + 1)) /log(2.)) + 1;
if ($PRIV(__m_size) >= $PRIV(__n_size)){
rwork = (float *) malloc ((10*$PRIV(__n_size) + 2 * $PRIV(__n_size) * smlsiz + 8 * $PRIV(__n_size) * size_i + 3 * smlsiz * $PRIV(__q_size) + pow((smlsiz+1),2)) * sizeof(float));
}
else{
rwork = (float *) malloc ((10*$PRIV(__m_size) + 2 * $PRIV(__m_size) * smlsiz + 8 * $PRIV(__m_size) * size_i + 3 * smlsiz * $PRIV(__q_size) + pow((smlsiz+1),2)) * sizeof(float));
}
%}
types(D) %{
smlsiz = ilaenv_(&c_nine, "ZGELSD", " ", &c_zero, &c_zero, &c_zero, &c_zero, (ftnlen)6, (ftnlen)1);
size_i = (integer) (log((double) minmn / (double) (smlsiz + 1)) /log(2.)) + 1;
if ($PRIV(__m_size) >= $PRIV(__n_size)){
rwork = (double *) malloc ((10*$PRIV(__n_size) + 2 * $PRIV(__n_size) * smlsiz + 8 * $PRIV(__n_size) * size_i + 3 * smlsiz * $PRIV(__q_size) + pow((smlsiz+1),2)) * sizeof(double));
}
else{
rwork = (double *) malloc ((10*$PRIV(__m_size) + 2 * $PRIV(__m_size) * smlsiz + 8 * $PRIV(__m_size) * size_i + 3 * smlsiz * $PRIV(__q_size) + pow((smlsiz+1),2)) * sizeof(double));
}
%}
nlvl = max(size_i, 0);
iwork = (integer *)malloc((3 * minmn * nlvl + 11 * minmn) * sizeof(integer));
$TFD(cgelsd_,zgelsd_)(
&$PRIV(__m_size),
&$PRIV(__n_size),
&$PRIV(__q_size),
$P(A),
&$PRIV(__m_size),
$P(B),
&$PRIV(__p_size),
$P(s),
$P(rcond),
$P(rank),
&tmp_work[0],
&lwork,
rwork,
iwork,
$P(info));
lwork = (integer )tmp_work[0];
{
types(F) %{
float *work = (float *)malloc(2 * lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2 * lwork * sizeof(double));
%}
$TFD(cgelsd_,zgelsd_)(
&$PRIV(__m_size),
&$PRIV(__n_size),
&$PRIV(__q_size),
$P(A),
&$PRIV(__m_size),
$P(B),
&$PRIV(__p_size),
$P(s),
$P(rcond),
$P(rank),
work,
&lwork,
rwork,
iwork,
$P(info));
free(work);
}
free (iwork);
free (rwork);
');
pp_defc("gglse", HandleBad => 0, Pars => '[phys]A(2,m,n); [phys]B(2,p,n);[io,phys]c(2,m);[phys]d(2,p);[o,phys]x(2,n);int [o,phys]info()', GenericTypes => [F,D], Code => generate_code '
integer lwork = -1;
types(F) %{
extern int cgglse_(integer *m, integer *n, integer *p, float *
a, integer *lda, float *b, integer *ldb, float *c__,
float *d__, float *x, float *work, integer *lwork,
integer *info);
float tmp_work[2];
%}
types(D) %{
extern int zgglse_(integer *m, integer *n, integer *p, double *
a, integer *lda, double *b, integer *ldb, double *c__,
double *d__, double *x, double *work, integer *lwork,
integer *info);
double tmp_work[2];
%}
$TFD(cgglse_,zgglse_)(
&$PRIV(__m_size),
&$PRIV(__n_size),
&$PRIV(__p_size),
$P(A),
&$PRIV(__m_size),
$P(B),
&$PRIV(__p_size),
$P(c),
$P(d),
$P(x),
&tmp_work[0],
&lwork,
$P(info));
lwork = (integer )tmp_work[0];
{
types(F) %{
float *work = (float *)malloc(2 * lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2 * lwork * sizeof(double));
%}
$TFD(cgglse_,zgglse_)(
&$PRIV(__m_size),
&$PRIV(__n_size),
&$PRIV(__p_size),
$P(A),
&$PRIV(__m_size),
$P(B),
&$PRIV(__p_size),
$P(c),
$P(d),
$P(x),
work,
&lwork,
$P(info));
free(work);
}
');
pp_defc("ggglm", HandleBad => 0, Pars => '[phys]A(2,n,m); [phys]B(2,n,p);[phys]d(2,n);[o,phys]x(2,m);[o,phys]y(2,p);int [o,phys]info()', GenericTypes => [F,D], Code => generate_code ' integer lwork = -1;
types(F) %{
extern int cggglm_(integer *n, integer *m, integer *p, float *
a, integer *lda, float *b, integer *ldb, float *d__,
float *x, float *y, float *work, integer *lwork,
integer *info);
float tmp_work[2];
%}
types(D) %{
extern int zggglm_(integer *n, integer *m, integer *p, double *
a, integer *lda, double *b, integer *ldb, double *d__,
double *x, double *y, double *work, integer *lwork,
integer *info);
double tmp_work[2];
%}
$TFD(cggglm_,zggglm_)(
&$PRIV(__n_size),
&$PRIV(__m_size),
&$PRIV(__p_size),
$P(A),
&$PRIV(__n_size),
$P(B),
&$PRIV(__n_size),
$P(d),
$P(x),
$P(y),
&tmp_work[0],
&lwork,
$P(info));
lwork = (integer )tmp_work[0];
{
types(F) %{
float *work = (float *)malloc(2 * lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2 * lwork * sizeof(double));
%}
$TFD(cggglm_,zggglm_)(
&$PRIV(__n_size),
&$PRIV(__m_size),
&$PRIV(__p_size),
$P(A),
&$PRIV(__n_size),
$P(B),
&$PRIV(__n_size),
$P(d),
$P(x),
$P(y),
work,
&lwork,
$P(info));
free(work);
}
');
################################################################################ # # COMPUTATIONAL LEVEL ROUTINES # ################################################################################ # TODO IPIV = min(m,n) pp_defc("getrf", HandleBad => 0, RedoDimsCode => '$SIZE(p) = $PDL(A)->ndims > 2 ? min($PDL(A)->dims[1], $PDL(A)->dims[2]) : 1;', Pars => '[io,phys]A(2,m,n); int [o,phys]ipiv(p); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code ' types(F) %{
extern int cgetrf_(integer *m, integer *n, float *a, integer *
lda, integer *ipiv, integer *info);
%}
types(D) %{
extern int zgetrf_(integer *m, integer *n, double *a, integer *
lda, integer *ipiv, integer *info);
%}
$TFD(cgetrf_,zgetrf_)(
&$PRIV(__m_size),
&$PRIV(__n_size),
$P(A),
&$PRIV(__m_size),
$P(ipiv),
$P(info));
');
pp_defc("getf2", HandleBad => 0, RedoDimsCode => '$SIZE(p) = $PDL(A)->ndims > 2 ? min($PDL(A)->dims[1], $PDL(A)->dims[2]) : 1;', Pars => '[io,phys]A(2,m,n); int [o,phys]ipiv(p); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code ' types(F) %{
extern int cgetf2_(integer *m, integer *n, float *a, integer *
lda, integer *ipiv, integer *info);
%}
types(D) %{
extern int zgetf2_(integer *m, integer *n, double *a, integer *
lda, integer *ipiv, integer *info);
%}
$TFD(cgetf2_,zgetf2_)(
&$PRIV(__m_size),
&$PRIV(__n_size),
$P(A),
&$PRIV(__m_size),
$P(ipiv),
$P(info));
');
pp_defc("sytrf", HandleBad => 0, Pars => '[io,phys]A(2,n,n); int uplo(); int [o,phys]ipiv(n); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code '
char puplo = \'U\';
integer lwork = -1;
types(F) %{
extern int csytrf_(char *uplo, integer *n, float *a, integer *
lda, integer *ipiv, float *work, integer *lwork, integer *info);
float tmp_work[2];
%}
types(D) %{
extern int zsytrf_(char *uplo, integer *n, double *a, integer *
lda, integer *ipiv, double *work, integer *lwork, integer *info);
double tmp_work[2];
%}
if ($uplo())
puplo = \'L\';
$TFD(csytrf_,zsytrf_)(
&puplo,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(ipiv),
&tmp_work[0],
&lwork,
$P(info));
lwork = (integer )tmp_work[0];
{
types(F) %{
float *work = (float *)malloc(2*lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2*lwork * sizeof(double));
%}
$TFD(csytrf_,zsytrf_)(
&puplo,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(ipiv),
work,
&lwork,
$P(info));
free (work);
}
');
pp_defc("sytf2", HandleBad => 0, Pars => '[io,phys]A(2,n,n); int uplo(); int [o,phys]ipiv(n); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code '
char puplo = \'U\';
types(F) %{
extern int csytf2_(char *uplo, integer *n, float *a, integer *
lda, integer *ipiv, integer *info);
%}
types(D) %{
extern int zsytf2_(char *uplo, integer *n, double *a, integer *
lda, integer *ipiv, integer *info);
%}
if ($uplo())
puplo = \'L\';
$TFD(csytf2_,zsytf2_)(
&puplo,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(ipiv),
$P(info));
');
pp_defc("chetrf", HandleBad => 0, Pars => '[io,phys]A(2,n,n); int uplo(); int [o,phys]ipiv(n); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code '
char puplo = \'U\';
integer lwork = -1;
integer blocksiz;
extern integer ilaenv_(integer *ispec, char *name__, char *opts, integer *n1,
integer *n2, integer *n3, integer *n4, ftnlen name_len, ftnlen
opts_len);
types(F) %{
extern int chetrf_(char *uplo, integer *n, float *a, integer *
lda, integer *ipiv, float *work, integer *lwork, integer *info);
float *work;
blocksiz = ilaenv_(&c_nine, "CHETRF", " ", &c_zero, &c_zero, &c_zero, &c_zero, (ftnlen)6, (ftnlen)1);
%}
types(D) %{
extern int zhetrf_(char *uplo, integer *n, double *a, integer *
lda, integer *ipiv, double *work, integer *lwork, integer *info);
double *work;
blocksiz = ilaenv_(&c_nine, "ZHETRF", " ", &c_zero, &c_zero, &c_zero, &c_zero, (ftnlen)6, (ftnlen)1);
%}
if ($uplo())
puplo = \'L\';
lwork = (integer ) $PRIV(__n_size) * blocksiz;
types(F) %{
work = (float *)malloc(2*lwork * sizeof(float));
%}
types(D) %{
work = (double *)malloc(2*lwork * sizeof(double));
%}
$TFD(chetrf_,zhetrf_)(
&puplo,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(ipiv),
work,
&lwork,
$P(info));
free (work);
',
Doc=>'
Complex version of sytrf for Hermitian matrix
');
pp_defc("hetf2", HandleBad => 0, Pars => '[io,phys]A(2,n,n); int uplo(); int [o,phys]ipiv(n); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code '
char puplo = \'U\';
types(F) %{
extern int chetf2_(char *uplo, integer *n, float *a, integer *
lda, integer *ipiv, integer *info);
%}
types(D) %{
extern int zhetf2_(char *uplo, integer *n, double *a, integer *
lda, integer *ipiv, integer *info);
%}
if ($uplo())
puplo = \'L\';
$TFD(chetf2_,zhetf2_)(
&puplo,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(ipiv),
$P(info));
',
Doc=>'
Complex version of sytf2 for Hermitian matrix
');
pp_defc("potrf", HandleBad => 0, Pars => '[io,phys]A(2,n,n); int uplo(); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code '
char puplo = \'U\';
types(F) %{
extern int cpotrf_(char *uplo, integer *n, float *a, integer *
lda, integer *info);
%}
types(D) %{
extern int zpotrf_(char *uplo, integer *n, double *a, integer *
lda, integer *info);
%}
if ($uplo())
puplo = \'L\';
$TFD(cpotrf_,zpotrf_)(
&puplo,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(info));
',
Doc=>'
Complex version of potrf for Hermitian positive definite matrix
');
pp_defc("potf2", HandleBad => 0, Pars => '[io,phys]A(2,n,n); int uplo(); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code '
char puplo = \'U\';
types(F) %{
extern int cpotf2_(char *uplo, integer *n, float *a, integer *
lda, integer *info);
%}
types(D) %{
extern int zpotf2_(char *uplo, integer *n, double *a, integer *
lda, integer *info);
%}
if ($uplo())
puplo = \'L\';
$TFD(cpotf2_,zpotf2_)(
&puplo,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(info));
',
Doc => '
Complex version of potf2 for Hermitian positive definite matrix
');
pp_defc("getri", HandleBad => 0, Pars => '[io,phys]A(2,n,n); int [phys]ipiv(n); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code '
integer lwork = -1;
types(F) %{
extern int cgetri_(integer *n, float *a, integer *lda, integer
*ipiv, float *work, integer *lwork, integer *info);
float tmp_work[2];
%}
types(D) %{
extern int zgetri_(integer *n, double *a, integer *lda, integer
*ipiv, double *work, integer *lwork, integer *info);
double tmp_work[2];
%}
$TFD(cgetri_,zgetri_)(
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(ipiv),
&tmp_work[0],
&lwork,
$P(info));
lwork = (integer )tmp_work[0];
{
types(F) %{
float *work = (float *)malloc(2*lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2*lwork * sizeof(double));
%}
$TFD(cgetri_,zgetri_)(
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(ipiv),
work,
&lwork,
$P(info));
free(work);
}
');
pp_defc("sytri", HandleBad => 0, Pars => '[io,phys]A(2,n,n); int uplo(); int [phys]ipiv(n); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code '
char puplo = \'U\';
types(F) %{
extern int csytri_(char *uplo, integer *n, float *a, integer *
lda, integer *ipiv, float *work, integer *info);
float *work = (float *)malloc(2*$PRIV(__n_size) * sizeof(float));
%}
types(D) %{
extern int zsytri_(char *uplo, integer *n, double *a, integer *
lda, integer *ipiv, double *work, integer *info);
double *work = (double *)malloc(2*$PRIV(__n_size) * sizeof(double));
%}
if ($uplo())
puplo = \'L\';
$TFD(csytri_, zsytri_)(
&puplo,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(ipiv),
work,
$P(info));
free(work);
');
pp_defc("hetri", HandleBad => 0, Pars => '[io,phys]A(2,n,n); int uplo(); int [phys]ipiv(n); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code ' char puplo = \'U\'; types(F) %{
extern int chetri_(char *uplo, integer *n, float *a, integer *
lda, integer *ipiv, float *work, integer *info);
float *work = (float *)malloc(2*$PRIV(__n_size) * sizeof(float));
%}
types(D) %{
extern int zhetri_(char *uplo, integer *n, double *a, integer *
lda, integer *ipiv, double *work, integer *info);
double *work = (double *)malloc(2*$PRIV(__n_size) * sizeof(double));
%}
if ($uplo())
puplo = \'L\';
$TFD(chetri_, zhetri_)(
&puplo,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(ipiv),
work,
$P(info));
free(work);
',
Doc => '
Complex version of sytri for Hermitian matrix
');
pp_defc("potri", HandleBad => 0, Pars => '[io,phys]A(2,n,n); int uplo(); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code '
char puplo = \'U\';
types(F) %{
extern int cpotri_(char *uplo, integer *n, float *a, integer *
lda, integer *info);
%}
types(D) %{
extern int zpotri_(char *uplo, integer *n, double *a, integer *
lda, integer *info);
%}
if ($uplo())
puplo = \'L\';
$TFD(cpotri_,zpotri_)(
&puplo,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(info));
');
pp_defc("trtri", HandleBad => 0, Pars => '[io,phys]A(2,n,n); int uplo(); int diag(); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code ' char puplo = \'U\'; char pdiag = \'N\'; types(F) %{
extern int ctrtri_(char *uplo, char *diag, integer *n, float *a, integer *
lda, integer *info);
%}
types(D) %{
extern int ztrtri_(char *uplo, char *diag, integer *n, double *a, integer *
lda, integer *info);
%}
if ($uplo())
puplo = \'L\';
if ($diag())
pdiag = \'U\';
$TFD(ctrtri_, ztrtri_)(
&puplo,
&pdiag,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(info));
');
pp_defc("trti2", HandleBad => 0, Pars => '[io,phys]A(2,n,n); int uplo(); int diag(); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code ' char puplo = \'U\'; char pdiag = \'N\'; types(F) %{
extern int ctrti2_(char *uplo, char *diag, integer *n, float *a, integer *
lda, integer *info);
%}
types(D) %{
extern int ztrti2_(char *uplo, char *diag, integer *n, double *a, integer *
lda, integer *info);
%}
if ($uplo())
puplo = \'L\';
if ($diag())
pdiag = \'U\';
$TFD(ctrti2_, ztrti2_)(
&puplo,
&pdiag,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(info));
');
pp_defc("getrs", HandleBad => 0, Pars => '[phys]A(2,n,n); int trans(); [io,phys]B(2,n,m); int [phys]ipiv(n); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code ' char transp = \'N\'; types(F) %{ extern int cgetrs_(char *trans, integer *n, integer *nrhs, float *a, integer *lda, integer *ipiv, float *b, integer * ldb, integer *info); %} types(D) %{ extern int zgetrs_(char *trans, integer *n, integer *nrhs, double *a, integer *lda, integer *ipiv, double *b, integer * ldb, integer *info); %} if($trans() == 1) transp = \'T\'; else if($trans() == 2) transp = \'C\';
$TFD(cgetrs_,zgetrs_)(
&transp,
&$PRIV(__n_size),
&$PRIV(__m_size),
$P(A),
&$PRIV(__n_size),
$P(ipiv),
$P(B),
&$PRIV(__n_size),
$P(info));
',
Doc=>'
Complex version of getrs
Arguments
=========
trans: = 0: No transpose;
= 1: Transpose;
= 2: Conjugate transpose;
');
pp_defc("sytrs", HandleBad => 0, Pars => '[phys]A(2,n,n); int uplo();[io,phys]B(2,n,m); int [phys]ipiv(n); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code ' char puplo = \'U\'; types(F) %{ extern int csytrs_(char *uplo, integer *n, integer *nrhs, float *a, integer *lda, integer *ipiv, float *b, integer * ldb, integer *info); %} types(D) %{ extern int zsytrs_(char *uplo, integer *n, integer *nrhs, double *a, integer *lda, integer *ipiv, double *b, integer * ldb, integer *info); %} if($uplo()) puplo = \'L\';
$TFD(csytrs_,zsytrs_)(
&puplo,
&$PRIV(__n_size),
&$PRIV(__m_size),
$P(A),
&$PRIV(__n_size),
$P(ipiv),
$P(B),
&$PRIV(__n_size),
$P(info));
');
pp_defc("hetrs", HandleBad => 0, Pars => '[phys]A(2,n,n); int uplo();[io,phys]B(2,n,m); int [phys]ipiv(n); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code ' char puplo = \'U\'; types(F) %{ extern int chetrs_(char *uplo, integer *n, integer *nrhs, float *a, integer *lda, integer *ipiv, float *b, integer * ldb, integer *info); %} types(D) %{ extern int zhetrs_(char *uplo, integer *n, integer *nrhs, double *a, integer *lda, integer *ipiv, double *b, integer * ldb, integer *info); %} if($uplo()) puplo = \'L\';
$TFD(chetrs_,zhetrs_)(
&puplo,
&$PRIV(__n_size),
&$PRIV(__m_size),
$P(A),
&$PRIV(__n_size),
$P(ipiv),
$P(B),
&$PRIV(__n_size),
$P(info));
',
Doc => '
Complex version of sytrs for Hermitian matrix
');
pp_defc("potrs", HandleBad => 0, Pars => '[phys]A(2,n,n); int uplo(); [io,phys]B(2,n,m); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code ' char puplo = \'U\'; types(F) %{ extern int cpotrs_(char *uplo, integer *n, integer *nrhs, float *a, integer *lda, float *b, integer *ldb, integer * info); %} types(D) %{ extern int zpotrs_(char *uplo, integer *n, integer *nrhs, double *a, integer *lda, double *b, integer *ldb, integer * info); %} if($uplo()) puplo = \'L\';
$TFD(cpotrs_,zpotrs_)(
&puplo,
&$PRIV(__n_size),
&$PRIV(__m_size),
$P(A),
&$PRIV(__n_size),
$P(B),
&$PRIV(__n_size),
$P(info));
',
Doc=>'
Complex version of potrs for Hermitian positive definite matrix
');
pp_defc("trtrs", HandleBad => 0, Pars => '[phys]A(2,n,n); int uplo(); int trans(); int diag();[io,phys]B(2,n,m); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code ' char puplo = \'U\'; char ptrans = \'N\'; char pdiag = \'N\'; types(F) %{ extern int ctrtrs_(char *uplo, char *trans, char *diag, integer *n, integer *nrhs, float *a, integer *lda, float *b, integer * ldb, integer *info); %} types(D) %{ extern int ztrtrs_(char *uplo, char *trans, char *diag, integer *n, integer *nrhs, double *a, integer *lda, double *b, integer * ldb, integer *info); %} if($uplo()) puplo = \'L\'; if($trans() == 1) ptrans = \'T\'; else if($trans() == 2) ptrans = \'C\'; if($diag()) pdiag = \'U\';
$TFD(ctrtrs_,ztrtrs_)(
&puplo,
&ptrans,
&pdiag,
&$PRIV(__n_size),
&$PRIV(__m_size),
$P(A),
&$PRIV(__n_size),
$P(B),
&$PRIV(__n_size),
$P(info));
',
Doc=>'
Complex version of trtrs
Arguments
=========
trans: = 0: No transpose;
= 1: Transpose;
= 2: Conjugate transpose;
');
pp_defc("latrs", HandleBad => 0, Pars => '[phys]A(2,n,n); int uplo(); int trans(); int diag(); int normin();[io,phys]x(2,n); [o,phys]scale();[io,phys]cnorm(n);int [o,phys]info()', GenericTypes => [F,D], Code => generate_code ' char puplo = \'U\'; char ptrans = \'N\'; char pdiag = \'N\'; char pnormin = \'N\';
types(F) %{
extern int clatrs_(char *uplo, char *trans, char *diag, char *
normin, integer *n, float *a, integer *lda, float *x,
float *scale, float *cnorm, integer *info);
%}
types(D) %{
extern int zlatrs_(char *uplo, char *trans, char *diag, char *
normin, integer *n, double *a, integer *lda, double *x,
double *scale, double *cnorm, integer *info);
%}
if($uplo())
puplo = \'L\';
if($trans())
ptrans = \'T\';
else if($trans() == 2)
ptrans = \'C\';
if($diag())
pdiag = \'U\';
if($normin())
pnormin = \'Y\';
$TFD(clatrs_,zlatrs_)(
&puplo,
&ptrans,
&pdiag,
&pnormin,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(x),
$P(scale),
$P(cnorm),
$P(info));
',
Doc=>'
Complex version of latrs
Arguments
=========
trans: = 0: No transpose;
= 1: Transpose;
= 2: Conjugate transpose;
');
pp_defc("gecon", HandleBad => 0, Pars => '[phys]A(2,n,n); int norm(); [phys]anorm(); [o,phys]rcond();int [o,phys]info()', GenericTypes => [F,D], Code => generate_code '
char pnorm = \'I\';
types(F) %{
extern int sgecon_(char *norm, integer *n, float *a, integer *
lda, float *anorm, float *rcond, float *work, float *
rwork, integer *info);
float *work = (float *) malloc(($PRIV(__n_size) * 4) * sizeof(float));
float *rwork = (float *) malloc(($PRIV(__n_size) * 2)* sizeof(integer));
%}
types(D) %{
extern int dgecon_(char *norm, integer *n, double *a, integer *
lda, double *anorm, double *rcond, double *work, double *
rwork, integer *info);
double *work = (double *) malloc(($PRIV(__n_size)*4) * sizeof(double));
double *rwork = (double *) malloc(($PRIV(__n_size)*2 )* sizeof(integer));
%}
if($norm())
pnorm = \'O\';
$TFD(sgecon_,dgecon_)(
&pnorm,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(anorm),
$P(rcond),
work,
rwork,
$P(info));
free (work);
free(rwork);
');
pp_defc("sycon", HandleBad => 0, Pars => '[phys]A(2,n,n); int uplo(); int ipiv(n); [phys]anorm(); [o,phys]rcond();int [o,phys]info()', GenericTypes => [F,D], Code => generate_code '
char puplo = \'U\';
types(F) %{
extern int csycon_(char *uplo, integer *n, float *a, integer *
lda, integer *ipiv, float *anorm, float *rcond, float *
work, integer *info);
float *work = (float *) malloc(($PRIV(__n_size) * 4) * sizeof(float));
%}
types(D) %{
extern int zsycon_(char *uplo, integer *n, double *a, integer *
lda, integer *ipiv, double *anorm, double *rcond, double *
work, integer *info);
double *work = (double *) malloc(($PRIV(__n_size)*4) * sizeof(double));
%}
if($uplo())
puplo = \'L\';
$TFD(csycon_,zsycon_)(
&puplo,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(ipiv),
$P(anorm),
$P(rcond),
work,
$P(info));
free (work);
');
pp_defc("hecon", HandleBad => 0, Pars => '[phys]A(2,n,n); int uplo(); int ipiv(n); [phys]anorm(); [o,phys]rcond();int [o,phys]info()', GenericTypes => [F,D], Code => generate_code '
char puplo = \'U\';
types(F) %{
extern int checon_(char *uplo, integer *n, float *a, integer *
lda, integer *ipiv, float *anorm, float *rcond, float *
work, integer *info);
float *work = (float *) malloc(($PRIV(__n_size) * 4) * sizeof(float));
%}
types(D) %{
extern int zhecon_(char *uplo, integer *n, double *a, integer *
lda, integer *ipiv, double *anorm, double *rcond, double *
work, integer *info);
double *work = (double *) malloc(($PRIV(__n_size)*4) * sizeof(double));
%}
if($uplo())
puplo = \'L\';
$TFD(checon_,zhecon_)(
&puplo,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(ipiv),
$P(anorm),
$P(rcond),
work,
$P(info));
free (work);
', Doc => '
Complex version of sycon for Hermitian matrix
');
pp_defc("pocon", HandleBad => 0, Pars => '[phys]A(2,n,n); int uplo(); [phys]anorm(); [o,phys]rcond();int [o,phys]info()', GenericTypes => [F,D], Code => generate_code '
char puplo = \'U\';
types(F) %{
extern int cpocon_(char *uplo, integer *n, float *a, integer *
lda, float *anorm, float *rcond, float *work, float *
rwork, integer *info);
float *work = (float *) malloc(($PRIV(__n_size) * 4) * sizeof(float));
float *rwork = (float *) malloc($PRIV(__n_size) * sizeof(integer));
%}
types(D) %{
extern int zpocon_(char *uplo, integer *n, double *a, integer *
lda, double *anorm, double *rcond, double *work, double *
rwork, integer *info);
double *work = (double *) malloc(($PRIV(__n_size)* 4) * sizeof(double));
double *rwork = (double *) malloc($PRIV(__n_size) * sizeof(integer));
%}
if($uplo())
puplo = \'L\';
$TFD(cpocon_, zpocon_)(
&puplo,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(anorm),
$P(rcond),
work,
rwork,
$P(info));
free (work);
free(rwork);
', Doc => '
Complex version of pocon for Hermitian positive definite matrix
');
pp_defc("trcon", HandleBad => 0, Pars => '[phys]A(2,n,n); int norm();int uplo();int diag(); [o,phys]rcond();int [o,phys]info()', GenericTypes => [F,D], Code => generate_code ' char puplo = \'U\'; char pdiag = \'N\'; char pnorm = \'I\'; types(F) %{ extern int strcon_(char *norm, char *uplo, char *diag,integer *n, float *a, integer * lda, float *rcond, float *work, float *rwork, integer *info); float *work = (float *) malloc(($PRIV(__n_size) * 4) * sizeof(float)); float *rwork = (float *) malloc($PRIV(__n_size) * sizeof(integer)); %} types(D) %{ extern int dtrcon_(char *norm, char *uplo, char *diag, integer *n, double *a, integer * lda, double *rcond, double * work, double *rwork, integer *info); double *work = (double *) malloc(($PRIV(__n_size)* 4) * sizeof(double)); double *rwork = (double *) malloc($PRIV(__n_size) * sizeof(integer)); %}
if($uplo())
puplo = \'L\';
if($diag())
pdiag = \'U\';
if($norm())
pnorm = \'O\';
$TFD(strcon_,dtrcon_)(
&pnorm,
&puplo,
&pdiag,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(rcond),
work,
rwork,
$P(info));
free (work);
free(rwork);
');
pp_defc("geqp3", HandleBad => 0, Pars => '[io,phys]A(2,m,n); int [io,phys]jpvt(n); [o,phys]tau(2,k); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code '
integer lwork = -1;
types(F) %{
extern int cgeqp3_(integer *m, integer *n, float *a, integer *
lda, integer *jpvt, float *tau, float *work, integer *lwork,
float *rwork, integer *info);
float tmp_work[2], *rwork;
rwork = (float *) malloc ($PRIV(__n_size) * 2 * sizeof(float));
%}
types(D) %{
extern int zgeqp3_(integer *m, integer *n, double *a, integer *
lda, integer *jpvt, double *tau, double *work, integer *lwork,
double *rwork, integer *info);
double tmp_work[2], *rwork;
rwork = (double *) malloc ($PRIV(__n_size) * 2 * sizeof(double));
%}
$TFD(cgeqp3_,zgeqp3_)(
&$PRIV(__m_size),
&$PRIV(__n_size),
$P(A),
&$PRIV(__m_size),
$P(jpvt),
$P(tau),
&tmp_work[0],
&lwork,
rwork,
$P(info));
lwork = (integer )tmp_work[0];
{
types(F) %{
float *work = (float *)malloc(2 * lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2 * lwork * sizeof(double));
%}
$TFD(cgeqp3_,zgeqp3_)(
&$PRIV(__m_size),
&$PRIV(__n_size),
$P(A),
&$PRIV(__m_size),
$P(jpvt),
$P(tau),
work,
&lwork,
rwork,
$P(info));
free(work);
}
free(rwork);
'
);
pp_defc("geqrf", HandleBad => 0, Pars => '[io,phys]A(2,m,n); [o,phys]tau(2,k); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code ' integer lwork = -1;
types(F) %{
extern int cgeqrf_(integer *m, integer *n, float *a, integer *
lda, float *tau, float *work, integer *lwork,
integer *info);
float tmp_work[2];
%}
types(D) %{
extern int zgeqrf_(integer *m, integer *n, double *a, integer *
lda, double *tau, double *work, integer *lwork,
integer *info);
double tmp_work[2];
%}
$TFD(cgeqrf_,zgeqrf_)(
&$PRIV(__m_size),
&$PRIV(__n_size),
$P(A),
&$PRIV(__m_size),
$P(tau),
&tmp_work[0],
&lwork,
$P(info));
lwork = (integer )tmp_work[0];
{
types(F) %{
float *work = (float *)malloc(2 * lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2 * lwork * sizeof(double));
%}
$TFD(cgeqrf_,zgeqrf_)(
&$PRIV(__m_size),
&$PRIV(__n_size),
$P(A),
&$PRIV(__m_size),
$P(tau),
work,
&lwork,
$P(info));
free(work);
}
');
pp_defc("ungqr", HandleBad => 0, Pars => '[io,phys]A(2,m,n); [phys]tau(2,k); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code ' integer lwork = -1;
types(F) %{
extern int cungqr_(integer *m, integer *n, integer *k, float *
a, integer *lda, float *tau, float *work, integer *lwork,
integer *info);
float tmp_work[2];
%}
types(D) %{
extern int zungqr_(integer *m, integer *n, integer *k, double *
a, integer *lda, double *tau, double *work, integer *lwork,
integer *info);
double tmp_work[2];
%}
$TFD(cungqr_, zungqr_)(
&$PRIV(__m_size),
&$PRIV(__n_size),
&$PRIV(__k_size),
$P(A),
&$PRIV(__m_size),
$P(tau),
&tmp_work[0],
&lwork,
$P(info));
lwork = (integer )tmp_work[0];
{
types(F) %{
float *work = (float *)malloc(2 * lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2 * lwork * sizeof(double));
%}
$TFD(cungqr_,zungqr_)(
&$PRIV(__m_size),
&$PRIV(__n_size),
&$PRIV(__k_size),
$P(A),
&$PRIV(__m_size),
$P(tau),
work,
&lwork,
$P(info));
free(work);
}
', Doc=>'
Complex version of orgqr
');
pp_defc("unmqr", HandleBad => 0, Pars => '[phys]A(2,p,k); int side(); int trans(); [phys]tau(2,k); [io,phys]C(2,m,n);int [o,phys]info()', GenericTypes => [F,D], Code => generate_code ' char ptrans = \'N\', pside = \'L\'; integer lwork = -1;
types(F) %{
extern int cunmqr_(char *side, char *trans, integer *m, integer *n,
integer *k, float *a, integer *lda, float *tau, float *
c__, integer *ldc, float *work, integer *lwork, integer *info);
float tmp_work[2];
%}
types(D) %{
extern int zunmqr_(char *side, char *trans, integer *m, integer *n,
integer *k, double *a, integer *lda, double *tau, double *
c__, integer *ldc, double *work, integer *lwork, integer *info);
double tmp_work[2];
%}
if($trans())
ptrans = \'C\';
if($side())
pside = \'R\';
$TFD(cunmqr_,zunmqr_)(
&pside,
&ptrans,
&$PRIV(__m_size),
&$PRIV(__n_size),
&$PRIV(__k_size),
$P(A),
&$PRIV(__p_size),
$P(tau),
$P(C),
&$PRIV(__m_size),
&tmp_work[0],
&lwork,
$P(info));
lwork = (integer )tmp_work[0];
{
types(F) %{
float *work = (float *)malloc(2 * lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2 * lwork * sizeof(double));
%}
$TFD(cunmqr_,zunmqr_)(
&pside,
&ptrans,
&$PRIV(__m_size),
&$PRIV(__n_size),
&$PRIV(__k_size),
$P(A),
&$PRIV(__p_size),
$P(tau),
$P(C),
&$PRIV(__m_size),
work,
&lwork,
$P(info));
free(work);
}
', Doc=>'
Complex version of ormqr. Here trans = 1 means conjugate transpose.
');
pp_defc("gelqf", HandleBad => 0, Pars => '[io,phys]A(2,m,n); [o,phys]tau(2,k); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code ' integer lwork = -1;
types(F) %{
extern int cgelqf_(integer *m, integer *n, float *a, integer *
lda, float *tau, float *work, integer *lwork, integer *info);
float tmp_work[2];
%}
types(D) %{
extern int zgelqf_(integer *m, integer *n, double *a, integer *
lda, double *tau, double *work, integer *lwork, integer *info);
double tmp_work[2];
%}
$TFD(cgelqf_,zgelqf_)(
&$PRIV(__m_size),
&$PRIV(__n_size),
$P(A),
&$PRIV(__m_size),
$P(tau),
&tmp_work[0],
&lwork,
$P(info));
lwork = (integer )tmp_work[0];
{
types(F) %{
float *work = (float *)malloc(2 * lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2 * lwork * sizeof(double));
%}
$TFD(cgelqf_,zgelqf_)(
&$PRIV(__m_size),
&$PRIV(__n_size),
$P(A),
&$PRIV(__m_size),
$P(tau),
work,
&lwork,
$P(info));
free(work);
}
');
pp_defc("unglq", HandleBad => 0, Pars => '[io,phys]A(2,m,n); [phys]tau(2,k); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code ' integer lwork = -1;
types(F) %{
extern int cunglq_(integer *m, integer *n, integer *k, float *
a, integer *lda, float *tau, float *work, integer *lwork,
integer *info);
float tmp_work[2];
%}
types(D) %{
extern int zunglq_(integer *m, integer *n, integer *k, double *
a, integer *lda, double *tau, double *work, integer *lwork,
integer *info);
double tmp_work[2];
%}
$TFD(cunglq_,zunglq_)(
&$PRIV(__m_size),
&$PRIV(__n_size),
&$PRIV(__k_size),
$P(A),
&$PRIV(__m_size),
$P(tau),
&tmp_work[0],
&lwork,
$P(info));
lwork = (integer )tmp_work[0];
{
types(F) %{
float *work = (float *)malloc(2 * lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2 * lwork * sizeof(double));
%}
$TFD(cunglq_,zunglq_)(
&$PRIV(__m_size),
&$PRIV(__n_size),
&$PRIV(__k_size),
$P(A),
&$PRIV(__m_size),
$P(tau),
work,
&lwork,
$P(info));
free(work);
}
', Doc=>'
Complex version of orglq
');
pp_defc("unmlq", HandleBad => 0, Pars => '[phys]A(2,k,p); int side(); int trans(); [phys]tau(2,k); [io,phys]C(2,m,n);int [o,phys]info()', GenericTypes => [F,D], Code => generate_code ' char ptrans = \'N\', pside = \'L\'; integer lwork = -1;
types(F) %{
extern int cunmlq_(char *side, char *trans, integer *m, integer *n,
integer *k, float *a, integer *lda, float *tau, float *
c__, integer *ldc, float *work, integer *lwork, integer *info);
float tmp_work[2];
%}
types(D) %{
extern int zunmlq_(char *side, char *trans, integer *m, integer *n,
integer *k, double *a, integer *lda, double *tau, double *
c__, integer *ldc, double *work, integer *lwork, integer *info);
double tmp_work[2];
%}
if($trans())
ptrans = \'C\';
if($side())
pside = \'R\';
$TFD(cunmlq_,zunmlq_)(
&pside,
&ptrans,
&$PRIV(__m_size),
&$PRIV(__n_size),
&$PRIV(__k_size),
$P(A),
&$PRIV(__k_size),
$P(tau),
$P(C),
&$PRIV(__m_size),
&tmp_work[0],
&lwork,
$P(info));
lwork = (integer )tmp_work[0];
{
types(F) %{
float *work = (float *)malloc(2 * lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2 * lwork * sizeof(double));
%}
$TFD(cunmlq_,zunmlq_)(
&pside,
&ptrans,
&$PRIV(__m_size),
&$PRIV(__n_size),
&$PRIV(__k_size),
$P(A),
&$PRIV(__k_size),
$P(tau),
$P(C),
&$PRIV(__m_size),
work,
&lwork,
$P(info));
free(work);
}
', Doc=>'
Complex version of ormlq. Here trans = 1 means conjugate transpose.
');
pp_defc("geqlf", HandleBad => 0, Pars => '[io,phys]A(2,m,n); [o,phys]tau(2,k); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code ' integer lwork = -1;
types(F) %{
extern int cgeqlf_(integer *m, integer *n, float *a, integer *
lda, float *tau, float *work, integer *lwork, integer *info);
float tmp_work[2];
%}
types(D) %{
extern int zgeqlf_(integer *m, integer *n, double *a, integer *
lda, double *tau, double *work, integer *lwork, integer *info);
double tmp_work[2];
%}
$TFD(cgeqlf_,zgeqlf_)(
&$PRIV(__m_size),
&$PRIV(__n_size),
$P(A),
&$PRIV(__m_size),
$P(tau),
&tmp_work[0],
&lwork,
$P(info));
lwork = (integer )tmp_work[0];
{
types(F) %{
float *work = (float *)malloc(2 * lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2 * lwork * sizeof(double));
%}
$TFD(cgeqlf_,zgeqlf_)(
&$PRIV(__m_size),
&$PRIV(__n_size),
$P(A),
&$PRIV(__m_size),
$P(tau),
work,
&lwork,
$P(info));
free(work);
}
');
pp_defc("ungql", HandleBad => 0, Pars => '[io,phys]A(2,m,n); [phys]tau(2,k); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code ' integer lwork = -1;
types(F) %{
extern int cungql_(integer *m, integer *n, integer *k, float *
a, integer *lda, float *tau, float *work, integer *lwork,
integer *info);
float tmp_work[2];
%}
types(D) %{
extern int zungql_(integer *m, integer *n, integer *k, double *
a, integer *lda, double *tau, double *work, integer *lwork,
integer *info);
double tmp_work[2];
%}
$TFD(cungql_,zungql_)(
&$PRIV(__m_size),
&$PRIV(__n_size),
&$PRIV(__k_size),
$P(A),
&$PRIV(__m_size),
$P(tau),
&tmp_work[0],
&lwork,
$P(info));
lwork = (integer )tmp_work[0];
{
types(F) %{
float *work = (float *)malloc(2 * lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2 * lwork * sizeof(double));
%}
$TFD(cungql_,zungql_)(
&$PRIV(__m_size),
&$PRIV(__n_size),
&$PRIV(__k_size),
$P(A),
&$PRIV(__m_size),
$P(tau),
work,
&lwork,
$P(info));
free(work);
}
', Doc=>'
Complex version of orgql.
');
pp_defc("unmql", HandleBad => 0, Pars => '[phys]A(2,p,k); int side(); int trans(); [phys]tau(2,k); [io,phys]C(2,m,n);int [o,phys]info()', GenericTypes => [F,D], Code => generate_code ' char ptrans = \'N\', pside = \'L\'; integer lwork = -1;
types(F) %{
extern int cunmql_(char *side, char *trans, integer *m, integer *n,
integer *k, float *a, integer *lda, float *tau, float *
c__, integer *ldc, float *work, integer *lwork, integer *info);
float tmp_work[2];
%}
types(D) %{
extern int zunmql_(char *side, char *trans, integer *m, integer *n,
integer *k, double *a, integer *lda, double *tau, double *
c__, integer *ldc, double *work, integer *lwork, integer *info);
double tmp_work[2];
%}
if($trans())
ptrans = \'C\';
if($side())
pside = \'R\';
$TFD(cunmql_,zunmql_)(
&pside,
&ptrans,
&$PRIV(__m_size),
&$PRIV(__n_size),
&$PRIV(__k_size),
$P(A),
&$PRIV(__p_size),
$P(tau),
$P(C),
&$PRIV(__m_size),
&tmp_work[0],
&lwork,
$P(info));
lwork = (integer )tmp_work[0];
{
types(F) %{
float *work = (float *)malloc(2 * lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2 * lwork * sizeof(double));
%}
$TFD(cunmql_,zunmql_)(
&pside,
&ptrans,
&$PRIV(__m_size),
&$PRIV(__n_size),
&$PRIV(__k_size),
$P(A),
&$PRIV(__p_size),
$P(tau),
$P(C),
&$PRIV(__m_size),
work,
&lwork,
$P(info));
free(work);
}
', Doc=>'
Complex version of ormql. Here trans = 1 means conjugate transpose.
');
pp_defc("gerqf", HandleBad => 0, Pars => '[io,phys]A(2,m,n); [o,phys]tau(2,k); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code ' integer lwork = -1;
types(F) %{
extern int cgerqf_(integer *m, integer *n, float *a, integer *
lda, float *tau, float *work, integer *lwork, integer *info);
float tmp_work[2];
%}
types(D) %{
extern int zgerqf_(integer *m, integer *n, double *a, integer *
lda, double *tau, double *work, integer *lwork, integer *info);
double tmp_work[2];
%}
$TFD(cgerqf_,zgerqf_)(
&$PRIV(__m_size),
&$PRIV(__n_size),
$P(A),
&$PRIV(__m_size),
$P(tau),
&tmp_work[0],
&lwork,
$P(info));
lwork = (integer )tmp_work[0];
{
types(F) %{
float *work = (float *)malloc(2 * lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2 * lwork * sizeof(double));
%}
$TFD(cgerqf_,zgerqf_)(
&$PRIV(__m_size),
&$PRIV(__n_size),
$P(A),
&$PRIV(__m_size),
$P(tau),
work,
&lwork,
$P(info));
free(work);
}
');
pp_defc("ungrq", HandleBad => 0, Pars => '[io,phys]A(2,m,n); [phys]tau(2,k); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code ' integer lwork = -1;
types(F) %{
extern int cungrq_(integer *m, integer *n, integer *k, float *
a, integer *lda, float *tau, float *work, integer *lwork,
integer *info);
float tmp_work[2];
%}
types(D) %{
extern int zungrq_(integer *m, integer *n, integer *k, double *
a, integer *lda, double *tau, double *work, integer *lwork,
integer *info);
double tmp_work[2];
%}
$TFD(cungrq_,zungrq_)(
&$PRIV(__m_size),
&$PRIV(__n_size),
&$PRIV(__k_size),
$P(A),
&$PRIV(__m_size),
$P(tau),
&tmp_work[0],
&lwork,
$P(info));
lwork = (integer )tmp_work[0];
{
types(F) %{
float *work = (float *)malloc(2 * lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2 * lwork * sizeof(double));
%}
$TFD(cungrq_,zungrq_)(
&$PRIV(__m_size),
&$PRIV(__n_size),
&$PRIV(__k_size),
$P(A),
&$PRIV(__m_size),
$P(tau),
work,
&lwork,
$P(info));
free(work);
}
', Doc=>'
Complex version of orgrq.
');
pp_defc("unmrq", HandleBad => 0, Pars => '[phys]A(2,k,p); int side(); int trans(); [phys]tau(2,k); [io,phys]C(2,m,n);int [o,phys]info()', GenericTypes => [F,D], Code => generate_code ' char ptrans = \'N\', pside = \'L\'; integer lwork = -1;
types(F) %{
extern int cunmrq_(char *side, char *trans, integer *m, integer *n,
integer *k, float *a, integer *lda, float *tau, float *
c__, integer *ldc, float *work, integer *lwork, integer *info);
float tmp_work[2];
%}
types(D) %{
extern int zunmrq_(char *side, char *trans, integer *m, integer *n,
integer *k, double *a, integer *lda, double *tau, double *
c__, integer *ldc, double *work, integer *lwork, integer *info);
double tmp_work[2];
%}
if($trans())
ptrans = \'C\';
if($side())
pside = \'R\';
$TFD(cunmrq_,zunmrq_)(
&pside,
&ptrans,
&$PRIV(__m_size),
&$PRIV(__n_size),
&$PRIV(__k_size),
$P(A),
&$PRIV(__k_size),
$P(tau),
$P(C),
&$PRIV(__m_size),
&tmp_work[0],
&lwork,
$P(info));
lwork = (integer )tmp_work[0];
{
types(F) %{
float *work = (float *)malloc(2 * lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2 * lwork * sizeof(double));
%}
$TFD(cunmrq_,zunmrq_)(
&pside,
&ptrans,
&$PRIV(__m_size),
&$PRIV(__n_size),
&$PRIV(__k_size),
$P(A),
&$PRIV(__k_size),
$P(tau),
$P(C),
&$PRIV(__m_size),
work,
&lwork,
$P(info));
free(work);
}
', Doc=>'
Complex version of ormrq. Here trans = 1 means conjugate transpose.
');
pp_defc("tzrzf", HandleBad => 0, Pars => '[io,phys]A(2,m,n); [o,phys]tau(2,k); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code ' integer lwork = -1;
types(F) %{
extern int ctzrzf_(integer *m, integer *n, float *a, integer *
lda, float *tau, float *work, integer *lwork, integer *info);
float tmp_work[2];
%}
types(D) %{
extern int ztzrzf_(integer *m, integer *n, double *a, integer *
lda, double *tau, double *work, integer *lwork, integer *info);
double tmp_work[2];
%}
$TFD(ctzrzf_,ztzrzf_)(
&$PRIV(__m_size),
&$PRIV(__n_size),
$P(A),
&$PRIV(__m_size),
$P(tau),
&tmp_work[0],
&lwork,
$P(info));
lwork = (integer )tmp_work[0];
{
types(F) %{
float *work = (float *)malloc(2 * lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2 * lwork * sizeof(double));
%}
$TFD(ctzrzf_,ztzrzf_)(
&$PRIV(__m_size),
&$PRIV(__n_size),
$P(A),
&$PRIV(__m_size),
$P(tau),
work,
&lwork,
$P(info));
free(work);
}
');
pp_defc("unmrz", HandleBad => 0, Pars => '[phys]A(2,k,p); int side(); int trans(); [phys]tau(2,k); [io,phys]C(2,m,n);int [o,phys]info()', GenericTypes => [F,D], Code => generate_code ' char ptrans = \'N\', pside = \'L\'; integer lwork = -1; integer kk = $SIZE(p) - $SIZE(k);
types(F) %{
extern int cunmrz_(char *side, char *trans, integer *m, integer *n,
integer *k, integer *l, float *a, integer *lda, float *tau, float *
c__, integer *ldc, float *work, integer *lwork, integer *info);
float tmp_work[2];
%}
types(D) %{
extern int zunmrz_(char *side, char *trans, integer *m, integer *n,
integer *k, integer *l, double *a, integer *lda, double *tau, double *
c__, integer *ldc, double *work, integer *lwork, integer *info);
double tmp_work[2];
%}
if($trans())
ptrans = \'C\';
if($side())
pside = \'R\';
$TFD(cunmrz_,zunmrz_)(
&pside,
&ptrans,
&$PRIV(__m_size),
&$PRIV(__n_size),
&$PRIV(__k_size),
&kk,
$P(A),
&$PRIV(__k_size),
$P(tau),
$P(C),
&$PRIV(__m_size),
&tmp_work[0],
&lwork,
$P(info));
lwork = (integer )tmp_work[0];
{
types(F) %{
float *work = (float *)malloc(2 * lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2 * lwork * sizeof(double));
%}
$TFD(cunmrz_,zunmrz_)(
&pside,
&ptrans,
&$PRIV(__m_size),
&$PRIV(__n_size),
&$PRIV(__k_size),
&kk,
$P(A),
&$PRIV(__k_size),
$P(tau),
$P(C),
&$PRIV(__m_size),
work,
&lwork,
$P(info));
free(work);
}
', Doc=>'
Complex version of ormrz. Here trans = 1 means conjugate transpose.
');
pp_defc("gehrd", HandleBad => 0, Pars => '[io,phys]A(2,n,n); int [phys]ilo();int [phys]ihi();[o,phys]tau(2,k); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code '
integer lwork = -1;
types(F) %{
extern int cgehrd_(integer *n, integer *ilo, integer *ihi,
float *a, integer *lda, float *tau, float *work,
integer *lwork, integer *info);
float tmp_work[2];
%}
types(D) %{
extern int zgehrd_(integer *n, integer *ilo, integer *ihi,
double *a, integer *lda, double *tau, double *work,
integer *lwork, integer *info);
double tmp_work[2];
%}
$TFD(cgehrd_,zgehrd_)(
&$PRIV(__n_size),
$P(ilo),
$P(ihi),
$P(A),
&$PRIV(__n_size),
$P(tau),
&tmp_work[0],
&lwork,
$P(info));
lwork = (integer )tmp_work[0];
{
types(F) %{
float *work = (float *)malloc(2*lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2*lwork * sizeof(double));
%}
$TFD(cgehrd_,zgehrd_)(
&$PRIV(__n_size),
$P(ilo),
$P(ihi),
$P(A),
&$PRIV(__n_size),
$P(tau),
work,
&lwork,
$P(info));
free(work);
}
');
pp_defc("unghr", HandleBad => 0, Pars => '[io,phys]A(2,n,n); int [phys]ilo();int [phys]ihi();[phys]tau(2,k); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code '
integer lwork = -1;
types(F) %{
extern int cunghr_(integer *n, integer *ilo, integer *ihi,
float *a, integer *lda, float *tau, float *work,
integer *lwork, integer *info);
float tmp_work[2];
%}
types(D) %{
extern int zunghr_(integer *n, integer *ilo, integer *ihi,
double *a, integer *lda, double *tau, double *work,
integer *lwork, integer *info);
double tmp_work[2];
%}
$TFD(cunghr_,zunghr_)(
&$PRIV(__n_size),
$P(ilo),
$P(ihi),
$P(A),
&$PRIV(__n_size),
$P(tau),
&tmp_work[0],
&lwork,
$P(info));
lwork = (integer )tmp_work[0];
{
types(F) %{
float *work = (float *)malloc(2*lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2*lwork * sizeof(double));
%}
$TFD(cunghr_,zunghr_)(
&$PRIV(__n_size),
$P(ilo),
$P(ihi),
$P(A),
&$PRIV(__n_size),
$P(tau),
work,
&lwork,
$P(info));
free(work);
}
', Doc=>'
Complex version of orghr
');
pp_defc("hseqr", HandleBad => 0, Pars => '[io,phys]H(2,n,n); int job();int compz();int [phys]ilo();int [phys]ihi();[o,phys]w(2,n); [o,phys]Z(2,m,m); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code '
char pcompz;
char pjob = \'E\';
integer lwork = -1;
types(F) %{
extern int chseqr_(char *job, char *compz, integer *n, integer *ilo,
integer *ihi, float *h__, integer *ldh, float *w,
float *z__, integer *ldz, float *work,
integer *lwork, integer *info);
float tmp_work[2];
%}
types(D) %{
extern int zhseqr_(char *job, char *compz, integer *n, integer *ilo,
integer *ihi, double *h__, integer *ldh, double *w,
double *z__, integer *ldz, double *work,
integer *lwork, integer *info);
double tmp_work[2];
%}
if($job())
pjob = \'S\';
switch ($compz())
{
case 1: pcompz = \'I\';
break;
case 2: pcompz = \'V\';
break;
default: pcompz = \'N\';
}
$TFD(chseqr_,zhseqr_)(
&pjob,
&pcompz,
&$PRIV(__n_size),
$P(ilo),
$P(ihi),
$P(H),
&$PRIV(__n_size),
$P(w),
$P(Z),
&$PRIV(__m_size),
&tmp_work[0],
&lwork,
$P(info));
lwork = (integer )tmp_work[0];
{
types(F) %{
float *work = (float *)malloc(2 * lwork * sizeof(float));
%}
types(D) %{
double *work = (double *)malloc(2 * lwork * sizeof(double));
%}
$TFD(chseqr_,zhseqr_)(
&pjob,
&pcompz,
&$PRIV(__n_size),
$P(ilo),
$P(ihi),
$P(H),
&$PRIV(__n_size),
$P(w),
$P(Z),
&$PRIV(__m_size),
work,
&lwork,
$P(info));
free(work);
}
');
pp_defc("trevc", HandleBad => 0, Pars => '[io,phys]T(2,n,n); int side();int howmny();int [phys]select(q);[io,phys]VL(2,m,r); [io,phys]VR(2,p,s);int [o,phys]m(); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code '
char pside,phowmny;
integer mm = 0;
types(F) %{
extern int ctrevc_(char *side, char *howmny, logical *select,
integer *n, float *t, integer *ldt, float *vl, integer *
ldvl, float *vr, integer *ldvr, integer *mm, integer *m,
float *work, float *rwork, integer *info);
float *work = (float *) malloc(5 * $SIZE(n) *sizeof(float));
%}
types(D) %{
extern int ztrevc_(char *side, char *howmny, logical *select,
integer *n, double *t, integer *ldt, double *vl, integer *
ldvl, double *vr, integer *ldvr, integer *mm, integer *m,
double *work, double *rwork, integer *info);
double *work = (double *) malloc (5 * $SIZE(n) * sizeof(double));
%}
switch ($howmny())
{
case 1: phowmny = \'B\';
break;
case 2: phowmny = \'S\';
break;
default: phowmny = \'A\';
}
switch ($side())
{
case 1: pside = \'R\';
mm = $SIZE(s);
break;
case 2: pside = \'L\';
mm = $SIZE(r);
break;
default:pside = \'B\';
mm = $SIZE(s);
}
$TFD(ctrevc_,ztrevc_)(
&pside,
&phowmny,
$P(select),
&$PRIV(__n_size),
$P(T),
&$PRIV(__n_size),
$P(VL),
&$PRIV(__m_size),
$P(VR),
&$PRIV(__p_size),
&mm,
$P(m),
&work[$SIZE(n)],
work,
$P(info));
free(work);
');
pp_defc("tgevc", HandleBad => 0, Pars => '[io,phys]A(2,n,n); int side();int howmny(); [io,phys]B(2,n,n);int [phys]select(q);[io,phys]VL(2,m,r); [io,phys]VR(2,p,s);int [o,phys]m(); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code '
char pside,phowmny;
integer mm = 0;
types(F) %{
extern int ctgevc_(char *side, char *howmny, logical *select,
integer *n, float *a, integer *lda, float *b, integer *ldb,
float *vl, integer *ldvl, float *vr, integer *ldvr,
integer *mm, integer *m, float *work, float *rwork, integer *info);
float *work = (float *) malloc(6 * $SIZE(n) *sizeof(float));
%}
types(D) %{
extern int ztgevc_(char *side, char *howmny, logical *select,
integer *n, double *a, integer *lda, double *b, integer *ldb,
double *vl, integer *ldvl, double *vr, integer *ldvr,
integer *mm, integer *m, double *work, double *rwork, integer *info);
double *work = (double *) malloc (6 * $SIZE(n) * sizeof(double));
%}
switch ($howmny())
{
case 1: phowmny = \'B\';
break;
case 2: phowmny = \'S\';
break;
default: phowmny = \'A\';
}
switch ($side())
{
case 1: pside = \'R\';
mm = $SIZE(s);
break;
case 2: pside = \'L\';
mm = $SIZE(r);
break;
default:pside = \'B\';
mm = $SIZE(s);
}
$TFD(ctgevc_,ztgevc_)(
&pside,
&phowmny,
$P(select),
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(B),
&$PRIV(__n_size),
$P(VL),
&$PRIV(__m_size),
$P(VR),
&$PRIV(__p_size),
&mm,
$P(m),
&work[2*$SIZE(n)],
work,
$P(info));
free(work);
');
pp_defc("gebal", HandleBad => 0, Pars => '[io,phys]A(2,n,n); int job(); int [o,phys]ilo();int [o,phys]ihi();[o,phys]scale(n); int [o,phys]info()', GenericTypes => [F,D], Code => generate_code '
char pjob;
types(F) %{
extern int cgebal_(char *job, integer *n, float *a, integer *
lda, integer *ilo, integer *ihi, float *scale, integer *info);
%}
types(D) %{
extern int zgebal_(char *job, integer *n, double *a, integer *
lda, integer *ilo, integer *ihi, double *scale, integer *info);
%}
switch ($job())
{
case 1: pjob = \'P\';
break;
case 2: pjob = \'S\';
break;
case 3: pjob = \'B\';
break;
default: pjob = \'N\';
}
$TFD(cgebal_,zgebal_)(
&pjob,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
$P(ilo),
$P(ihi),
$P(scale),
$P(info));
');
################################################################################# pp_defc("lange", HandleBad => 0, Pars => '[phys]A(2,n,m); int norm(); [o]b()', GenericTypes => [F,D], Code => ' char pnorm;
types(F) %{
extern float clange_(char *norm, integer *m, integer *n, float *a, integer
*lda, float *work);
float *work;
%}
types(D) %{
extern double zlange_(char *norm, integer *m, integer *n, double *a, integer
*lda, double *work);
double *work;
%}
switch ($norm())
{
case 1: pnorm = \'O\';
break;
case 2: pnorm = \'I\';
types(F) %{
work = (float *)malloc($PRIV(__n_size) * sizeof(float));
%}
types(D) %{
work = (double *)malloc($PRIV(__n_size) * sizeof(double));
%}
break;
case 3: pnorm = \'F\';
break;
default: pnorm = \'M\';
}
$b() = $TFD(clange_,zlange_)(
&pnorm,
&$PRIV(__n_size),
&$PRIV(__m_size),
$P(A),
&$PRIV(__n_size),
work);
if ($norm() == 2)
free (work);
');
pp_defc("lansy", HandleBad => 0, Pars => '[phys]A(2, n,n); int uplo(); int norm(); [o]b()', GenericTypes => [F,D], Code => '
char pnorm, puplo = \'U\';
types(F) %{
extern float clansy_(char *norm, char *uplo, integer *n, float *a, integer
*lda, float *work);
float *work;
%}
types(D) %{
extern double zlansy_(char *norm, char *uplo, integer *n, double *a, integer
*lda, double *work);
double *work;
%}
switch ($norm())
{
case 1: pnorm = \'O\';
types(F) %{
work = (float *)malloc($PRIV(__n_size) * sizeof(float));
%}
types(D) %{
work = (double *)malloc($PRIV(__n_size) * sizeof(double));
%}
break;
case 2: pnorm = \'I\';
types(F) %{
work = (float *)malloc($PRIV(__n_size) * sizeof(float));
%}
types(D) %{
work = (double *)malloc($PRIV(__n_size) * sizeof(double));
%}
break;
case 3: pnorm = \'F\';
break;
default: pnorm = \'M\';
}
if($uplo())
puplo = \'L\';
$b() = $TFD(clansy_,zlansy_)(
&pnorm,
&puplo,
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
work);
if ($norm() == 2 || $norm() == 1)
free (work);
');
pp_defc("lantr", HandleBad => 0, Pars => '[phys]A(2,m,n);int uplo();int norm();int diag();[o]b()', GenericTypes => [F,D], Code => '
char pnorm, puplo = \'U\';
char pdiag = \'N\';
types(F) %{
extern float clantr_(char *norm, char *uplo, char *diag, integer *m, integer *n, float *a, integer
*lda, float *work);
float *work;
%}
types(D) %{
extern double zlantr_(char *norm, char *uplo, char *diag, integer *m, integer *n, double *a, integer
*lda, double *work);
double *work;
%}
switch ($norm())
{
case 1: pnorm = \'O\';
break;
case 2: pnorm = \'I\';
types(F) %{
work = (float *)malloc($PRIV(__m_size) * sizeof(float));
%}
types(D) %{
work = (double *)malloc($PRIV(__m_size) * sizeof(double));
%}
break;
case 3: pnorm = \'F\';
break;
default: pnorm = \'M\';
}
if($uplo())
puplo = \'L\';
if($diag())
pdiag = \'U\';
$b() = $TFD(clantr_,zlantr_)(
&pnorm,
&puplo,
&pdiag,
&$PRIV(__m_size),
&$PRIV(__n_size),
$P(A),
&$PRIV(__n_size),
work);
if ($norm() == 2)
free (work);
');
################################################################################ # # BLAS ROUTINES # ################################################################################
pp_defc("gemm", HandleBad => 0, Pars => '[phys]A(2,m,n); int transa(); int transb(); [phys]B(2,p,q);[phys]alpha(2); [phys]beta(2); [io,phys]C(2,r,s)', GenericTypes => [F,D], Code => ' char ptransa = \'N\'; char ptransb = \'N\';
types(F) %{
extern int cgemm_(char *transa, char *transb, integer *m, integer *
n, integer *k, float *alpha, float *a, integer *lda,
float *b, integer *ldb, float *beta, float *c__,
integer *ldc);
%}
types(D) %{
extern int zgemm_(char *transa, char *transb, integer *m, integer *
n, integer *k, double *alpha, double *a, integer *lda,
double *b, integer *ldb, double *beta, double *c__,
integer *ldc);
%}
integer kk = $transa() ? $SIZE(m) : $SIZE(n);
if ($transa() == 1)
ptransa = \'T\';
else if ($transa() == 2)
ptransa = \'C\';
if ($transb())
ptransb = \'T\';
else if ($transb() == 2)
ptransb = \'C\';
$TFD(cgemm_,zgemm_)(
&ptransa,
&ptransb,
&$PRIV(__r_size),
&$PRIV(__s_size),
&kk,
$P(alpha),
$P(A),
&$PRIV(__m_size),
$P(B),
&$PRIV(__p_size),
$P(beta),
$P(C),
&$PRIV(__r_size));
',
Doc=>'
Complex version of gemm.
Arguments
=========
transa: = 0: No transpose;
= 1: Transpose;
= 2: Conjugate transpose;
transb: = 0: No transpose;
= 1: Transpose;
= 2: Conjugate transpose;
');
if ($config{CBLAS}){ pp_def("rmcgemm", HandleBad => 0, Pars => '[phys]A(2,m,n); int transa(); int transb(); [phys]B(2,p,q);[phys]alpha(2); [phys]beta(2); [io,phys]C(2,r,s)', GenericTypes => [F,D], Code => ' int ptransa, ptransb; types(F) %{ extern void cblas_cgemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, const int M, const int N, const int K, const void *alpha, const void *A, const int lda, const void *B, const int ldb, const void *beta, void *C, const int ldc); %} types(D) %{ extern void cblas_zgemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, const int M, const int N, const int K, const void *alpha, const void *A, const int lda, const void *B, const int ldb, const void *beta, void *C, const int ldc); %} integer kk = $transa() ? $SIZE(n) : $SIZE(m);
switch($transa()){
case 1: ptransa = CblasTrans;
break;
case 2: ptransa = CblasConjTrans;
break;
default:ptransa = CblasNoTrans;
}
switch($transb()){
case 1: ptransb = CblasTrans;
break;
case 2: ptransb = CblasConjTrans;
break;
default:ptransb = CblasNoTrans;
}
$TFD(cblas_cgemm,cblas_zgemm)(
CblasRowMajor,
ptransa,
ptransb,
$PRIV(__s_size),
$PRIV(__r_size),
kk,
$P(alpha),
$P(A),
$PRIV(__m_size),
$P(B),
$PRIV(__p_size),
$P(beta),
$P(C),
$PRIV(__r_size));
',
Doc=>'
Complex version of rmgemm.
Arguments
=========
transa: = 0: No transpose;
= 1: Transpose;
= 2: Conjugate transpose;
transb: = 0: No transpose;
= 1: Transpose;
= 2: Conjugate transpose;
'); }
pp_defc("mmult", HandleBad => 0, Pars => '[phys]A(2,m,n); [phys]B(2,p,m); [o,phys]C(2,p,n)', GenericTypes => [F,D], Code => ' char ptrans = \'N\'; types(F) %{ extern int cgemm_(char *transa, char *transb, integer *m, integer * n, integer *k, float *alpha, float *a, integer *lda, float *b, integer *ldb, float *beta, float *c__, integer *ldc); float alpha[2] = {1,0}; float beta[2] = {0,0}; %} types(D) %{ extern int zgemm_(char *transa, char *transb, integer *m, integer * n, integer *k, double *alpha, double *a, integer *lda, double *b, integer *ldb, double *beta, double *c__, integer *ldc); double alpha[2] = {1,0}; double beta[2] = {0,0}; %}
$TFD(cgemm_,zgemm_)(
&ptrans,
&ptrans,
&$PRIV(__p_size),
&$PRIV(__n_size),
&$PRIV(__m_size),
&alpha[0],
$P(B),
&$PRIV(__p_size),
$P(A),
&$PRIV(__m_size),
&beta[0],
$P(C),
&$PRIV(__p_size));
');
if ($config{STRASSEN}){
pp_defc("smmult",
HandleBad => 0,
Pars => '[phys]A(2,m,n); [phys]B(2,p,m); [o,phys]C(2,p,n)',
GenericTypes => [F,D],
Code => '
char ptrans = \'N\';
types(F) %{
extern int cgemmb_(char *transa, char *transb, integer *m, integer *
n, integer *k, float *alpha, float *a, integer *lda,
float *b, integer *ldb, float *beta, float *c__,
integer *ldc);
float alpha[2] = {1,0};
float beta[2] = {0,0};
%}
types(D) %{
extern int zgemmb_(char *transa, char *transb, integer *m, integer *
n, integer *k, double *alpha, double *a, integer *lda,
double *b, integer *ldb, double *beta, double *c__,
integer *ldc);
double alpha[2] = {1,0};
double beta[2] = {0,0};
%}
$TFD(cgemmb_,zgemmb_)(
&ptrans,
&ptrans,
&$PRIV(__p_size),
&$PRIV(__n_size),
&$PRIV(__m_size),
&alpha[0],
$P(B),
&$PRIV(__p_size),
$P(A),
&$PRIV(__m_size),
&beta[0],
$P(C),
&$PRIV(__p_size));
');
}
pp_defc("crossprod", HandleBad => 0, Pars => '[phys]A(2,n,m); [phys]B(2,p,m); [o,phys]C(2,p,n)', GenericTypes => [F,D], Code => ' char btrans = \'N\'; char atrans = \'C\'; types(F) %{ extern int cgemm_(char *transa, char *transb, integer *m, integer * n, integer *k, float *alpha, float *a, integer *lda, float *b, integer *ldb, float *beta, float *c__, integer *ldc); float alpha[2] = {1,0}; float beta[2] = {0,0}; %} types(D) %{ extern int zgemm_(char *transa, char *transb, integer *m, integer * n, integer *k, double *alpha, double *a, integer *lda, double *b, integer *ldb, double *beta, double *c__, integer *ldc); double alpha[2] = {1,0}; double beta[2] = {0,0}; %}
$TFD(cgemm_,zgemm_)(
&btrans,
&atrans,
&$PRIV(__p_size),
&$PRIV(__n_size),
&$PRIV(__m_size),
&alpha[0],
$P(B),
&$PRIV(__p_size),
$P(A),
&$PRIV(__n_size),
&beta[0],
$P(C),
&$PRIV(__p_size));
');
pp_defc("syrk", HandleBad => 0, Pars => '[phys]A(2,m,n); int uplo(); int trans(); [phys]alpha(2); [phys]beta(2); [io,phys]C(2,p,p)', RedoDimsCode => '$SIZE(p) = $trans() ? $SIZE(n) : $SIZE(m);', GenericTypes => [F,D], Code => '
char puplo = \'U\';
char ptrans = \'N\';
types(F) %{
extern int csyrk_(char *uplo, char *trans, integer *n, integer *k,
float *alpha, float *a, integer *lda, float *beta,
float *c__, integer *ldc);
%}
types(D) %{
extern int zsyrk_(char *uplo, char *trans, integer *n, integer *k,
double *alpha, double *a, integer *lda, double *beta,
double *c__, integer *ldc);
%}
integer kk = $trans() ? $SIZE(m) : $SIZE(n);
if ($uplo())
puplo = \'L\';
if ($trans())
ptrans = \'T\';
$TFD(csyrk_,zsyrk_)(
&puplo,
&ptrans,
&$PRIV(__p_size),
&kk,
$P(alpha),
$P(A),
&$PRIV(__m_size),
$P(beta),
$P(C),
&$PRIV(__p_size));
');
if ($config{CBLAS}){ pp_def("rmcsyrk", HandleBad => 0, Pars => '[phys]A(2,m,n); int uplo(); int trans(); [phys]alpha(2); [phys]beta(2); [io,phys]C(2,p,p)', GenericTypes => [F,D], Code => ' int puplo = CblasUpper; int ptrans;
types(F) %{
extern void cblas_csyrk(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo,
const enum CBLAS_TRANSPOSE Trans, const int N, const int K,
const void *alpha, const void *A, const int lda,
const void *beta, void *C, const int ldc);
%}
types(D) %{
extern void cblas_zsyrk(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo,
const enum CBLAS_TRANSPOSE Trans, const int N, const int K,
const void *alpha, const void *A, const int lda,
const void *beta, void *C, const int ldc);
%}
integer kk = $trans() ? $SIZE(n) : $SIZE(m);
if ($uplo())
puplo = CblasLower;
switch($trans()){
case 1: ptrans = CblasTrans;
break;
case 2: ptrans = CblasConjTrans;
break;
default:ptrans = CblasNoTrans;
}
$TFD(cblas_csyrk,cblas_zsyrk)(
CblasRowMajor,
puplo,
ptrans,
$PRIV(__p_size),
kk,
$P(alpha),
$P(A),
$PRIV(__m_size),
$P(beta),
$P(C),
$PRIV(__p_size));
',
Doc=>'
Complex version of rmsyrk
'); } pp_defc("dot", HandleBad => 0, Pars => '[phys]a(2,n);int [phys]inca();[phys]b(2,n);int [phys]incb();[o,phys]c(2)', GenericTypes => [F,D], Code => ' types(F) %{ extern float cdotu_(float *ret, integer *n, float *dx, integer *incx, float *dy, integer *incy); %} types(D) %{ extern double zdotu_(double *ret, integer *n, double *dx, integer *incx, double *dy, integer *incy); %} integer n = (integer ) $PRIV(__n_size)/$inca();
$TFD(cdotu_,zdotu_)(
$P(c),
&n,
$P(a),
$P(inca),
$P(b),
$P(incb));
');
pp_def("cdotc", HandleBad => 0, Pars => '[phys]a(2,n);int [phys]inca();[phys]b(2,n);int [phys]incb();[o,phys]c(2)', GenericTypes => [F,D], Code => ' types(F) %{ extern float cdotc_(float *ret, integer *n, float *dx, integer *incx, float *dy, integer *incy); %} types(D) %{ extern double zdotc_(double *ret, integer *n, double *dx, integer *incx, double *dy, integer *incy); %} integer n = (integer ) $PRIV(__n_size)/$inca();
$TFD(cdotc_,zdotc_)(
$P(c),
&n,
$P(a),
$P(inca),
$P(b),
$P(incb));
',
Doc=>'
Forms the dot product of two vectors, conjugating the first vector.
');
pp_defc("axpy", HandleBad => 0, Pars => '[phys]a(2,n);int [phys]inca();[phys] alpha(2);[io,phys]b(2,n);int [phys]incb()', GenericTypes => [F,D], Code => '
types(F) %{
extern int caxpy_(integer *n, float *da, float *dx,
integer *incx, float *dy, integer *incy);
%}
types(D) %{
extern int zaxpy_(integer *n, double *da, double *dx,
integer *incx, double *dy, integer *incy);
%}
integer n = (integer ) $PRIV(__n_size)/$inca();
$TFD(caxpy_,zaxpy_)(
&n,
$P(alpha),
$P(a),
$P(inca),
$P(b),
$P(incb));
');
pp_defc("nrm2", HandleBad => 0, Pars => '[phys]a(2,n);int [phys]inca();[o,phys]b()', GenericTypes => [F,D], Code => '
types(F) %{
extern float scnrm2_(integer *n, float *dx,
integer *incx);
%}
types(D) %{
extern double dznrm2_(integer *n, double *dx,
integer *incx);
%}
integer n = (integer ) $PRIV(__n_size)/$inca();
$b() = $TFD(scnrm2_,dznrm2_)(
&n,
$P(a),
$P(inca));
');
pp_defc("asum", HandleBad => 0, Pars => '[phys]a(2,n);int [phys]inca();[o,phys]b()', GenericTypes => [F,D], Code => '
types(F) %{
extern float scasum_(integer *n, float *dx,
integer *incx);
%}
types(D) %{
extern double dzasum_(integer *n, double *dx,
integer *incx);
%}
integer n = (integer ) $PRIV(__n_size)/$inca();
$b() = $TFD(scasum_,dzasum_)(
&n,
$P(a),
$P(inca));
');
pp_defc("scal", HandleBad => 0, Pars => '[io,phys]a(2,n);int [phys]inca();[phys]scale(2)', GenericTypes => [F,D], Code => '
types(F) %{
extern int cscal_(integer *n, float *sa,
float *dx, integer *incx);
%}
types(D) %{
extern int zscal_(integer *n, double *sa,
double *dx,integer *incx);
%}
integer n = (integer ) $PRIV(__n_size)/$inca();
$TFD(cscal_,zscal_)(
&n,
$P(scale),
$P(a),
$P(inca));
');
pp_def("sscal", HandleBad => 0, Pars => '[io,phys]a(2,n);int [phys]inca();[phys]scale()', GenericTypes => [F], Code => '
extern int csscal_(integer *n, float *sa,
float *dx, integer *incx);
integer n = (integer ) $PRIV(__n_size)/$inca();
csscal_( &n, $P(scale), $P(a), $P(inca));
',
Doc=>'
Scales a complex vector by a real constant.
');
pp_defc("rotg", HandleBad => 0, Pars => '[io,phys]a(2);[phys]b(2);[o,phys]c(); [o,phys]s(2)', GenericTypes => [F,D], Code => '
types(F) %{
extern int crotg_(float *dx,
float *dy, float *c, float *s);
%}
types(D) %{
extern int zrotg_(double *dx,
double *dy, double *c, double *s);
%}
$TFD(crotg_,zrotg_)(
$P(a),
$P(b),
$P(c),
$P(s)
);
');
################################################################################ # # LAPACK AUXILIARY ROUTINES # ################################################################################ pp_defc("lacpy", HandleBad => 0, Pars => '[phys]A(2,m,n); int uplo(); [o,phys]B(2,p,n)', GenericTypes => [F,D], Code => ' char puplo;
types(F) %{
extern int clacpy_(char *uplo, integer *m, integer *n, float *
a, integer *lda, float *b, integer *ldb);
%}
types(D) %{
extern int zlacpy_(char *uplo, integer *m, integer *n, double *
a, integer *lda, double *b, integer *ldb);
%}
switch ($uplo())
{
case 0: puplo = \'U\';
break;
case 1: puplo = \'L\';
break;
default: puplo = \'A\';
}
$TFD(clacpy_,zlacpy_)(
&puplo,
&$PRIV(__m_size),
&$PRIV(__n_size),
$P(A),
&$PRIV(__m_size),
$P(B),
&$PRIV(__p_size));
');
pp_defc("laswp", HandleBad => 0, Pars => '[io,phys]A(2,m,n); int [phys]k1(); int [phys]k2(); int [phys]ipiv(p);int [phys]inc()', GenericTypes => [F,D], Code => ' types(F) %{ extern int claswp_(integer *n, float *a, integer *lda, integer *k1, integer *k2, integer *ipiv, integer *incx); %} types(D) %{ extern int zlaswp_(integer *n, double *a, integer *lda, integer *k1, integer *k2, integer *ipiv, integer *incx); %}
$TFD(claswp_,zlaswp_)(
&$PRIV(__n_size),
$P(A),
&$PRIV(__m_size),
$P(k1),
$P(k2),
$P(ipiv),
$P(inc));
');
################################################################################ # # OTHER AUXILIARY ROUTINES # ################################################################################ pp_def( 'ctricpy', Pars => 'A(c=2,m,n);int uplo();[o] C(c=2,m,n)', Code => ' PDL_Long i, j, k;
if ($uplo())
{
for (i = 0; i < $SIZE(n);i++)
{
k = min(i,($SIZE(m)-1));
for (j = 0; j <= k; j++)
{
$C(c=>0,m=>j,n=>i) = $A(c=>0,m=>j,n=>i);
$C(c=>1,m=>j,n=>i) = $A(c=>1,m=>j,n=>i);
}
}
}
else
{
for (i = 0; i < $SIZE(n);i++)
{
for (j = i; j < $SIZE(m); j++)
{
$C(c=>0,m=>j,n=>i) = $A(c=>0,m=>j,n=>i);
$C(c=>1,m=>j,n=>i) = $A(c=>1,m=>j,n=>i);
}
if (i >= $SIZE(m))
break;
}
}
',
Doc => <<EOT
Copy triangular part to another matrix. If uplo == 0 copy upper triangular part.
Combine two 3D piddles into a single piddle. This routine does backward and forward dataflow automatically.
Copyright (C) Grégory Vanuxem 2005-2007.
This library is free software; you can redistribute it and/or modify it under the terms of the artistic license as specified in the Artistic file.
