lib/cfft.c

   1 /* from fftpkg package in cmlib and netlib -- translated by f2c and modified */
   2
   3 #include "xmath.h"
   4
   5 /*
   6   Public Routine
   7 */
   8
   9 /*
  10  * cfft computes the forward or backward complex discrete fourier transform.
  11  *
  12  * Input Parameters:
  13  *
  14  * n      The length of the complex sequence c. The method is
  15  *        more efficient when n is the product of small primes.
  16  *
  17  * c      A complex array of length n which contains the sequence
  18  *
  19  * wsave  a real work array which must be dimensioned at least 4n+15
  20  *        in the program that calls cfft.
  21  * isign  1 for transform, -1 for inverse transform.
  22  *        A call of cfft with isign = 1 followed by a call of cfft with
  23  *        isign = -1 will multiply the sequence by n.
  24  *
  25  * Output Parameters:
  26  *
  27  * c      For j=1,...,n
  28  *
  29  *             c(j)=the sum from k=1,...,n of
  30  *
  31  *                   c(k)*exp(-isign*i*(j-1)*(k-1)*2*pi/n)
  32  *
  33  *                         where i=sqrt(-1)
  34  */
  35
  36 int cfft(n, c, wsave, isign)
  37      int n;
  38      double *c, *wsave;
  39      int isign;
  40 {
  41   int iw1, iw2;
  42
  43   /* Parameter adjustments */
  44   --c;
  45   --wsave;
  46
  47   /* Function Body */
  48   if (n != 1) {
  49     iw1 = n + n + 1;
  50     iw2 = iw1 + n + n;
  51     cffti1_(&n, &wsave[iw1], &wsave[iw2]);
  52     cfft1_(&n, &c[1], &wsave[1], &wsave[iw1], &wsave[iw2], &isign);
  53   }
  54   return 0;
  55 }
  56
  57 /*
  58   Internal Routines
  59 */
  60
  61 static int cfft1_(n, c, ch, wa, ifac, isign)
  62      int *n;
  63      double *c, *ch, *wa;
  64      int *ifac;
  65      int *isign;
  66 {
  67   /* System generated locals */
  68   int i_1;
  69
  70   /* Local variables */
  71   int idot;
  72   int i, k1, l1, l2, n2, na, nf, ip, iw, ix2, ix3, ix4, nac, ido, idl1;
  73
  74   /* Parameter adjustments */
  75   --c;
  76   --ch;
  77   --wa;
  78   --ifac;
  79
  80   /* Function Body */
  81   nf = ifac[2];
  82   na = 0;
  83   l1 = 1;
  84   iw = 1;
  85   i_1 = nf;
  86   for (k1 = 1; k1 <= i_1; ++k1) {
  87     ip = ifac[k1 + 2];
  88     l2 = ip * l1;
  89     ido = *n / l2;
  90     idot = ido + ido;
  91     idl1 = idot * l1;
  92     if (ip == 4) {
  93       ix2 = iw + idot;
  94       ix3 = ix2 + idot;
  95       if (na == 0) {
  96         pass4_(&idot, &l1, &c[1], &ch[1], &wa[iw], &wa[ix2], &wa[ix3], isign);
  97       }
  98       else {
  99         pass4_(&idot, &l1, &ch[1], &c[1], &wa[iw], &wa[ix2], &wa[ix3], isign);
 100       }
 101       na = 1 - na;
 102     }
 103     else if (ip == 2) {
 104       if (na == 0) {
 105         pass2_(&idot, &l1, &c[1], &ch[1], &wa[iw], isign);
 106       }
 107       else {
 108         pass2_(&idot, &l1, &ch[1], &c[1], &wa[iw], isign);
 109       }
 110       na = 1 - na;
 111     }
 112     else if (ip == 3) {
 113       ix2 = iw + idot;
 114       if (na == 0) {
 115         pass3_(&idot, &l1, &c[1], &ch[1], &wa[iw], &wa[ix2], isign);
 116       }
 117       else {
 118         pass3_(&idot, &l1, &ch[1], &c[1], &wa[iw], &wa[ix2], isign);
 119       }
 120       na = 1 - na;
 121     }
 122     else if (ip == 5) {
 123       ix2 = iw + idot;
 124       ix3 = ix2 + idot;
 125       ix4 = ix3 + idot;
 126       if (na == 0) {
 127         pass5_(&idot, &l1, &c[1], &ch[1], &wa[iw], &wa[ix2], &wa[ix3],
 128                &wa[ix4], isign);
 129       }
 130       else {
 131         pass5_(&idot, &l1, &ch[1], &c[1], &wa[iw], &wa[ix2], &wa[ix3],
 132                &wa[ix4], isign);
 133       }
 134       na = 1 - na;
 135     }
 136     else {
 137       if (na == 0) {
 138         pass_(&nac, &idot, &ip, &l1, &idl1, &c[1], &c[1], &c[1], &ch[1],
 139               &ch[1], &wa[iw], isign);
 140       }
 141       else {
 142         pass_(&nac, &idot, &ip, &l1, &idl1, &ch[1], &ch[1], &ch[1], &c[1],
 143               &c[1], &wa[iw], isign);
 144       }
 145       if (nac != 0) {
 146         na = 1 - na;
 147       }
 148     }
 149     l1 = l2;
 150     iw += (ip - 1) * idot;
 151   }
 152   if (na != 0) {
 153     n2 = *n + *n;
 154     i_1 = n2;
 155     for (i = 1; i <= i_1; ++i) {
 156       c[i] = ch[i];
 157     }
 158   }
 159   return 0;
 160 }
 161
 162 static int cffti1_(n, wa, ifac)
 163      int *n;
 164      double *wa;
 165      int *ifac;
 166 {
 167   /* Initialized data */
 168   static int ntryh[4] = { 3,4,2,5 };
 169
 170   /* System generated locals */
 171   int i_1, i_2, i_3;
 172
 173   /* Local variables */
 174   double argh;
 175   int idot, ntry, i, j;
 176   double argld;
 177   int i1, k1, l1, l2, ib;
 178   double fi;
 179   int ld, ii, nf, ip, nl, nq, nr;
 180   double arg;
 181   int ido, ipm;
 182   double tpi;
 183
 184   /* Parameter adjustments */
 185   --wa;
 186   --ifac;
 187
 188   /* Function Body */
 189   nl = *n;
 190   nf = 0;
 191   j = 0;
 192
 193  L101:
 194   ++j;
 195   if (j - 4 <= 0) ntry = ntryh[j - 1];
 196   else ntry += 2;
 197  L104:
 198   nq = nl / ntry;
 199   nr = nl - ntry * nq;
 200   if (nr != 0) goto L101;
 201   ++nf;
 202   ifac[nf + 2] = ntry;
 203   nl = nq;
 204   if (ntry == 2 && nf != 1) {
 205     i_1 = nf;
 206     for (i = 2; i <= i_1; ++i) {
 207       ib = nf - i + 2;
 208       ifac[ib + 2] = ifac[ib + 1];
 209     }
 210     ifac[3] = 2;
 211   }
 212   if (nl != 1) goto L104;
 213
 214   ifac[1] = *n;
 215   ifac[2] = nf;
 216   tpi = 6.28318530717959;
 217   argh = tpi / (double) (*n);
 218   i = 2;
 219   l1 = 1;
 220   i_1 = nf;
 221   for (k1 = 1; k1 <= i_1; ++k1) {
 222     ip = ifac[k1 + 2];
 223     ld = 0;
 224     l2 = l1 * ip;
 225     ido = *n / l2;
 226     idot = ido + ido + 2;
 227     ipm = ip - 1;
 228     i_2 = ipm;
 229     for (j = 1; j <= i_2; ++j) {
 230       i1 = i;
 231       wa[i - 1] = 1.0;
 232       wa[i] = 0.0;
 233       ld += l1;
 234       fi = 0.0;
 235       argld = (double) ld * argh;
 236       i_3 = idot;
 237       for (ii = 4; ii <= i_3; ii += 2) {
 238         i += 2;
 239         fi += 1.0;
 240         arg = fi * argld;
 241         wa[i - 1] = cos(arg);
 242         wa[i] = sin(arg);
 243       }
 244       if (ip > 5) {
 245         wa[i1 - 1] = wa[i - 1];
 246         wa[i1] = wa[i];
 247       }
 248     }
 249     l1 = l2;
 250   }
 251   return 0;
 252 }
 253
 254 static int pass_(nac, ido, ip, l1, idl1, cc, c1, c2, ch, ch2, wa, isign)
 255      int *nac;
 256      int *ido, *ip, *l1, *idl1;
 257      double *cc, *c1, *c2, *ch, *ch2, *wa;
 258      int *isign;
 259 {
 260   /* System generated locals */
 261   int ch_dim1, ch_dim2, ch_offset, cc_dim1, cc_dim2, cc_offset, c1_dim1,
 262       c1_dim2, c1_offset, c2_dim1, c2_offset, ch2_dim1, ch2_offset,
 263       i_1, i_2, i_3;
 264
 265   /* Local variables */
 266   int idij, idlj, idot, ipph, i, j, k, l, jc, lc, ik, idj, idl, inc, idp;
 267   double wai, war;
 268   int ipp2;
 269
 270   /* Parameter adjustments */
 271   cc_dim1 = *ido;
 272   cc_dim2 = *ip;
 273   cc_offset = cc_dim1 * (cc_dim2 + 1) + 1;
 274   cc -= cc_offset;
 275   c1_dim1 = *ido;
 276   c1_dim2 = *l1;
 277   c1_offset = c1_dim1 * (c1_dim2 + 1) + 1;
 278   c1 -= c1_offset;
 279   c2_dim1 = *idl1;
 280   c2_offset = c2_dim1 + 1;
 281   c2 -= c2_offset;
 282   ch_dim1 = *ido;
 283   ch_dim2 = *l1;
 284   ch_offset = ch_dim1 * (ch_dim2 + 1) + 1;
 285   ch -= ch_offset;
 286   ch2_dim1 = *idl1;
 287   ch2_offset = ch2_dim1 + 1;
 288   ch2 -= ch2_offset;
 289   --wa;
 290
 291   /* Function Body */
 292   idot = *ido / 2;
 293   ipp2 = *ip + 2;
 294   ipph = (*ip + 1) / 2;
 295   idp = *ip * *ido;
 296
 297   if (*ido >= *l1) {
 298     i_1 = ipph;
 299     for (j = 2; j <= i_1; ++j) {
 300       jc = ipp2 - j;
 301       i_2 = *l1;
 302       for (k = 1; k <= i_2; ++k) {
 303         i_3 = *ido;
 304         for (i = 1; i <= i_3; ++i) {
 305           ch[i + (k + j * ch_dim2) * ch_dim1] =
 306             cc[i + (j + k * cc_dim2) * cc_dim1]
 307               + cc[i + (jc + k * cc_dim2) * cc_dim1];
 308           ch[i + (k + jc * ch_dim2) * ch_dim1] =
 309             cc[i + (j + k * cc_dim2) * cc_dim1]
 310               - cc[i + (jc + k * cc_dim2) * cc_dim1];
 311         }
 312       }
 313     }
 314     i_1 = *l1;
 315     for (k = 1; k <= i_1; ++k) {
 316       i_2 = *ido;
 317       for (i = 1; i <= i_2; ++i) {
 318         ch[i + (k + ch_dim2) * ch_dim1] = cc[i + (k * cc_dim2 + 1) * cc_dim1];
 319       }
 320     }
 321   }
 322   else {
 323     i_1 = ipph;
 324     for (j = 2; j <= i_1; ++j) {
 325       jc = ipp2 - j;
 326       i_2 = *ido;
 327       for (i = 1; i <= i_2; ++i) {
 328         i_3 = *l1;
 329         for (k = 1; k <= i_3; ++k) {
 330           ch[i + (k + j * ch_dim2) * ch_dim1] =
 331             cc[i + (j + k * cc_dim2) * cc_dim1]
 332               + cc[i + (jc + k * cc_dim2) * cc_dim1];
 333           ch[i + (k + jc * ch_dim2) * ch_dim1] =
 334             cc[i + (j + k * cc_dim2) * cc_dim1]
 335               - cc[i + (jc + k * cc_dim2) * cc_dim1];
 336         }
 337       }
 338     }
 339     i_1 = *ido;
 340     for (i = 1; i <= i_1; ++i) {
 341       i_2 = *l1;
 342       for (k = 1; k <= i_2; ++k) {
 343         ch[i + (k + ch_dim2) * ch_dim1] = cc[i + (k * cc_dim2 + 1) * cc_dim1];
 344       }
 345     }
 346   }
 347   idl = 2 - *ido;
 348   inc = 0;
 349   i_1 = ipph;
 350   for (l = 2; l <= i_1; ++l) {
 351     lc = ipp2 - l;
 352     idl += *ido;
 353     i_2 = *idl1;
 354     for (ik = 1; ik <= i_2; ++ik) {
 355       c2[ik + l * c2_dim1] = ch2[ik + ch2_dim1] + wa[idl - 1]
 356         * ch2[ik + (ch2_dim1 << 1)];
 357       c2[ik + lc * c2_dim1] = - *isign * wa[idl] * ch2[ik + *ip * ch2_dim1];
 358     }
 359     idlj = idl;
 360     inc += *ido;
 361     i_2 = ipph;
 362     for (j = 3; j <= i_2; ++j) {
 363       jc = ipp2 - j;
 364       idlj += inc;
 365       if (idlj > idp) {
 366         idlj -= idp;
 367       }
 368       war = wa[idlj - 1];
 369       wai = wa[idlj];
 370       i_3 = *idl1;
 371       for (ik = 1; ik <= i_3; ++ik) {
 372         c2[ik + l * c2_dim1] += war * ch2[ik + j * ch2_dim1];
 373         c2[ik + lc * c2_dim1] -= *isign * wai * ch2[ik + jc * ch2_dim1];
 374       }
 375     }
 376   }
 377   i_1 = ipph;
 378   for (j = 2; j <= i_1; ++j) {
 379     i_2 = *idl1;
 380     for (ik = 1; ik <= i_2; ++ik) {
 381       ch2[ik + ch2_dim1] += ch2[ik + j * ch2_dim1];
 382     }
 383   }
 384   i_1 = ipph;
 385   for (j = 2; j <= i_1; ++j) {
 386     jc = ipp2 - j;
 387     i_2 = *idl1;
 388     for (ik = 2; ik <= i_2; ik += 2) {
 389       ch2[ik - 1 + j * ch2_dim1] = c2[ik - 1 + j * c2_dim1]
 390         - c2[ik + jc * c2_dim1];
 391       ch2[ik - 1 + jc * ch2_dim1] = c2[ik - 1 + j * c2_dim1]
 392         + c2[ik + jc * c2_dim1];
 393       ch2[ik + j * ch2_dim1] = c2[ik + j * c2_dim1]
 394         + c2[ik - 1 + jc * c2_dim1];
 395       ch2[ik + jc * ch2_dim1] = c2[ik + j * c2_dim1]
 396         - c2[ik - 1 + jc * c2_dim1];
 397     }
 398   }
 399   *nac = 1;
 400   if (*ido != 2) {
 401     *nac = 0;
 402     i_1 = *idl1;
 403     for (ik = 1; ik <= i_1; ++ik) {
 404       c2[ik + c2_dim1] = ch2[ik + ch2_dim1];
 405     }
 406     i_1 = *ip;
 407     for (j = 2; j <= i_1; ++j) {
 408       i_2 = *l1;
 409       for (k = 1; k <= i_2; ++k) {
 410         c1[(k + j * c1_dim2) * c1_dim1 + 1] =
 411           ch[(k + j * ch_dim2) * ch_dim1 + 1];
 412         c1[(k + j * c1_dim2) * c1_dim1 + 2] =
 413           ch[(k + j * ch_dim2) * ch_dim1 + 2];
 414       }
 415     }
 416     if (idot <= *l1) {
 417       idij = 0;
 418       i_1 = *ip;
 419       for (j = 2; j <= i_1; ++j) {
 420         idij += 2;
 421         i_2 = *ido;
 422         for (i = 4; i <= i_2; i += 2) {
 423           idij += 2;
 424           i_3 = *l1;
 425           for (k = 1; k <= i_3; ++k) {
 426             c1[i - 1 + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1]
 427               * ch[i - 1 + (k + j * ch_dim2) * ch_dim1] + *isign * wa[idij]
 428                 * ch[i + (k + j * ch_dim2) * ch_dim1];
 429             c1[i + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1]
 430               * ch[i + (k + j * ch_dim2) * ch_dim1] - *isign * wa[idij]
 431                 * ch[i - 1 + (k + j * ch_dim2) * ch_dim1];
 432           }
 433         }
 434       }
 435     }
 436     else {
 437       idj = 2 - *ido;
 438       i_1 = *ip;
 439       for (j = 2; j <= i_1; ++j) {
 440         idj += *ido;
 441         i_2 = *l1;
 442         for (k = 1; k <= i_2; ++k) {
 443           idij = idj;
 444           i_3 = *ido;
 445           for (i = 4; i <= i_3; i += 2) {
 446             idij += 2;
 447             c1[i - 1 + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1]
 448               * ch[i - 1 + (k + j * ch_dim2) * ch_dim1]
 449                 + *isign * wa[idij] * ch[i + (k + j * ch_dim2) * ch_dim1];
 450             c1[i + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1]
 451               * ch[i + (k + j * ch_dim2) * ch_dim1] - *isign * wa[idij]
 452                 * ch[i - 1 + (k + j * ch_dim2) * ch_dim1];
 453           }
 454         }
 455       }
 456     }
 457   }
 458   return 0;
 459 }
 460
 461 static int pass2_(ido, l1, cc, ch, wa1, isign)
 462      int *ido, *l1;
 463      double *cc, *ch, *wa1;
 464      int *isign;
 465 {
 466   /* System generated locals */
 467   int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i_1, i_2;
 468
 469   /* Local variables */
 470   int i, k;
 471   double ti2, tr2;
 472
 473   /* Parameter adjustments */
 474   cc_dim1 = *ido;
 475   cc_offset = cc_dim1 * 3 + 1;
 476   cc -= cc_offset;
 477   ch_dim1 = *ido;
 478   ch_dim2 = *l1;
 479   ch_offset = ch_dim1 * (ch_dim2 + 1) + 1;
 480   ch -= ch_offset;
 481   --wa1;
 482
 483   /* Function Body */
 484   if (*ido <= 2) {
 485     i_1 = *l1;
 486     for (k = 1; k <= i_1; ++k) {
 487       ch[(k + ch_dim2) * ch_dim1 + 1] = cc[((k << 1) + 1) * cc_dim1 + 1]
 488         + cc[((k << 1) + 2) * cc_dim1 + 1];
 489       ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = cc[((k << 1) + 1) * cc_dim1 + 1]
 490         - cc[((k << 1) + 2) * cc_dim1 + 1];
 491       ch[(k + ch_dim2) * ch_dim1 + 2] = cc[((k << 1) + 1) * cc_dim1 + 2]
 492         + cc[((k << 1) + 2) * cc_dim1 + 2];
 493       ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = cc[((k << 1) + 1) * cc_dim1 + 2]
 494         - cc[((k << 1) + 2) * cc_dim1 + 2];
 495     }
 496   }
 497   else {
 498     i_1 = *l1;
 499     for (k = 1; k <= i_1; ++k) {
 500       i_2 = *ido;
 501       for (i = 2; i <= i_2; i += 2) {
 502         ch[i - 1 + (k + ch_dim2) * ch_dim1]
 503           = cc[i - 1 + ((k << 1) + 1) * cc_dim1]
 504             + cc[i - 1 + ((k << 1) + 2) * cc_dim1];
 505         tr2 = cc[i - 1 + ((k << 1) + 1) * cc_dim1]
 506           - cc[i - 1 + ((k << 1) + 2) * cc_dim1];
 507         ch[i + (k + ch_dim2) * ch_dim1] = cc[i + ((k << 1) + 1) * cc_dim1]
 508           + cc[i + ((k << 1) + 2) * cc_dim1];
 509         ti2 = cc[i + ((k << 1) + 1) * cc_dim1]
 510           - cc[i + ((k << 1) + 2) * cc_dim1];
 511         ch[i + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * ti2
 512           - *isign * wa1[i] * tr2;
 513         ch[i - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * tr2
 514           + *isign * wa1[i] * ti2;
 515       }
 516     }
 517   }
 518   return 0;
 519 }
 520
 521 static int pass3_(ido, l1, cc, ch, wa1, wa2, isign)
 522      int *ido, *l1;
 523      double *cc, *ch, *wa1, *wa2;
 524      int *isign;
 525 {
 526   /* System generated locals */
 527   int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i_1, i_2;
 528
 529   /* Local variables */
 530   double taui, taur;
 531   int i, k;
 532   double ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
 533
 534   /* Parameter adjustments */
 535   cc_dim1 = *ido;
 536   cc_offset = (cc_dim1 << 2) + 1;
 537   cc -= cc_offset;
 538   ch_dim1 = *ido;
 539   ch_dim2 = *l1;
 540   ch_offset = ch_dim1 * (ch_dim2 + 1) + 1;
 541   ch -= ch_offset;
 542   --wa1;
 543   --wa2;
 544
 545   /* Function Body */
 546   taur = -.5;
 547   taui = -(*isign) * .866025403784439;
 548   if (*ido == 2) {
 549     i_1 = *l1;
 550     for (k = 1; k <= i_1; ++k) {
 551       tr2 = cc[(k * 3 + 2) * cc_dim1 + 1] + cc[(k * 3 + 3) * cc_dim1 + 1];
 552       cr2 = cc[(k * 3 + 1) * cc_dim1 + 1] + taur * tr2;
 553       ch[(k + ch_dim2) * ch_dim1 + 1] = cc[(k * 3 + 1) * cc_dim1 + 1] + tr2;
 554       ti2 = cc[(k * 3 + 2) * cc_dim1 + 2] + cc[(k * 3 + 3) * cc_dim1 + 2];
 555       ci2 = cc[(k * 3 + 1) * cc_dim1 + 2] + taur * ti2;
 556       ch[(k + ch_dim2) * ch_dim1 + 2] = cc[(k * 3 + 1) * cc_dim1 + 2] + ti2;
 557       cr3 = taui * (cc[(k * 3 + 2) * cc_dim1 + 1]
 558                     - cc[(k * 3 + 3) * cc_dim1 + 1]);
 559       ci3 = taui * (cc[(k * 3 + 2) * cc_dim1 + 2]
 560                     - cc[(k * 3 + 3) * cc_dim1 + 2]);
 561       ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = cr2 - ci3;
 562       ch[(k + ch_dim2 * 3) * ch_dim1 + 1] = cr2 + ci3;
 563       ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = ci2 + cr3;
 564       ch[(k + ch_dim2 * 3) * ch_dim1 + 2] = ci2 - cr3;
 565     }
 566   }
 567   else {
 568     i_1 = *l1;
 569     for (k = 1; k <= i_1; ++k) {
 570       i_2 = *ido;
 571       for (i = 2; i <= i_2; i += 2) {
 572         tr2 = cc[i - 1 + (k * 3 + 2) * cc_dim1]
 573           + cc[i - 1 + (k * 3 + 3) * cc_dim1];
 574         cr2 = cc[i - 1 + (k * 3 + 1) * cc_dim1] + taur * tr2;
 575         ch[i - 1 + (k + ch_dim2) * ch_dim1] = cc[i - 1 + (k * 3 + 1)
 576                                                  * cc_dim1] + tr2;
 577         ti2 = cc[i + (k * 3 + 2) * cc_dim1] + cc[i + (k * 3 + 3) * cc_dim1];
 578         ci2 = cc[i + (k * 3 + 1) * cc_dim1] + taur * ti2;
 579         ch[i + (k + ch_dim2) * ch_dim1] = cc[i + (k * 3 + 1) * cc_dim1] + ti2;
 580         cr3 = taui * (cc[i - 1 + (k * 3 + 2) * cc_dim1]
 581                       - cc[i - 1 + (k * 3 + 3) * cc_dim1]);
 582         ci3 = taui * (cc[i + (k * 3 + 2) * cc_dim1]
 583                       - cc[i + (k * 3 + 3) * cc_dim1]);
 584         dr2 = cr2 - ci3;
 585         dr3 = cr2 + ci3;
 586         di2 = ci2 + cr3;
 587         di3 = ci2 - cr3;
 588         ch[i + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * di2
 589           - *isign * wa1[i] * dr2;
 590         ch[i - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * dr2
 591           + *isign * wa1[i] * di2;
 592         ch[i + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * di3
 593           - *isign * wa2[i] * dr3;
 594         ch[i - 1 + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * dr3
 595           + *isign * wa2[i] * di3;
 596       }
 597     }
 598   }
 599   return 0;
 600 }
 601
 602 static int pass4_(ido, l1, cc, ch, wa1, wa2, wa3, isign)
 603      int *ido, *l1;
 604      double *cc, *ch, *wa1, *wa2, *wa3;
 605      int *isign;
 606 {
 607   /* System generated locals */
 608   int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i_1, i_2;
 609
 610   /* Local variables */
 611   int i, k;
 612   double ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
 613
 614   /* Parameter adjustments */
 615   cc_dim1 = *ido;
 616   cc_offset = cc_dim1 * 5 + 1;
 617   cc -= cc_offset;
 618   ch_dim1 = *ido;
 619   ch_dim2 = *l1;
 620   ch_offset = ch_dim1 * (ch_dim2 + 1) + 1;
 621   ch -= ch_offset;
 622   --wa1;
 623   --wa2;
 624   --wa3;
 625
 626   /* Function Body */
 627   if (*ido == 2) {
 628     i_1 = *l1;
 629     for (k = 1; k <= i_1; ++k) {
 630       ti1 = cc[((k << 2) + 1) * cc_dim1 + 2]
 631         - cc[((k << 2) + 3) * cc_dim1 + 2];
 632       ti2 = cc[((k << 2) + 1) * cc_dim1 + 2]
 633         + cc[((k << 2) + 3) * cc_dim1 + 2];
 634       tr4 = *isign * (cc[((k << 2) + 2) * cc_dim1 + 2]
 635                       - cc[((k << 2) + 4) * cc_dim1 + 2]);
 636       ti3 = cc[((k << 2) + 2) * cc_dim1 + 2]
 637         + cc[((k << 2) + 4) * cc_dim1 + 2];
 638       tr1 = cc[((k << 2) + 1) * cc_dim1 + 1]
 639         - cc[((k << 2) + 3) * cc_dim1 + 1];
 640       tr2 = cc[((k << 2) + 1) * cc_dim1 + 1]
 641         + cc[((k << 2) + 3) * cc_dim1 + 1];
 642       ti4 = *isign * (cc[((k << 2) + 4) * cc_dim1 + 1]
 643                       - cc[((k << 2) + 2) * cc_dim1 + 1]);
 644       tr3 = cc[((k << 2) + 2) * cc_dim1 + 1]
 645         + cc[((k << 2) + 4) * cc_dim1 + 1];
 646       ch[(k + ch_dim2) * ch_dim1 + 1] = tr2 + tr3;
 647       ch[(k + ch_dim2 * 3) * ch_dim1 + 1] = tr2 - tr3;
 648       ch[(k + ch_dim2) * ch_dim1 + 2] = ti2 + ti3;
 649       ch[(k + ch_dim2 * 3) * ch_dim1 + 2] = ti2 - ti3;
 650       ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = tr1 + tr4;
 651       ch[(k + (ch_dim2 << 2)) * ch_dim1 + 1] = tr1 - tr4;
 652       ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = ti1 + ti4;
 653       ch[(k + (ch_dim2 << 2)) * ch_dim1 + 2] = ti1 - ti4;
 654     }
 655   }
 656   else {
 657     i_1 = *l1;
 658     for (k = 1; k <= i_1; ++k) {
 659       i_2 = *ido;
 660       for (i = 2; i <= i_2; i += 2) {
 661         ti1 = cc[i + ((k << 2) + 1) * cc_dim1]
 662           - cc[i + ((k << 2) + 3) * cc_dim1];
 663         ti2 = cc[i + ((k << 2) + 1) * cc_dim1]
 664           + cc[i + ((k << 2) + 3) * cc_dim1];
 665         ti3 = cc[i + ((k << 2) + 2) * cc_dim1]
 666           + cc[i + ((k << 2) + 4) * cc_dim1];
 667         tr4 = *isign * (cc[i + ((k << 2) + 2) * cc_dim1]
 668                         - cc[i + ((k << 2) + 4) * cc_dim1]);
 669         tr1 = cc[i - 1 + ((k << 2) + 1) * cc_dim1]
 670           - cc[i - 1 + ((k << 2) + 3) * cc_dim1];
 671         tr2 = cc[i - 1 + ((k << 2) + 1) * cc_dim1]
 672           + cc[i - 1 + ((k << 2) + 3) * cc_dim1];
 673         ti4 = *isign * (cc[i - 1 + ((k << 2) + 4) * cc_dim1]
 674                         - cc[i - 1 + ((k << 2) + 2) * cc_dim1]);
 675         tr3 = cc[i - 1 + ((k << 2) + 2) * cc_dim1]
 676           + cc[i - 1 + ((k << 2) + 4) * cc_dim1];
 677         ch[i - 1 + (k + ch_dim2) * ch_dim1] = tr2 + tr3;
 678         cr3 = tr2 - tr3;
 679         ch[i + (k + ch_dim2) * ch_dim1] = ti2 + ti3;
 680         ci3 = ti2 - ti3;
 681         cr2 = tr1 + tr4;
 682         cr4 = tr1 - tr4;
 683         ci2 = ti1 + ti4;
 684         ci4 = ti1 - ti4;
 685         ch[i - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * cr2
 686           + *isign * wa1[i] * ci2;
 687         ch[i + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * ci2
 688           - *isign * wa1[i] * cr2;
 689         ch[i - 1 + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * cr3
 690           + *isign * wa2[i] * ci3;
 691         ch[i + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * ci3
 692           - *isign * wa2[i] * cr3;
 693         ch[i - 1 + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i - 1] * cr4
 694           + *isign * wa3[i] * ci4;
 695         ch[i + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i - 1] * ci4
 696           - *isign * wa3[i] * cr4;
 697       }
 698     }
 699   }
 700   return 0;
 701 }
 702
 703 static int pass5_(ido, l1, cc, ch, wa1, wa2, wa3, wa4, isign)
 704      int *ido, *l1;
 705      double *cc, *ch, *wa1, *wa2, *wa3, *wa4;
 706      int *isign;
 707 {
 708   /* System generated locals */
 709   int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i_1, i_2;
 710
 711   /* Local variables */
 712   int i, k;
 713   double ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4, ti2, ti3,
 714          ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5, ti11, ti12, tr11,
 715          tr12;
 716
 717   /* Parameter adjustments */
 718   cc_dim1 = *ido;
 719   cc_offset = cc_dim1 * 6 + 1;
 720   cc -= cc_offset;
 721   ch_dim1 = *ido;
 722   ch_dim2 = *l1;
 723   ch_offset = ch_dim1 * (ch_dim2 + 1) + 1;
 724   ch -= ch_offset;
 725   --wa1;
 726   --wa2;
 727   --wa3;
 728   --wa4;
 729
 730   /* Function Body */
 731   tr11 = .309016994374947;
 732   ti11 = -(*isign) * .951056516295154;
 733   tr12 = -.809016994374947;
 734   ti12 = -(*isign) * .587785252292473;
 735   if (*ido == 2) {
 736     i_1 = *l1;
 737     for (k = 1; k <= i_1; ++k) {
 738       ti5 = cc[(k * 5 + 2) * cc_dim1 + 2] - cc[(k * 5 + 5) * cc_dim1 + 2];
 739       ti2 = cc[(k * 5 + 2) * cc_dim1 + 2] + cc[(k * 5 + 5) * cc_dim1 + 2];
 740       ti4 = cc[(k * 5 + 3) * cc_dim1 + 2] - cc[(k * 5 + 4) * cc_dim1 + 2];
 741       ti3 = cc[(k * 5 + 3) * cc_dim1 + 2] + cc[(k * 5 + 4) * cc_dim1 + 2];
 742       tr5 = cc[(k * 5 + 2) * cc_dim1 + 1] - cc[(k * 5 + 5) * cc_dim1 + 1];
 743       tr2 = cc[(k * 5 + 2) * cc_dim1 + 1] + cc[(k * 5 + 5) * cc_dim1 + 1];
 744       tr4 = cc[(k * 5 + 3) * cc_dim1 + 1] - cc[(k * 5 + 4) * cc_dim1 + 1];
 745       tr3 = cc[(k * 5 + 3) * cc_dim1 + 1] + cc[(k * 5 + 4) * cc_dim1 + 1];
 746       ch[(k + ch_dim2) * ch_dim1 + 1] = cc[(k * 5 + 1) * cc_dim1 + 1]
 747         + tr2 + tr3;
 748       ch[(k + ch_dim2) * ch_dim1 + 2] = cc[(k * 5 + 1) * cc_dim1 + 2]
 749         + ti2 + ti3;
 750       cr2 = cc[(k * 5 + 1) * cc_dim1 + 1] + tr11 * tr2 + tr12 * tr3;
 751       ci2 = cc[(k * 5 + 1) * cc_dim1 + 2] + tr11 * ti2 + tr12 * ti3;
 752       cr3 = cc[(k * 5 + 1) * cc_dim1 + 1] + tr12 * tr2 + tr11 * tr3;
 753       ci3 = cc[(k * 5 + 1) * cc_dim1 + 2] + tr12 * ti2 + tr11 * ti3;
 754       cr5 = ti11 * tr5 + ti12 * tr4;
 755       ci5 = ti11 * ti5 + ti12 * ti4;
 756       cr4 = ti12 * tr5 - ti11 * tr4;
 757       ci4 = ti12 * ti5 - ti11 * ti4;
 758       ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = cr2 - ci5;
 759       ch[(k + ch_dim2 * 5) * ch_dim1 + 1] = cr2 + ci5;
 760       ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = ci2 + cr5;
 761       ch[(k + ch_dim2 * 3) * ch_dim1 + 2] = ci3 + cr4;
 762       ch[(k + ch_dim2 * 3) * ch_dim1 + 1] = cr3 - ci4;
 763       ch[(k + (ch_dim2 << 2)) * ch_dim1 + 1] = cr3 + ci4;
 764       ch[(k + (ch_dim2 << 2)) * ch_dim1 + 2] = ci3 - cr4;
 765       ch[(k + ch_dim2 * 5) * ch_dim1 + 2] = ci2 - cr5;
 766     }
 767   }
 768   else {
 769     i_1 = *l1;
 770     for (k = 1; k <= i_1; ++k) {
 771       i_2 = *ido;
 772       for (i = 2; i <= i_2; i += 2) {
 773         ti5 = cc[i + (k * 5 + 2) * cc_dim1] - cc[i + (k * 5 + 5) * cc_dim1];
 774         ti2 = cc[i + (k * 5 + 2) * cc_dim1] + cc[i + (k * 5 + 5) * cc_dim1];
 775         ti4 = cc[i + (k * 5 + 3) * cc_dim1] - cc[i + (k * 5 + 4) * cc_dim1];
 776         ti3 = cc[i + (k * 5 + 3) * cc_dim1] + cc[i + (k * 5 + 4) * cc_dim1];
 777         tr5 = cc[i - 1 + (k * 5 + 2) * cc_dim1]
 778           - cc[i - 1 + (k * 5 + 5) * cc_dim1];
 779         tr2 = cc[i - 1 + (k * 5 + 2) * cc_dim1]
 780           + cc[i - 1 + (k * 5 + 5) * cc_dim1];
 781         tr4 = cc[i - 1 + (k * 5 + 3) * cc_dim1]
 782           - cc[i - 1 + (k * 5 + 4) * cc_dim1];
 783         tr3 = cc[i - 1 + (k * 5 + 3) * cc_dim1]
 784           + cc[i - 1 + (k * 5 + 4) * cc_dim1];
 785         ch[i - 1 + (k + ch_dim2) * ch_dim1] = cc[i - 1 + (k * 5 + 1)
 786                                                  * cc_dim1] + tr2 + tr3;
 787         ch[i + (k + ch_dim2) * ch_dim1] = cc[i + (k * 5 + 1)
 788                                              * cc_dim1] + ti2 + ti3;
 789         cr2 = cc[i - 1 + (k * 5 + 1) * cc_dim1] + tr11 * tr2 + tr12 * tr3;
 790         ci2 = cc[i + (k * 5 + 1) * cc_dim1] + tr11 * ti2 + tr12 * ti3;
 791
 792         cr3 = cc[i - 1 + (k * 5 + 1) * cc_dim1] + tr12 * tr2 + tr11 * tr3;
 793         ci3 = cc[i + (k * 5 + 1) * cc_dim1] + tr12 * ti2 + tr11 * ti3;
 794
 795         cr5 = ti11 * tr5 + ti12 * tr4;
 796         ci5 = ti11 * ti5 + ti12 * ti4;
 797         cr4 = ti12 * tr5 - ti11 * tr4;
 798         ci4 = ti12 * ti5 - ti11 * ti4;
 799         dr3 = cr3 - ci4;
 800         dr4 = cr3 + ci4;
 801         di3 = ci3 + cr4;
 802         di4 = ci3 - cr4;
 803         dr5 = cr2 + ci5;
 804         dr2 = cr2 - ci5;
 805         di5 = ci2 - cr5;
 806         di2 = ci2 + cr5;
 807         ch[i - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * dr2
 808           + *isign * wa1[i] * di2;
 809         ch[i + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * di2
 810           - *isign * wa1[i] * dr2;
 811         ch[i - 1 + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * dr3
 812           + *isign * wa2[i] * di3;
 813         ch[i + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * di3
 814           - *isign * wa2[i] * dr3;
 815         ch[i - 1 + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i - 1] * dr4
 816           + *isign * wa3[i] * di4;
 817         ch[i + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i - 1] * di4
 818           - *isign * wa3[i] * dr4;
 819         ch[i - 1 + (k + ch_dim2 * 5) * ch_dim1] = wa4[i - 1] * dr5
 820           + *isign * wa4[i] * di5;
 821         ch[i + (k + ch_dim2 * 5) * ch_dim1] = wa4[i - 1] * di5
 822           - *isign * wa4[i] * dr5;
 823       }
 824     }
 825   }
 826   return 0;
 827 }