{ # required because PAUSE will reject undefined $VERSION
require ExtUtils::MM;
my $MM = bless { NAME => 'Fake' }, 'MM';
my $global_version = $MM->parse_version('lib/PDL/LinearAlgebra.pm');
pp_setversion($global_version);
}
{ no warnings 'once'; # pass info back to Makefile.PL
$PDL::Core::Dev::EXTRAS{$::PDLMOD}{OBJECT} .= join '', map " $::PDLBASE-$_\$(OBJ_EXT)", qw(selectfunc);
}
if ($^O =~ /MSWin/) {
pp_addhdr('
#include <float.h>
');
}
pp_addhdr('
#include <math.h>
#define CONCAT_(A, B) A ## B
#define CONCAT(A, B) CONCAT_(A, B)
#define FORTRAN(name) CONCAT(name, F77_USCORE)
typedef PDL_Long logical;
typedef PDL_Long integer;
typedef PDL_Long ftnlen;
#ifndef abs
#define abs(x) ((x) >= 0 ? (x) : -(x))
#endif
#ifndef max
#define max(a,b) ((a) >= (b) ? (a) : (b))
#endif
');
pp_addpm({At=>'Top'},<<'EOD');
use PDL::Core;
use PDL::Slices;
use PDL::Ops qw//;
use PDL::Math qw/floor/;
use PDL::MatrixOps qw/identity/;
use PDL::LinearAlgebra;
use PDL::LinearAlgebra::Real qw //;
use PDL::LinearAlgebra::Complex qw //;
use strict;
=encoding utf8
=head1 NAME
PDL::LinearAlgebra::Trans - Linear Algebra based transcendental functions for PDL
=head1 SYNOPSIS
use PDL::LinearAlgebra::Trans;
$a = random (100,100);
$sqrt = msqrt($a);
=head1 DESCRIPTION
This module provides some transcendental functions for matrices.
Moreover it provides sec, asec, sech, asech, cot, acot, acoth, coth, csc,
acsc, csch, acsch. Beware, importing this module will overwrite the hidden
PDL routine sec. If you need to call it specify its origin module : PDL::Basic::sec(args)
=cut
EOD
pp_add_exported('', ' mexp mexpts mlog msqrt mpow
mcos msin mtan msec mcsc mcot
mcosh msinh mtanh msech mcsch mcoth
macos masin matan masec macsc macot
macosh masinh matanh masech macsch macoth
sec asec sech asech
cot acot acoth coth mfun
csc acsc csch acsch toreal pi');
pp_def('geexp',
Pars => '[io]A(n,n);int deg();scale();[io]trace();int [o]ns();int [o]info(); int [t]ipiv(n); [t]wsp(wspn=CALC(3*$SIZE(n)*$SIZE(n)));',
GenericTypes => [D],
Code => '
/* ----------------------------------------------------------------------|
int dgpadm_(integer *ideg, integer *m, double *t,
double *h__, integer *ldh, double *wsp, integer *lwsp,
integer *ipiv, integer *iexph, integer *ns, integer *iflag)
-----Purpose----------------------------------------------------------|
Computes exp(t*H), the matrix exponential of a general matrix in
full, using the irreducible rational Pade approximation to the
exponential function exp(x) = r(x) = (+/-)( I + 2*(q(x)/p(x)) ),
combined with scaling-and-squaring.
-----Arguments--------------------------------------------------------|
ideg|deg : (input) the degre of the diagonal Pade to be used.
a value of 6 is generally satisfactory.
m : (input) order of H.
H|A(ldh,m): (input) argument matrix.
t|scale : (input) time-scale (can be < 0).
wsp|exp(lwsp) : (workspace/output) lwsp .ge. 4*m*m+ideg+1.
ipiv(m) : (workspace)
>>>> iexph|index : (output) number such that wsp(iexph) points to exp(tH)
i.e., exp(tH) is located at wsp(iexph ... iexph+m*m-1)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
NOTE: if the routine was called with wsp(iptr),
then exp(tH) will start at wsp(iptr+iexph-1).
ns : (output) number of scaling-squaring used.
iflag|info : (output) exit flag.
0 - no problem
> 0 - Singularity in LU factorization when solving
Pade approximation
----------------------------------------------------------------------|
Roger B. Sidje (rbs@maths.uq.edu.au)
EXPOKIT: Software Package for Computing Matrix Exponentials.
ACM - Transactions On Mathematical Software, 24(1):130-156, 1998
Gregory Vanuxem (vanuxemg@yahoo.fr)
Minor modifications (Fortran to C with help of f2c, memory allocation,
PDL adaptation, null entry matrix, error return,
exp(tH) in entry matrix, trace normalization).
----------------------------------------------------------------------|
*/
extern int FORTRAN(dscal)(integer *, double *, double *,
integer *);
extern int FORTRAN(dgemm)(char *, char *, integer *, integer *, integer *,
double *, double *, integer *, double *,
integer *, double *, double *, integer *,
ftnlen, ftnlen);
extern int FORTRAN(dgesv)(integer *, integer *, double *,
integer *, integer *, double *, integer *, integer *);
extern int FORTRAN(daxpy)(integer *, double *, double *,
integer *, double *, integer *);
extern double FORTRAN(dlange)(char *norm, integer *m, integer *n, double *a, integer
*lda, double *work);
integer i__1, i__2, i__, j, k;
integer ip, mm, iq, ih2, ifree, iused, iodd, iget, iput;
double cp, cq, scale, scale2, hnorm, tracen;
double *coef, *wsp = $P(wsp);
integer c__1 = 1;
integer c__2 = 2;
double c_b7 = 0.;
double c_b11 = 1.;
double c_b19 = -1.;
double c_b23 = 2.;
double *h__ = $P(A);
$info() = 0;
/////////////
mm = $SIZE(n) * $SIZE(n);
coef = (double *) malloc(($deg() + 1 + mm) * sizeof(double));
//////////
// --- initialise pointers ...
ih2 = $deg() + 1;
ip = 0;
iq = mm;
ifree = iq + mm;
tracen = 0;
if ($trace())
{
loop(n0) %{ tracen += h__[n0 + n0 * $SIZE(n)]; %}
tracen /= $SIZE(n);
if (tracen > 0)
{
loop(n0) %{ h__[n0 + n0 * $SIZE(n)] -= tracen; %}
}
}
// --- scaling: seek ns such that ||t*H/2^ns|| < 1/2;
// and set scale = t/2^ns ...
// Compute Infinity norm
hnorm = FORTRAN(dlange)("I", &(integer){$SIZE(n)}, &(integer){$SIZE(n)}, $P(A),
&(integer){$SIZE(n)}, wsp);
hnorm = abs($scale() * hnorm);
if (hnorm == 0.)
{
free(coef);
coef = $P(A);
$ns() = 0;
loop(n0,n1) %{
coef[n1 + n0 * $SIZE(n)] = (n0 == n1) ? 1 : 0;
%}
}
else
{
i__2 = (integer ) (log(hnorm) / log(2.0) + 2);
$ns() = max(0, i__2);
scale = $scale() / pow(2.0, (double )$ns());
scale2 = scale * scale;
// --- compute Pade coefficients ...
i__ = $deg() + 1;
j = ($deg() << 1) + 1;
coef[0] = 1.;
i__1 = $deg();
for (k = 1; k <= i__1; ++k) {
coef[k] = coef[k - 1] * (double) (i__ - k) / (double) (k * (j - k));
}
// --- H2 = scale2*H*H ...
FORTRAN(dgemm)("n", "n", &(integer){$SIZE(n)}, &(integer){$SIZE(n)}, &(integer){$SIZE(n)}, &scale2, h__, &(integer){$SIZE(n)}, h__,
&(integer){$SIZE(n)}, &c_b7, &coef[ih2], &(integer){$SIZE(n)}, (ftnlen)1, (ftnlen)1);
// --- initialize p (numerator) and q (denominator) ...
cp = coef[$deg() - 1];
cq = coef[$deg()];
loop(n0) %{
loop(n1) %{
wsp[ip + n0 * $SIZE(n) + n1] = 0.;
wsp[iq + n0 * $SIZE(n) + n1] = 0.;
%}
wsp[ip + n0 * ($SIZE(n) + 1)] = cp;
wsp[iq + n0 * ($SIZE(n) + 1)] = cq;
%}
// --- Apply Horner rule ...
iodd = 1;
k = $deg() - 2;
do{
iused = iodd * iq + (1 - iodd) * ip;
FORTRAN(dgemm)("n", "n", &(integer){$SIZE(n)}, &(integer){$SIZE(n)}, &(integer){$SIZE(n)}, &c_b11, &wsp[iused], &(integer){$SIZE(n)}, &coef[ih2], &(integer){$SIZE(n)}, &c_b7, &
wsp[ifree], &(integer){$SIZE(n)}, (ftnlen)1, (ftnlen)1);
loop(n0) %{
wsp[ifree + n0 * ($SIZE(n) + 1)] += coef[k];
%}
ip = (1 - iodd) * ifree + iodd * ip;
iq = iodd * ifree + (1 - iodd) * iq;
ifree = iused;
iodd = 1 - iodd;
--k;
}
while (k >= 0);
// --- Obtain (+/-)(I + 2*(p\q)) ...
if (iodd == 1)
{
FORTRAN(dgemm)("n", "n", &(integer){$SIZE(n)}, &(integer){$SIZE(n)}, &(integer){$SIZE(n)}, &scale, &wsp[iq], &(integer){$SIZE(n)}, h__, &(integer){$SIZE(n)}, &
c_b7, &wsp[ifree], &(integer){$SIZE(n)}, (ftnlen)1, (ftnlen)1);
iq = ifree;
}
else
{
FORTRAN(dgemm)("n", "n", &(integer){$SIZE(n)}, &(integer){$SIZE(n)}, &(integer){$SIZE(n)}, &scale, &wsp[ip], &(integer){$SIZE(n)}, h__, &(integer){$SIZE(n)}, &
c_b7, &wsp[ifree], &(integer){$SIZE(n)}, (ftnlen)1, (ftnlen)1);
ip = ifree;
}
FORTRAN(daxpy)(&mm, &c_b19, &wsp[ip], &c__1, &wsp[iq], &c__1);
FORTRAN(dgesv)(&(integer){$SIZE(n)}, &(integer){$SIZE(n)}, &wsp[iq], &(integer){$SIZE(n)}, $P(ipiv), &wsp[ip], &(integer){$SIZE(n)}, $P(info));
if ($info() != 0)
{
free(coef);
}
else
{
FORTRAN(dscal)(&mm, &c_b23, &wsp[ip], &c__1);
loop(n0) %{
wsp[ip + n0 * ($SIZE(n) + 1)] += 1.;
%}
iput = ip;
if ($ns() == 0 && iodd == 1)
FORTRAN(dscal)(&mm, &c_b19, &wsp[ip], &c__1);
else{
// -- squaring : exp(t*H) = (exp(t*H))^(2^ns) ...
iodd = 1;
i__1 = $ns();
for (k = 0; k < i__1; k++) {
iget = iodd * ip + (1 - iodd) * iq;
iput = (1 - iodd) * ip + iodd * iq;
FORTRAN(dgemm)("n", "n", &(integer){$SIZE(n)}, &(integer){$SIZE(n)}, &(integer){$SIZE(n)}, &c_b11, &wsp[iget], &(integer){$SIZE(n)}, &wsp[iget], &(integer){$SIZE(n)}, &c_b7,
&wsp[iput], &(integer){$SIZE(n)}, (ftnlen)1, (ftnlen)1);
iodd = 1 - iodd;
}
}
free(coef);
coef = $P(A);
i__2 = iput + mm;
i__ = 0;
$trace() = tracen;
if (tracen > 0)
{
scale = exp(tracen);
for (i__1 = iput; i__1 < i__2 ; i__1++)
coef[i__++] = scale * wsp[i__1];
}
else
{
for (i__1 = iput; i__1 < i__2 ; i__1++)
coef[i__++] = wsp[i__1];
}
}
}
',
Doc => <<EOT
=for ref
Computes exp(t*A), the matrix exponential of a general matrix,
using the irreducible rational Pade approximation to the
exponential function exp(x) = r(x) = (+/-)( I + 2*(q(x)/p(x)) ),
combined with scaling-and-squaring and optionally normalization of the trace.
The algorithm is described in Roger B. Sidje (rbs@maths.uq.edu.au)
"EXPOKIT: Software Package for Computing Matrix Exponentials".
ACM - Transactions On Mathematical Software, 24(1):130-156, 1998
A: On input argument matrix. On output exp(t*A).
Use Fortran storage type.
deg: the degre of the diagonal Pade to be used.
a value of 6 is generally satisfactory.
scale: time-scale (can be < 0).
trace: on input, boolean value indicating whether or not perform
a trace normalization. On output value used.
ns: on output number of scaling-squaring used.
info: exit flag.
0 - no problem
> 0 - Singularity in LU factorization when solving
Pade approximation
=for example
$a = random(5,5);
$trace = pdl(1);
$a->t->geexp(6,1,$trace, ($ns = null), ($info = null));
=cut
EOT
);
do './pp_defc.pl'; die if $@;
pp_defc('geexp',
Pars => '[io]A(2,n,n);int deg();scale();int trace();int [o]ns();int [o]info(); int [t] ipiv(n)',
Doc => '
=for ref
Complex version of geexp. The value used for trace normalization is not returned.
The algorithm is described in Roger B. Sidje (rbs@maths.uq.edu.au)
"EXPOKIT: Software Package for Computing Matrix Exponentials".
ACM - Transactions On Mathematical Software, 24(1):130-156, 1998
',
GenericTypes => [D],
Code => '
/* ----------------------------------------------------------------------|
int dgpadm_(integer *ideg, integer *m, double *t,
PDL_CDouble *h__, integer *ldh, PDL_CDouble *wsp, integer *lwsp,
integer *ipiv, integer *iexph, integer *ns, integer *iflag)
-----Purpose----------------------------------------------------------|
Computes exp(t*H), the matrix exponential of a general matrix in
full, using the irreducible rational Pade approximation to the
exponential function exp(x) = r(x) = (+/-)( I + 2*(q(x)/p(x)) ),
combined with scaling-and-squaring.
-----Arguments--------------------------------------------------------|
ideg|deg : (input) the degre of the diagonal Pade to be used.
a value of 6 is generally satisfactory.
m : (input) order of H.
H|A(ldh,m): (input) argument matrix.
t|scale : (input) time-scale (can be < 0).
wsp|exp(lwsp) : (workspace/output) lwsp .ge. 4*m*m+ideg+1.
ipiv(m) : (workspace)
>>>> iexph : (output) number such that wsp(iexph) points to exp(tH)
i.e., exp(tH) is located at wsp(iexph ... iexph+m*m-1)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
NOTE: if the routine was called with wsp(iptr),
then exp(tH) will start at wsp(iptr+iexph-1).
ns : (output) number of scaling-squaring used.
iflag|info : (output) exit flag.
0 - no problem
> 0 - Singularity in LU factorization when solving
Pade approximation
----------------------------------------------------------------------|
Roger B. Sidje (rbs@maths.uq.edu.au)
EXPOKIT: Software Package for Computing Matrix Exponentials.
ACM - Transactions On Mathematical Software, 24(1):130-156, 1998
Gregory Vanuxem (vanuxemg@yahoo.fr)
Minor modifications (Fortran to C with help of f2c, memory allocation,
PDL adaptation, null entry matrix, error return,
exp(tH) in entry matrix), trace normalization.
----------------------------------------------------------------------|
*/
static PDL_CDouble c_b1 = 0;
static PDL_CDouble c_b2 = 1;
static integer c__2 = 2;
static integer c__1 = 1;
static double c_b19 = 2.;
static double c_b21 = -1.;
integer i__1, i__2, i__3, i__, j, k;
integer ip, mm, iq, ih2, iodd, iget, iput, ifree, iused;
double hnorm,d__1, d__2;
PDL_CDouble scale, tracen, cp, cq, z__1;
PDL_CDouble scale2;
PDL_CDouble *wsp, *h__;
extern int FORTRAN(zgemm)(char *, char *, integer *, integer *,
integer *, PDL_CDouble *, PDL_CDouble *, integer *,
PDL_CDouble *, integer *, PDL_CDouble *, PDL_CDouble *,
integer *, ftnlen, ftnlen);
extern int FORTRAN(zgesv)(integer *, integer *, PDL_CDouble *,
integer *, integer *, PDL_CDouble *, integer *, integer *);
extern int FORTRAN(zaxpy)(integer *, PDL_CDouble *,
PDL_CDouble *, integer *, PDL_CDouble *, integer *), FORTRAN(zdscal)(
integer *, double *, PDL_CDouble *, integer *);
extern double FORTRAN(zlange)(char *norm, integer *m, integer *n, PDL_CDouble *a,
integer *lda, PDL_CDouble *work);
h__ = (PDL_CDouble *) $P(A);
wsp = malloc (( $deg() + 1 + 4 * $SIZE(n) * $SIZE(n)) * sizeof (PDL_CDouble));
mm = $SIZE(n) * $SIZE(n);
$info() = 0;
ih2 = $deg() + 1;
ip = ih2 + mm;
iq = ip + mm;
ifree = iq + mm;
tracen = 0;
if ($trace())
{
loop(n0) %{
tracen += h__[n0 + n0 * $SIZE(n)];
%}
tracen /= $SIZE(n);
if (creal(tracen) < 0){
tracen = cimag(tracen);
}
loop(n0) %{
h__[n0 + n0 * $SIZE(n)] -= tracen;
%}
}
// --- scaling: seek ns such that ||t*H/2^ns|| < 1/2;
// and set scale = t/2^ns ...
// Compute Infinity norm
hnorm = FORTRAN(zlange)("I", &(integer){$SIZE(n)}, &(integer){$SIZE(n)}, h__,
&(integer){$SIZE(n)}, wsp);
hnorm = abs($scale() * hnorm);
if (hnorm == 0.)
{
$ns() = 0;
loop(n0,n1) %{
h__[n1 + n0 * $SIZE(n) ] = (n0 == n1) ? 1 : 0;
%}
}
else
{
i__2 = (integer) ((log(hnorm) / log(2.)) + 2);
$ns() = max(0,i__2);
scale = $scale() / pow(2, (double )$ns());
scale2 = scale * scale;
// --- compute Pade coefficients ...
i__ = $deg() + 1;
j = ($deg() << 1) + 1;
wsp[0] = 1;
for (k = 1; k <= $deg(); ++k) {
d__1 = (double) (i__ - k);
d__2 = (double) (k * (j - k));
wsp[k] = (d__1 * wsp[k - 1]) / d__2;
}
// --- H2 = scale2*H*H ...
FORTRAN(zgemm)("n", "n", &(integer){$SIZE(n)}, &(integer){$SIZE(n)}, &(integer){$SIZE(n)}, &scale2, &h__[0], &(integer){$SIZE(n)}, &h__[0],
&(integer){$SIZE(n)}, &c_b1, &wsp[ih2], &(integer){$SIZE(n)}, (ftnlen)1, (ftnlen)1);
// --- initialise p (numerator) and q (denominator) ...
cp = wsp[$deg() - 1];
cq = wsp[$deg()];
loop(n0) %{
loop(n1) %{
i__3 = ip + n0 * $SIZE(n) + n1;
wsp[i__3] = 0.;
i__3 = iq + n0 * $SIZE(n) + n1;
wsp[i__3] = 0.;
%}
i__2 = ip + n0 * ($SIZE(n) + 1);
wsp[i__2] = cp;
i__2 = iq + n0 * ($SIZE(n) + 1);
wsp[i__2] = cq;
%}
// --- Apply Horner rule ...
iodd = 1;
k = $deg() - 2;
do
{
iused = iodd * iq + (1 - iodd) * ip;
FORTRAN(zgemm)("n", "n", &(integer){$SIZE(n)}, &(integer){$SIZE(n)}, &(integer){$SIZE(n)}, &c_b2, &wsp[iused], &(integer){$SIZE(n)}, &wsp[ih2], &(integer){$SIZE(n)}, &c_b1, &
wsp[ifree], &(integer){$SIZE(n)}, (ftnlen)1, (ftnlen)1);
loop(n0) %{
i__1 = ifree + n0 * ($SIZE(n) + 1);
i__2 = ifree + n0 * ($SIZE(n) + 1);
i__3 = k;
wsp[i__1] = wsp[i__2] + wsp[i__3];
%}
ip = (1 - iodd) * ifree + iodd * ip;
iq = iodd * ifree + (1 - iodd) * iq;
ifree = iused;
iodd = 1 - iodd;
--k;
}
while(k >= 0);
// --- Obtain (+/-)(I + 2*(p\q)) ...
if (iodd != 0)
{
FORTRAN(zgemm)("n", "n", &(integer){$SIZE(n)}, &(integer){$SIZE(n)}, &(integer){$SIZE(n)}, &scale, &wsp[iq], &(integer){$SIZE(n)}, &h__[0], &(integer){$SIZE(n)}, &
c_b1, &wsp[ifree], &(integer){$SIZE(n)}, (ftnlen)1, (ftnlen)1);
iq = ifree;
}else
{
FORTRAN(zgemm)("n", "n", &(integer){$SIZE(n)}, &(integer){$SIZE(n)}, &(integer){$SIZE(n)}, &scale, &wsp[ip], &(integer){$SIZE(n)}, &h__[0], &(integer){$SIZE(n)}, &
c_b1, &wsp[ifree], &(integer){$SIZE(n)}, (ftnlen)1, (ftnlen)1);
ip = ifree;
}
z__1 = -1.;
FORTRAN(zaxpy)(&mm, &z__1, &wsp[ip], &c__1, &wsp[iq], &c__1);
FORTRAN(zgesv)(&(integer){$SIZE(n)}, &(integer){$SIZE(n)}, &wsp[iq], &(integer){$SIZE(n)}, $P(ipiv), &wsp[ip], &(integer){$SIZE(n)}, $P(info));
if ($info() == 0)
{
FORTRAN(zdscal)(&mm, &c_b19, &wsp[ip], &c__1);
loop(n0) %{
i__2 = ip + n0 * ($SIZE(n) + 1);
i__3 = ip + n0 * ($SIZE(n) + 1);
wsp[i__2] = wsp[i__3] + 1.;
%}
iput = ip;
if ($ns() == 0 && iodd != 0)
FORTRAN(zdscal)(&mm, &c_b21, &wsp[ip], &c__1);
else
{
// -- squaring : exp(t*H) = (exp(t*H))^(2^ns) ...
iodd = 1;
for (k = 0; k < $ns(); k++)
{
iget = iodd * ip + (1 - iodd) * iq;
iput = (1 - iodd) * ip + iodd * iq;
FORTRAN(zgemm)("n", "n", &(integer){$SIZE(n)}, &(integer){$SIZE(n)}, &(integer){$SIZE(n)}, &c_b2, &wsp[iget], &(integer){$SIZE(n)}, &wsp[iget], &(integer){$SIZE(n)}, &c_b1,
&wsp[iput], &(integer){$SIZE(n)}, (ftnlen)1, (ftnlen)1);
iodd = 1 - iodd;
}
}
}
i__2 = iput + mm;
i__ = 0;
if ($trace())
{
tracen = cexp(tracen);
for (i__1 = iput; i__1 < i__2 ; i__1++)
{
h__[i__++] = tracen * wsp[i__1];
}
}
else
{
for (i__1 = iput; i__1 < i__2 ; i__1++)
{
h__[i__++] = wsp[i__1];
}
}
}
free(wsp);
',
);
pp_defc('trsqrt',
Pars => '[io]A(2,n,n);int uplo();[o] B(2,n,n);int [o]info()',
Doc => '
=for ref
Root square of complex triangular matrix. Uses a recurrence of Björck and Hammarling.
(See Nicholas J. Higham. A new sqrtm for MATLAB. Numerical Analysis
Report No. 336, Manchester Centre for Computational Mathematics,
If uplo is true, A is lower triangular.
',
GenericTypes => [D],
Code => '
PDL_CDouble *cb, *ca;
PDL_CDouble s, snum, sdenum;
double *b;
double tt, dn;
integer i, j, k, l, ind, ind1, ind2;
ca = (PDL_CDouble *) $P(A);
b = (void *)$P(B);
cb = (PDL_CDouble *) $P(B);
#define subscr(a_1,a_2)\
($uplo()) ? (a_2) * $SIZE(n) + (a_1) : (a_1) * $SIZE(n) + (a_2)
$info() = 0;
ind = $SIZE(n) * $SIZE(n) * 2;
for (i = 0; i < ind;i++)
b[i] = 0;
loop(n0) %{
ind = subscr(n0,n0);
cb[ind] = csqrt(ca[ind]);
%}
for (k = 0; k < $SIZE(n)-1; k++)
{
for (i = 0; i < $SIZE(n)-(k+1); i++)
{
j = i + k + 1;
ind = subscr(i,j);
s = 0;
for (l = i+1; l < j; l++)
{
ind1 = subscr(i,l);
ind2 = subscr(l,j);
s += cb[ind1] * cb[ind2];
}
ind1 = subscr(i,i);
ind2 = subscr(j,j);
sdenum = cb[ind1] + cb[ind2];
snum = ca[ind] - s;
if (fabs (creal(sdenum)) > fabs (cimag(sdenum)))
{
if (creal(sdenum) == 0){
$info() = -1;
goto ENDCTRSQRT;
}
tt = cimag(sdenum) / creal(sdenum);
dn = creal(sdenum) + tt * cimag(sdenum);
cb[ind] = (snum + tt * (-I * snum)) / dn;
}
else
{
if (cimag(sdenum) == 0){
$info() = -1;
goto ENDCTRSQRT;
}
tt = creal(sdenum) / cimag(sdenum);
dn = creal(sdenum) * tt + cimag(sdenum);
cb[ind] = (snum * tt + (-I * snum)) / dn;
}
}
}
ENDCTRSQRT:
;
#undef subscr
',
);
pp_addhdr('
void dfunc_wrapper(PDL_CDouble *p, integer n, SV* dfunc);
');
pp_defc('trfun',
Pars => '[io]A(2,n,n);int uplo();[o] B(2,n,n);int [o]info(); complex [t]diag(n);',
OtherPars => "SV* func" ,
Doc => '
=for ref
Apply an arbitrary function to a complex triangular matrix. Uses a recurrence of Parlett.
If uplo is true, A is lower triangular.
',
GenericTypes => [D],
Code => '
PDL_CDouble *cb, *ca;
PDL_CDouble s, snum, sdenum;
double *b;
double tt, dn;
integer i, j, k, l, p, ind, ind1, ind2, tmp, c__1;
extern double FORTRAN(zdotu)(PDL_CDouble *ret, integer *n, PDL_CDouble *dx, integer *incx, PDL_CDouble *dy,
integer *incy);
ca = (PDL_CDouble *) $P(A);
b = (double *)$P(B);
cb = (PDL_CDouble *) $P(B);
c__1 = 1;
#define subscr(a_1,a_2) \
( $uplo() ) ? (a_2)+ $SIZE(n) * (a_1) : (a_1)+ $SIZE(n) * (a_2)
$info() = 0;
ind = $SIZE(n) * $SIZE(n) * 2;
for (i = 0; i < ind;i++)
b[i] = 0;
loop (n) %{
$diag() = ca[subscr(n,n)];
%}
dfunc_wrapper($P(diag), $SIZE(n), $COMP(func));
loop (n) %{
cb[subscr(n,n)] = $diag();
%}
for (p = 1; p < $SIZE(n); p++)
{
tmp = $SIZE(n) - p;
for (i = 0; i < tmp; i++)
{
j = i + p;
//$s = $T(,($j),($i))->Cmul($F(,($j),($j))->Csub($F(,($i),($i))));
ind1 = subscr(i,i);
ind2 = subscr(j,j);
ind = subscr(j,i);
s = snum = ca[ind] * (cb[ind2] - cb[ind1]);
if (i < (j-1))
{
//$s = $s + $T(,$i+1:$j-1,($i))->cdot(1, $F(,($j), $i+1:$j-1),1)->Csub($F(,$i+1:$j-1,($i))->cdot(1,$T(,($j), $i+1:$j-1),1));
// $T(,$i+1:$j-1,($i))->cdot(1, $F(,($j), $i+1:$j-1),1)
l = 0;
for(k = i+1; k < j; k++)
$diag(n=>l++) = cb[subscr(j,k)];
ind = subscr(i+1,i);
// TODO : Humm
FORTRAN(zdotu)(&s, &l, &ca[ind], &c__1, $P(diag), &c__1);
snum += s;
// $F(,$i+1:$j-1,($i))->cdot(1,$T(,($j), $i+1:$j-1),1)
l = 0;
for(k = i+1; k < j; k++)
$diag(n=>l++) = ca[subscr(j,k)];
ind = subscr(i+1,i);
// TODO : Humm
FORTRAN(zdotu)(&s, &l, &cb[ind], &c__1, $P(diag), &c__1);
snum -= s;
ind = subscr(j,i);
}
//$sdenum = $T(,($j),($j))->Csub($T(,($i),($i)));
sdenum = ca[ind2] - ca[ind1];
//$s = $s / $sdenum;
if (fabs (creal(sdenum)) > fabs (cimag(sdenum)))
{
if (creal(sdenum) == 0)
{
$info() = -1;
goto ENDCTRFUN;
}
tt = cimag(sdenum) / creal(sdenum);
dn = creal(sdenum) + tt * cimag(sdenum);
cb[ind] = (snum + tt * (-I * snum)) / dn;
}
else
{
if (cimag(sdenum) == 0)
{
$info() = -1;
goto ENDCTRFUN;
}
tt = creal(sdenum) / cimag(sdenum);
dn = creal(sdenum) * tt + cimag(sdenum);
cb[ind] = (snum * tt + (-I * snum)) / dn;
}
}
}
ENDCTRFUN:
;
#undef subscr
',
);
pp_addpm(<<'EOD');
my $pi;
BEGIN { $pi = pdl(3.1415926535897932384626433832795029) }
sub pi () { $pi->copy };
*sec = \&PDL::sec;
sub PDL::sec{1/cos($_[0])}
*csc = \&PDL::csc;
sub PDL::csc($) {1/sin($_[0])}
*cot = \&PDL::cot;
sub PDL::cot($) {1/(sin($_[0])/cos($_[0]))}
*sech = \&PDL::sech;
sub PDL::sech($){1/pdl($_[0])->cosh}
*csch = \&PDL::csch;
sub PDL::csch($) {1/pdl($_[0])->sinh}
*coth = \&PDL::coth;
sub PDL::coth($) {1/pdl($_[0])->tanh}
*asec = \&PDL::asec;
sub PDL::asec($) {my $tmp = 1/pdl($_[0]) ; $tmp->acos}
*acsc = \&PDL::acsc;
sub PDL::acsc($) {my $tmp = 1/pdl($_[0]) ; $tmp->asin}
*acot = \&PDL::acot;
sub PDL::acot($) {my $tmp = 1/pdl($_[0]) ; $tmp->atan}
*asech = \&PDL::asech;
sub PDL::asech($) {my $tmp = 1/pdl($_[0]) ; $tmp->acosh}
*acsch = \&PDL::acsch;
sub PDL::acsch($) {my $tmp = 1/pdl($_[0]) ; $tmp->asinh}
*acoth = \&PDL::acoth;
sub PDL::acoth($) {my $tmp = 1/pdl($_[0]) ; $tmp->atanh}
my $_tol = 9.99999999999999e-15;
sub toreal {
$_tol = $_[1] if defined $_[1];
return $_[0] if $_[0]->isempty or $_[0]->type->real;
my ($min, $max, $tmp);
($min, $max) = $_[0]->im->minmax;
return $_[0]->re->sever unless (abs($min) > $_tol || abs($max) > $_tol);
$_[0];
}
=head2 mlog
=for ref
Return matrix logarithm of a square matrix.
=for usage
PDL = mlog(PDL(A))
=for example
my $a = random(10,10);
my $log = mlog($a);
=cut
*mlog = \&PDL::mlog;
sub PDL::mlog {
&PDL::LinearAlgebra::_square;
my ($m, $tol) = @_;
mfun($m, sub{$_[0].=log $_[0]} , 0, $tol);
}
=head2 msqrt
=for ref
Return matrix square root (principal) of a square matrix.
=for usage
PDL = msqrt(PDL(A))
=for example
my $a = random(10,10);
my $sqrt = msqrt($a);
=cut
*msqrt = \&PDL::msqrt;
sub PDL::msqrt {
&PDL::LinearAlgebra::_square;
my ($m, $tol) = @_;
$m = $m->r2C unless $m->_is_complex;
my ($t, undef, $z, undef, $info) = $m->mschur(1);
if ($info){
warn "msqrt: Can't compute Schur form\n";
return;
}
($t, $info) = $t->ctrsqrt(0);
if($info){
warn "msqrt: can't compute square root\n";
return;
}
$m = $z x $t x $z->t(1);
return $m->_is_complex ? $m : toreal($m, $tol);
}
=head2 mexp
=for ref
Return matrix exponential of a square matrix.
=for usage
PDL = mexp(PDL(A))
=for example
my $a = random(10,10);
my $exp = mexp($a);
=cut
*mexp = \&PDL::mexp;
sub PDL::mexp {
&PDL::LinearAlgebra::_square;
my ($m, $order, $trace) = @_;
$trace = 1 unless defined $trace;
$order = 6 unless defined $order;
$m = $m->copy;
$m->t->_call_method('geexp', $order, 1, $trace, my $ns = PDL->null, my $info = PDL->null);
if ($info){
warn "mexp: Error $info";
}
else{
return $m;
}
}
*mexpts = \&PDL::mexpts;
sub PDL::mexpts {
&PDL::LinearAlgebra::_square;
my ($m, $order, $tol) = @_;
my @dims = $m->dims;
my ($em, $trm);
$order = 20 unless defined $order;
$em = $m->_is_complex ? diag(r2C(ones($dims[1]))) : diag(ones($dims[1]));
$trm = $em->copy;
for (1..($order - 1)){
$trm = $trm x ($m / $_);
$em += $trm;
}
return $m->_is_complex ? $em : toreal($em, $tol);
}
=head2 mpow
=for ref
Return matrix power of a square matrix.
=for usage
PDL = mpow(PDL(A), SCALAR(exponent))
=for example
my $a = random(10,10);
my $powered = mpow($a,2.5);
=cut
*mpow = \&PDL::mpow;
sub PDL::mpow {
&PDL::LinearAlgebra::_square;
my $di = $_[0]->dims_internal;
my ($m, $power, $tol, $eigen) = @_;
my @dims = $m->dims;
my $ret;
if (UNIVERSAL::isa($power,'PDL') and $power->dims > 1){
my ($e, $v) = $m->meigen(0,1);
$ret = $v * ($e**$power) x $v->minv;
}
elsif( 1/$dims[$di] * 1000 > abs($power) and !$eigen){
$ret = identity($dims[$di]);
$ret = $ret->r2C if $m->_is_complex;
my $pow = floor($power);
$pow++ if ($power < 0 and $power != $pow);
# TODO: what a beautiful thing (is it a game ?)
for(my $i = 0; $i < abs($pow); $i++){$ret x= $m;}
$ret = $ret->minv if $power < 0;
if ($power = $power - $pow){
if($power == 0.5){
my $v = $m->msqrt;
$ret = ($pow == 0) ? $v : $ret x $v;
}
else{
my ($e, $v) = $m->meigen(0,1);
$ret = ($pow == 0) ? ($v * $e**$power x $v->minv) : $ret->r2C x ($v * $e**$power x $v->minv);
}
}
}
else{
my ($e, $v) = $m->meigen(0,1);
$ret = $v * $e**$power x $v->minv;
}
return $m->_is_complex ? $ret : toreal($ret, $tol);
}
=head2 mcos
=for ref
Return matrix cosine of a square matrix.
=for usage
PDL = mcos(PDL(A))
=for example
my $a = random(10,10);
my $cos = mcos($a);
=cut
sub _i {
i();
}
*mcos = \&PDL::mcos;
sub PDL::mcos {
&PDL::LinearAlgebra::_square;
my $m = shift;
my $i = _i();
return $m->_is_complex ? (mexp($i*$m) + mexp(- $i*$m)) / 2 : mexp($i*$m)->re->sever;
}
=head2 macos
=for ref
Return matrix inverse cosine of a square matrix.
=for usage
PDL = macos(PDL(A))
=for example
my $a = random(10,10);
my $acos = macos($a);
=cut
*macos = \&PDL::macos;
sub PDL::macos {
&PDL::LinearAlgebra::_square;
my $di = $_[0]->dims_internal;
my ($m, $tol) = @_;
my @dims = $m->dims;
my $id = identity($dims[$di]); $id = $id->r2C if $m->_is_complex;
my $i = _i();
my $ret = $i * mlog( ($m->r2C - $i * msqrt( ($id - $m x $m), $tol)));
return $m->_is_complex ? $ret : toreal($ret, $tol);
}
=head2 msin
=for ref
Return matrix sine of a square matrix.
=for usage
PDL = msin(PDL(A))
=for example
my $a = random(10,10);
my $sin = msin($a);
=cut
*msin = \&PDL::msin;
sub PDL::msin {
&PDL::LinearAlgebra::_square;
my $m = shift;
my $i = _i();
$m->_is_complex ? (mexp($i*$m) - mexp(- $i*$m))/(2*$i) : mexp($i*$m)->im->sever;
}
=head2 masin
=for ref
Return matrix inverse sine of a square matrix.
=for usage
PDL = masin(PDL(A))
=for example
my $a = random(10,10);
my $asin = masin($a);
=cut
*masin = \&PDL::masin;
sub PDL::masin {
&PDL::LinearAlgebra::_square;
my $di = $_[0]->dims_internal;
my ($m, $tol) = @_;
my @dims = $m->dims;
my $i = _i();
my $id = identity($dims[$di]); $id = $id->r2C if $m->_is_complex;
my $ret = (-$i) * mlog((($i*$m) + msqrt($id - $m x $m, $tol)));
return $m->_is_complex ? $ret : toreal($ret, $tol);
}
=head2 mtan
=for ref
Return matrix tangent of a square matrix.
=for usage
PDL = mtan(PDL(A))
=for example
my $a = random(10,10);
my $tan = mtan($a);
=cut
*mtan = \&PDL::mtan;
sub PDL::mtan {
&PDL::LinearAlgebra::_square;
my ($m, $id) = @_;
my @dims = $m->dims;
return scalar msolvex(mcos($m), msin($m),equilibrate=>1) unless $id;
my $i = _i();
if ($m->_is_complex){
my $di = $_[0]->dims_internal;
$id = identity($dims[$di])->r2C;
$m = mexp(-2*$i*$m);
return scalar msolvex( ($id + $m ),( (- $i) * ($id - $m)),equilibrate=>1);
}
else{
$m = mexp($i * $m);
return scalar $m->re->msolvex($m->im,equilibrate=>1);
}
}
=head2 matan
=for ref
Return matrix inverse tangent of a square matrix.
=for usage
PDL = matan(PDL(A))
=for example
my $a = random(10,10);
my $atan = matan($a);
=cut
*matan = \&PDL::matan;
sub PDL::matan {
&PDL::LinearAlgebra::_square;
my $di = $_[0]->dims_internal;
my ($m, $tol) = @_;
my @dims = $m->dims;
my $i = _i();
my $id = identity($dims[$di]); $id = $id->r2C if $m->_is_complex;
my $ret = -$i/2 * mlog( scalar PDL::msolvex( ($id - $i*$m) ,($id + $i*$m),equilibrate=>1 ));
return $m->_is_complex ? $ret : toreal($ret, $tol);
}
=head2 mcot
=for ref
Return matrix cotangent of a square matrix.
=for usage
PDL = mcot(PDL(A))
=for example
my $a = random(10,10);
my $cot = mcot($a);
=cut
*mcot = \&PDL::mcot;
sub PDL::mcot {
&PDL::LinearAlgebra::_square;
my ($m, $id) = @_;
my @dims = $m->dims;
return scalar msolvex(msin($m),mcos($m),equilibrate=>1) unless $id;
my $i = _i();
if ($m->_is_complex){
my $di = $_[0]->dims_internal;
$id = identity($dims[$di])->r2C;
$m = mexp(-2*$i*$m);
return scalar msolvex( ($id - $m), ($i * ($id + $m)), equilibrate=>1);
}
else{
$m = mexp($i * $m);
return scalar $m->im->msolvex($m->re,equilibrate=>1);
}
}
=head2 macot
=for ref
Return matrix inverse cotangent of a square matrix.
=for usage
PDL = macot(PDL(A))
=for example
my $a = random(10,10);
my $acot = macot($a);
=cut
*macot = \&PDL::macot;
sub PDL::macot {
&PDL::LinearAlgebra::_square;
my ($m, $tol, $id) = @_;
my ($inv, $info) = $m->minv;
if ($info){
warn "macot: singular matrix";
return;
}
return matan($inv,$tol);
}
=head2 msec
=for ref
Return matrix secant of a square matrix.
=for usage
PDL = msec(PDL(A))
=for example
my $a = random(10,10);
my $sec = msec($a);
=cut
*msec = \&PDL::msec;
sub PDL::msec {
&PDL::LinearAlgebra::_square;
my $m = shift;
my $i = _i();
return $m->_is_complex ? PDL::minv(mexp($i*$m) + mexp(- $i*$m)) * 2 : scalar PDL::minv(re(mexp($i*$m)));
}
=head2 masec
=for ref
Return matrix inverse secant of a square matrix.
=for usage
PDL = masec(PDL(A))
=for example
my $a = random(10,10);
my $asec = masec($a);
=cut
*masec = \&PDL::masec;
sub PDL::masec {
&PDL::LinearAlgebra::_square;
my ($m, $tol) = @_;
my ($inv, $info) = $m->minv;
if ($info){
warn "masec: singular matrix";
return;
}
return macos($inv,$tol);
}
=head2 mcsc
=for ref
Return matrix cosecant of a square matrix.
=for usage
PDL = mcsc(PDL(A))
=for example
my $a = random(10,10);
my $csc = mcsc($a);
=cut
*mcsc = \&PDL::mcsc;
sub PDL::mcsc {
&PDL::LinearAlgebra::_square;
my $m = shift;
my $i = _i();
$m->_is_complex ? PDL::minv(mexp($i*$m) - mexp(- $i*$m)) * 2*$i : scalar PDL::minv(im(mexp($i*$m)));
}
=head2 macsc
=for ref
Return matrix inverse cosecant of a square matrix.
=for usage
PDL = macsc(PDL(A))
=for example
my $a = random(10,10);
my $acsc = macsc($a);
=cut
*macsc = \&PDL::macsc;
sub PDL::macsc {
&PDL::LinearAlgebra::_square;
my ($m, $tol) = @_;
my ($inv, $info) = $m->minv;
if ($info){
warn "macsc: singular matrix";
return;
}
return masin($inv,$tol);
}
=head2 mcosh
=for ref
Return matrix hyperbolic cosine of a square matrix.
=for usage
PDL = mcosh(PDL(A))
=for example
my $a = random(10,10);
my $cos = mcosh($a);
=cut
*mcosh = \&PDL::mcosh;
sub PDL::mcosh {
&PDL::LinearAlgebra::_square;
my $m = shift;
( $m->mexp + mexp(-$m) )/2;
}
=head2 macosh
=for ref
Return matrix hyperbolic inverse cosine of a square matrix.
=for usage
PDL = macosh(PDL(A))
=for example
my $a = random(10,10);
my $acos = macosh($a);
=cut
*macosh = \&PDL::macosh;
sub PDL::macosh {
&PDL::LinearAlgebra::_square;
my ($m, $tol) = @_;
my @dims = $m->dims;
my $di = $_[0]->dims_internal;
my $id = identity($dims[$di]); $id = $id->r2C if $m->_is_complex;
my $ret = msqrt($m x $m - $id);
$m = $m->r2C if $ret->getndims > @dims;
mlog($m + $ret, $tol);
}
=head2 msinh
=for ref
Return matrix hyperbolic sine of a square matrix.
=for usage
PDL = msinh(PDL(A))
=for example
my $a = random(10,10);
my $sinh = msinh($a);
=cut
*msinh = \&PDL::msinh;
sub PDL::msinh {
&PDL::LinearAlgebra::_square;
my $m = shift;
( $m->mexp - mexp(-$m) )/2;
}
=head2 masinh
=for ref
Return matrix hyperbolic inverse sine of a square matrix.
=for usage
PDL = masinh(PDL(A))
=for example
my $a = random(10,10);
my $asinh = masinh($a);
=cut
*masinh = \&PDL::masinh;
sub PDL::masinh {
&PDL::LinearAlgebra::_square;
my ($m, $tol) = @_;
my @dims = $m->dims;
my $di = $_[0]->dims_internal;
my $id = identity($dims[$di]); $id = $id->r2C if $m->_is_complex;
my $ret = msqrt($m x $m + $id);
$m = $m->r2C if $ret->getndims > @dims;
mlog(($m + $ret), $tol);
}
=head2 mtanh
=for ref
Return matrix hyperbolic tangent of a square matrix.
=for usage
PDL = mtanh(PDL(A))
=for example
my $a = random(10,10);
my $tanh = mtanh($a);
=cut
*mtanh = \&PDL::mtanh;
sub PDL::mtanh {
&PDL::LinearAlgebra::_square;
my ($m, $id) = @_;
my @dims = $m->dims;
return scalar msolvex(mcosh($m), msinh($m),equilibrate=>1) unless $id;
my $di = $_[0]->dims_internal;
$id = identity($dims[$di]); $id = $id->r2C if $m->_is_complex;
$m = mexp(-2*$m);
return scalar msolvex( ($id + $m ),($id - $m), equilibrate=>1);
}
=head2 matanh
=for ref
Return matrix hyperbolic inverse tangent of a square matrix.
=for usage
PDL = matanh(PDL(A))
=for example
my $a = random(10,10);
my $atanh = matanh($a);
=cut
*matanh = \&PDL::matanh;
sub PDL::matanh {
&PDL::LinearAlgebra::_square;
my ($m, $tol) = @_;
my @dims = $m->dims;
my $di = $_[0]->dims_internal;
my $id = identity($dims[$di]); $id = $id->r2C if $m->_is_complex;
mlog( scalar msolvex( ($id - $m ),($id + $m),equilibrate=>1), $tol ) / 2;
}
=head2 mcoth
=for ref
Return matrix hyperbolic cotangent of a square matrix.
=for usage
PDL = mcoth(PDL(A))
=for example
my $a = random(10,10);
my $coth = mcoth($a);
=cut
*mcoth = \&PDL::mcoth;
sub PDL::mcoth {
&PDL::LinearAlgebra::_square;
my ($m, $id) = @_;
my @dims = $m->dims;
scalar msolvex(msinh($m), mcosh($m),equilibrate=>1) unless $id;
my $di = $_[0]->dims_internal;
$id = identity($dims[$di]); $id = $id->r2C if $m->_is_complex;
$m = mexp(-2*$m);
return scalar msolvex( ($id - $m ),($id + $m),equilibrate=>1);
}
=head2 macoth
=for ref
Return matrix hyperbolic inverse cotangent of a square matrix.
=for usage
PDL = macoth(PDL(A))
=for example
my $a = random(10,10);
my $acoth = macoth($a);
=cut
*macoth = \&PDL::macoth;
sub PDL::macoth {
&PDL::LinearAlgebra::_square;
my ($m, $tol) = @_;
my ($inv, $info) = $m->minv;
if ($info){
warn "macoth: singular matrix";
return;
}
return matanh($inv,$tol);
}
=head2 msech
=for ref
Return matrix hyperbolic secant of a square matrix.
=for usage
PDL = msech(PDL(A))
=for example
my $a = random(10,10);
my $sech = msech($a);
=cut
*msech = \&PDL::msech;
sub PDL::msech {
&PDL::LinearAlgebra::_square;
my $m = shift;
PDL::minv( $m->mexp + mexp(-$m) ) * 2;
}
=head2 masech
=for ref
Return matrix hyperbolic inverse secant of a square matrix.
=for usage
PDL = masech(PDL(A))
=for example
my $a = random(10,10);
my $asech = masech($a);
=cut
*masech = \&PDL::masech;
sub PDL::masech {
&PDL::LinearAlgebra::_square;
my ($m, $tol) = @_;
my ($inv, $info) = $m->minv;
if ($info){
warn "masech: singular matrix";
return;
}
return macosh($inv,$tol);
}
=head2 mcsch
=for ref
Return matrix hyperbolic cosecant of a square matrix.
=for usage
PDL = mcsch(PDL(A))
=for example
my $a = random(10,10);
my $csch = mcsch($a);
=cut
*mcsch = \&PDL::mcsch;
sub PDL::mcsch {
&PDL::LinearAlgebra::_square;
my $m = shift;
PDL::minv( $m->mexp - mexp(-$m) ) * 2;
}
=head2 macsch
=for ref
Return matrix hyperbolic inverse cosecant of a square matrix.
=for usage
PDL = macsch(PDL(A))
=for example
my $a = random(10,10);
my $acsch = macsch($a);
=cut
*macsch = \&PDL::macsch;
sub PDL::macsch {
&PDL::LinearAlgebra::_square;
my ($m, $tol) = @_;
my ($inv, $info) = $m->minv;
if ($info){
warn "macsch: singular matrix";
return;
}
return masinh($inv,$tol);
}
=head2 mfun
=for ref
Return matrix function of second argument of a square matrix.
Function will be applied on a complex ndarray.
=for usage
PDL = mfun(PDL(A),'cos')
=for example
my $a = random(10,10);
my $fun = mfun($a,'cos');
sub sinbycos2{
$_[0]->set_inplace(0);
$_[0] .= $_[0]->Csin/$_[0]->Ccos**2;
}
# Try diagonalization
$fun = mfun($a, \&sinbycos2,1);
# Now try Schur/Parlett
$fun = mfun($a, \&sinbycos2);
# Now with function.
scalar msolve($a->mcos->mpow(2), $a->msin);
=cut
*mfun = \&PDL::mfun;
sub PDL::mfun {
&PDL::LinearAlgebra::_square;
my ($m, $method, $diag, $tol) = @_;
my @dims = $m->dims;
if ($diag){
my ($e, $v) = $m->meigen(0,1);
my ($inv, $info) = $v->minv;
unless ($info){
eval {$v = ($v * $e->$method) x $v->minv;};
if ($@){
warn "mfun: Error $@\n";
return;
}
}
else{
warn "mfun: Non invertible matrix in computation of $method\n";
return;
}
return $m->_is_complex ? $v : toreal($v, $tol);
}
else{
$m = $m->r2C unless $m->_is_complex;
my ($t, undef, $z, undef, $info) = $m->mschur(1);
if ($info){
warn "mfun: Can't compute Schur form\n";
return;
}
($t, $info) = $t->ctrfun(0,$method);
if($info){
warn "mfun: Can't compute $method\n";
return;
}
$m = $z x $t x $z->t(1);
return $m->_is_complex ? $m : toreal($m, $tol);
}
}
=head1 TODO
Improve error return and check singularity.
Improve (msqrt,mlog) / r2C
=head1 AUTHOR
Copyright (C) Grégory Vanuxem 2005-2018.
This library is free software; you can redistribute it and/or modify
it under the terms of the Perl Artistic License as in the file Artistic_2
in this distribution.
=cut
EOD
pp_done();
1;