PLaSK library
Loading...
Searching...
No Matches
mradfg.c
Go to the documentation of this file.
1/* mradfg.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/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
16/* * * */
17/* * copyright (c) 2011 by UCAR * */
18/* * * */
19/* * University Corporation for Atmospheric Research * */
20/* * * */
21/* * all rights reserved * */
22/* * * */
23/* * FFTPACK version 5.1 * */
24/* * * */
25/* * A Fortran Package of Fast Fourier * */
26/* * * */
27/* * Subroutines and Example Programs * */
28/* * * */
29/* * by * */
30/* * * */
31/* * Paul Swarztrauber and Dick Valent * */
32/* * * */
33/* * of * */
34/* * * */
35/* * the National Center for Atmospheric Research * */
36/* * * */
37/* * Boulder, Colorado (80307) U.S.A. * */
38/* * * */
39/* * which is sponsored by * */
40/* * * */
41/* * the National Science Foundation * */
42/* * * */
43/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
44
45/* Subroutine */ int mradfg_(integer *m, integer *ido, integer *ip, integer *
46 l1, integer *idl1, doublereal *cc, doublereal *c1, doublereal *c2,
47 integer *im1, integer *in1, doublereal *ch, doublereal *ch2, integer *
48 im2, integer *in2, doublereal *wa)
49{
50 /* System generated locals */
51 integer ch_dim1, ch_dim2, ch_dim3, ch_offset, cc_dim1, cc_dim2, cc_dim3,
52 cc_offset, c1_dim1, c1_dim2, c1_dim3, c1_offset, c2_dim1, c2_dim2,
53 c2_offset, ch2_dim1, ch2_dim2, ch2_offset, i__1, i__2, i__3,
54 i__4, i__5;
55
56 /* Builtin functions */
57 double atan(doublereal), cos(doublereal), sin(doublereal);
58
59 /* Local variables */
60 integer i__, j, k, l, j2, m1, m2, ic, jc, lc, ik, is;
61 doublereal dc2, ai1, ai2;
62 integer m1d;
63 doublereal ar1, ar2, ds2;
64 integer m2s, nbd;
65 doublereal dcp, arg, dsp, tpi, ar1h, ar2h;
66 integer idp2, ipp2, idij, ipph;
67
68
69 /* Parameter adjustments */
70 --wa;
71 c2_dim1 = *in1;
72 c2_dim2 = *idl1;
73 c2_offset = 1 + c2_dim1 * (1 + c2_dim2);
74 c2 -= c2_offset;
75 c1_dim1 = *in1;
76 c1_dim2 = *ido;
77 c1_dim3 = *l1;
78 c1_offset = 1 + c1_dim1 * (1 + c1_dim2 * (1 + c1_dim3));
79 c1 -= c1_offset;
80 cc_dim1 = *in1;
81 cc_dim2 = *ido;
82 cc_dim3 = *ip;
83 cc_offset = 1 + cc_dim1 * (1 + cc_dim2 * (1 + cc_dim3));
84 cc -= cc_offset;
85 ch2_dim1 = *in2;
86 ch2_dim2 = *idl1;
87 ch2_offset = 1 + ch2_dim1 * (1 + ch2_dim2);
88 ch2 -= ch2_offset;
89 ch_dim1 = *in2;
90 ch_dim2 = *ido;
91 ch_dim3 = *l1;
92 ch_offset = 1 + ch_dim1 * (1 + ch_dim2 * (1 + ch_dim3));
93 ch -= ch_offset;
94
95 /* Function Body */
96 m1d = (*m - 1) * *im1 + 1;
97 m2s = 1 - *im2;
98 tpi = atan(1.) * 8.;
99 arg = tpi / (doublereal) (*ip);
100 dcp = cos(arg);
101 dsp = sin(arg);
102 ipph = (*ip + 1) / 2;
103 ipp2 = *ip + 2;
104 idp2 = *ido + 2;
105 nbd = (*ido - 1) / 2;
106 if (*ido == 1) {
107 goto L119;
108 }
109 i__1 = *idl1;
110 for (ik = 1; ik <= i__1; ++ik) {
111 m2 = m2s;
112 i__2 = m1d;
113 i__3 = *im1;
114 for (m1 = 1; i__3 < 0 ? m1 >= i__2 : m1 <= i__2; m1 += i__3) {
115 m2 += *im2;
116 ch2[m2 + (ik + ch2_dim2) * ch2_dim1] = c2[m1 + (ik + c2_dim2) *
117 c2_dim1];
118/* L1001: */
119 }
120/* L101: */
121 }
122 i__1 = *ip;
123 for (j = 2; j <= i__1; ++j) {
124 i__3 = *l1;
125 for (k = 1; k <= i__3; ++k) {
126 m2 = m2s;
127 i__2 = m1d;
128 i__4 = *im1;
129 for (m1 = 1; i__4 < 0 ? m1 >= i__2 : m1 <= i__2; m1 += i__4) {
130 m2 += *im2;
131 ch[m2 + ((k + j * ch_dim3) * ch_dim2 + 1) * ch_dim1] = c1[m1
132 + ((k + j * c1_dim3) * c1_dim2 + 1) * c1_dim1];
133/* L1002: */
134 }
135/* L102: */
136 }
137/* L103: */
138 }
139 if (nbd > *l1) {
140 goto L107;
141 }
142 is = -(*ido);
143 i__1 = *ip;
144 for (j = 2; j <= i__1; ++j) {
145 is += *ido;
146 idij = is;
147 i__3 = *ido;
148 for (i__ = 3; i__ <= i__3; i__ += 2) {
149 idij += 2;
150 i__4 = *l1;
151 for (k = 1; k <= i__4; ++k) {
152 m2 = m2s;
153 i__2 = m1d;
154 i__5 = *im1;
155 for (m1 = 1; i__5 < 0 ? m1 >= i__2 : m1 <= i__2; m1 += i__5) {
156 m2 += *im2;
157 ch[m2 + (i__ - 1 + (k + j * ch_dim3) * ch_dim2) * ch_dim1]
158 = wa[idij - 1] * c1[m1 + (i__ - 1 + (k + j *
159 c1_dim3) * c1_dim2) * c1_dim1] + wa[idij] * c1[m1
160 + (i__ + (k + j * c1_dim3) * c1_dim2) * c1_dim1];
161 ch[m2 + (i__ + (k + j * ch_dim3) * ch_dim2) * ch_dim1] =
162 wa[idij - 1] * c1[m1 + (i__ + (k + j * c1_dim3) *
163 c1_dim2) * c1_dim1] - wa[idij] * c1[m1 + (i__ - 1
164 + (k + j * c1_dim3) * c1_dim2) * c1_dim1];
165/* L1004: */
166 }
167/* L104: */
168 }
169/* L105: */
170 }
171/* L106: */
172 }
173 goto L111;
174L107:
175 is = -(*ido);
176 i__1 = *ip;
177 for (j = 2; j <= i__1; ++j) {
178 is += *ido;
179 i__3 = *l1;
180 for (k = 1; k <= i__3; ++k) {
181 idij = is;
182 i__4 = *ido;
183 for (i__ = 3; i__ <= i__4; i__ += 2) {
184 idij += 2;
185 m2 = m2s;
186 i__5 = m1d;
187 i__2 = *im1;
188 for (m1 = 1; i__2 < 0 ? m1 >= i__5 : m1 <= i__5; m1 += i__2) {
189 m2 += *im2;
190 ch[m2 + (i__ - 1 + (k + j * ch_dim3) * ch_dim2) * ch_dim1]
191 = wa[idij - 1] * c1[m1 + (i__ - 1 + (k + j *
192 c1_dim3) * c1_dim2) * c1_dim1] + wa[idij] * c1[m1
193 + (i__ + (k + j * c1_dim3) * c1_dim2) * c1_dim1];
194 ch[m2 + (i__ + (k + j * ch_dim3) * ch_dim2) * ch_dim1] =
195 wa[idij - 1] * c1[m1 + (i__ + (k + j * c1_dim3) *
196 c1_dim2) * c1_dim1] - wa[idij] * c1[m1 + (i__ - 1
197 + (k + j * c1_dim3) * c1_dim2) * c1_dim1];
198/* L1008: */
199 }
200/* L108: */
201 }
202/* L109: */
203 }
204/* L110: */
205 }
206L111:
207 if (nbd < *l1) {
208 goto L115;
209 }
210 i__1 = ipph;
211 for (j = 2; j <= i__1; ++j) {
212 jc = ipp2 - j;
213 i__3 = *l1;
214 for (k = 1; k <= i__3; ++k) {
215 i__4 = *ido;
216 for (i__ = 3; i__ <= i__4; i__ += 2) {
217 m2 = m2s;
218 i__2 = m1d;
219 i__5 = *im1;
220 for (m1 = 1; i__5 < 0 ? m1 >= i__2 : m1 <= i__2; m1 += i__5) {
221 m2 += *im2;
222 c1[m1 + (i__ - 1 + (k + j * c1_dim3) * c1_dim2) * c1_dim1]
223 = ch[m2 + (i__ - 1 + (k + j * ch_dim3) * ch_dim2)
224 * ch_dim1] + ch[m2 + (i__ - 1 + (k + jc *
225 ch_dim3) * ch_dim2) * ch_dim1];
226 c1[m1 + (i__ - 1 + (k + jc * c1_dim3) * c1_dim2) *
227 c1_dim1] = ch[m2 + (i__ + (k + j * ch_dim3) *
228 ch_dim2) * ch_dim1] - ch[m2 + (i__ + (k + jc *
229 ch_dim3) * ch_dim2) * ch_dim1];
230 c1[m1 + (i__ + (k + j * c1_dim3) * c1_dim2) * c1_dim1] =
231 ch[m2 + (i__ + (k + j * ch_dim3) * ch_dim2) *
232 ch_dim1] + ch[m2 + (i__ + (k + jc * ch_dim3) *
233 ch_dim2) * ch_dim1];
234 c1[m1 + (i__ + (k + jc * c1_dim3) * c1_dim2) * c1_dim1] =
235 ch[m2 + (i__ - 1 + (k + jc * ch_dim3) * ch_dim2) *
236 ch_dim1] - ch[m2 + (i__ - 1 + (k + j * ch_dim3) *
237 ch_dim2) * ch_dim1];
238/* L1012: */
239 }
240/* L112: */
241 }
242/* L113: */
243 }
244/* L114: */
245 }
246 goto L121;
247L115:
248 i__1 = ipph;
249 for (j = 2; j <= i__1; ++j) {
250 jc = ipp2 - j;
251 i__3 = *ido;
252 for (i__ = 3; i__ <= i__3; i__ += 2) {
253 i__4 = *l1;
254 for (k = 1; k <= i__4; ++k) {
255 m2 = m2s;
256 i__5 = m1d;
257 i__2 = *im1;
258 for (m1 = 1; i__2 < 0 ? m1 >= i__5 : m1 <= i__5; m1 += i__2) {
259 m2 += *im2;
260 c1[m1 + (i__ - 1 + (k + j * c1_dim3) * c1_dim2) * c1_dim1]
261 = ch[m2 + (i__ - 1 + (k + j * ch_dim3) * ch_dim2)
262 * ch_dim1] + ch[m2 + (i__ - 1 + (k + jc *
263 ch_dim3) * ch_dim2) * ch_dim1];
264 c1[m1 + (i__ - 1 + (k + jc * c1_dim3) * c1_dim2) *
265 c1_dim1] = ch[m2 + (i__ + (k + j * ch_dim3) *
266 ch_dim2) * ch_dim1] - ch[m2 + (i__ + (k + jc *
267 ch_dim3) * ch_dim2) * ch_dim1];
268 c1[m1 + (i__ + (k + j * c1_dim3) * c1_dim2) * c1_dim1] =
269 ch[m2 + (i__ + (k + j * ch_dim3) * ch_dim2) *
270 ch_dim1] + ch[m2 + (i__ + (k + jc * ch_dim3) *
271 ch_dim2) * ch_dim1];
272 c1[m1 + (i__ + (k + jc * c1_dim3) * c1_dim2) * c1_dim1] =
273 ch[m2 + (i__ - 1 + (k + jc * ch_dim3) * ch_dim2) *
274 ch_dim1] - ch[m2 + (i__ - 1 + (k + j * ch_dim3) *
275 ch_dim2) * ch_dim1];
276/* L1016: */
277 }
278/* L116: */
279 }
280/* L117: */
281 }
282/* L118: */
283 }
284 goto L121;
285L119:
286 i__1 = *idl1;
287 for (ik = 1; ik <= i__1; ++ik) {
288 m2 = m2s;
289 i__3 = m1d;
290 i__4 = *im1;
291 for (m1 = 1; i__4 < 0 ? m1 >= i__3 : m1 <= i__3; m1 += i__4) {
292 m2 += *im2;
293 c2[m1 + (ik + c2_dim2) * c2_dim1] = ch2[m2 + (ik + ch2_dim2) *
294 ch2_dim1];
295/* L1020: */
296 }
297/* L120: */
298 }
299L121:
300 i__1 = ipph;
301 for (j = 2; j <= i__1; ++j) {
302 jc = ipp2 - j;
303 i__4 = *l1;
304 for (k = 1; k <= i__4; ++k) {
305 m2 = m2s;
306 i__3 = m1d;
307 i__2 = *im1;
308 for (m1 = 1; i__2 < 0 ? m1 >= i__3 : m1 <= i__3; m1 += i__2) {
309 m2 += *im2;
310 c1[m1 + ((k + j * c1_dim3) * c1_dim2 + 1) * c1_dim1] = ch[m2
311 + ((k + j * ch_dim3) * ch_dim2 + 1) * ch_dim1] + ch[
312 m2 + ((k + jc * ch_dim3) * ch_dim2 + 1) * ch_dim1];
313 c1[m1 + ((k + jc * c1_dim3) * c1_dim2 + 1) * c1_dim1] = ch[m2
314 + ((k + jc * ch_dim3) * ch_dim2 + 1) * ch_dim1] - ch[
315 m2 + ((k + j * ch_dim3) * ch_dim2 + 1) * ch_dim1];
316/* L1022: */
317 }
318/* L122: */
319 }
320/* L123: */
321 }
322
323 ar1 = 1.;
324 ai1 = 0.;
325 i__1 = ipph;
326 for (l = 2; l <= i__1; ++l) {
327 lc = ipp2 - l;
328 ar1h = dcp * ar1 - dsp * ai1;
329 ai1 = dcp * ai1 + dsp * ar1;
330 ar1 = ar1h;
331 i__4 = *idl1;
332 for (ik = 1; ik <= i__4; ++ik) {
333 m2 = m2s;
334 i__2 = m1d;
335 i__3 = *im1;
336 for (m1 = 1; i__3 < 0 ? m1 >= i__2 : m1 <= i__2; m1 += i__3) {
337 m2 += *im2;
338 ch2[m2 + (ik + l * ch2_dim2) * ch2_dim1] = c2[m1 + (ik +
339 c2_dim2) * c2_dim1] + ar1 * c2[m1 + (ik + (c2_dim2 <<
340 1)) * c2_dim1];
341 ch2[m2 + (ik + lc * ch2_dim2) * ch2_dim1] = ai1 * c2[m1 + (ik
342 + *ip * c2_dim2) * c2_dim1];
343/* L1024: */
344 }
345/* L124: */
346 }
347 dc2 = ar1;
348 ds2 = ai1;
349 ar2 = ar1;
350 ai2 = ai1;
351 i__4 = ipph;
352 for (j = 3; j <= i__4; ++j) {
353 jc = ipp2 - j;
354 ar2h = dc2 * ar2 - ds2 * ai2;
355 ai2 = dc2 * ai2 + ds2 * ar2;
356 ar2 = ar2h;
357 i__3 = *idl1;
358 for (ik = 1; ik <= i__3; ++ik) {
359 m2 = m2s;
360 i__2 = m1d;
361 i__5 = *im1;
362 for (m1 = 1; i__5 < 0 ? m1 >= i__2 : m1 <= i__2; m1 += i__5) {
363 m2 += *im2;
364 ch2[m2 + (ik + l * ch2_dim2) * ch2_dim1] += ar2 * c2[m1 +
365 (ik + j * c2_dim2) * c2_dim1];
366 ch2[m2 + (ik + lc * ch2_dim2) * ch2_dim1] += ai2 * c2[m1
367 + (ik + jc * c2_dim2) * c2_dim1];
368/* L1025: */
369 }
370/* L125: */
371 }
372/* L126: */
373 }
374/* L127: */
375 }
376 i__1 = ipph;
377 for (j = 2; j <= i__1; ++j) {
378 i__4 = *idl1;
379 for (ik = 1; ik <= i__4; ++ik) {
380 m2 = m2s;
381 i__3 = m1d;
382 i__5 = *im1;
383 for (m1 = 1; i__5 < 0 ? m1 >= i__3 : m1 <= i__3; m1 += i__5) {
384 m2 += *im2;
385 ch2[m2 + (ik + ch2_dim2) * ch2_dim1] += c2[m1 + (ik + j *
386 c2_dim2) * c2_dim1];
387/* L1028: */
388 }
389/* L128: */
390 }
391/* L129: */
392 }
393
394 if (*ido < *l1) {
395 goto L132;
396 }
397 i__1 = *l1;
398 for (k = 1; k <= i__1; ++k) {
399 i__4 = *ido;
400 for (i__ = 1; i__ <= i__4; ++i__) {
401 m2 = m2s;
402 i__5 = m1d;
403 i__3 = *im1;
404 for (m1 = 1; i__3 < 0 ? m1 >= i__5 : m1 <= i__5; m1 += i__3) {
405 m2 += *im2;
406 cc[m1 + (i__ + (k * cc_dim3 + 1) * cc_dim2) * cc_dim1] = ch[
407 m2 + (i__ + (k + ch_dim3) * ch_dim2) * ch_dim1];
408/* L1030: */
409 }
410/* L130: */
411 }
412/* L131: */
413 }
414 goto L135;
415L132:
416 i__1 = *ido;
417 for (i__ = 1; i__ <= i__1; ++i__) {
418 i__4 = *l1;
419 for (k = 1; k <= i__4; ++k) {
420 m2 = m2s;
421 i__3 = m1d;
422 i__5 = *im1;
423 for (m1 = 1; i__5 < 0 ? m1 >= i__3 : m1 <= i__3; m1 += i__5) {
424 m2 += *im2;
425 cc[m1 + (i__ + (k * cc_dim3 + 1) * cc_dim2) * cc_dim1] = ch[
426 m2 + (i__ + (k + ch_dim3) * ch_dim2) * ch_dim1];
427/* L1033: */
428 }
429/* L133: */
430 }
431/* L134: */
432 }
433L135:
434 i__1 = ipph;
435 for (j = 2; j <= i__1; ++j) {
436 jc = ipp2 - j;
437 j2 = j + j;
438 i__4 = *l1;
439 for (k = 1; k <= i__4; ++k) {
440 m2 = m2s;
441 i__5 = m1d;
442 i__3 = *im1;
443 for (m1 = 1; i__3 < 0 ? m1 >= i__5 : m1 <= i__5; m1 += i__3) {
444 m2 += *im2;
445 cc[m1 + (*ido + (j2 - 2 + k * cc_dim3) * cc_dim2) * cc_dim1] =
446 ch[m2 + ((k + j * ch_dim3) * ch_dim2 + 1) * ch_dim1];
447 cc[m1 + ((j2 - 1 + k * cc_dim3) * cc_dim2 + 1) * cc_dim1] =
448 ch[m2 + ((k + jc * ch_dim3) * ch_dim2 + 1) * ch_dim1];
449/* L1036: */
450 }
451/* L136: */
452 }
453/* L137: */
454 }
455 if (*ido == 1) {
456 return 0;
457 }
458 if (nbd < *l1) {
459 goto L141;
460 }
461 i__1 = ipph;
462 for (j = 2; j <= i__1; ++j) {
463 jc = ipp2 - j;
464 j2 = j + j;
465 i__4 = *l1;
466 for (k = 1; k <= i__4; ++k) {
467 i__3 = *ido;
468 for (i__ = 3; i__ <= i__3; i__ += 2) {
469 ic = idp2 - i__;
470 m2 = m2s;
471 i__5 = m1d;
472 i__2 = *im1;
473 for (m1 = 1; i__2 < 0 ? m1 >= i__5 : m1 <= i__5; m1 += i__2) {
474 m2 += *im2;
475 cc[m1 + (i__ - 1 + (j2 - 1 + k * cc_dim3) * cc_dim2) *
476 cc_dim1] = ch[m2 + (i__ - 1 + (k + j * ch_dim3) *
477 ch_dim2) * ch_dim1] + ch[m2 + (i__ - 1 + (k + jc *
478 ch_dim3) * ch_dim2) * ch_dim1];
479 cc[m1 + (ic - 1 + (j2 - 2 + k * cc_dim3) * cc_dim2) *
480 cc_dim1] = ch[m2 + (i__ - 1 + (k + j * ch_dim3) *
481 ch_dim2) * ch_dim1] - ch[m2 + (i__ - 1 + (k + jc *
482 ch_dim3) * ch_dim2) * ch_dim1];
483 cc[m1 + (i__ + (j2 - 1 + k * cc_dim3) * cc_dim2) *
484 cc_dim1] = ch[m2 + (i__ + (k + j * ch_dim3) *
485 ch_dim2) * ch_dim1] + ch[m2 + (i__ + (k + jc *
486 ch_dim3) * ch_dim2) * ch_dim1];
487 cc[m1 + (ic + (j2 - 2 + k * cc_dim3) * cc_dim2) * cc_dim1]
488 = ch[m2 + (i__ + (k + jc * ch_dim3) * ch_dim2) *
489 ch_dim1] - ch[m2 + (i__ + (k + j * ch_dim3) *
490 ch_dim2) * ch_dim1];
491/* L1038: */
492 }
493/* L138: */
494 }
495/* L139: */
496 }
497/* L140: */
498 }
499 return 0;
500L141:
501 i__1 = ipph;
502 for (j = 2; j <= i__1; ++j) {
503 jc = ipp2 - j;
504 j2 = j + j;
505 i__4 = *ido;
506 for (i__ = 3; i__ <= i__4; i__ += 2) {
507 ic = idp2 - i__;
508 i__3 = *l1;
509 for (k = 1; k <= i__3; ++k) {
510 m2 = m2s;
511 i__2 = m1d;
512 i__5 = *im1;
513 for (m1 = 1; i__5 < 0 ? m1 >= i__2 : m1 <= i__2; m1 += i__5) {
514 m2 += *im2;
515 cc[m1 + (i__ - 1 + (j2 - 1 + k * cc_dim3) * cc_dim2) *
516 cc_dim1] = ch[m2 + (i__ - 1 + (k + j * ch_dim3) *
517 ch_dim2) * ch_dim1] + ch[m2 + (i__ - 1 + (k + jc *
518 ch_dim3) * ch_dim2) * ch_dim1];
519 cc[m1 + (ic - 1 + (j2 - 2 + k * cc_dim3) * cc_dim2) *
520 cc_dim1] = ch[m2 + (i__ - 1 + (k + j * ch_dim3) *
521 ch_dim2) * ch_dim1] - ch[m2 + (i__ - 1 + (k + jc *
522 ch_dim3) * ch_dim2) * ch_dim1];
523 cc[m1 + (i__ + (j2 - 1 + k * cc_dim3) * cc_dim2) *
524 cc_dim1] = ch[m2 + (i__ + (k + j * ch_dim3) *
525 ch_dim2) * ch_dim1] + ch[m2 + (i__ + (k + jc *
526 ch_dim3) * ch_dim2) * ch_dim1];
527 cc[m1 + (ic + (j2 - 2 + k * cc_dim3) * cc_dim2) * cc_dim1]
528 = ch[m2 + (i__ + (k + jc * ch_dim3) * ch_dim2) *
529 ch_dim1] - ch[m2 + (i__ + (k + j * ch_dim3) *
530 ch_dim2) * ch_dim1];
531/* L1042: */
532 }
533/* L142: */
534 }
535/* L143: */
536 }
537/* L144: */
538 }
539 return 0;
540} /* mradfg_ */
541