Mercurial > hg > octave-nkf

comparison liboctave/dMatrix.cc @ 4773:ccfbd6047a54

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2004-02-16 19:02:32 by jwe]
author jwe
date Mon, 16 Feb 2004 19:02:33 +0000
parents 334a27c8f453
children 44046bbaa52c
comparison
equal deleted inserted replaced
4772:9eed17b2c8d1 4773:ccfbd6047a54
50 #include "mx-dm-m.h" 50 #include "mx-dm-m.h"
51 #include "mx-inlines.cc" 51 #include "mx-inlines.cc"
52 #include "oct-cmplx.h" 52 #include "oct-cmplx.h"
53 #include "quit.h" 53 #include "quit.h"
54 54
55 #ifdef HAVE_FFTW 55 #if defined (HAVE_FFTW3)
56 #include "oct-fftw.h" 56 #include "oct-fftw.h"
57 #endif 57 #endif
58 58
59 // Fortran functions we call. 59 // Fortran functions we call.
60 60
764 Matrix Vr = V.extract (0, 0, nc-1, r); 764 Matrix Vr = V.extract (0, 0, nc-1, r);
765 return Vr * D * Ur.transpose (); 765 return Vr * D * Ur.transpose ();
766 } 766 }
767 } 767 }
768 768
769 #ifdef HAVE_FFTW 769 #if defined (HAVE_FFTW3)
770 770
771 ComplexMatrix 771 ComplexMatrix
772 Matrix::fourier (void) const 772 Matrix::fourier (void) const
773 {
774 size_t nr = rows ();
775 size_t nc = cols ();
776
777 ComplexMatrix retval (nr, nc);
778
779 size_t npts, nsamples;
780
781 if (nr == 1 || nc == 1)
782 {
783 npts = nr > nc ? nr : nc;
784 nsamples = 1;
785 }
786 else
787 {
788 npts = nr;
789 nsamples = nc;
790 }
791
792 const double *in (fortran_vec ());
793 Complex *out (retval.fortran_vec ());
794
795 octave_fftw::fft (in, out, npts, nsamples);
796
797 return retval;
798 }
799
800 ComplexMatrix
801 Matrix::ifourier (void) const
773 { 802 {
774 size_t nr = rows (); 803 size_t nr = rows ();
775 size_t nc = cols (); 804 size_t nc = cols ();
776 805
777 ComplexMatrix retval (nr, nc); 806 ComplexMatrix retval (nr, nc);
791 820
792 ComplexMatrix tmp (*this); 821 ComplexMatrix tmp (*this);
793 Complex *in (tmp.fortran_vec ()); 822 Complex *in (tmp.fortran_vec ());
794 Complex *out (retval.fortran_vec ()); 823 Complex *out (retval.fortran_vec ());
795 824
796 for (size_t i = 0; i < nsamples; i++) 825 octave_fftw::ifft (in, out, npts, nsamples);
797 {
798 OCTAVE_QUIT;
799
800 octave_fftw::fft (&in[npts * i], &out[npts * i], npts);
801 }
802
803 return retval;
804 }
805
806 ComplexMatrix
807 Matrix::ifourier (void) const
808 {
809 size_t nr = rows ();
810 size_t nc = cols ();
811
812 ComplexMatrix retval (nr, nc);
813
814 size_t npts, nsamples;
815
816 if (nr == 1 || nc == 1)
817 {
818 npts = nr > nc ? nr : nc;
819 nsamples = 1;
820 }
821 else
822 {
823 npts = nr;
824 nsamples = nc;
825 }
826
827 ComplexMatrix tmp (*this);
828 Complex *in (tmp.fortran_vec ());
829 Complex *out (retval.fortran_vec ());
830
831 for (size_t i = 0; i < nsamples; i++)
832 {
833 OCTAVE_QUIT;
834
835 octave_fftw::ifft (&in[npts * i], &out[npts * i], npts);
836 }
837 826
838 return retval; 827 return retval;
839 } 828 }
840 829
841 ComplexMatrix 830 ComplexMatrix
842 Matrix::fourier2d (void) const 831 Matrix::fourier2d (void) const
843 { 832 {
844 int nr = rows (); 833 dim_vector dv(rows (), cols ());
845 int nc = cols (); 834
846 835 const double *in = fortran_vec ();
847 ComplexMatrix retval (*this); 836 ComplexMatrix retval (rows (), cols ());
848 // Note the order of passing the rows and columns to account for 837 octave_fftw::fftNd (in, retval.fortran_vec (), 2, dv);
849 // column-major storage.
850 octave_fftw::fft2d (retval.fortran_vec (), nc, nr);
851 838
852 return retval; 839 return retval;
853 } 840 }
854 841
855 ComplexMatrix 842 ComplexMatrix
856 Matrix::ifourier2d (void) const 843 Matrix::ifourier2d (void) const
857 { 844 {
858 int nr = rows (); 845 dim_vector dv(rows (), cols ());
859 int nc = cols ();
860 846
861 ComplexMatrix retval (*this); 847 ComplexMatrix retval (*this);
862 // Note the order of passing the rows and columns to account for 848 Complex *out (retval.fortran_vec ());
863 // column-major storage. 849
864 octave_fftw::ifft2d (retval.fortran_vec (), nc, nr); 850 octave_fftw::ifftNd (out, out, 2, dv);
865 851
866 return retval; 852 return retval;
867 } 853 }
868 854
869 #else 855 #else
996 nn = 4*npts+15; 982 nn = 4*npts+15;
997 983
998 wsave.resize (nn); 984 wsave.resize (nn);
999 pwsave = wsave.fortran_vec (); 985 pwsave = wsave.fortran_vec ();
1000 986
1001 Array<Complex> row (npts); 987 Array<Complex> tmp (npts);
1002 Complex *prow = row.fortran_vec (); 988 Complex *prow = tmp.fortran_vec ();
1003 989
1004 F77_FUNC (cffti, CFFTI) (npts, pwsave); 990 F77_FUNC (cffti, CFFTI) (npts, pwsave);
1005 991
1006 for (int j = 0; j < nsamples; j++) 992 for (int j = 0; j < nsamples; j++)
1007 { 993 {
1065 nn = 4*npts+15; 1051 nn = 4*npts+15;
1066 1052
1067 wsave.resize (nn); 1053 wsave.resize (nn);
1068 pwsave = wsave.fortran_vec (); 1054 pwsave = wsave.fortran_vec ();
1069 1055
1070 Array<Complex> row (npts); 1056 Array<Complex> tmp (npts);
1071 Complex *prow = row.fortran_vec (); 1057 Complex *prow = tmp.fortran_vec ();
1072 1058
1073 F77_FUNC (cffti, CFFTI) (npts, pwsave); 1059 F77_FUNC (cffti, CFFTI) (npts, pwsave);
1074 1060
1075 for (int j = 0; j < nsamples; j++) 1061 for (int j = 0; j < nsamples; j++)
1076 { 1062 {