On Wed, Dec 09, 2015 at 09:52:22AM +0100, Bryan Jurish wrote:
> morning all,
>
> I've just released a new module PDL::SVDSLEPc
> ...
> <https://metacpan.org/release/PDL-SVDLIBC> on my machine here (I believe
> it's using ARPACK under the hood, which was mentioned on this list a few
> years back, but my searches have failed to turn up any PDL-ARPACK interface
> code).
I used arpack in May 2014 and posted a message (below) and an example
(attached) to the list:
On Mon, 12 May 2014 14:16:09 -0500, Luis Mochan wrote:
> I once made some pdl code that interfaces with ARPACK to calculate
> photonic bands, but it may be used for eigenvalues of sparse matrices
> or of arbitrary functions which perform a linear transformation (i.e.,
> not necessarily using a matrix to respresent the transformation). I
> enclose an example that calculates the eigenvalues of the a second
> diference operator (discrete analogy to the second derivative, as in
> 1D phonon calculations) and compares with the analytical solution (I
> don't remember the details). Hope the interface proves useful.
Best regards,
Luis
--
o
W. Luis Mochán, | tel:(52)(777)329-1734 /<(*)
Instituto de Ciencias Físicas, UNAM | fax:(52)(777)317-5388 `>/ /\
Apdo. Postal 48-3, 62251 | (*)/\/ \
Cuernavaca, Morelos, México | [email protected] /\_/\__/
GPG: 791EB9EB, C949 3F81 6D9B 1191 9A16 C2DF 5F0A C52B 791E B9EB
#!/usr/bin/env perl
# arpack for photonic bands with complex eigenvalues
use strict;
use warnings;
use feature q(say);
use PDL;
use PDL::NiceSlice;
use PDL::Graphics::PGPLOT;
use PDL::Complex;
*znaupdPP=\&PDL::znaupdPP;
*zneupdPP=\&PDL::zneupdPP;
my $N=100; #size of system
my $nev=50; #eigenvalues
my $ncv=2*$nev; #number of Lanczos vectors
# Mode=1 A*x = lambda*x, A symmetric
#show analytical solution (divide by 2 to account for degeneracy)
my $PI=4*atan2(1,1);
my $seq=2*$PI*(sequence($N))->qsort/$N/2;
my $res=sqrt(2*(1-cos($seq)));
line $res;
hold;
# setup arguments to eigen solver
my @args=(
bmat=>'I', #standard problem
n=>$N, #size of system
nev=>$nev, #how many eigenvalues
tol=>0, #1e-8, #default tolerance 0=> machine precission
ncv=>$ncv, #Number of lanczos vectors
ishift=> 1, #kind of shift
mxiter=> 1250, #max number of iterations
info=>0, #random initial residual
function=>\&func #function to multiply matrix.vector
);
#find smallest (algebraic) eigenvalues
# Which may be LM SM LR SR LI SI
my ($values,$vectors)= areigens(@args, which=>'SR');
my $realpart=$values((0))->qsort;
printf "Eigenvalues=$realpart\n";
#printf "Eigenvectors=$vectors\n";
line $realpart->xvals+$N-$nev, sqrt(-$realpart(-1:0));
#find largest (algebraic) eigenvalues
my ($values1,$vectors1)= areigens(@args, which=>'LR');
my $realpart1=$values1((0))->qsort;
printf "Eigenvalues=$realpart1\n";
$realpart1->where($realpart1>0).=0;
#printf "Eigenvectors=$vectors1\n";
line $realpart1->xvals, sqrt(-$realpart1(-1:0));
print "listo?\n";
my $r=<>;
#interface
sub areigens {
my $mode=1; #solve only type 1 problems, for the time being
#First change hash argument into parameters
my %arg=@_;
my $ido=0;
my $bmat=uc $arg{bmat};
my $n=$arg{n}; #size of matrix
my $which=$arg{which};
my $nev=$arg{nev};
my $tol=$arg{tol};
my $resid=$arg{resid} || zeroes(double, 2, $n); #input/output residuals
my $ncv=$arg{ncv}; #desired lanczos vectors
my $v=zeroes(double, 2, $n, $ncv); #space for vectors
# my $ldv=$n; # leading dimension of $v
my $ishift=$arg{ishift}; #kind of shifts
my $mxiter=$arg{mxiter}; #Max num. of arnoldi iterations
my $iparam=pdl(long, $ishift, 0, $mxiter, 1, 0, 0, $mode, 0, 0, 0, 0);
my $ipntr=zeroes(long,14); #output
my $workd=zeroes(double, 2, 3*$n); # work array
my $lworkl =$ncv*(3*$ncv+5);
my $workl=zeroes(double, 2, $lworkl); #work space
my $rwork=zeroes(double, $ncv);
my $info=$arg{info}; #kind of initial residual, error flag
my $function=$arg{function};
# test inputs
die "BMAT must be I or G, not $bmat" unless $bmat eq 'I' or $bmat eq 'G';
my $t=0;
$t ||= $which eq $_ foreach qw(LM SM LR SR LI SI);
die "WHICH must be any of LM SM LR SR LI SI, not $which" unless $t;
die "0<NEV<N violated by NEV=$nev" unless 0<$nev and $nev <= $n;
$t=$resid->dim(1);
die "RESID must be of at least size N, not $t" unless $n<=$t;
die "NEV<NCV<=N violated by $ncv" unless $nev<$ncv and $ncv<=$n;
die "ISHIFT must be 0 or 1" unless $ishift==0 or $ishift == 1;
die "MXITER must be nonnegative" unless 0<=$mxiter;
die "MODE must be 1 for the time being" unless $mode == 1;
#convert whatever is possible to a pdl
($ido, $n, $nev, $info)=
map long($_), ($ido, $n, $nev, $info);
$tol=double($tol);
while(1){
znaupdPP($ido, $nev, $tol, $resid, $v, $iparam,
$ipntr, $workd, $workl, $rwork, $info, $bmat, $which);
last if $ido == 99;
die "Unexpected \$ido=$ido" unless $ido==1 or $ido==-1;
my $offx=$ipntr((1-1))-1; #careful with offsets
my $offy=$ipntr((2-1))-1;
my $x=$workd(:, $offx:$offx+$n-1); # vector to operate on
my $y=$workd(:, $offy:$offy+$n-1); # vector to leave result
$y.=&$function($x); #operate on $x and leave result in $y
#within workd;
}
die "Error in znaupd, \$info=$info" unless $info == 0;
#Obtain the eigenvalues and eigenvectors
my $rvec=long(1); #require eigenvectors
my $howmny='A'; #compute all nev eigenvectors
my $select=ones(long, $ncv); #workspace or selection criteria
my $d=zeroes(double, 2, $nev); #workspace for eigenvalues (nev+1???)
my $z=zeroes(double, 2, $n,$nev); #eigenvectors
my $sigma=zeroes(2); #shift. unused
my $workev=zeroes(double,2,2*$ncv);
zneupdPP($rvec, $select, $d, $z, $sigma, $workev, $tol, $resid, $v,
$iparam, $ipntr, $workd, $workl, $rwork, $info,
$howmny, $bmat, $which);
die "Error in zneupd, \$info=$info" unless $info == 0;
return ($d, $z);
}
sub func {
# function to test dsaupd. Eigenvectors are cartesian
# my $vector=shift @_;
# my $a=sequence($vector->dim(0))+1;
# my $y= $a*$vector; # scale along axes
# return $y;
#second derivative: periodic boundaries
my $vector=(shift @_)->mv(0,1); #mv r-i part to simplify rotations below
# might need to cast to complex value with cplx (with dataflow)
# before making complex operations. Unneded in simple cases.
return ($vector->rotate(1)-2*$vector+$vector->rotate(-1))->mv(1,0);
}
#call dsaupd
#c ( IDO, BMAT, N, WHICH, NEV, TOL, RESID, NCV, V, LDV, IPARAM,
#c IPNTR, WORKD, WORKL, LWORKL, INFO )
#c
no PDL::NiceSlice;
use Inline Pdlpp=>Config=>
#MYEXTLIB=>'/home/mochan/txt/cache/12/arpack/ARPACK/libarpack_Debian.a',
LIBS=>'-larpack -lgfortran';
use Inline Pdlpp => <<'FIN_PP';
pp_def('znaupdPP',
Pars=>'int [o] ido(); int nev(); tol(); [o] resid(ri=2, n);'
.'[o] v(ri=2, n,ncv); int [o] iparam(l=11); '
.'int [o] ipntr(lp=14); [o] workd(ri=2,tn); [o] workl(ri=2,lw);'
.'[o] rwork(ncv); int [o] info();',
OtherPars=>'char *bmat; char *which;',
GenericTypes=>['D'],
Code=> q{
int N=$SIZE(n);
int LDV=N;
int NCV=$SIZE(ncv);
int LWORKL=$SIZE(lw);
znaupd_($P(ido), $COMP(bmat), &N, $COMP(which), $P(nev),
$P(tol), $P(resid), &NCV, $P(v), &LDV, $P(iparam),
$P(ipntr), $P(workd), $P(workl), &LWORKL, $P(rwork),
$P(info));
},
);
pp_def('zneupdPP',
Pars=>'int rvec(); int select(ncv); [o] d(ri=2, nev); '
.'[o] z(ri=2, n,nev); sigma(ri=2);'
.'workev(ri=2,ncv2); tol(); resid(ri=2,n); '
.'v(ri=2,n,ncv); int iparam(l=11);'
.'int ipntr(lp=14); workd(ri=2,tn); workl(ri=2,lw);'
.'rwork(ncv); int [o] info();',
OtherPars=>'char *howmny; char *bmat; char *which; ',
GenericTypes=>['D'],
Code=> q{
int N=$SIZE(n);
int LDZ=N;
int NCV=$SIZE(ncv);
int NEV=$SIZE(nev);
int LDV=N;
int LWORKL=$SIZE(lw);
zneupd_($P(rvec), $COMP(howmny), $P(select), $P(d), $P(z),
&LDZ, $P(sigma), $P(workev), $COMP(bmat), &N, $COMP(which),
&NEV, $P(tol),$P(resid), &NCV, $P(v), &LDV,
$P(iparam), $P(ipntr), $P(workd), $P(workl),
&LWORKL, $P(rwork), $P(info));
},
);
FIN_PP
------------------------------------------------------------------------------
_______________________________________________
pdl-general mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/pdl-general
