PLaSK library
Loading...
Searching...
No Matches
mcstf1.c
Go to the documentation of this file.
1/* mcstf1.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_n5 = -5;
18
19/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
20/* * * */
21/* * copyright (c) 2011 by UCAR * */
22/* * * */
23/* * University Corporation for Atmospheric Research * */
24/* * * */
25/* * all rights reserved * */
26/* * * */
27/* * FFTPACK version 5.1 * */
28/* * * */
29/* * A Fortran Package of Fast Fourier * */
30/* * * */
31/* * Subroutines and Example Programs * */
32/* * * */
33/* * by * */
34/* * * */
35/* * Paul Swarztrauber and Dick Valent * */
36/* * * */
37/* * of * */
38/* * * */
39/* * the National Center for Atmospheric Research * */
40/* * * */
41/* * Boulder, Colorado (80307) U.S.A. * */
42/* * * */
43/* * which is sponsored by * */
44/* * * */
45/* * the National Science Foundation * */
46/* * * */
47/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
48
49/* Subroutine */ int mcstf1_(integer *lot, integer *jump, integer *n, integer
50 *inc, doublereal *x, doublereal *wsave, doublereal *dsum, doublereal *
51 work, integer *ier)
52{
53 /* System generated locals */
54 integer x_dim1, x_offset, i__1, i__2, i__3;
55
56 /* Builtin functions */
57 double log(doublereal);
58
59 /* Local variables */
60 integer i__, k, m, m1;
61 doublereal t1, t2;
62 integer kc, lj;
63 doublereal xi;
64 integer nm1, np1;
65 doublereal x1h;
66 integer ns2;
67 doublereal tx2, x1p3;
68 integer ier1;
69 doublereal snm1;
70 integer modn, lenx, lnwk, lnsv;
71 extern /* Subroutine */ int rfftmf_(integer *, integer *, integer *,
72 integer *, doublereal *, integer *, doublereal *, integer *,
73 doublereal *, integer *, integer *), xerfft_(char *, integer *,
74 ftnlen);
75
76 /* Parameter adjustments */
77 x_dim1 = *inc;
78 x_offset = 1 + x_dim1;
79 x -= x_offset;
80 --wsave;
81 --dsum;
82
83 /* Function Body */
84 *ier = 0;
85 nm1 = *n - 1;
86 np1 = *n + 1;
87 ns2 = *n / 2;
88 lj = (*lot - 1) * *jump + 1;
89 if ((i__1 = *n - 2) < 0) {
90 goto L200;
91 } else if (i__1 == 0) {
92 goto L101;
93 } else {
94 goto L102;
95 }
96L101:
97 i__1 = lj;
98 i__2 = *jump;
99 for (m = 1; i__2 < 0 ? m >= i__1 : m <= i__1; m += i__2) {
100 x1h = x[m + x_dim1] + x[m + (x_dim1 << 1)];
101 x[m + (x_dim1 << 1)] = (x[m + x_dim1] - x[m + (x_dim1 << 1)]) * .5;
102 x[m + x_dim1] = x1h * .5;
103/* L111: */
104 }
105 goto L200;
106L102:
107 if (*n > 3) {
108 goto L103;
109 }
110 i__2 = lj;
111 i__1 = *jump;
112 for (m = 1; i__1 < 0 ? m >= i__2 : m <= i__2; m += i__1) {
113 x1p3 = x[m + x_dim1] + x[m + x_dim1 * 3];
114 tx2 = x[m + (x_dim1 << 1)] + x[m + (x_dim1 << 1)];
115 x[m + (x_dim1 << 1)] = (x[m + x_dim1] - x[m + x_dim1 * 3]) * .5;
116 x[m + x_dim1] = (x1p3 + tx2) * .25;
117 x[m + x_dim1 * 3] = (x1p3 - tx2) * .25;
118/* L112: */
119 }
120 goto L200;
121L103:
122 m1 = 0;
123 i__1 = lj;
124 i__2 = *jump;
125 for (m = 1; i__2 < 0 ? m >= i__1 : m <= i__1; m += i__2) {
126 ++m1;
127 dsum[m1] = x[m + x_dim1] - x[m + *n * x_dim1];
128 x[m + x_dim1] += x[m + *n * x_dim1];
129/* L113: */
130 }
131 i__2 = ns2;
132 for (k = 2; k <= i__2; ++k) {
133 m1 = 0;
134 i__1 = lj;
135 i__3 = *jump;
136 for (m = 1; i__3 < 0 ? m >= i__1 : m <= i__1; m += i__3) {
137 ++m1;
138 kc = np1 - k;
139 t1 = x[m + k * x_dim1] + x[m + kc * x_dim1];
140 t2 = x[m + k * x_dim1] - x[m + kc * x_dim1];
141 dsum[m1] += wsave[kc] * t2;
142 t2 = wsave[k] * t2;
143 x[m + k * x_dim1] = t1 - t2;
144 x[m + kc * x_dim1] = t1 + t2;
145/* L114: */
146 }
147/* L104: */
148 }
149 modn = *n % 2;
150 if (modn == 0) {
151 goto L124;
152 }
153 i__2 = lj;
154 i__3 = *jump;
155 for (m = 1; i__3 < 0 ? m >= i__2 : m <= i__2; m += i__3) {
156 x[m + (ns2 + 1) * x_dim1] += x[m + (ns2 + 1) * x_dim1];
157/* L123: */
158 }
159L124:
160 lenx = (*lot - 1) * *jump + *inc * (nm1 - 1) + 1;
161 lnsv = nm1 + (integer) (log((doublereal) nm1) / log(2.)) + 4;
162 lnwk = *lot * nm1;
163
164 rfftmf_(lot, jump, &nm1, inc, &x[x_offset], &lenx, &wsave[*n + 1], &lnsv,
165 work, &lnwk, &ier1);
166 if (ier1 != 0) {
167 *ier = 20;
168 xerfft_("MCSTF1", &c_n5, (ftnlen)6);
169 goto L200;
170 }
171
172 snm1 = 1. / (doublereal) nm1;
173 i__3 = *lot;
174 for (m = 1; m <= i__3; ++m) {
175 dsum[m] = snm1 * dsum[m];
176/* L10: */
177 }
178 if (nm1 % 2 != 0) {
179 goto L30;
180 }
181 i__3 = lj;
182 i__2 = *jump;
183 for (m = 1; i__2 < 0 ? m >= i__3 : m <= i__3; m += i__2) {
184 x[m + nm1 * x_dim1] += x[m + nm1 * x_dim1];
185/* L20: */
186 }
187L30:
188 i__2 = *n;
189 for (i__ = 3; i__ <= i__2; i__ += 2) {
190 m1 = 0;
191 i__3 = lj;
192 i__1 = *jump;
193 for (m = 1; i__1 < 0 ? m >= i__3 : m <= i__3; m += i__1) {
194 ++m1;
195 xi = x[m + i__ * x_dim1] * .5;
196 x[m + i__ * x_dim1] = x[m + (i__ - 1) * x_dim1] * .5;
197 x[m + (i__ - 1) * x_dim1] = dsum[m1];
198 dsum[m1] += xi;
199/* L115: */
200 }
201/* L105: */
202 }
203 if (modn != 0) {
204 goto L117;
205 }
206 m1 = 0;
207 i__2 = lj;
208 i__1 = *jump;
209 for (m = 1; i__1 < 0 ? m >= i__2 : m <= i__2; m += i__1) {
210 ++m1;
211 x[m + *n * x_dim1] = dsum[m1];
212/* L116: */
213 }
214L117:
215 i__1 = lj;
216 i__2 = *jump;
217 for (m = 1; i__2 < 0 ? m >= i__1 : m <= i__1; m += i__2) {
218 x[m + x_dim1] *= .5;
219 x[m + *n * x_dim1] *= .5;
220/* L118: */
221 }
222
223L200:
224 return 0;
225} /* mcstf1_ */
226