PLaSK library
Loading...
Searching...
No Matches
rfft2f.c
Go to the documentation of this file.
1/* rfft2f.f -- translated by f2c (version 20100827).
2 You must link the resulting object file with libf2c:
3 on Microsoft Windows system, link with libf2c.lib;
4 on Linux or Unix systems, link with .../path/to/libf2c.a -lm
5 or, if you install libf2c.a in a standard place, with -lf2c -lm
6 -- in that order, at the end of the command line, as in
7 cc *.o -lf2c -lm
8 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
9
10 http://www.netlib.org/f2c/libf2c.zip
11*/
12
13#include "f2c.h"
14
15/* Table of constant values */
16
17static integer c__6 = 6;
18static integer c__8 = 8;
19static integer c_n6 = -6;
20static integer c__1 = 1;
21static integer c_n5 = -5;
22
23/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
24/* * * */
25/* * copyright (c) 2011 by UCAR * */
26/* * * */
27/* * University Corporation for Atmospheric Research * */
28/* * * */
29/* * all rights reserved * */
30/* * * */
31/* * FFTPACK version 5.1 * */
32/* * * */
33/* * A Fortran Package of Fast Fourier * */
34/* * * */
35/* * Subroutines and Example Programs * */
36/* * * */
37/* * by * */
38/* * * */
39/* * Paul Swarztrauber and Dick Valent * */
40/* * * */
41/* * of * */
42/* * * */
43/* * the National Center for Atmospheric Research * */
44/* * * */
45/* * Boulder, Colorado (80307) U.S.A. * */
46/* * * */
47/* * which is sponsored by * */
48/* * * */
49/* * the National Science Foundation * */
50/* * * */
51/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
52
53/* Subroutine */ int rfft2f_(integer *ldim, integer *l, integer *m,
54 doublereal *r__, doublereal *wsave, integer *lensav, doublereal *work,
55 integer *lenwrk, integer *ier)
56{
57 /* System generated locals */
58 integer r_dim1, r_offset, i__1, i__2, i__3;
59
60 /* Builtin functions */
61 double log(doublereal);
62
63 /* Local variables */
64 integer i__, j;
65 extern /* Subroutine */ int r2w_(integer *, integer *, integer *, integer
66 *, doublereal *, doublereal *), w2r_(integer *, integer *,
67 integer *, integer *, doublereal *, doublereal *);
68 integer ldh, ldw, ldx, ier1, modl, modm, mmsav, lwsav, mwsav;
69 extern /* Subroutine */ int cfftmf_(integer *, integer *, integer *,
70 integer *, doublereal *, integer *, doublereal *, integer *,
71 doublereal *, integer *, integer *), rfftmf_(integer *, integer *,
72 integer *, integer *, doublereal *, integer *, doublereal *,
73 integer *, doublereal *, integer *, integer *), xerfft_(char *,
74 integer *, ftnlen);
75
76
77
78/* INITIALIZE IER */
79
80 /* Parameter adjustments */
81 r_dim1 = *ldim;
82 r_offset = 1 + r_dim1;
83 r__ -= r_offset;
84 --wsave;
85 --work;
86
87 /* Function Body */
88 *ier = 0;
89
90/* VERIFY LENSAV */
91
92 lwsav = *l + (integer) (log((doublereal) (*l)) / log(2.)) + 4;
93 mwsav = (*m << 1) + (integer) (log((doublereal) (*m)) / log(2.)) + 4;
94 mmsav = *m + (integer) (log((doublereal) (*m)) / log(2.)) + 4;
95
96 if (*lensav < lwsav + mwsav + mmsav) {
97 *ier = 2;
98 xerfft_("RFFT2F", &c__6, (ftnlen)6);
99 goto L100;
100 }
101
102/* VERIFY LENWRK */
103
104 if (*lenwrk < (*l + 1) * *m) {
105 *ier = 3;
106 xerfft_("RFFT2F", &c__8, (ftnlen)6);
107 goto L100;
108 }
109
110/* VERIFY LDIM IS AS BIG AS L */
111
112 if (*ldim < *l) {
113 *ier = 5;
114 xerfft_("RFFT2F", &c_n6, (ftnlen)6);
115 goto L100;
116 }
117
118/* TRANSFORM FIRST DIMENSION OF ARRAY */
119
120 i__1 = *m * *ldim;
121 i__2 = *l + (integer) (log((doublereal) (*l)) / log(2.)) + 4;
122 rfftmf_(m, ldim, l, &c__1, &r__[r_offset], &i__1, &wsave[1], &i__2, &work[
123 1], lenwrk, &ier1);
124
125 if (ier1 != 0) {
126 *ier = 20;
127 xerfft_("RFFT2F", &c_n5, (ftnlen)6);
128 goto L100;
129 }
130
131 ldx = ((*l + 1) / 2 << 1) - 1;
132 i__1 = ldx;
133 for (i__ = 2; i__ <= i__1; ++i__) {
134 i__2 = *m;
135 for (j = 1; j <= i__2; ++j) {
136 r__[i__ + j * r_dim1] *= .5;
137 }
138 }
139 i__1 = *m;
140 for (j = 1; j <= i__1; ++j) {
141 i__2 = ldx;
142 for (i__ = 3; i__ <= i__2; i__ += 2) {
143 r__[i__ + j * r_dim1] = -r__[i__ + j * r_dim1];
144 }
145 }
146
147/* PRINT*, 'FORWARD TRANSFORM IN THE I DIRECTION' */
148/* DO I=1,L */
149/* PRINT*, (R(I,J),J=1,M) */
150/* END DO */
151
152/* RESHUFFLE TO ADD IN NYQUIST IMAGINARY COMPONENTS */
153
154 modl = *l % 2;
155 modm = *m % 2;
156
157/* TRANSFORM SECOND DIMENSION OF ARRAY */
158
159 i__1 = *m * *ldim;
160 rfftmf_(&c__1, &c__1, m, ldim, &r__[r_offset], &i__1, &wsave[lwsav +
161 mwsav + 1], &mmsav, &work[1], lenwrk, &ier1);
162 i__1 = ((*m + 1) / 2 << 1) - 1;
163 for (j = 2; j <= i__1; ++j) {
164 r__[j * r_dim1 + 1] *= .5;
165 }
166 i__1 = *m;
167 for (j = 3; j <= i__1; j += 2) {
168 r__[j * r_dim1 + 1] = -r__[j * r_dim1 + 1];
169 }
170 ldh = (*l + 1) / 2;
171 if (ldh > 1) {
172 ldw = ldh + ldh;
173
174/* R AND WORK ARE SWITCHED BECAUSE THE THE FIRST DIMENSION */
175/* OF THE INPUT TO COMPLEX CFFTMF MUST BE EVEN. */
176
177 r2w_(ldim, &ldw, l, m, &r__[r_offset], &work[1]);
178 i__1 = ldh - 1;
179 i__2 = ldh * *m;
180 i__3 = *l * *m;
181 cfftmf_(&i__1, &c__1, m, &ldh, &work[2], &i__2, &wsave[lwsav + 1], &
182 mwsav, &r__[r_offset], &i__3, &ier1);
183 if (ier1 != 0) {
184 *ier = 20;
185 xerfft_("RFFT2F", &c_n5, (ftnlen)6);
186 goto L100;
187 }
188 w2r_(ldim, &ldw, l, m, &r__[r_offset], &work[1]);
189 }
190
191 if (modl == 0) {
192 i__1 = *m * *ldim;
193 rfftmf_(&c__1, &c__1, m, ldim, &r__[*l + r_dim1], &i__1, &wsave[lwsav
194 + mwsav + 1], &mmsav, &work[1], lenwrk, &ier1);
195 i__1 = ((*m + 1) / 2 << 1) - 1;
196 for (j = 2; j <= i__1; ++j) {
197 r__[*l + j * r_dim1] *= .5;
198 }
199 i__1 = *m;
200 for (j = 3; j <= i__1; j += 2) {
201 r__[*l + j * r_dim1] = -r__[*l + j * r_dim1];
202 }
203 }
204
205/* PRINT*, 'FORWARD TRANSFORM IN THE J DIRECTION' */
206/* DO I=1,L */
207/* PRINT*, (R(I,J),J=1,M) */
208/* END DO */
209
210 if (ier1 != 0) {
211 *ier = 20;
212 xerfft_("RFFT2F", &c_n5, (ftnlen)6);
213 goto L100;
214 }
215
216
217L100:
218 return 0;
219} /* rfft2f_ */
220