Mercurial > hg > octave-nkf

comparison liboctave/dMatrix.cc @ 4153:6b96ce9f5743

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2002-11-06 20:38:49 by jwe]
author jwe
date Wed, 06 Nov 2002 20:38:50 +0000
parents 7d9bda865012
children 5719210fff4c
comparison
equal deleted inserted replaced
4152:f14251d33b01 4153:6b96ce9f5743
47 #include "mx-base.h" 47 #include "mx-base.h"
48 #include "mx-m-dm.h" 48 #include "mx-m-dm.h"
49 #include "mx-dm-m.h" 49 #include "mx-dm-m.h"
50 #include "mx-inlines.cc" 50 #include "mx-inlines.cc"
51 #include "oct-cmplx.h" 51 #include "oct-cmplx.h"
52 #include "quit.h"
52 53
53 #ifdef HAVE_FFTW 54 #ifdef HAVE_FFTW
54 #include "oct-fftw.h" 55 #include "oct-fftw.h"
55 #endif 56 #endif
56 57
679 Complex *in (tmp.fortran_vec ()); 680 Complex *in (tmp.fortran_vec ());
680 Complex *out (retval.fortran_vec ()); 681 Complex *out (retval.fortran_vec ());
681 682
682 for (size_t i = 0; i < nsamples; i++) 683 for (size_t i = 0; i < nsamples; i++)
683 { 684 {
685 OCTAVE_QUIT;
686
684 octave_fftw::fft (&in[npts * i], &out[npts * i], npts); 687 octave_fftw::fft (&in[npts * i], &out[npts * i], npts);
685 } 688 }
686 689
687 return retval; 690 return retval;
688 } 691 }
712 Complex *in (tmp.fortran_vec ()); 715 Complex *in (tmp.fortran_vec ());
713 Complex *out (retval.fortran_vec ()); 716 Complex *out (retval.fortran_vec ());
714 717
715 for (size_t i = 0; i < nsamples; i++) 718 for (size_t i = 0; i < nsamples; i++)
716 { 719 {
720 OCTAVE_QUIT;
721
717 octave_fftw::ifft (&in[npts * i], &out[npts * i], npts); 722 octave_fftw::ifft (&in[npts * i], &out[npts * i], npts);
718 } 723 }
719 724
720 return retval; 725 return retval;
721 } 726 }
780 Complex *tmp_data = retval.fortran_vec (); 785 Complex *tmp_data = retval.fortran_vec ();
781 786
782 F77_FUNC (cffti, CFFTI) (npts, pwsave); 787 F77_FUNC (cffti, CFFTI) (npts, pwsave);
783 788
784 for (int j = 0; j < nsamples; j++) 789 for (int j = 0; j < nsamples; j++)
785 F77_FUNC (cfftf, CFFTF) (npts, &tmp_data[npts*j], pwsave); 790 {
791 OCTAVE_QUIT;
792
793 F77_FUNC (cfftf, CFFTF) (npts, &tmp_data[npts*j], pwsave);
794 }
786 795
787 return retval; 796 return retval;
788 } 797 }
789 798
790 ComplexMatrix 799 ComplexMatrix
817 Complex *tmp_data = retval.fortran_vec (); 826 Complex *tmp_data = retval.fortran_vec ();
818 827
819 F77_FUNC (cffti, CFFTI) (npts, pwsave); 828 F77_FUNC (cffti, CFFTI) (npts, pwsave);
820 829
821 for (int j = 0; j < nsamples; j++) 830 for (int j = 0; j < nsamples; j++)
822 F77_FUNC (cfftb, CFFTB) (npts, &tmp_data[npts*j], pwsave); 831 {
832 OCTAVE_QUIT;
833
834 F77_FUNC (cfftb, CFFTB) (npts, &tmp_data[npts*j], pwsave);
835 }
823 836
824 for (int j = 0; j < npts*nsamples; j++) 837 for (int j = 0; j < npts*nsamples; j++)
825 tmp_data[j] = tmp_data[j] / static_cast<double> (npts); 838 tmp_data[j] = tmp_data[j] / static_cast<double> (npts);
826 839
827 return retval; 840 return retval;
857 Complex *tmp_data = retval.fortran_vec (); 870 Complex *tmp_data = retval.fortran_vec ();
858 871
859 F77_FUNC (cffti, CFFTI) (npts, pwsave); 872 F77_FUNC (cffti, CFFTI) (npts, pwsave);
860 873
861 for (int j = 0; j < nsamples; j++) 874 for (int j = 0; j < nsamples; j++)
862 F77_FUNC (cfftf, CFFTF) (npts, &tmp_data[npts*j], pwsave); 875 {
876 OCTAVE_QUIT;
877
878 F77_FUNC (cfftf, CFFTF) (npts, &tmp_data[npts*j], pwsave);
879 }
863 880
864 npts = nc; 881 npts = nc;
865 nsamples = nr; 882 nsamples = nr;
866 nn = 4*npts+15; 883 nn = 4*npts+15;
867 884
873 890
874 F77_FUNC (cffti, CFFTI) (npts, pwsave); 891 F77_FUNC (cffti, CFFTI) (npts, pwsave);
875 892
876 for (int j = 0; j < nsamples; j++) 893 for (int j = 0; j < nsamples; j++)
877 { 894 {
895 OCTAVE_QUIT;
896
878 for (int i = 0; i < npts; i++) 897 for (int i = 0; i < npts; i++)
879 prow[i] = tmp_data[i*nr + j]; 898 prow[i] = tmp_data[i*nr + j];
880 899
881 F77_FUNC (cfftf, CFFTF) (npts, prow, pwsave); 900 F77_FUNC (cfftf, CFFTF) (npts, prow, pwsave);
882 901
917 Complex *tmp_data = retval.fortran_vec (); 936 Complex *tmp_data = retval.fortran_vec ();
918 937
919 F77_FUNC (cffti, CFFTI) (npts, pwsave); 938 F77_FUNC (cffti, CFFTI) (npts, pwsave);
920 939
921 for (int j = 0; j < nsamples; j++) 940 for (int j = 0; j < nsamples; j++)
922 F77_FUNC (cfftb, CFFTB) (npts, &tmp_data[npts*j], pwsave); 941 {
942 OCTAVE_QUIT;
943
944 F77_FUNC (cfftb, CFFTB) (npts, &tmp_data[npts*j], pwsave);
945 }
923 946
924 for (int j = 0; j < npts*nsamples; j++) 947 for (int j = 0; j < npts*nsamples; j++)
925 tmp_data[j] = tmp_data[j] / static_cast<double> (npts); 948 tmp_data[j] = tmp_data[j] / static_cast<double> (npts);
926 949
927 npts = nc; 950 npts = nc;
936 959
937 F77_FUNC (cffti, CFFTI) (npts, pwsave); 960 F77_FUNC (cffti, CFFTI) (npts, pwsave);
938 961
939 for (int j = 0; j < nsamples; j++) 962 for (int j = 0; j < nsamples; j++)
940 { 963 {
964 OCTAVE_QUIT;
965
941 for (int i = 0; i < npts; i++) 966 for (int i = 0; i < npts; i++)
942 prow[i] = tmp_data[i*nr + j]; 967 prow[i] = tmp_data[i*nr + j];
943 968
944 F77_FUNC (cfftb, CFFTB) (npts, prow, pwsave); 969 F77_FUNC (cfftb, CFFTB) (npts, prow, pwsave);
945 970
1608 // apply inverse scaling to computed exponential 1633 // apply inverse scaling to computed exponential
1609 for (int i = 0; i < nc; i++) 1634 for (int i = 0; i < nc; i++)
1610 for (int j = 0; j < nc; j++) 1635 for (int j = 0; j < nc; j++)
1611 retval(i,j) *= dscale(i) / dscale(j); 1636 retval(i,j) *= dscale(i) / dscale(j);
1612 1637
1638 OCTAVE_QUIT;
1639
1613 // construct balancing permutation vector 1640 // construct balancing permutation vector
1614 Array<int> ipermute (nc); 1641 Array<int> ipermute (nc);
1615 for (int i = 0; i < nc; i++) 1642 for (int i = 0; i < nc; i++)
1616 ipermute(i) = i; // identity permutation 1643 ipermute(i) = i; // identity permutation
1617 1644
1635 1662
1636 // construct inverse balancing permutation vector 1663 // construct inverse balancing permutation vector
1637 Array<int> invpvec (nc); 1664 Array<int> invpvec (nc);
1638 for (int i = 0; i < nc; i++) 1665 for (int i = 0; i < nc; i++)
1639 invpvec(ipermute(i)) = i; // Thanks to R. A. Lippert for this method 1666 invpvec(ipermute(i)) = i; // Thanks to R. A. Lippert for this method
1667
1668 OCTAVE_QUIT;
1640 1669
1641 Matrix tmpMat = retval; 1670 Matrix tmpMat = retval;
1642 for (int i = 0; i < nc; i++) 1671 for (int i = 0; i < nc; i++)
1643 for (int j = 0; j < nc; j++) 1672 for (int j = 0; j < nc; j++)
1644 retval(i,j) = tmpMat(invpvec(i),invpvec(j)); 1673 retval(i,j) = tmpMat(invpvec(i),invpvec(j));