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