diff --git a/Changes b/Changes index 3e12fb0..7232ccb 100644 --- a/Changes +++ b/Changes @@ -1,3 +1,5 @@ +- remove support for PDL::Complex entirely + 0.40 2024-09-06 - moved augment, mstack, tricpy, t to PDL 2.091 diff --git a/Complex/complex.pd b/Complex/complex.pd index 3a635c2..645ffe0 100644 --- a/Complex/complex.pd +++ b/Complex/complex.pd @@ -80,22 +80,6 @@ pp_addpm({At=>'Top'},<<'EOD'); use strict; use PDL::LinearAlgebra::Real; -{ - package # hide from CPAN - 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 ; } -} - =encoding utf8 =head1 NAME @@ -325,7 +309,7 @@ EOF rwork = ($GENERIC() *)malloc( 7 * lwork * sizeof($GENERIC())); break; - } + } lwork = -1; FORTRAN($TFD(c,z)gesdd)( &tra, @@ -383,9 +367,9 @@ The SVD is written pp_defc("ggsvd", _decl => <<'EOF', extern int FORTRAN($TFD(c,z)ggsvd3)(char *jobu, char *jobv, char *jobq, integer *m, - integer *n, integer *p, integer *k, integer *l, $GENERIC() *a, - integer *lda, $GENERIC() *b, integer *ldb, $GENERIC() *alpha, - $GENERIC() *beta, $GENERIC() *u, integer *ldu, $GENERIC() *v, integer + integer *n, integer *p, integer *k, integer *l, $GENERIC() *a, + integer *lda, $GENERIC() *b, integer *ldb, $GENERIC() *alpha, + $GENERIC() *beta, $GENERIC() *u, integer *ldu, $GENERIC() *v, integer *ldv, $GENERIC() *q, integer *ldq, $GENERIC() *work, integer *lwork, $GENERIC() *rwork, integer *iwork, integer *info); EOF @@ -398,7 +382,7 @@ EOF ', Code => generate_code ' char pjobu = \'N\'; - char pjobv = \'N\'; + char pjobv = \'N\'; char pjobq = \'N\'; $GENERIC() *work, twork[2]; integer lwork = -1; @@ -439,7 +423,7 @@ EOF lwork = (integer) twork[0]; work = ($GENERIC() *)malloc(2*(lwork * sizeof($GENERIC()))); - + FORTRAN($TFD(c,z)ggsvd3)( &pjobu, &pjobv, @@ -870,7 +854,7 @@ EOF &lwork, $P(rwork), $P(bwork), - $P(info)); + $P(info)); lwork = (integer )tmp_work[0]; @@ -903,13 +887,13 @@ Complex version of L 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; + select_func(complex(w)) is true; Note that a selected complex eigenvalue may no longer - satisfy select_func(PDL::Complex(w)) = 1 after ordering, since + satisfy select_func(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. - + '); @@ -1009,9 +993,9 @@ Complex version of L 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; + select_func(complex(w)) is true; Note that a selected complex eigenvalue may no longer - satisfy select_func(PDL::Complex(w)) = 1 after ordering, since + satisfy select_func(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. @@ -1114,9 +1098,9 @@ Complex version of L 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; + select_func(complex(w), 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 + satisfy select_func(complex(w),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. @@ -1232,9 +1216,9 @@ Complex version of L 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; + select_func(complex(w), 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 + satisfy select_func(complex(w),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. @@ -1670,8 +1654,8 @@ EOF $P(w), &tmp_work[0], &lwork, - &tmp_rwork, - &lrwork, + &tmp_rwork, + &lrwork, &tmp_iwork, &liwork, $P(info)); @@ -1924,9 +1908,9 @@ Doc => ' Complex version of L. trans: Specifies the form of the system of equations: - = 0: A * X = B (No transpose) + = 0: A * X = B (No transpose) = 1: A\' * X = B (Transpose) - = 2: A**H * X = B (Conjugate transpose) + = 2: A**H * X = B (Conjugate transpose) '); pp_defc("sysv", @@ -2321,7 +2305,7 @@ EOF Doc=>' =for ref -Solves overdetermined or underdetermined complex linear systems +Solves overdetermined or underdetermined complex linear systems involving an M-by-N matrix A, or its conjugate-transpose. Complex version of L. @@ -2981,8 +2965,8 @@ EOF Complex version of L - Arguments - ========= + Arguments + ========= trans: = 0: No transpose; = 1: Transpose; = 2: Conjugate transpose; @@ -3112,15 +3096,15 @@ EOF $P(B), &(integer){$SIZE(n)}, $P(info)); -', +', Doc=>' =for ref Complex version of L - Arguments - ========= + Arguments + ========= trans: = 0: No transpose; = 1: Transpose; = 2: Conjugate transpose; @@ -3131,7 +3115,7 @@ Complex version of L pp_defc("latrs", _decl => <<'EOF', extern int FORTRAN($TFD(c,z)latrs)(char *uplo, char *trans, char *diag, char * - normin, integer *n, $GENERIC() *a, integer *lda, $GENERIC() *x, + normin, integer *n, $GENERIC() *a, integer *lda, $GENERIC() *x, $GENERIC() *scale, $GENERIC() *cnorm, integer *info); EOF HandleBad => 0, @@ -3172,8 +3156,8 @@ EOF Complex version of L - Arguments - ========= + Arguments + ========= trans: = 0: No transpose; = 1: Transpose; = 2: Conjugate transpose; @@ -3397,7 +3381,7 @@ EOF pp_defc("ungqr", _decl => <<'EOF', extern int FORTRAN($TFD(c,z)ungqr)(integer *m, integer *n, integer *k, $GENERIC() * - a, integer *lda, $GENERIC() *tau, $GENERIC() *work, integer *lwork, + a, integer *lda, $GENERIC() *tau, $GENERIC() *work, integer *lwork, integer *info); EOF HandleBad => 0, @@ -3540,7 +3524,7 @@ EOF pp_defc("unglq", _decl => <<'EOF', extern int FORTRAN($TFD(c,z)unglq)(integer *m, integer *n, integer *k, $GENERIC() * - a, integer *lda, $GENERIC() *tau, $GENERIC() *work, integer *lwork, + a, integer *lda, $GENERIC() *tau, $GENERIC() *work, integer *lwork, integer *info); EOF HandleBad => 0, @@ -3682,7 +3666,7 @@ EOF pp_defc("ungql", _decl => <<'EOF', extern int FORTRAN($TFD(c,z)ungql)(integer *m, integer *n, integer *k, $GENERIC() * - a, integer *lda, $GENERIC() *tau, $GENERIC() *work, integer *lwork, + a, integer *lda, $GENERIC() *tau, $GENERIC() *work, integer *lwork, integer *info); EOF HandleBad => 0, @@ -3822,7 +3806,7 @@ EOF pp_defc("ungrq", _decl => <<'EOF', extern int FORTRAN($TFD(c,z)ungrq)(integer *m, integer *n, integer *k, $GENERIC() * - a, integer *lda, $GENERIC() *tau, $GENERIC() *work, integer *lwork, + a, integer *lda, $GENERIC() *tau, $GENERIC() *work, integer *lwork, integer *info); EOF HandleBad => 0, @@ -4025,7 +4009,7 @@ Complex version of L. Here trans = 1 means conju pp_defc("gehrd", _decl => <<'EOF', extern int FORTRAN($TFD(c,z)gehrd)(integer *n, integer *ilo, integer *ihi, - $GENERIC() *a, integer *lda, $GENERIC() *tau, $GENERIC() *work, + $GENERIC() *a, integer *lda, $GENERIC() *tau, $GENERIC() *work, integer *lwork, integer *info); EOF HandleBad => 0, @@ -4064,7 +4048,7 @@ EOF pp_defc("unghr", _decl => <<'EOF', extern int FORTRAN($TFD(c,z)unghr)(integer *n, integer *ilo, integer *ihi, - $GENERIC() *a, integer *lda, $GENERIC() *tau, $GENERIC() *work, + $GENERIC() *a, integer *lda, $GENERIC() *tau, $GENERIC() *work, integer *lwork, integer *info); EOF HandleBad => 0, @@ -4110,7 +4094,7 @@ pp_defc("hseqr", _decl => <<'EOF', extern int FORTRAN($TFD(c,z)hseqr)(char *job, char *compz, integer *n, integer *ilo, integer *ihi, $GENERIC() *h__, integer *ldh, $GENERIC() *w, - $GENERIC() *z__, integer *ldz, $GENERIC() *work, + $GENERIC() *z__, integer *ldz, $GENERIC() *work, integer *lwork, integer *info); EOF HandleBad => 0, @@ -4172,7 +4156,7 @@ pp_defc("trevc", _decl => <<'EOF', extern int FORTRAN($TFD(c,z)trevc)(char *side, char *howmny, logical *select, integer *n, $GENERIC() *t, integer *ldt, $GENERIC() *vl, integer * - ldvl, $GENERIC() *vr, integer *ldvr, integer *mm, integer *m, + ldvl, $GENERIC() *vr, integer *ldvr, integer *mm, integer *m, $GENERIC() *work, $GENERIC() *rwork, integer *info); EOF HandleBad => 0, @@ -4214,7 +4198,7 @@ EOF $P(VR), &(integer){$SIZE(p)}, &mm, - $P(m), + $P(m), $P(work) + $SIZE(n), $P(work), $P(info)); @@ -4224,7 +4208,7 @@ pp_defc("tgevc", _decl => <<'EOF', extern int FORTRAN($TFD(c,z)tgevc)(char *side, char *howmny, logical *select, integer *n, $GENERIC() *a, integer *lda, $GENERIC() *b, integer *ldb, - $GENERIC() *vl, integer *ldvl, $GENERIC() *vr, integer *ldvr, + $GENERIC() *vl, integer *ldvl, $GENERIC() *vr, integer *ldvr, integer *mm, integer *m, $GENERIC() *work, $GENERIC() *rwork, integer *info); EOF HandleBad => 0, @@ -4266,7 +4250,7 @@ EOF $P(VR), &(integer){$SIZE(p)}, &mm, - $P(m), + $P(m), $P(work)+2*$SIZE(n), $P(work), $P(info)); @@ -4455,8 +4439,8 @@ EOF Complex version of L. - Arguments - ========= + Arguments + ========= transa: = 0: No transpose; = 1: Transpose; = 2: Conjugate transpose; @@ -4517,8 +4501,8 @@ Doc=>' Complex version of L. - Arguments - ========= + Arguments + ========= transa: = 0: No transpose; = 1: Transpose; = 2: Conjugate transpose; @@ -4736,7 +4720,7 @@ EOF Doc=>' =for ref -Forms the dot product of two vectors, conjugating the first +Forms the dot product of two vectors, conjugating the first vector. '); @@ -4835,7 +4819,7 @@ EOF $P(a), $P(b), $P(c), - $P(s) + $P(s) ); '); @@ -4906,8 +4890,8 @@ EOF Code => ' int i,j,k; $GENERIC() tr[2], b[2]; - //$GENERIC() *tmp; - char ptrans = \'N\'; + //$GENERIC() *tmp; + char ptrans = \'N\'; $GENERIC() alpha[2] = {1,0}; $GENERIC() beta[2] = {0,0}; loop(n0,n1) %{ diff --git a/Complex/selectfunc.c b/Complex/selectfunc.c index aeb45f6..5f1e1cd 100644 --- a/Complex/selectfunc.c +++ b/Complex/selectfunc.c @@ -15,7 +15,7 @@ extern Core *PDL; ENTER; SAVETMPS; PUSHMARK(sp); \ SV *svpdl = sv_newmortal(); \ PDL->SetSV_PDL(svpdl, pdlvar); \ - svpdl = sv_bless(svpdl, bless_stash); /* bless in PDL::Complex */ \ + svpdl = sv_bless(svpdl, bless_stash); \ XPUSHs(svpdl); \ PUTBACK; @@ -37,13 +37,10 @@ int xerbla_(char *sub, int *info) { return 0; } { \ dSP; \ PDL_Indx odims[] = {0}; \ - PDL_Indx pc_dims[] = {2}; \ - SV *pcv = perl_get_sv("PDL::Complex::VERSION", 0); \ - char use_native = !pcv || !SvOK(pcv); \ - PDL_Indx *dims = use_native ? NULL : pc_dims; \ - PDL_Indx ndims = use_native ? 0 : sizeof(pc_dims)/sizeof(pc_dims[0]); \ - int type_add = use_native ? PDL_CF - PDL_F : 0; \ - HV *bless_stash = gv_stashpv(use_native ? "PDL" : "PDL::Complex", 0); \ + PDL_Indx *dims = NULL; \ + PDL_Indx ndims = 0; \ + int type_add = PDL_CF - PDL_F; \ + HV *bless_stash = gv_stashpv("PDL", 0); \ init \ int count = perl_call_sv(letter ## select_func, G_SCALAR); \ SPAGAIN; \ diff --git a/MANIFEST b/MANIFEST index 3d0a8ea..dde5159 100644 --- a/MANIFEST +++ b/MANIFEST @@ -17,7 +17,6 @@ t/1.t t/cgtsv.t t/common.pl t/gtsv.t -t/legacy.t Trans/Makefile.PL Trans/selectfunc.c Trans/trans.pd diff --git a/Real/real.pd b/Real/real.pd index d67516f..bd9e1eb 100644 --- a/Real/real.pd +++ b/Real/real.pd @@ -10509,11 +10509,7 @@ sub PDL::cplx_eigen { my ($eigre, $eigim, $eigvec, $fortran, @outs) = @_; my ($outval, $outvec) = @outs ? @outs : (null, null); PDL::_cplx_eigen_int($eigre, $eigim, $eigvec, $fortran, $outval, $outvec); - if (defined $PDL::Complex::VERSION) { - $_ = bless($_, 'PDL::Complex') for $outval, $outvec; - } else { - $_ = $_->slice('(0)')->czip($_->slice('(1)')) for $outval, $outvec; - } + $_ = $_->slice('(0)')->czip($_->slice('(1)')) for $outval, $outvec; ($outval, $outvec); } EOF diff --git a/Trans/selectfunc.c b/Trans/selectfunc.c index adcc0ca..cc12725 100644 --- a/Trans/selectfunc.c +++ b/Trans/selectfunc.c @@ -15,19 +15,16 @@ void dfunc_wrapper(void *p, integer n, SV* dfunc) { dSP ; PDL_Indx odims[] = {0}; - PDL_Indx pc_dims[] = {2,n}; PDL_Indx nat_dims[] = {n}; - SV *pcv = perl_get_sv("PDL::Complex::VERSION", 0); - char use_native = !pcv || !SvOK(pcv); - PDL_Indx *dims = use_native ? nat_dims : pc_dims; - PDL_Indx ndims = use_native ? 1 : 2; - int type_add = use_native ? PDL_CF - PDL_F : 0; + PDL_Indx *dims = nat_dims; + PDL_Indx ndims = 1; + int type_add = PDL_CF - PDL_F; pdl *pdl = PDL->pdlnew(); PDL->setdims(pdl, dims, ndims); pdl->datatype = PDL_D + type_add; pdl->data = p; pdl->state |= PDL_DONTTOUCHDATA | PDL_ALLOCATED; - HV *bless_stash = gv_stashpv(use_native ? "PDL" : "PDL::Complex", 0); + HV *bless_stash = gv_stashpv("PDL", 0); ENTER ; SAVETMPS ; PUSHMARK(sp) ; SV *pdl1 = sv_newmortal(); PDL->SetSV_PDL(pdl1, pdl); diff --git a/Trans/trans.pd b/Trans/trans.pd index 72edc9d..1dbcd0b 100644 --- a/Trans/trans.pd +++ b/Trans/trans.pd @@ -1035,7 +1035,7 @@ Return matrix cosine of a square matrix. =cut sub _i { - defined $PDL::Complex::VERSION ? PDL::Complex::i() : i(); + i(); } *mcos = \&PDL::mcos; @@ -1743,7 +1743,6 @@ sub PDL::mfun { my ($e, $v) = $m->meigen(0,1); my ($inv, $info) = $v->minv; unless ($info){ - $method = 'PDL::Complex::'.$method unless ref($method); eval {$v = ($v * $e->$method) x $v->minv;}; if ($@){ warn "mfun: Error $@\n"; @@ -1763,7 +1762,6 @@ sub PDL::mfun { warn "mfun: Can't compute Schur form\n"; return; } - $method = 'PDL::Complex::'.$method unless ref($method); ($t, $info) = $t->ctrfun(0,$method); if($info){ warn "mfun: Can't compute $method\n"; diff --git a/lib/PDL/LinearAlgebra.pm b/lib/PDL/LinearAlgebra.pm index c6e9c5f..563d6e5 100644 --- a/lib/PDL/LinearAlgebra.pm +++ b/lib/PDL/LinearAlgebra.pm @@ -34,14 +34,6 @@ $VERSION = eval $VERSION; my $_laerror = BARF; -{ -package # hide from CPAN indexer - PDL::Complex; - -our $floatformat = "%4.4g"; # Default print format for long numbers -our $doubleformat = "%6.6g"; -our @ISA = @ISA ? @ISA : 'PDL'; # so still operates when no PDL::Complex -} ######################################################################## =encoding utf8 @@ -141,8 +133,6 @@ Supports threading. sub PDL::dims_internal {0} sub PDL::dims_internal_values {()} -sub PDL::Complex::dims_internal {1} -sub PDL::Complex::dims_internal_values {(2)} sub PDL::_similar { my @di_vals = $_[0]->dims_internal_values; my ($m, @vdims) = @_; @@ -150,22 +140,17 @@ sub PDL::_similar { } sub PDL::_similar_null { ref($_[0])->null } sub _complex_null { - (defined &PDL::Complex::new_from_specification ? 'PDL::Complex' : 'PDL')->null + PDL->null } sub _ecplx { my ($re, $im) = @_; $re = PDL->topdl($re); - return $re if UNIVERSAL::isa($re,'PDL::Complex') or !$re->type->real; + return $re if !$re->type->real; Carp::confess("Usage: _ecplx(re,im) or (complex)") if !defined $im; $im = PDL->topdl($im); - return $re->czip($im) if !defined &PDL::Complex::new_from_specification; - my $ret = PDL::Complex->new_from_specification($re->type, 2, $re->dims); - $ret->slice('(0),') .= $re; - $ret->slice('(1),') .= $im; - return $ret; + $re->czip($im); } sub PDL::_is_complex { !$_[0]->type->real } -sub PDL::Complex::_is_complex {1} sub PDL::_norm { my ($m, $real, $trans) = @_; # If trans == true => transpose output matrix @@ -284,7 +269,6 @@ sub PDL::issym { my ($m, $tol, $conj) = @_; $tol //= ($m->type >= double) ? 1e-8 : 1e-5; $m = $m - $m->t($conj); - $m = $m->clump(2) if $m->isa('PDL::Complex'); my ($min,$max) = PDL::Ufunc::minmaximum($m); $min = $min->minimum; $max = $max->maximum; @@ -341,7 +325,7 @@ sub PDL::diag { $z = $a->slice("$slice_prefix$diag:, :-@{[$diag+1]}")->diagonal(@diag_args); } else{$z = $a->diagonal(@diag_args);} - $a->isa('PDL::Complex') ? $z->complex : $z; + $z; } use attributes 'PDL', \&PDL::diag, 'lvalue'; @@ -430,12 +414,6 @@ sub PDL::_call_method { $method = $method->[!$m->type->real ? 1 : 0]; $m->$method(@args); } -sub PDL::Complex::_call_method { - my ($m, $method, @args) = @_; - $method = [$method, "c$method"] if !ref $method; - $method = $method->[1]; - $m->$method(@args); -} *mcrossprod = \&PDL::mcrossprod; sub PDL::mcrossprod { &_2d_array; @@ -548,7 +526,6 @@ sub PDL::mdet { my $m_orig = my $m = shift->copy; $m->_call_method('getrf', my $ipiv = null, my $info = null); $m = $m->diagonal($di,$di+1); - $m = $m->complex if $m_orig->isa('PDL::Complex'); $m = $m->prodover; $m = $m * ((PDL::Ufunc::sumover(sequence($ipiv->dim(0))->plus(1,0) != $ipiv)%2)*(-2)+1); $info = which($info != 0); @@ -957,13 +934,12 @@ sub PDL::mlu { my $l = $m->mtri(1); my $slice_prefix = ',' x $di; my $smallerm1 = ($dims[$di] < $dims[$di+1] ? $dims[$di] : $dims[$di+1]) - 1; - my $one = $m->isa('PDL::Complex') ? PDL::Complex::r2C(1) : 1; if ($dims[$di+1] > $dims[$di]) { $u = $u->slice("$slice_prefix,:$smallerm1")->sever; - $l->slice("$slice_prefix :$smallerm1, :$smallerm1")->diagonal($di,$di+1) .= $one; + $l->slice("$slice_prefix :$smallerm1, :$smallerm1")->diagonal($di,$di+1) .= 1; } else { $l = $l->slice("$slice_prefix:$smallerm1")->sever if $dims[$di+1] < $dims[$di]; - $l->diagonal($di,$di+1) .= $one; + $l->diagonal($di,$di+1) .= 1; } $l, $u, $ipiv, $info; } diff --git a/pp_defc.pl b/pp_defc.pl index e00c6ce..ecd2695 100644 --- a/pp_defc.pl +++ b/pp_defc.pl @@ -6,7 +6,6 @@ sub pp_defc { my $decl = delete $hash{_decl}; $decl =~ s/\$GENERIC\(\)\s*\*/void */g; # dodge float vs float complex ptr problem $hash{Code} = "$decl\n$hash{Code}"; - pp_def("__Cc$function", %hash); my %hash2 = %hash; $hash2{Pars} = join ';', map s/\(2(?:,|(?=\)))/(/ ? "complex $_" : $_, split /;/, $hash2{Pars}; if ($hash2{RedoDimsCode}) { @@ -32,10 +31,6 @@ =head2 c$function =cut sub PDL::c$function { - barf "Cannot mix PDL::Complex and native-complex" if - (grep ref(\$_) eq 'PDL::Complex', \@_) and - (grep UNIVERSAL::isa(\$_, 'PDL') && !\$_->type->real, \@_); - goto &PDL::__Cc$function if grep ref(\$_) eq 'PDL::Complex', \@_; goto &PDL::__Nc$function; } *c$function = \\&PDL::c$function; diff --git a/t/cgtsv.t b/t/cgtsv.t index 18309e2..22429c8 100644 --- a/t/cgtsv.t +++ b/t/cgtsv.t @@ -3,48 +3,35 @@ use strict; use warnings; use PDL; -use PDL::Complex (); use PDL::NiceSlice; use PDL::LinearAlgebra::Complex; use constant N=>10; use Test::More; -my $PCi = PDL::Complex::i(); - for my $D (3..N+2) { #first differences #solve (1+i)(b_{n+1}-b_n)=1-i with homogeneous BCs - my $c=zeroes(2, $D)->complex; - my $d=-ones($D)*(1+$PCi); - my $e=ones($D)*(1+$PCi); $e->(,(-1)).=0; - my $b=ones($D)*(1-$PCi); $b->(,(-1)).=(1-$D)*(1-$PCi); + my $c=zeroes(cdouble, $D); + my $d=-ones($D)*czip(1, 1); + my $e=ones($D)*czip(1, 1); $e->((-1)).=0; + my $b=ones($D)*czip(1, -1); $b->((-1)).=(1-$D)*czip(1, -1); my $info=pdl(short,0); cgtsv($c, $d, $e, $b, $info); - my $r=sequence($D)*(1-$PCi)/(1+$PCi); - ok($b->complex->approx($r)->all, "1st diff. cgtsv in D=$D") - or diag "info: ", $info, "\nGot: ", $b, "\nExpected: ", $r; - - $c=zeroes(cdouble, $D); - $d=-ones($D)*czip(1, 1); - $e=ones($D)*czip(1, 1); $e->((-1)).=0; - $b=ones($D)*czip(1, -1); $b->((-1)).=(1-$D)*czip(1, -1); - $info=pdl(short,0); - cgtsv($c, $d, $e, $b, $info); - $r=sequence($D)*czip(1, -1)/czip(1, 1); + my $r=sequence($D)*czip(1, -1)/czip(1, 1); ok($b->approx($r)->all, "1st diff. native cgtsv in D=$D") or diag "info: ", $info, "\nGot: ", $b, "\nExpected: ", $r; } for my $D (3..N+2){ #second differences #solve b_{n+1}-2{b_n}+b_{n-1}=1 with kinda homogeneous BCs - my $c=ones($D)*(1+$PCi); $c->(,(-1)).=0; - my $d=-2*ones($D)*(1+$PCi); - my $e=ones($D)*(1+$PCi); $e->(,(-1)).=0; - my $b=ones($D)*(1-$PCi); + my $c=ones($D)*(1+i()); $c->(,(-1)).=0; + my $d=-2*ones($D)*(1+i()); + my $e=ones($D)*(1+i()); $e->(,(-1)).=0; + my $b=ones($D)*(1-i()); my $info=pdl(short,0); cgtsv($c, $d, $e, $b, $info); - my $x=sequence($D)+0*$PCi; - my $r=(-$D/2-($D-1)/2*$x+1/2*$x*$x)*(1-$PCi)/(1+$PCi); - ok($b->complex->approx($r)->all, "2nd diff. cgtsv in D=$D") + my $x=sequence(cdouble, $D); + my $r=(-$D/2-($D-1)/2*$x+1/2*$x*$x)*(1-i())/(1+i()); + ok($b->approx($r)->all, "2nd diff. cgtsv in D=$D") or diag "info: ", $info, "\nGot: ", $b, "\nExpected: ", $r; } diff --git a/t/legacy.t b/t/legacy.t deleted file mode 100644 index e5a3e42..0000000 --- a/t/legacy.t +++ /dev/null @@ -1,46 +0,0 @@ -use strict; -use warnings; -use PDL::LiteF; -use PDL::MatrixOps qw(identity); -use PDL::Complex; -use PDL::LinearAlgebra; -use PDL::LinearAlgebra::Trans qw //; -use PDL::LinearAlgebra::Complex; -use Test::More; - -sub fapprox { - my($a,$b) = @_; - ($a-$b)->abs->max < 0.001; -} -# PDL::Complex only -sub runtest { - local $Test::Builder::Level = $Test::Builder::Level + 1; - my ($in, $method, $expected_cplx, $extra) = @_; - $expected_cplx = $expected_cplx->[1] if ref $expected_cplx eq 'ARRAY'; - my @cplx = ref($expected_cplx) eq 'ARRAY' ? @$expected_cplx : $expected_cplx; - $_ = PDL::Complex::r2C($_) for $in; - my ($got) = $in->$method(map ref() && ref() ne 'CODE' ? PDL::Complex::r2C($_) : $_, @{$extra||[]}); - my $ok = grep fapprox($got, PDL::Complex::r2C($_)), @cplx; - ok $ok, "PDL::Complex $method" or diag "got(".ref($got)."):$got\nexpected:@cplx"; -} - -my $aa = cplx random(2,2,2); -runtest($aa, 't', $aa->xchg(1,2)); - -$aa = sequence(2,2,2)->cplx + 1; -runtest($aa, '_norm', my $aa_exp = PDL::Complex->from_native(pdl <<'EOF')); -[ - [0.223606+0.223606i 0.670820+0.670820i] - [0.410997+0.410997i 0.575396+0.575396i] -] -EOF -runtest($aa, '_norm', $aa_exp->abs, [1]); -runtest($aa, '_norm', $aa_exp->t, [0,1]); - -$aa = pdl('[[[0 1] [2 3] [4 5]] [[6 7] [8 9] [10 11]] [[12 13] [14 15] [16 17]]] ')->cplx; -my $up = pdl('[[[0 1] [2 3] [4 5]] [[0 0] [8 9] [10 11]] [[0 0] [0 0] [16 17]]]')->cplx; -my $lo = pdl('[[[0 1] [0 0] [0 0]] [[6 7] [8 9] [0 0]] [[12 13] [14 15] [16 17]]]')->cplx; - -do './t/common.pl'; die if $@; - -done_testing;