Diff of /Trans/trans.pd [e99cbd] .. [591b21] Maximize Restore

  Switch to unified view

a/Trans/trans.pd b/Trans/trans.pd
1
our $VERSION = '0.08_01';
1
our $VERSION = 0.04;
2
pp_setversion(qq{'$VERSION'});
3
$VERSION = eval $VERSION;
4
5
use PDL::Exporter;
2
use PDL::Exporter;
6
3
7
if ($^O =~ /MSWin/) {
4
if ($^O =~ /MSWin/) {
8
pp_addhdr('
5
pp_addhdr('
9
#include <float.h>
6
#include <float.h>
10
');
7
');
11
}
8
}
12
9
13
pp_addhdr('
10
pp_addhdr('
14
#include <math.h>
11
#include <math.h>
15
16
#if defined(PDL_CORE_VERSION) && PDL_CORE_VERSION < 10
17
typedef PDL_Long PDL_Indx;
18
#endif
19
20
typedef PDL_Long logical;
12
typedef PDL_Long logical;
21
typedef PDL_Long integer;
13
typedef PDL_Long integer;
22
typedef PDL_Long ftnlen;
14
typedef PDL_Long ftnlen;
23
typedef struct { double r, i; } dcomplex;
15
typedef struct { double r, i; } dcomplex;
24
16
...
...
90
}
82
}
91
83
92
');
84
');
93
85
94
86
87
pp_setversion(0.04);
95
88
96
pp_addpm({At=>'Top'},<<'EOD');
89
pp_addpm({At=>'Top'},<<'EOD');
97
use PDL::Func;
90
use PDL::Func;
98
use PDL::Core;
91
use PDL::Core;
99
use PDL::Slices;
92
use PDL::Slices;
...
...
103
use PDL::NiceSlice;
96
use PDL::NiceSlice;
104
use PDL::LinearAlgebra;
97
use PDL::LinearAlgebra;
105
use PDL::LinearAlgebra::Real qw //;
98
use PDL::LinearAlgebra::Real qw //;
106
use PDL::LinearAlgebra::Complex qw //;
99
use PDL::LinearAlgebra::Complex qw //;
107
use strict;
100
use strict;
108
109
=encoding Latin-1
110
101
111
=head1 NAME
102
=head1 NAME
112
103
113
PDL::LinearAlgebra::Trans - Linear Algebra based transcendental functions for PDL
104
PDL::LinearAlgebra::Trans - Linear Algebra based transcendental functions for PDL
114
105
...
...
257
248
258
    // ---  scaling: seek ns such that ||t*H/2^ns|| < 1/2; 
249
    // ---  scaling: seek ns such that ||t*H/2^ns|| < 1/2; 
259
    //     and set scale = t/2^ns ... 
250
    //     and set scale = t/2^ns ... 
260
251
261
    // Compute Infinity norm
252
    // Compute Infinity norm
262
    hnorm = dlange_("I", &(integer){$PRIV(__n_size)}, &(integer){$PRIV(__n_size)}, $P(A), 
253
    hnorm = dlange_("I", &$PRIV(__n_size), &$PRIV(__n_size), $P(A), 
263
            &(integer){$PRIV(__n_size)}, wsp);
254
            &$PRIV(__n_size), wsp);
264
    
255
    
265
    hnorm = abs($scale() * hnorm);
256
    hnorm = abs($scale() * hnorm);
266
    if (hnorm == 0.)
257
    if (hnorm == 0.)
267
    {
258
    {
268
        free(ipiv);
259
        free(ipiv);
...
...
294
            coef[k] = coef[k - 1] * (double) (i__ - k) / (double) (k * (j - k));
285
            coef[k] = coef[k - 1] * (double) (i__ - k) / (double) (k * (j - k));
295
        }
286
        }
296
        
287
        
297
        // ---  H2 = scale2*H*H ... 
288
        // ---  H2 = scale2*H*H ... 
298
        
289
        
299
        dgemm_("n", "n", &(integer){$PRIV(__n_size)}, &(integer){$PRIV(__n_size)}, &(integer){$PRIV(__n_size)}, &scale2, h__, &(integer){$PRIV(__n_size)}, h__, 
290
        dgemm_("n", "n", &$PRIV(__n_size), &$PRIV(__n_size), &$PRIV(__n_size), &scale2, h__, &$PRIV(__n_size), h__, 
300
            &(integer){$PRIV(__n_size)}, &c_b7, &coef[ih2], &(integer){$PRIV(__n_size)}, (ftnlen)1, (ftnlen)1);
291
            &$PRIV(__n_size), &c_b7, &coef[ih2], &$PRIV(__n_size), (ftnlen)1, (ftnlen)1);
301
        
292
        
302
        // ---  initialize p (numerator) and q (denominator) ...
293
        // ---  initialize p (numerator) and q (denominator) ...
303
        cp = coef[$deg() - 1];
294
        cp = coef[$deg() - 1];
304
        cq = coef[$deg()];
295
        cq = coef[$deg()];
305
        for (j = 0; j < $PRIV(__n_size); j++) {
296
        for (j = 0; j < $PRIV(__n_size); j++) {
...
...
314
        // ---  Apply Horner rule ... 
305
        // ---  Apply Horner rule ... 
315
        iodd = 1;
306
        iodd = 1;
316
        k = $deg() - 2;
307
        k = $deg() - 2;
317
        do{
308
        do{
318
            iused = iodd * iq + (1 - iodd) * ip;
309
            iused = iodd * iq + (1 - iodd) * ip;
319
            dgemm_("n", "n", &(integer){$PRIV(__n_size)}, &(integer){$PRIV(__n_size)}, &(integer){$PRIV(__n_size)}, &c_b11, &wsp[iused], &(integer){$PRIV(__n_size)}, &coef[ih2], &(integer){$PRIV(__n_size)}, &c_b7, &
310
            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, &
320
                wsp[ifree], &(integer){$PRIV(__n_size)}, (ftnlen)1, (ftnlen)1);
311
                wsp[ifree], &$PRIV(__n_size), (ftnlen)1, (ftnlen)1);
321
            for (j = 0; j < $PRIV(__n_size); j++) {
312
            for (j = 0; j < $PRIV(__n_size); j++) {
322
                wsp[ifree + j  * ($PRIV(__n_size) + 1)] += coef[k];
313
                wsp[ifree + j  * ($PRIV(__n_size) + 1)] += coef[k];
323
            }
314
            }
324
            ip = (1 - iodd) * ifree + iodd * ip;
315
            ip = (1 - iodd) * ifree + iodd * ip;
325
            iq = iodd * ifree + (1 - iodd) * iq;
316
            iq = iodd * ifree + (1 - iodd) * iq;
...
...
332
    
323
    
333
        // ---  Obtain (+/-)(I + 2*(p\q)) ... 
324
        // ---  Obtain (+/-)(I + 2*(p\q)) ... 
334
        
325
        
335
        if (iodd == 1)
326
        if (iodd == 1)
336
        {
327
        {
337
            dgemm_("n", "n", &(integer){$PRIV(__n_size)}, &(integer){$PRIV(__n_size)}, &(integer){$PRIV(__n_size)}, &scale, &wsp[iq], &(integer){$PRIV(__n_size)}, h__, &(integer){$PRIV(__n_size)}, &
328
            dgemm_("n", "n", &$PRIV(__n_size), &$PRIV(__n_size), &$PRIV(__n_size), &scale, &wsp[iq], &$PRIV(__n_size), h__, &$PRIV(__n_size), &
338
                c_b7, &wsp[ifree], &(integer){$PRIV(__n_size)}, (ftnlen)1, (ftnlen)1);
329
                c_b7, &wsp[ifree], &$PRIV(__n_size), (ftnlen)1, (ftnlen)1);
339
            iq = ifree;
330
            iq = ifree;
340
        }
331
        }
341
        else
332
        else
342
        {
333
        {
343
            dgemm_("n", "n", &(integer){$PRIV(__n_size)}, &(integer){$PRIV(__n_size)}, &(integer){$PRIV(__n_size)}, &scale, &wsp[ip], &(integer){$PRIV(__n_size)}, h__, &(integer){$PRIV(__n_size)}, &
334
            dgemm_("n", "n", &$PRIV(__n_size), &$PRIV(__n_size), &$PRIV(__n_size), &scale, &wsp[ip], &$PRIV(__n_size), h__, &$PRIV(__n_size), &
344
                c_b7, &wsp[ifree], &(integer){$PRIV(__n_size)}, (ftnlen)1, (ftnlen)1);
335
                c_b7, &wsp[ifree], &$PRIV(__n_size), (ftnlen)1, (ftnlen)1);
345
            ip = ifree;
336
            ip = ifree;
346
        }
337
        }
347
        daxpy_(&mm, &c_b19, &wsp[ip], &c__1, &wsp[iq], &c__1);
338
        daxpy_(&mm, &c_b19, &wsp[ip], &c__1, &wsp[iq], &c__1);
348
        dgesv_(&(integer){$PRIV(__n_size)}, &(integer){$PRIV(__n_size)}, &wsp[iq], &(integer){$PRIV(__n_size)}, ipiv, &wsp[ip], &(integer){$PRIV(__n_size)}, $P(info));
339
        dgesv_(&$PRIV(__n_size), &$PRIV(__n_size), &wsp[iq], &$PRIV(__n_size), ipiv, &wsp[ip], &$PRIV(__n_size), $P(info));
349
    
340
    
350
        if ($info() != 0)
341
        if ($info() != 0)
351
        {
342
        {
352
            free(ipiv);
343
            free(ipiv);
353
            free(coef);
344
            free(coef);
...
...
367
                iodd = 1;
358
                iodd = 1;
368
                i__1 = $ns();
359
                i__1 = $ns();
369
                for (k = 0; k < i__1; k++) {
360
                for (k = 0; k < i__1; k++) {
370
                    iget = iodd * ip + (1 - iodd) * iq;
361
                    iget = iodd * ip + (1 - iodd) * iq;
371
                    iput = (1 - iodd) * ip + iodd * iq;
362
                    iput = (1 - iodd) * ip + iodd * iq;
372
                    dgemm_("n", "n", &(integer){$PRIV(__n_size)}, &(integer){$PRIV(__n_size)}, &(integer){$PRIV(__n_size)}, &c_b11, &wsp[iget], &(integer){$PRIV(__n_size)}, &wsp[iget], &(integer){$PRIV(__n_size)}, &c_b7,
363
                    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,
373
                        &wsp[iput], &(integer){$PRIV(__n_size)}, (ftnlen)1, (ftnlen)1);
364
                        &wsp[iput], &$PRIV(__n_size), (ftnlen)1, (ftnlen)1);
374
                    iodd = 1 - iodd;
365
                    iodd = 1 - iodd;
375
                }
366
                }
376
            }
367
            }
377
            
368
            
378
            free(coef);
369
            free(coef);
...
...
575
    }
566
    }
576
567
577
    // ---  scaling: seek ns such that ||t*H/2^ns|| < 1/2; 
568
    // ---  scaling: seek ns such that ||t*H/2^ns|| < 1/2; 
578
    //     and set scale = t/2^ns ... 
569
    //     and set scale = t/2^ns ... 
579
    // Compute Infinity norm
570
    // Compute Infinity norm
580
    hnorm = zlange_("I", &(integer){$PRIV(__n_size)}, &(integer){$PRIV(__n_size)}, &h__[0], 
571
    hnorm = zlange_("I", &$PRIV(__n_size), &$PRIV(__n_size), &h__[0], 
581
            &(integer){$PRIV(__n_size)}, &wsp[0].r);
572
            &$PRIV(__n_size), &wsp[0].r);
582
573
583
        hnorm = abs($scale() * hnorm);
574
        hnorm = abs($scale() * hnorm);
584
        if (hnorm == 0.)
575
        if (hnorm == 0.)
585
        {
576
        {
586
        $ns() = 0;
577
        $ns() = 0;
...
...
614
            wsp[i__2].r = (d__1 * wsp[i__3].r) / d__2;
605
            wsp[i__2].r = (d__1 * wsp[i__3].r) / d__2;
615
            wsp[i__2].i = (d__1 * wsp[i__3].i) / d__2;
606
            wsp[i__2].i = (d__1 * wsp[i__3].i) / d__2;
616
            }
607
            }
617
        
608
        
618
        // ---  H2 = scale2*H*H ...
609
        // ---  H2 = scale2*H*H ...
619
            zgemm_("n", "n", &(integer){$PRIV(__n_size)}, &(integer){$PRIV(__n_size)}, &(integer){$PRIV(__n_size)}, &scale2, &h__[0], &(integer){$PRIV(__n_size)}, &h__[0], 
610
            zgemm_("n", "n", &$PRIV(__n_size), &$PRIV(__n_size), &$PRIV(__n_size), &scale2, &h__[0], &$PRIV(__n_size), &h__[0], 
620
                &(integer){$PRIV(__n_size)}, &c_b1, &wsp[ih2], &(integer){$PRIV(__n_size)}, (ftnlen)1, (ftnlen)1);
611
                &$PRIV(__n_size), &c_b1, &wsp[ih2], &$PRIV(__n_size), (ftnlen)1, (ftnlen)1);
621
        
612
        
622
        // ---  initialise p (numerator) and q (denominator) ...
613
        // ---  initialise p (numerator) and q (denominator) ...
623
        
614
        
624
            i__1 = $deg() - 1;
615
            i__1 = $deg() - 1;
625
            cp.r = wsp[i__1].r;
616
            cp.r = wsp[i__1].r;
...
...
647
            iodd = 1;
638
            iodd = 1;
648
            k = $deg() - 2;
639
            k = $deg() - 2;
649
            do
640
            do
650
            {
641
            {
651
                iused = iodd * iq + (1 - iodd) * ip;
642
                iused = iodd * iq + (1 - iodd) * ip;
652
                zgemm_("n", "n", &(integer){$PRIV(__n_size)}, &(integer){$PRIV(__n_size)}, &(integer){$PRIV(__n_size)}, &c_b2, &wsp[iused], &(integer){$PRIV(__n_size)}, &wsp[ih2], &(integer){$PRIV(__n_size)}, &c_b1, &
643
                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, &
653
                    wsp[ifree], &(integer){$PRIV(__n_size)}, (ftnlen)1, (ftnlen)1);
644
                    wsp[ifree], &$PRIV(__n_size), (ftnlen)1, (ftnlen)1);
654
                for (j = 0; j < $PRIV(__n_size); j++) {
645
                for (j = 0; j < $PRIV(__n_size); j++) {
655
                i__1 = ifree + j  * ($PRIV(__n_size) + 1);
646
                i__1 = ifree + j  * ($PRIV(__n_size) + 1);
656
                i__2 = ifree + j  * ($PRIV(__n_size) + 1);
647
                i__2 = ifree + j  * ($PRIV(__n_size) + 1);
657
                i__3 = k;
648
                i__3 = k;
658
                
649
                
...
...
671
662
672
        // ---  Obtain (+/-)(I + 2*(p\q)) ...
663
        // ---  Obtain (+/-)(I + 2*(p\q)) ...
673
        
664
        
674
            if (iodd != 0)
665
            if (iodd != 0)
675
            {
666
            {
676
            zgemm_("n", "n", &(integer){$PRIV(__n_size)}, &(integer){$PRIV(__n_size)}, &(integer){$PRIV(__n_size)}, &scale, &wsp[iq], &(integer){$PRIV(__n_size)}, &h__[0], &(integer){$PRIV(__n_size)}, &
667
            zgemm_("n", "n", &$PRIV(__n_size), &$PRIV(__n_size), &$PRIV(__n_size), &scale, &wsp[iq], &$PRIV(__n_size), &h__[0], &$PRIV(__n_size), &
677
                c_b1, &wsp[ifree], &(integer){$PRIV(__n_size)}, (ftnlen)1, (ftnlen)1);
668
                c_b1, &wsp[ifree], &$PRIV(__n_size), (ftnlen)1, (ftnlen)1);
678
            iq = ifree;
669
            iq = ifree;
679
            }else
670
            }else
680
            {
671
            {
681
            zgemm_("n", "n", &(integer){$PRIV(__n_size)}, &(integer){$PRIV(__n_size)}, &(integer){$PRIV(__n_size)}, &scale, &wsp[ip], &(integer){$PRIV(__n_size)}, &h__[0], &(integer){$PRIV(__n_size)}, &
672
            zgemm_("n", "n", &$PRIV(__n_size), &$PRIV(__n_size), &$PRIV(__n_size), &scale, &wsp[ip], &$PRIV(__n_size), &h__[0], &$PRIV(__n_size), &
682
                c_b1, &wsp[ifree], &(integer){$PRIV(__n_size)}, (ftnlen)1, (ftnlen)1);
673
                c_b1, &wsp[ifree], &$PRIV(__n_size), (ftnlen)1, (ftnlen)1);
683
            ip = ifree;
674
            ip = ifree;
684
            }
675
            }
685
            z__1.r = -1.;
676
            z__1.r = -1.;
686
            z__1.i = -0.;
677
            z__1.i = -0.;
687
            zaxpy_(&mm, &z__1, &wsp[ip], &c__1, &wsp[iq], &c__1);
678
            zaxpy_(&mm, &z__1, &wsp[ip], &c__1, &wsp[iq], &c__1);
688
            zgesv_(&(integer){$PRIV(__n_size)}, &(integer){$PRIV(__n_size)}, &wsp[iq], &(integer){$PRIV(__n_size)}, ipiv, &wsp[ip], &(integer){$PRIV(__n_size)}, $P(info));
679
            zgesv_(&$PRIV(__n_size), &$PRIV(__n_size), &wsp[iq], &$PRIV(__n_size), ipiv, &wsp[ip], &$PRIV(__n_size), $P(info));
689
680
690
            if ($info() == 0) 
681
            if ($info() == 0) 
691
            {
682
            {
692
                zdscal_(&mm, &c_b19, &wsp[ip], &c__1);
683
                zdscal_(&mm, &c_b19, &wsp[ip], &c__1);
693
                for (j = 0; j < $PRIV(__n_size); j++)
684
                for (j = 0; j < $PRIV(__n_size); j++)
...
...
708
                    iodd = 1;
699
                    iodd = 1;
709
                    for (k = 0; k < $ns(); k++)
700
                    for (k = 0; k < $ns(); k++)
710
                    {
701
                    {
711
                    iget = iodd * ip + (1 - iodd) * iq;
702
                    iget = iodd * ip + (1 - iodd) * iq;
712
                    iput = (1 - iodd) * ip + iodd * iq;
703
                    iput = (1 - iodd) * ip + iodd * iq;
713
                    zgemm_("n", "n", &(integer){$PRIV(__n_size)}, &(integer){$PRIV(__n_size)}, &(integer){$PRIV(__n_size)}, &c_b2, &wsp[iget], &(integer){$PRIV(__n_size)}, &wsp[iget], &(integer){$PRIV(__n_size)}, &c_b1, 
704
                    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, 
714
                        &wsp[iput], &(integer){$PRIV(__n_size)}, (ftnlen)1, (ftnlen)1);
705
                        &wsp[iput], &$PRIV(__n_size), (ftnlen)1, (ftnlen)1);
715
                    iodd = 1 - iodd;
706
                    iodd = 1 - iodd;
716
                    }
707
                    }
717
                }
708
                }
718
            }
709
            }
719
710
...
...
849
   int count;
840
   int count;
850
   SV *pdl1;
841
   SV *pdl1;
851
   HV *bless_stash;
842
   HV *bless_stash;
852
843
853
   pdl *pdl;
844
   pdl *pdl;
854
   PDL_Indx odims[1];
845
   PDL_Long odims[1];
855
846
856
   PDL_Indx dims[] = {2,n};
847
   PDL_Long dims[] = {2,n};
857
   pdl = PDL->pdlnew();
848
   pdl = PDL->pdlnew();
858
   PDL->setdims (pdl, dims, 2);
849
   PDL->setdims (pdl, dims, 2);
859
   pdl->datatype = PDL_D;
850
   pdl->datatype = PDL_D;
860
   pdl->data = (double *) &p[0].r;
851
   pdl->data = (double *) &p[0].r;
861
   pdl->state |= PDL_DONTTOUCHDATA | PDL_ALLOCATED;
852
   pdl->state |= PDL_DONTTOUCHDATA | PDL_ALLOCATED;