PLaSK library
Loading...
Searching...
No Matches
mcstb1.c
Go to the documentation of this file.
1/* mcstb1.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 mcstb1_(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, x2;
62 integer kc, lj;
63 doublereal xi;
64 integer nm1, np1;
65 doublereal x1h;
66 integer ns2;
67 doublereal x1p3;
68 integer ier1, modn, lenx, lnwk, lnsv;
69 doublereal fnm1s2, fnm1s4;
70 extern /* Subroutine */ int rfftmf_(integer *, integer *, integer *,
71 integer *, doublereal *, integer *, doublereal *, integer *,
72 doublereal *, integer *, integer *), xerfft_(char *, integer *,
73 ftnlen);
74
75 /* Parameter adjustments */
76 x_dim1 = *inc;
77 x_offset = 1 + x_dim1;
78 x -= x_offset;
79 --wsave;
80 --dsum;
81
82 /* Function Body */
83 *ier = 0;
84 nm1 = *n - 1;
85 np1 = *n + 1;
86 ns2 = *n / 2;
87 lj = (*lot - 1) * *jump + 1;
88 if ((i__1 = *n - 2) < 0) {
89 goto L106;
90 } else if (i__1 == 0) {
91 goto L101;
92 } else {
93 goto L102;
94 }
95L101:
96 i__1 = lj;
97 i__2 = *jump;
98 for (m = 1; i__2 < 0 ? m >= i__1 : m <= i__1; m += i__2) {
99 x1h = x[m + x_dim1] + x[m + (x_dim1 << 1)];
100 x[m + (x_dim1 << 1)] = x[m + x_dim1] - x[m + (x_dim1 << 1)];
101 x[m + x_dim1] = x1h;
102/* L111: */
103 }
104 return 0;
105L102:
106 if (*n > 3) {
107 goto L103;
108 }
109 i__2 = lj;
110 i__1 = *jump;
111 for (m = 1; i__1 < 0 ? m >= i__2 : m <= i__2; m += i__1) {
112 x1p3 = x[m + x_dim1] + x[m + x_dim1 * 3];
113 x2 = x[m + (x_dim1 << 1)];
114 x[m + (x_dim1 << 1)] = x[m + x_dim1] - x[m + x_dim1 * 3];
115 x[m + x_dim1] = x1p3 + x2;
116 x[m + x_dim1 * 3] = x1p3 - x2;
117/* L112: */
118 }
119 return 0;
120L103:
121 i__1 = lj;
122 i__2 = *jump;
123 for (m = 1; i__2 < 0 ? m >= i__1 : m <= i__1; m += i__2) {
124 x[m + x_dim1] += x[m + x_dim1];
125 x[m + *n * x_dim1] += x[m + *n * x_dim1];
126/* L118: */
127 }
128 m1 = 0;
129 i__2 = lj;
130 i__1 = *jump;
131 for (m = 1; i__1 < 0 ? m >= i__2 : m <= i__2; m += i__1) {
132 ++m1;
133 dsum[m1] = x[m + x_dim1] - x[m + *n * x_dim1];
134 x[m + x_dim1] += x[m + *n * x_dim1];
135/* L113: */
136 }
137 i__1 = ns2;
138 for (k = 2; k <= i__1; ++k) {
139 m1 = 0;
140 i__2 = lj;
141 i__3 = *jump;
142 for (m = 1; i__3 < 0 ? m >= i__2 : m <= i__2; m += i__3) {
143 ++m1;
144 kc = np1 - k;
145 t1 = x[m + k * x_dim1] + x[m + kc * x_dim1];
146 t2 = x[m + k * x_dim1] - x[m + kc * x_dim1];
147 dsum[m1] += wsave[kc] * t2;
148 t2 = wsave[k] * t2;
149 x[m + k * x_dim1] = t1 - t2;
150 x[m + kc * x_dim1] = t1 + t2;
151/* L114: */
152 }
153/* L104: */
154 }
155 modn = *n % 2;
156 if (modn == 0) {
157 goto L124;
158 }
159 i__1 = lj;
160 i__3 = *jump;
161 for (m = 1; i__3 < 0 ? m >= i__1 : m <= i__1; m += i__3) {
162 x[m + (ns2 + 1) * x_dim1] += x[m + (ns2 + 1) * x_dim1];
163/* L123: */
164 }
165L124:
166 lenx = (*lot - 1) * *jump + *inc * (nm1 - 1) + 1;
167 lnsv = nm1 + (integer) (log((doublereal) nm1) / log(2.)) + 4;
168 lnwk = *lot * nm1;
169
170 rfftmf_(lot, jump, &nm1, inc, &x[x_offset], &lenx, &wsave[*n + 1], &lnsv,
171 work, &lnwk, &ier1);
172 if (ier1 != 0) {
173 *ier = 20;
174 xerfft_("MCSTB1", &c_n5, (ftnlen)6);
175 goto L106;
176 }
177
178 fnm1s2 = (doublereal) nm1 / 2.;
179 m1 = 0;
180 i__3 = lj;
181 i__1 = *jump;
182 for (m = 1; i__1 < 0 ? m >= i__3 : m <= i__3; m += i__1) {
183 ++m1;
184 dsum[m1] *= .5;
185 x[m + x_dim1] = fnm1s2 * x[m + x_dim1];
186/* L10: */
187 }
188 if (nm1 % 2 != 0) {
189 goto L30;
190 }
191 i__1 = lj;
192 i__3 = *jump;
193 for (m = 1; i__3 < 0 ? m >= i__1 : m <= i__1; m += i__3) {
194 x[m + nm1 * x_dim1] += x[m + nm1 * x_dim1];
195/* L20: */
196 }
197L30:
198 fnm1s4 = (doublereal) nm1 / 4.;
199 i__3 = *n;
200 for (i__ = 3; i__ <= i__3; i__ += 2) {
201 m1 = 0;
202 i__1 = lj;
203 i__2 = *jump;
204 for (m = 1; i__2 < 0 ? m >= i__1 : m <= i__1; m += i__2) {
205 ++m1;
206 xi = fnm1s4 * x[m + i__ * x_dim1];
207 x[m + i__ * x_dim1] = fnm1s4 * x[m + (i__ - 1) * x_dim1];
208 x[m + (i__ - 1) * x_dim1] = dsum[m1];
209 dsum[m1] += xi;
210/* L115: */
211 }
212/* L105: */
213 }
214 if (modn != 0) {
215 return 0;
216 }
217 m1 = 0;
218 i__3 = lj;
219 i__2 = *jump;
220 for (m = 1; i__2 < 0 ? m >= i__3 : m <= i__3; m += i__2) {
221 ++m1;
222 x[m + *n * x_dim1] = dsum[m1];
223/* L116: */
224 }
225L106:
226 return 0;
227} /* mcstb1_ */
228