PLaSK library
Loading...
Searching...
No Matches
msntb1.c
Go to the documentation of this file.
1/* msntb1.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 msntb1_(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 fnp1s4, xhold, srt3s2;
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 L200;
85 } else if (i__1 == 0) {
86 goto L102;
87 } else {
88 goto L103;
89 }
90L102:
91 srt3s2 = sqrt(3.) / 2.;
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 = srt3s2 * (x[m + x_dim1] + x[m + (x_dim1 << 1)]);
96 x[m + (x_dim1 << 1)] = srt3s2 * (x[m + x_dim1] - x[m + (x_dim1 << 1)])
97 ;
98 x[m + x_dim1] = xhold;
99/* L112: */
100 }
101 goto L200;
102L103:
103 np1 = *n + 1;
104 ns2 = *n / 2;
105 i__2 = ns2;
106 for (k = 1; k <= i__2; ++k) {
107 kc = np1 - k;
108 m1 = 0;
109 i__1 = lj;
110 i__3 = *jump;
111 for (m = 1; i__3 < 0 ? m >= i__1 : m <= i__1; m += i__3) {
112 ++m1;
113 t1 = x[m + k * x_dim1] - x[m + kc * x_dim1];
114 t2 = wsave[k] * (x[m + k * x_dim1] + x[m + kc * x_dim1]);
115 xh[m1 + (k + 1) * xh_dim1] = t1 + t2;
116 xh[m1 + (kc + 1) * xh_dim1] = t2 - t1;
117/* L114: */
118 }
119/* L104: */
120 }
121 modn = *n % 2;
122 if (modn == 0) {
123 goto L124;
124 }
125 m1 = 0;
126 i__2 = lj;
127 i__3 = *jump;
128 for (m = 1; i__3 < 0 ? m >= i__2 : m <= i__2; m += i__3) {
129 ++m1;
130 xh[m1 + (ns2 + 2) * xh_dim1] = x[m + (ns2 + 1) * x_dim1] * 4.;
131/* L123: */
132 }
133L124:
134 i__3 = *lot;
135 for (m = 1; m <= i__3; ++m) {
136 xh[m + xh_dim1] = 0.;
137/* L127: */
138 }
139 lnxh = *lot - 1 + *lot * (np1 - 1) + 1;
140 lnsv = np1 + (integer) (log((doublereal) np1) / log(2.)) + 4;
141 lnwk = *lot * np1;
142
143 rfftmf_(lot, &c__1, &np1, lot, &xh[xh_offset], &lnxh, &wsave[ns2 + 1], &
144 lnsv, work, &lnwk, &ier1);
145 if (ier1 != 0) {
146 *ier = 20;
147 xerfft_("MSNTB1", &c_n5, (ftnlen)6);
148 goto L200;
149 }
150
151 if (np1 % 2 != 0) {
152 goto L30;
153 }
154 i__3 = *lot;
155 for (m = 1; m <= i__3; ++m) {
156 xh[m + np1 * xh_dim1] += xh[m + np1 * xh_dim1];
157/* L20: */
158 }
159L30:
160 fnp1s4 = (doublereal) np1 / 4.;
161 m1 = 0;
162 i__3 = lj;
163 i__2 = *jump;
164 for (m = 1; i__2 < 0 ? m >= i__3 : m <= i__3; m += i__2) {
165 ++m1;
166 x[m + x_dim1] = fnp1s4 * xh[m1 + xh_dim1];
167 dsum[m1] = x[m + x_dim1];
168/* L125: */
169 }
170 i__2 = *n;
171 for (i__ = 3; i__ <= i__2; i__ += 2) {
172 m1 = 0;
173 i__3 = lj;
174 i__1 = *jump;
175 for (m = 1; i__1 < 0 ? m >= i__3 : m <= i__3; m += i__1) {
176 ++m1;
177 x[m + (i__ - 1) * x_dim1] = fnp1s4 * xh[m1 + i__ * xh_dim1];
178 dsum[m1] += fnp1s4 * xh[m1 + (i__ - 1) * xh_dim1];
179 x[m + i__ * x_dim1] = dsum[m1];
180/* L115: */
181 }
182/* L105: */
183 }
184 if (modn != 0) {
185 goto L200;
186 }
187 m1 = 0;
188 i__2 = lj;
189 i__1 = *jump;
190 for (m = 1; i__1 < 0 ? m >= i__2 : m <= i__2; m += i__1) {
191 ++m1;
192 x[m + *n * x_dim1] = fnp1s4 * xh[m1 + (*n + 1) * xh_dim1];
193/* L116: */
194 }
195
196L200:
197 return 0;
198} /* msntb1_ */
199