PDL::LinearAlgebra::Trans - Linear Algebra based transcendental functions for PDL
use PDL::LinearAlgebra::Trans; $a = random (100,100); $sqrt = msqrt($a);
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)
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,phys]A(n,n);int deg();scale();[io]trace();int [o]ns();int [o]info()', 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 dscal_(integer *, double *, double *,
integer *);
extern int dgemm_(char *, char *, integer *, integer *, integer *,
double *, double *, integer *, double *,
integer *, double *, double *, integer *,
ftnlen, ftnlen);
extern int dgesv_(integer *, integer *, double *,
integer *, integer *, double *, integer *, integer *);
extern int daxpy_(integer *, double *, double *,
integer *, double *, integer *);
extern double 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;
integer *ipiv;
double cp, cq, scale, scale2, hnorm, tracen;
double *coef;
double *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 = $PRIV(__n_size) * $PRIV(__n_size); ipiv = (integer *) malloc ( $PRIV(__n_size) * sizeof (integer)); coef = (double *) malloc ( (($deg() + 1) + $PRIV(__n_size) * $PRIV(__n_size)) * sizeof (double)); wsp = (double *) malloc ( 3 * mm * sizeof (double)); //////////
// --- initialise pointers ...
ih2 = $deg() + 1;
ip = 0;
iq = mm;
ifree = iq + mm;
tracen = 0;
if ($trace())
{
for (i__ = 0; i__ < $PRIV(__n_size); i__++)
tracen += h__[i__ + i__ * $PRIV(__n_size)];
tracen /= $PRIV(__n_size);
if (tracen > 0)
{
for (i__ = 0; i__ < $PRIV(__n_size); i__++)
h__[i__ + i__ * $PRIV(__n_size)] -= tracen;
}
}
// --- scaling: seek ns such that ||t*H/2^ns|| < 1/2;
// and set scale = t/2^ns ...
// Compute Infinity norm
hnorm = dlange_("I", &$PRIV(__n_size), &$PRIV(__n_size), $P(A),
&$PRIV(__n_size), wsp);
hnorm = abs($scale() * hnorm);
if (hnorm == 0.)
{
free(ipiv);
free(wsp);
free(coef);
coef = $P(A);
$ns() = 0;
for(i__ = 0 ; i__ < $PRIV(__n_size) ; i__++ ){
for(i__1 = 0 ; i__1 < $PRIV(__n_size) ; i__1++ ){
coef[i__1 + i__ * $PRIV(__n_size) ] = (i__ == i__1) ? 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 ...
dgemm_("n", "n", &$PRIV(__n_size), &$PRIV(__n_size), &$PRIV(__n_size), &scale2, h__, &$PRIV(__n_size), h__,
&$PRIV(__n_size), &c_b7, &coef[ih2], &$PRIV(__n_size), (ftnlen)1, (ftnlen)1);
// --- initialize p (numerator) and q (denominator) ...
cp = coef[$deg() - 1];
cq = coef[$deg()];
for (j = 0; j < $PRIV(__n_size); j++) {
for (i__ = 1; i__ <= $PRIV(__n_size); ++i__) {
wsp[ip + j * $PRIV(__n_size) + i__ - 1] = 0.;
wsp[iq + j * $PRIV(__n_size) + i__ - 1] = 0.;
}
wsp[ip + j * ($PRIV(__n_size) + 1)] = cp;
wsp[iq + j * ($PRIV(__n_size) + 1)] = cq;
}
// --- Apply Horner rule ...
iodd = 1;
k = $deg() - 2;
do{
iused = iodd * iq + (1 - iodd) * ip;
dgemm_("n", "n", &$PRIV(__n_size), &$PRIV(__n_size), &$PRIV(__n_size), &c_b11, &wsp[iused], &$PRIV(__n_size), &coef[ih2], &$PRIV(__n_size), &c_b7, &
wsp[ifree], &$PRIV(__n_size), (ftnlen)1, (ftnlen)1);
for (j = 0; j < $PRIV(__n_size); j++) {
wsp[ifree + j * ($PRIV(__n_size) + 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)
{
dgemm_("n", "n", &$PRIV(__n_size), &$PRIV(__n_size), &$PRIV(__n_size), &scale, &wsp[iq], &$PRIV(__n_size), h__, &$PRIV(__n_size), &
c_b7, &wsp[ifree], &$PRIV(__n_size), (ftnlen)1, (ftnlen)1);
iq = ifree;
}
else
{
dgemm_("n", "n", &$PRIV(__n_size), &$PRIV(__n_size), &$PRIV(__n_size), &scale, &wsp[ip], &$PRIV(__n_size), h__, &$PRIV(__n_size), &
c_b7, &wsp[ifree], &$PRIV(__n_size), (ftnlen)1, (ftnlen)1);
ip = ifree;
}
daxpy_(&mm, &c_b19, &wsp[ip], &c__1, &wsp[iq], &c__1);
dgesv_(&$PRIV(__n_size), &$PRIV(__n_size), &wsp[iq], &$PRIV(__n_size), ipiv, &wsp[ip], &$PRIV(__n_size), $P(info));
if ($info() != 0)
{
free(ipiv);
free(coef);
free(wsp);
}
else
{
dscal_(&mm, &c_b23, &wsp[ip], &c__1);
for (j = 0; j < $PRIV(__n_size); j++) {
wsp[ip + j * ($PRIV(__n_size) + 1)] += 1.;
}
iput = ip;
if ($ns() == 0 && iodd == 1)
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;
dgemm_("n", "n", &$PRIV(__n_size), &$PRIV(__n_size), &$PRIV(__n_size), &c_b11, &wsp[iget], &$PRIV(__n_size), &wsp[iget], &$PRIV(__n_size), &c_b7,
&wsp[iput], &$PRIV(__n_size), (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];
}
free(wsp);
free(ipiv);
}
}
',
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 optionaly 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
$a = random(5,5); $trace = pdl(1); $a->xchg(0,1)->geexp(6,1,$trace, ($ns = null), ($info = null));
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, dcomplex *h__, integer *ldh, dcomplex *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 dcomplex c_b1 = {0.,0.};
static dcomplex c_b2 = {1.,0.};
static integer c__2 = 2;
static integer c__1 = 1;
static double c_b19 = 2.;
static double c_b21 = -1.;
integer *ipiv;
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;
dcomplex scale, tracen, cp, cq, z__1;
dcomplex scale2;
dcomplex *wsp, *h__;
extern int zgemm_(char *, char *, integer *, integer *,
integer *, dcomplex *, dcomplex *, integer *,
dcomplex *, integer *, dcomplex *, dcomplex *,
integer *, ftnlen, ftnlen);
extern int zgesv_(integer *, integer *, dcomplex *,
integer *, integer *, dcomplex *, integer *, integer *);
extern int zaxpy_(integer *, dcomplex *,
dcomplex *, integer *, dcomplex *, integer *), zdscal_(
integer *, double *, dcomplex *, integer *);
extern double zlange_(char *norm, integer *m, integer *n, dcomplex *a,
integer *lda, double *work);
h__ = (dcomplex *) $P(A);
wsp = (dcomplex *) malloc (( $deg() + 1 + 4 * $PRIV(__n_size) * $PRIV(__n_size)) * sizeof (dcomplex));
ipiv = (integer *) malloc ( $PRIV(__n_size) * sizeof (integer));
mm = $PRIV(__n_size) * $PRIV(__n_size);
$info() = 0;
ih2 = $deg() + 1;
ip = ih2 + mm;
iq = ip + mm;
ifree = iq + mm;
tracen.r = 0;
tracen.i = 0;
if ($trace())
{
for (i__ = 0; i__ < $PRIV(__n_size); i__++)
{
tracen.r = tracen.r + h__[i__ + i__ * $PRIV(__n_size)].r;
tracen.i = tracen.i + h__[i__ + i__ * $PRIV(__n_size)].i;
}
tracen.r = tracen.r / $PRIV(__n_size);
tracen.i = tracen.i / $PRIV(__n_size);
if (tracen.r < 0){
tracen.r = tracen.i;
tracen.i = 0;
}
for (i__ = 0; i__ < $PRIV(__n_size); i__++)
{
h__[i__ + i__ * $PRIV(__n_size)].r = h__[i__ + i__ * $PRIV(__n_size)].r - tracen.r ;
h__[i__ + i__ * $PRIV(__n_size)].i = h__[i__ + i__ * $PRIV(__n_size)].i - tracen.i ;
}
}
// --- scaling: seek ns such that ||t*H/2^ns|| < 1/2;
// and set scale = t/2^ns ...
// Compute Infinity norm
hnorm = zlange_("I", &$PRIV(__n_size), &$PRIV(__n_size), &h__[0],
&$PRIV(__n_size), &wsp[0].r);
hnorm = abs($scale() * hnorm);
if (hnorm == 0.)
{
$ns() = 0;
for(i__ = 0 ; i__ < $PRIV(__n_size) ; i__++ ){
for(i__1 = 0 ; i__1 < $PRIV(__n_size) ; i__1++ ){
h__[i__1 + i__ * $PRIV(__n_size) ].r = (i__ == i__1) ? 1 : 0;
h__[i__1 + i__ * $PRIV(__n_size) ].i = 0;
}
}
}
else
{
i__2 = (integer) ((log(hnorm) / log(2.)) + 2);
$ns() = max(0,i__2);
scale.r = $scale() / pow(2, (double )$ns());
scale.i = 0;
scale2.r = scale.r * scale.r - scale.i * scale.i;
scale2.i = scale.r * scale.i + scale.i * scale.r;
// --- compute Pade coefficients ...
i__ = $deg() + 1;
j = ($deg() << 1) + 1;
wsp[0].r = 1;
wsp[0].i = 0;
for (k = 1; k <= $deg(); ++k) {
i__2 = k;
i__3 = k - 1;
d__1 = (double) (i__ - k);
d__2 = (double) (k * (j - k));
wsp[i__2].r = (d__1 * wsp[i__3].r) / d__2;
wsp[i__2].i = (d__1 * wsp[i__3].i) / d__2;
}
// --- H2 = scale2*H*H ...
zgemm_("n", "n", &$PRIV(__n_size), &$PRIV(__n_size), &$PRIV(__n_size), &scale2, &h__[0], &$PRIV(__n_size), &h__[0],
&$PRIV(__n_size), &c_b1, &wsp[ih2], &$PRIV(__n_size), (ftnlen)1, (ftnlen)1);
// --- initialise p (numerator) and q (denominator) ...
i__1 = $deg() - 1;
cp.r = wsp[i__1].r;
cp.i = wsp[i__1].i;
i__1 = $deg();
cq.r = wsp[i__1].r;
cq.i = wsp[i__1].i;
for (j = 0; j < $PRIV(__n_size); j++) {
for (i__ = 1; i__ <= $PRIV(__n_size); ++i__) {
i__3 = ip + j * $PRIV(__n_size) + i__ - 1;
wsp[i__3].r = 0., wsp[i__3].i = 0.;
i__3 = iq + j * $PRIV(__n_size) + i__ - 1;
wsp[i__3].r = 0., wsp[i__3].i = 0.;
}
i__2 = ip + j * ($PRIV(__n_size) + 1);
wsp[i__2].r = cp.r;
wsp[i__2].i = cp.i;
i__2 = iq + j * ($PRIV(__n_size) + 1);
wsp[i__2].r = cq.r;
wsp[i__2].i = cq.i;
}
// --- Apply Horner rule ...
iodd = 1;
k = $deg() - 2;
do
{
iused = iodd * iq + (1 - iodd) * ip;
zgemm_("n", "n", &$PRIV(__n_size), &$PRIV(__n_size), &$PRIV(__n_size), &c_b2, &wsp[iused], &$PRIV(__n_size), &wsp[ih2], &$PRIV(__n_size), &c_b1, &
wsp[ifree], &$PRIV(__n_size), (ftnlen)1, (ftnlen)1);
for (j = 0; j < $PRIV(__n_size); j++) {
i__1 = ifree + j * ($PRIV(__n_size) + 1);
i__2 = ifree + j * ($PRIV(__n_size) + 1);
i__3 = k;
z__1.r = wsp[i__2].r + wsp[i__3].r;
z__1.i = wsp[i__2].i + wsp[i__3].i;
wsp[i__1].r = z__1.r;
wsp[i__1].i = z__1.i;
}
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)
{
zgemm_("n", "n", &$PRIV(__n_size), &$PRIV(__n_size), &$PRIV(__n_size), &scale, &wsp[iq], &$PRIV(__n_size), &h__[0], &$PRIV(__n_size), &
c_b1, &wsp[ifree], &$PRIV(__n_size), (ftnlen)1, (ftnlen)1);
iq = ifree;
}else
{
zgemm_("n", "n", &$PRIV(__n_size), &$PRIV(__n_size), &$PRIV(__n_size), &scale, &wsp[ip], &$PRIV(__n_size), &h__[0], &$PRIV(__n_size), &
c_b1, &wsp[ifree], &$PRIV(__n_size), (ftnlen)1, (ftnlen)1);
ip = ifree;
}
z__1.r = -1.;
z__1.i = -0.;
zaxpy_(&mm, &z__1, &wsp[ip], &c__1, &wsp[iq], &c__1);
zgesv_(&$PRIV(__n_size), &$PRIV(__n_size), &wsp[iq], &$PRIV(__n_size), ipiv, &wsp[ip], &$PRIV(__n_size), $P(info));
if ($info() == 0)
{
zdscal_(&mm, &c_b19, &wsp[ip], &c__1);
for (j = 0; j < $PRIV(__n_size); j++)
{
i__2 = ip + j * ($PRIV(__n_size) + 1);
i__3 = ip + j * ($PRIV(__n_size) + 1);
wsp[i__2].r = wsp[i__3].r + 1.;
wsp[i__2].i = wsp[i__3].i + 0.;
}
iput = ip;
if ($ns() == 0 && iodd != 0)
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;
zgemm_("n", "n", &$PRIV(__n_size), &$PRIV(__n_size), &$PRIV(__n_size), &c_b2, &wsp[iget], &$PRIV(__n_size), &wsp[iget], &$PRIV(__n_size), &c_b1,
&wsp[iput], &$PRIV(__n_size), (ftnlen)1, (ftnlen)1);
iodd = 1 - iodd;
}
}
}
i__2 = iput + mm;
i__ = 0;
if ($trace())
{
// exp(tracen)
d__1 = exp(tracen.r);
tracen.r = d__1 * cos(tracen.i);
tracen.i = d__1 * sin(tracen.i);
for (i__1 = iput; i__1 < i__2 ; i__1++)
{
h__[i__].r = tracen.r * wsp[i__1].r - tracen.i * wsp[i__1].i;
h__[i__++].i = tracen.r * wsp[i__1].i + tracen.i * wsp[i__1].r;
}
}
else
{
for (i__1 = iput; i__1 < i__2 ; i__1++)
{
h__[i__].r = wsp[i__1].r;
h__[i__++].i = wsp[i__1].i;
}
}
}
free(wsp);
free(ipiv);
',
);
pp_def('ctrsqrt', Pars => '[io,phys]A(2,n,n);int uplo();[phys,o] B(2,n,n);int [o]info()', Doc => '
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, Manchester, England, January 1999. It\'s available at http://www.ma.man.ac.uk/~higham/pap-mf.html) If uplo is true, A is lower triangular.
', GenericTypes => [D], Code => ' dcomplex *cb, *ca; dcomplex s, snum, sdenum; double *b; double tt, dn; integer i, j, k, l, ind, ind1, ind2;
ca = (dcomplex *) $P(A);
b = $P(B);
cb = (dcomplex *) $P(B);
#define subscr(a_1,a_2)\ ($uplo()) ? (a_2) * $PRIV(__n_size) + (a_1) : (a_1) * $PRIV(__n_size) + (a_2)
$info() = 0;
ind = $PRIV(__n_size) * $PRIV(__n_size) * 2;
for (i = 0; i < ind;i++)
b[i] = 0;
for (i = 0; i < $PRIV(__n_size); i++){
ind = subscr(i,i);
CSQRT(ca[ind].r, ca[ind].i, cb[ind].r, cb[ind].i)
}
for (k = 0; k < $PRIV(__n_size)-1; k++)
{
for (i = 0; i < $PRIV(__n_size)-(k+1); i++)
{
j = i + k + 1;
ind = subscr(i,j);
s.r = 0;
s.i = 0;
for (l = i+1; l < j; l++)
{
ind1 = subscr(i,l);
ind2 = subscr(l,j);
s.r += cb[ind1].r * cb[ind2].r - cb[ind1].i * cb[ind2].i;
s.i += cb[ind1].r * cb[ind2].i + cb[ind1].i * cb[ind2].r;
}
ind1 = subscr(i,i);
ind2 = subscr(j,j);
sdenum.r = cb[ind1].r + cb[ind2].r;
sdenum.i = cb[ind1].i + cb[ind2].i;
snum.r = ca[ind].r - s.r;
snum.i = ca[ind].i - s.i;
if (fabs (sdenum.r) > fabs (sdenum.i))
{
if (sdenum.r == 0){
$info() = -1;
goto ENDCTRSQRT;
}
tt = sdenum.i / sdenum.r;
dn = sdenum.r + tt * sdenum.i;
cb[ind].r = (snum.r + tt * snum.i) / dn;
cb[ind].i = (snum.i - tt * snum.r) / dn;
}
else
{
if (sdenum.i == 0){
$info() = -1;
goto ENDCTRSQRT;
}
tt = sdenum.r / sdenum.i;
dn = sdenum.r * tt + sdenum.i;
cb[ind].r = (snum.r * tt + snum.i) / dn;
cb[ind].i = (snum.i * tt - snum.r) / dn;
}
}
}
ENDCTRSQRT:
;
#undef subscr
',
);
pp_addhdr('
void dfunc_wrapper(dcomplex *p, integer n, SV* dfunc) { dSP ;
int count;
SV *pdl1;
HV *bless_stash;
pdl *pdl;
PDL_Long odims[1];
PDL_Long dims[] = {2,n};
pdl = PDL->pdlnew();
PDL->setdims (pdl, dims, 2);
pdl->datatype = PDL_D;
pdl->data = (double *) &p[0].r;
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(dfunc, G_SCALAR);
SPAGAIN;
// For pdl_free
odims[0] = 0;
PDL->setdims (pdl, odims, 0);
pdl->state &= ~ (PDL_ALLOCATED |PDL_DONTTOUCHDATA);
pdl->data=NULL;
if (count !=1)
croak("Error calling perl function\n");
PUTBACK ; FREETMPS ; LEAVE ;
}
');
pp_def('ctrfun', Pars => '[io,phys]A(2,n,n);int uplo();[phys,o] B(2,n,n);int [o]info()', OtherPars => "SV* func" , Doc => '
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 => ' dcomplex *cb, *ca, *diag; dcomplex s, snum, sdenum; double *b; double tt, dn; integer i, j, k, l, p, ind, ind1, ind2, tmp, c__1; extern double zdotu_(double *ret, integer *n, double *dx, integer *incx, double *dy, integer *incy); diag = (dcomplex *) malloc ($PRIV(__n_size) * sizeof(dcomplex));
ca = (dcomplex *) $P(A);
b = $P(B);
cb = (dcomplex *) $P(B);
c__1 = 1;
#define subscr(a_1,a_2) \ ( $uplo() ) ? (a_2)+ $PRIV(__n_size) * (a_1) : (a_1)+ $PRIV(__n_size) * (a_2)
$info() = 0;
ind = $PRIV(__n_size) * $PRIV(__n_size) * 2;
for (i = 0; i < ind;i++)
b[i] = 0;
for (i = 0; i < $PRIV(__n_size);i++)
{
ind = subscr(i,i);
diag[i].r = ca[ind].r;
diag[i].i = ca[ind].i;
}
dfunc_wrapper(diag, $PRIV(__n_size), $PRIV(func));
for (i = 0; i < $PRIV(__n_size); i++)
{
ind = subscr(i,i);
cb[ind].r = diag[i].r;
cb[ind].i = diag[i].i;
}
for (p = 1; p < $PRIV(__n_size); p++)
{
tmp = $PRIV(__n_size) - 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);
s.r = cb[ind2].r - cb[ind1].r;
s.i = cb[ind2].i - cb[ind1].i;
ind = subscr(j,i);
CMULT(ca[ind].r, ca[ind].i, s.r, s.i, snum.r, snum.i)
s.r = snum.r;
s.i = snum.i;
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++)
{
ind = subscr(j,k);
diag[l].r = cb[ind].r;
diag[l++].i = cb[ind].i;
}
ind = subscr(i+1,i);
// TODO : Humm
zdotu_(&s.r, &l, &ca[ind].r, &c__1, &diag[0].r, &c__1);
snum.r += s.r;
snum.i += s.i;
// $F(,$i+1:$j-1,($i))->cdot(1,$T(,($j), $i+1:$j-1),1)
l = 0;
for(k = i+1; k < j; k++)
{
ind = subscr(j,k);
diag[l].r = ca[ind].r;
diag[l++].i = ca[ind].i;
}
ind = subscr(i+1,i);
// TODO : Humm
zdotu_(&s.r, &l, &cb[ind].r, &c__1, &diag[0].r, &c__1);
snum.r -= s.r;
snum.i -= s.i;
ind = subscr(j,i);
}
//$sdenum = $T(,($j),($j))->Csub($T(,($i),($i)));
sdenum.r = ca[ind2].r - ca[ind1].r;
sdenum.i = ca[ind2].i - ca[ind1].i;
//$s = $s / $sdenum;
if (fabs (sdenum.r) > fabs (sdenum.i))
{
if (sdenum.r == 0)
{
$info() = -1;
goto ENDCTRFUN;
}
tt = sdenum.i / sdenum.r;
dn = sdenum.r + tt * sdenum.i;
cb[ind].r = (snum.r + tt * snum.i) / dn;
cb[ind].i = (snum.i - tt * snum.r) / dn;
}
else
{
if (sdenum.i == 0)
{
$info() = -1;
goto ENDCTRFUN;
}
tt = sdenum.r / sdenum.i;
dn = sdenum.r * tt + sdenum.i;
cb[ind].r = (snum.r * tt + snum.i) / dn;
cb[ind].i = (snum.i * tt - snum.r) / 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{ return $_[0] if $_[0]->isempty; $_tol = $_[1] if defined $_[1]; my ($min, $max, $tmp); ($min, $max) = $_[0]->slice('(1)')->minmax; return re($_[0])->sever unless (abs($min) > $_tol || abs($max) > $_tol); $_[0]; }
Return matrix logarithm of a square matrix.
PDL = mlog(PDL(A))
my $a = random(10,10); my $log = mlog($a);
Return matrix square root (principal) of a square matrix.
PDL = msqrt(PDL(A))
my $a = random(10,10); my $sqrt = msqrt($a);
Return matrix exponential of a square matrix.
PDL = mexp(PDL(A))
my $a = random(10,10); my $exp = mexp($a);
Return matrix power of a square matrix.
PDL = mpow(PDL(A), SCALAR(exponent))
my $a = random(10,10); my $powered = mpow($a,2.5);
Return matrix cosine of a square matrix.
PDL = mcos(PDL(A))
my $a = random(10,10); my $cos = mcos($a);
Return matrix inverse cosine of a square matrix.
PDL = macos(PDL(A))
my $a = random(10,10); my $acos = macos($a);
Return matrix sine of a square matrix.
PDL = msin(PDL(A))
my $a = random(10,10); my $sin = msin($a);
Return matrix inverse sine of a square matrix.
PDL = masin(PDL(A))
my $a = random(10,10); my $asin = masin($a);
Return matrix tangent of a square matrix.
PDL = mtan(PDL(A))
my $a = random(10,10); my $tan = mtan($a);
Return matrix inverse tangent of a square matrix.
PDL = matan(PDL(A))
my $a = random(10,10); my $atan = matan($a);
Return matrix cotangent of a square matrix.
PDL = mcot(PDL(A))
my $a = random(10,10); my $cot = mcot($a);
Return matrix inverse cotangent of a square matrix.
PDL = macot(PDL(A))
my $a = random(10,10); my $acot = macot($a);
Return matrix secant of a square matrix.
PDL = msec(PDL(A))
my $a = random(10,10); my $sec = msec($a);
Return matrix inverse secant of a square matrix.
PDL = masec(PDL(A))
my $a = random(10,10); my $asec = masec($a);
Return matrix cosecant of a square matrix.
PDL = mcsc(PDL(A))
my $a = random(10,10); my $csc = mcsc($a);
Return matrix inverse cosecant of a square matrix.
PDL = macsc(PDL(A))
my $a = random(10,10); my $acsc = macsc($a);
Return matrix hyperbolic cosine of a square matrix.
PDL = mcosh(PDL(A))
my $a = random(10,10); my $cos = mcosh($a);
Return matrix hyperbolic inverse cosine of a square matrix.
PDL = macosh(PDL(A))
my $a = random(10,10); my $acos = macosh($a);
Return matrix hyperbolic sine of a square matrix.
PDL = msinh(PDL(A))
my $a = random(10,10); my $sinh = msinh($a);
Return matrix hyperbolic inverse sine of a square matrix.
PDL = masinh(PDL(A))
my $a = random(10,10); my $asinh = masinh($a);
Return matrix hyperbolic tangent of a square matrix.
PDL = mtanh(PDL(A))
my $a = random(10,10); my $tanh = mtanh($a);
Return matrix hyperbolic inverse tangent of a square matrix.
PDL = matanh(PDL(A))
my $a = random(10,10); my $atanh = matanh($a);
Return matrix hyperbolic cotangent of a square matrix.
PDL = mcoth(PDL(A))
my $a = random(10,10); my $coth = mcoth($a);
Return matrix hyperbolic inverse cotangent of a square matrix.
PDL = macoth(PDL(A))
my $a = random(10,10); my $acoth = macoth($a);
Return matrix hyperbolic secant of a square matrix.
PDL = msech(PDL(A))
my $a = random(10,10); my $sech = msech($a);
Return matrix hyperbolic inverse secant of a square matrix.
PDL = masech(PDL(A))
my $a = random(10,10); my $asech = masech($a);
Return matrix hyperbolic cosecant of a square matrix.
PDL = mcsch(PDL(A))
my $a = random(10,10); my $csch = mcsch($a);
Return matrix hyperbolic inverse cosecant of a square matrix.
PDL = macsch(PDL(A))
my $a = random(10,10); my $acsch = macsch($a);
Return matrix function of second argument of a square matrix. Function will be applied on a PDL::Complex object.
PDL = mfun(PDL(A),'cos')
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);
Improve error return and check singularity. Improve (msqrt,mlog) / r2C
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.
