1 | /* --------------- apcxx_cstyle.icc -----------------
|
---|
2 |
|
---|
3 | Example illustrating C++ scientific programming
|
---|
4 | using SOPHYA library arrays and FFT computation
|
---|
5 | through FFTServerInterface.
|
---|
6 | This example, performs the same operations as in
|
---|
7 | apcxx.icc, through C-style programming with
|
---|
8 | explicit loops. This C_style should be avoided when
|
---|
9 | programming in C++.
|
---|
10 |
|
---|
11 | R. Ansari 08/2001
|
---|
12 |
|
---|
13 | ---- Computation steps :
|
---|
14 | > Create a matrix (NL x NC ) (mtx)
|
---|
15 | > fill it with a gaussian distributed random values
|
---|
16 | > make a copy of the matrix (mtxs)
|
---|
17 | > create a 1D filter in Fourier space
|
---|
18 | > Loop over matrix rows k
|
---|
19 | >> Extract row k (fline)
|
---|
20 | >> compute 1D Fourier transform FFTForward (fline)
|
---|
21 | >> apply filter in Fourier space
|
---|
22 | >> compute backward 1D FFT
|
---|
23 | >> Replace matrix row with the filtered values
|
---|
24 |
|
---|
25 |
|
---|
26 | this example code can be
|
---|
27 | - included in a main program
|
---|
28 | - executed using runcxx
|
---|
29 | csh> runcxx -tmpdir /tmp -f apcxx_cstyle.icc
|
---|
30 | - executed within spiapp
|
---|
31 | Cmd> c++execfrf apcxx_cstyle.icc
|
---|
32 | */
|
---|
33 |
|
---|
34 | // Select computation on float or double (r_4 r_8)
|
---|
35 | #define FTYP r_4
|
---|
36 |
|
---|
37 | // Number of matrix lines and colums
|
---|
38 | int NL = 1024;
|
---|
39 | int NC = 4096;
|
---|
40 | int i, j;
|
---|
41 | cout << " apc_cstyle : NL= " << NL << " NC= " << NC << endl;
|
---|
42 | PrtTim("apc_cstyle_Start");
|
---|
43 |
|
---|
44 | //BaseArray::SetDefaultMemoryMapping(BaseArray::CMemoryMapping);
|
---|
45 | //BaseArray::SetMaxPrint(10, 3);
|
---|
46 |
|
---|
47 |
|
---|
48 | // Creation of the initial matrix
|
---|
49 | TMatrix< FTYP > mtx(NL, NC);
|
---|
50 | // Filling matrix with gaussian random values
|
---|
51 | // mtx = RandomSequence(RandomSequence::Gaussian, 15., 3.);
|
---|
52 | for(i=0; i<NL; i++)
|
---|
53 | for(j=0; j<NC; j++) mtx(i, j) = 15.+NorRand()*3.;
|
---|
54 |
|
---|
55 | // Making a copy of the original matrix
|
---|
56 | //mtxs = mtx;
|
---|
57 | TMatrix< FTYP > mtxs(NL, NC);
|
---|
58 | for(i=0; i<NL; i++)
|
---|
59 | for(j=0; j<NC; j++) mtxs(i, j) = mtx(i, j);
|
---|
60 |
|
---|
61 |
|
---|
62 | // Creation and initialization of the Fourier filter filt(nu) = 1/(1+0.3*nu)
|
---|
63 | int LFFT = NC/2+1;
|
---|
64 | TVector< complex< FTYP > > filt(LFFT, BaseArray::RowVector);
|
---|
65 | filt(0) = 1.;
|
---|
66 | for(int i=1; i<filt.Size(); i++)
|
---|
67 | filt(i) = 1./(1+0.3*(double)i);
|
---|
68 |
|
---|
69 | // Creation of the FFTServer
|
---|
70 | FFTPackServer ffts;
|
---|
71 | ffts.setNormalize(true);
|
---|
72 |
|
---|
73 | PrtTim("apc_cstyle_AfterInit");
|
---|
74 |
|
---|
75 | // Vectors for FFT operations
|
---|
76 | TVector< FTYP > fline(NC, BaseArray::RowVector);
|
---|
77 | TVector< complex< FTYP > > vfft;
|
---|
78 |
|
---|
79 | for(int k=0; k<NL; k++) {
|
---|
80 | // fline = mtx.Row(k); Matrix row extraction
|
---|
81 | for(j=0; j<NC; j++) fline(j) = mtx(k, j);
|
---|
82 | ffts.FFTForward(fline, vfft); // Compute 1D forward FFT
|
---|
83 | // Applying filter in Fourier space f(nu) = f(nu)*filter(nu)
|
---|
84 | // vfft.MulElt(filt);
|
---|
85 | for(j=0; j<LFFT; j++) vfft(j) *= filt(j);
|
---|
86 | ffts.FFTBackward(vfft, fline); // backward FFT
|
---|
87 | // mtx.Row(k) = fline; // replace matrix row with filtered values
|
---|
88 | for(j=0; j<NC; j++) mtx(k, j) = fline(j);
|
---|
89 | }
|
---|
90 |
|
---|
91 | PrtTim("apc_cstyle_AfterFFTLoop");
|
---|
92 |
|
---|
93 | // Macro KeepObj can be used with runcxx or within (s)piapp
|
---|
94 | // KeepObj(mtx);
|
---|
95 | // KeepObj(mtxs);
|
---|
96 |
|
---|