PLaSK library
Loading...
Searching...
No Matches
mradbg.c
Go to the documentation of this file.
1/* mradbg.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 mradbg_(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 idp2 = *ido + 2;
103 nbd = (*ido - 1) / 2;
104 ipp2 = *ip + 2;
105 ipph = (*ip + 1) / 2;
106 if (*ido < *l1) {
107 goto L103;
108 }
109 i__1 = *l1;
110 for (k = 1; k <= i__1; ++k) {
111 i__2 = *ido;
112 for (i__ = 1; i__ <= i__2; ++i__) {
113 m2 = m2s;
114 i__3 = m1d;
115 i__4 = *im1;
116 for (m1 = 1; i__4 < 0 ? m1 >= i__3 : m1 <= i__3; m1 += i__4) {
117 m2 += *im2;
118 ch[m2 + (i__ + (k + ch_dim3) * ch_dim2) * ch_dim1] = cc[m1 + (
119 i__ + (k * cc_dim3 + 1) * cc_dim2) * cc_dim1];
120/* L1001: */
121 }
122/* L101: */
123 }
124/* L102: */
125 }
126 goto L106;
127L103:
128 i__1 = *ido;
129 for (i__ = 1; i__ <= i__1; ++i__) {
130 i__2 = *l1;
131 for (k = 1; k <= i__2; ++k) {
132 m2 = m2s;
133 i__4 = m1d;
134 i__3 = *im1;
135 for (m1 = 1; i__3 < 0 ? m1 >= i__4 : m1 <= i__4; m1 += i__3) {
136 m2 += *im2;
137 ch[m2 + (i__ + (k + ch_dim3) * ch_dim2) * ch_dim1] = cc[m1 + (
138 i__ + (k * cc_dim3 + 1) * cc_dim2) * cc_dim1];
139/* L1004: */
140 }
141/* L104: */
142 }
143/* L105: */
144 }
145L106:
146 i__1 = ipph;
147 for (j = 2; j <= i__1; ++j) {
148 jc = ipp2 - j;
149 j2 = j + j;
150 i__2 = *l1;
151 for (k = 1; k <= i__2; ++k) {
152 m2 = m2s;
153 i__3 = m1d;
154 i__4 = *im1;
155 for (m1 = 1; i__4 < 0 ? m1 >= i__3 : m1 <= i__3; m1 += i__4) {
156 m2 += *im2;
157 ch[m2 + ((k + j * ch_dim3) * ch_dim2 + 1) * ch_dim1] = cc[m1
158 + (*ido + (j2 - 2 + k * cc_dim3) * cc_dim2) * cc_dim1]
159 + cc[m1 + (*ido + (j2 - 2 + k * cc_dim3) * cc_dim2) *
160 cc_dim1];
161 ch[m2 + ((k + jc * ch_dim3) * ch_dim2 + 1) * ch_dim1] = cc[m1
162 + ((j2 - 1 + k * cc_dim3) * cc_dim2 + 1) * cc_dim1] +
163 cc[m1 + ((j2 - 1 + k * cc_dim3) * cc_dim2 + 1) *
164 cc_dim1];
165/* L1007: */
166 }
167/* L107: */
168 }
169/* L108: */
170 }
171 if (*ido == 1) {
172 goto L116;
173 }
174 if (nbd < *l1) {
175 goto L112;
176 }
177 i__1 = ipph;
178 for (j = 2; j <= i__1; ++j) {
179 jc = ipp2 - j;
180 i__2 = *l1;
181 for (k = 1; k <= i__2; ++k) {
182 i__4 = *ido;
183 for (i__ = 3; i__ <= i__4; i__ += 2) {
184 ic = idp2 - i__;
185 m2 = m2s;
186 i__3 = m1d;
187 i__5 = *im1;
188 for (m1 = 1; i__5 < 0 ? m1 >= i__3 : m1 <= i__3; m1 += i__5) {
189 m2 += *im2;
190 ch[m2 + (i__ - 1 + (k + j * ch_dim3) * ch_dim2) * ch_dim1]
191 = cc[m1 + (i__ - 1 + ((j << 1) - 1 + k * cc_dim3)
192 * cc_dim2) * cc_dim1] + cc[m1 + (ic - 1 + ((j <<
193 1) - 2 + k * cc_dim3) * cc_dim2) * cc_dim1];
194 ch[m2 + (i__ - 1 + (k + jc * ch_dim3) * ch_dim2) *
195 ch_dim1] = cc[m1 + (i__ - 1 + ((j << 1) - 1 + k *
196 cc_dim3) * cc_dim2) * cc_dim1] - cc[m1 + (ic - 1
197 + ((j << 1) - 2 + k * cc_dim3) * cc_dim2) *
198 cc_dim1];
199 ch[m2 + (i__ + (k + j * ch_dim3) * ch_dim2) * ch_dim1] =
200 cc[m1 + (i__ + ((j << 1) - 1 + k * cc_dim3) *
201 cc_dim2) * cc_dim1] - cc[m1 + (ic + ((j << 1) - 2
202 + k * cc_dim3) * cc_dim2) * cc_dim1];
203 ch[m2 + (i__ + (k + jc * ch_dim3) * ch_dim2) * ch_dim1] =
204 cc[m1 + (i__ + ((j << 1) - 1 + k * cc_dim3) *
205 cc_dim2) * cc_dim1] + cc[m1 + (ic + ((j << 1) - 2
206 + k * cc_dim3) * cc_dim2) * cc_dim1];
207/* L1009: */
208 }
209/* L109: */
210 }
211/* L110: */
212 }
213/* L111: */
214 }
215 goto L116;
216L112:
217 i__1 = ipph;
218 for (j = 2; j <= i__1; ++j) {
219 jc = ipp2 - j;
220 i__2 = *ido;
221 for (i__ = 3; i__ <= i__2; i__ += 2) {
222 ic = idp2 - i__;
223 i__4 = *l1;
224 for (k = 1; k <= i__4; ++k) {
225 m2 = m2s;
226 i__5 = m1d;
227 i__3 = *im1;
228 for (m1 = 1; i__3 < 0 ? m1 >= i__5 : m1 <= i__5; m1 += i__3) {
229 m2 += *im2;
230 ch[m2 + (i__ - 1 + (k + j * ch_dim3) * ch_dim2) * ch_dim1]
231 = cc[m1 + (i__ - 1 + ((j << 1) - 1 + k * cc_dim3)
232 * cc_dim2) * cc_dim1] + cc[m1 + (ic - 1 + ((j <<
233 1) - 2 + k * cc_dim3) * cc_dim2) * cc_dim1];
234 ch[m2 + (i__ - 1 + (k + jc * ch_dim3) * ch_dim2) *
235 ch_dim1] = cc[m1 + (i__ - 1 + ((j << 1) - 1 + k *
236 cc_dim3) * cc_dim2) * cc_dim1] - cc[m1 + (ic - 1
237 + ((j << 1) - 2 + k * cc_dim3) * cc_dim2) *
238 cc_dim1];
239 ch[m2 + (i__ + (k + j * ch_dim3) * ch_dim2) * ch_dim1] =
240 cc[m1 + (i__ + ((j << 1) - 1 + k * cc_dim3) *
241 cc_dim2) * cc_dim1] - cc[m1 + (ic + ((j << 1) - 2
242 + k * cc_dim3) * cc_dim2) * cc_dim1];
243 ch[m2 + (i__ + (k + jc * ch_dim3) * ch_dim2) * ch_dim1] =
244 cc[m1 + (i__ + ((j << 1) - 1 + k * cc_dim3) *
245 cc_dim2) * cc_dim1] + cc[m1 + (ic + ((j << 1) - 2
246 + k * cc_dim3) * cc_dim2) * cc_dim1];
247/* L1013: */
248 }
249/* L113: */
250 }
251/* L114: */
252 }
253/* L115: */
254 }
255L116:
256 ar1 = 1.;
257 ai1 = 0.;
258 i__1 = ipph;
259 for (l = 2; l <= i__1; ++l) {
260 lc = ipp2 - l;
261 ar1h = dcp * ar1 - dsp * ai1;
262 ai1 = dcp * ai1 + dsp * ar1;
263 ar1 = ar1h;
264 i__2 = *idl1;
265 for (ik = 1; ik <= i__2; ++ik) {
266 m2 = m2s;
267 i__4 = m1d;
268 i__3 = *im1;
269 for (m1 = 1; i__3 < 0 ? m1 >= i__4 : m1 <= i__4; m1 += i__3) {
270 m2 += *im2;
271 c2[m1 + (ik + l * c2_dim2) * c2_dim1] = ch2[m2 + (ik +
272 ch2_dim2) * ch2_dim1] + ar1 * ch2[m2 + (ik + (
273 ch2_dim2 << 1)) * ch2_dim1];
274 c2[m1 + (ik + lc * c2_dim2) * c2_dim1] = ai1 * ch2[m2 + (ik +
275 *ip * ch2_dim2) * ch2_dim1];
276/* L1017: */
277 }
278/* L117: */
279 }
280 dc2 = ar1;
281 ds2 = ai1;
282 ar2 = ar1;
283 ai2 = ai1;
284 i__2 = ipph;
285 for (j = 3; j <= i__2; ++j) {
286 jc = ipp2 - j;
287 ar2h = dc2 * ar2 - ds2 * ai2;
288 ai2 = dc2 * ai2 + ds2 * ar2;
289 ar2 = ar2h;
290 i__3 = *idl1;
291 for (ik = 1; ik <= i__3; ++ik) {
292 m2 = m2s;
293 i__4 = m1d;
294 i__5 = *im1;
295 for (m1 = 1; i__5 < 0 ? m1 >= i__4 : m1 <= i__4; m1 += i__5) {
296 m2 += *im2;
297 c2[m1 + (ik + l * c2_dim2) * c2_dim1] += ar2 * ch2[m2 + (
298 ik + j * ch2_dim2) * ch2_dim1];
299 c2[m1 + (ik + lc * c2_dim2) * c2_dim1] += ai2 * ch2[m2 + (
300 ik + jc * ch2_dim2) * ch2_dim1];
301/* L1018: */
302 }
303/* L118: */
304 }
305/* L119: */
306 }
307/* L120: */
308 }
309 i__1 = ipph;
310 for (j = 2; j <= i__1; ++j) {
311 i__2 = *idl1;
312 for (ik = 1; ik <= i__2; ++ik) {
313 m2 = m2s;
314 i__3 = m1d;
315 i__5 = *im1;
316 for (m1 = 1; i__5 < 0 ? m1 >= i__3 : m1 <= i__3; m1 += i__5) {
317 m2 += *im2;
318 ch2[m2 + (ik + ch2_dim2) * ch2_dim1] += ch2[m2 + (ik + j *
319 ch2_dim2) * ch2_dim1];
320/* L1021: */
321 }
322/* L121: */
323 }
324/* L122: */
325 }
326 i__1 = ipph;
327 for (j = 2; j <= i__1; ++j) {
328 jc = ipp2 - j;
329 i__2 = *l1;
330 for (k = 1; k <= i__2; ++k) {
331 m2 = m2s;
332 i__5 = m1d;
333 i__3 = *im1;
334 for (m1 = 1; i__3 < 0 ? m1 >= i__5 : m1 <= i__5; m1 += i__3) {
335 m2 += *im2;
336 ch[m2 + ((k + j * ch_dim3) * ch_dim2 + 1) * ch_dim1] = c1[m1
337 + ((k + j * c1_dim3) * c1_dim2 + 1) * c1_dim1] - c1[
338 m1 + ((k + jc * c1_dim3) * c1_dim2 + 1) * c1_dim1];
339 ch[m2 + ((k + jc * ch_dim3) * ch_dim2 + 1) * ch_dim1] = c1[m1
340 + ((k + j * c1_dim3) * c1_dim2 + 1) * c1_dim1] + c1[
341 m1 + ((k + jc * c1_dim3) * c1_dim2 + 1) * c1_dim1];
342/* L1023: */
343 }
344/* L123: */
345 }
346/* L124: */
347 }
348 if (*ido == 1) {
349 goto L132;
350 }
351 if (nbd < *l1) {
352 goto L128;
353 }
354 i__1 = ipph;
355 for (j = 2; j <= i__1; ++j) {
356 jc = ipp2 - j;
357 i__2 = *l1;
358 for (k = 1; k <= i__2; ++k) {
359 i__3 = *ido;
360 for (i__ = 3; i__ <= i__3; i__ += 2) {
361 m2 = m2s;
362 i__5 = m1d;
363 i__4 = *im1;
364 for (m1 = 1; i__4 < 0 ? m1 >= i__5 : m1 <= i__5; m1 += i__4) {
365 m2 += *im2;
366 ch[m2 + (i__ - 1 + (k + j * ch_dim3) * ch_dim2) * ch_dim1]
367 = c1[m1 + (i__ - 1 + (k + j * c1_dim3) * c1_dim2)
368 * c1_dim1] - c1[m1 + (i__ + (k + jc * c1_dim3) *
369 c1_dim2) * c1_dim1];
370 ch[m2 + (i__ - 1 + (k + jc * ch_dim3) * ch_dim2) *
371 ch_dim1] = c1[m1 + (i__ - 1 + (k + j * c1_dim3) *
372 c1_dim2) * c1_dim1] + c1[m1 + (i__ + (k + jc *
373 c1_dim3) * c1_dim2) * c1_dim1];
374 ch[m2 + (i__ + (k + j * ch_dim3) * ch_dim2) * ch_dim1] =
375 c1[m1 + (i__ + (k + j * c1_dim3) * c1_dim2) *
376 c1_dim1] + c1[m1 + (i__ - 1 + (k + jc * c1_dim3) *
377 c1_dim2) * c1_dim1];
378 ch[m2 + (i__ + (k + jc * ch_dim3) * ch_dim2) * ch_dim1] =
379 c1[m1 + (i__ + (k + j * c1_dim3) * c1_dim2) *
380 c1_dim1] - c1[m1 + (i__ - 1 + (k + jc * c1_dim3) *
381 c1_dim2) * c1_dim1];
382/* L1025: */
383 }
384/* L125: */
385 }
386/* L126: */
387 }
388/* L127: */
389 }
390 goto L132;
391L128:
392 i__1 = ipph;
393 for (j = 2; j <= i__1; ++j) {
394 jc = ipp2 - j;
395 i__2 = *ido;
396 for (i__ = 3; i__ <= i__2; i__ += 2) {
397 i__3 = *l1;
398 for (k = 1; k <= i__3; ++k) {
399 m2 = m2s;
400 i__4 = m1d;
401 i__5 = *im1;
402 for (m1 = 1; i__5 < 0 ? m1 >= i__4 : m1 <= i__4; m1 += i__5) {
403 m2 += *im2;
404 ch[m2 + (i__ - 1 + (k + j * ch_dim3) * ch_dim2) * ch_dim1]
405 = c1[m1 + (i__ - 1 + (k + j * c1_dim3) * c1_dim2)
406 * c1_dim1] - c1[m1 + (i__ + (k + jc * c1_dim3) *
407 c1_dim2) * c1_dim1];
408 ch[m2 + (i__ - 1 + (k + jc * ch_dim3) * ch_dim2) *
409 ch_dim1] = c1[m1 + (i__ - 1 + (k + j * c1_dim3) *
410 c1_dim2) * c1_dim1] + c1[m1 + (i__ + (k + jc *
411 c1_dim3) * c1_dim2) * c1_dim1];
412 ch[m2 + (i__ + (k + j * ch_dim3) * ch_dim2) * ch_dim1] =
413 c1[m1 + (i__ + (k + j * c1_dim3) * c1_dim2) *
414 c1_dim1] + c1[m1 + (i__ - 1 + (k + jc * c1_dim3) *
415 c1_dim2) * c1_dim1];
416 ch[m2 + (i__ + (k + jc * ch_dim3) * ch_dim2) * ch_dim1] =
417 c1[m1 + (i__ + (k + j * c1_dim3) * c1_dim2) *
418 c1_dim1] - c1[m1 + (i__ - 1 + (k + jc * c1_dim3) *
419 c1_dim2) * c1_dim1];
420/* L1029: */
421 }
422/* L129: */
423 }
424/* L130: */
425 }
426/* L131: */
427 }
428L132:
429 if (*ido == 1) {
430 return 0;
431 }
432 i__1 = *idl1;
433 for (ik = 1; ik <= i__1; ++ik) {
434 m2 = m2s;
435 i__2 = m1d;
436 i__3 = *im1;
437 for (m1 = 1; i__3 < 0 ? m1 >= i__2 : m1 <= i__2; m1 += i__3) {
438 m2 += *im2;
439 c2[m1 + (ik + c2_dim2) * c2_dim1] = ch2[m2 + (ik + ch2_dim2) *
440 ch2_dim1];
441/* L1033: */
442 }
443/* L133: */
444 }
445 i__1 = *ip;
446 for (j = 2; j <= i__1; ++j) {
447 i__3 = *l1;
448 for (k = 1; k <= i__3; ++k) {
449 m2 = m2s;
450 i__2 = m1d;
451 i__5 = *im1;
452 for (m1 = 1; i__5 < 0 ? m1 >= i__2 : m1 <= i__2; m1 += i__5) {
453 m2 += *im2;
454 c1[m1 + ((k + j * c1_dim3) * c1_dim2 + 1) * c1_dim1] = ch[m2
455 + ((k + j * ch_dim3) * ch_dim2 + 1) * ch_dim1];
456/* L1034: */
457 }
458/* L134: */
459 }
460/* L135: */
461 }
462 if (nbd > *l1) {
463 goto L139;
464 }
465 is = -(*ido);
466 i__1 = *ip;
467 for (j = 2; j <= i__1; ++j) {
468 is += *ido;
469 idij = is;
470 i__3 = *ido;
471 for (i__ = 3; i__ <= i__3; i__ += 2) {
472 idij += 2;
473 i__5 = *l1;
474 for (k = 1; k <= i__5; ++k) {
475 m2 = m2s;
476 i__2 = m1d;
477 i__4 = *im1;
478 for (m1 = 1; i__4 < 0 ? m1 >= i__2 : m1 <= i__2; m1 += i__4) {
479 m2 += *im2;
480 c1[m1 + (i__ - 1 + (k + j * c1_dim3) * c1_dim2) * c1_dim1]
481 = wa[idij - 1] * ch[m2 + (i__ - 1 + (k + j *
482 ch_dim3) * ch_dim2) * ch_dim1] - wa[idij] * ch[m2
483 + (i__ + (k + j * ch_dim3) * ch_dim2) * ch_dim1];
484 c1[m1 + (i__ + (k + j * c1_dim3) * c1_dim2) * c1_dim1] =
485 wa[idij - 1] * ch[m2 + (i__ + (k + j * ch_dim3) *
486 ch_dim2) * ch_dim1] + wa[idij] * ch[m2 + (i__ - 1
487 + (k + j * ch_dim3) * ch_dim2) * ch_dim1];
488/* L1036: */
489 }
490/* L136: */
491 }
492/* L137: */
493 }
494/* L138: */
495 }
496 goto L143;
497L139:
498 is = -(*ido);
499 i__1 = *ip;
500 for (j = 2; j <= i__1; ++j) {
501 is += *ido;
502 i__3 = *l1;
503 for (k = 1; k <= i__3; ++k) {
504 idij = is;
505 i__5 = *ido;
506 for (i__ = 3; i__ <= i__5; i__ += 2) {
507 idij += 2;
508 m2 = m2s;
509 i__4 = m1d;
510 i__2 = *im1;
511 for (m1 = 1; i__2 < 0 ? m1 >= i__4 : m1 <= i__4; m1 += i__2) {
512 m2 += *im2;
513 c1[m1 + (i__ - 1 + (k + j * c1_dim3) * c1_dim2) * c1_dim1]
514 = wa[idij - 1] * ch[m2 + (i__ - 1 + (k + j *
515 ch_dim3) * ch_dim2) * ch_dim1] - wa[idij] * ch[m2
516 + (i__ + (k + j * ch_dim3) * ch_dim2) * ch_dim1];
517 c1[m1 + (i__ + (k + j * c1_dim3) * c1_dim2) * c1_dim1] =
518 wa[idij - 1] * ch[m2 + (i__ + (k + j * ch_dim3) *
519 ch_dim2) * ch_dim1] + wa[idij] * ch[m2 + (i__ - 1
520 + (k + j * ch_dim3) * ch_dim2) * ch_dim1];
521/* L1040: */
522 }
523/* L140: */
524 }
525/* L141: */
526 }
527/* L142: */
528 }
529L143:
530 return 0;
531} /* mradbg_ */
532