/*
DFT++ is a density functional package developed by the research group
of Professor Tomas Arias
Copyright 1996-2003 Sohrab Ismail-Beigi
This file is part of DFT++.
DFT++ is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
DFT++ is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with DFT++; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
Please see the file CREDITS for a list of authors.
For academic users, we request that publications using results obtained with
this software reference
"New algebraic formulation of density functional calculation," by Sohrab Ismail-Beigi
and T.A. Arias, Computer Physics Communications 128:1-2, 1-45 (June 2000).
and, if using the wavelet basis, further reference
"Multiresolution analysis of electronic structure: semicardinal and wavelet bases,"
T.A. Arias, Reviews of Modern Physics 71:1, 267-311 (January 1999).
and
"Robust ab initio calculation of condensed matter: transparent convergence through
semicardinal multiresolution analysis,'' I.P. Daykov, T.A. Arias, and
Torkel D. Engeness, Physical Review Letters, 90:21, 216402 (May 2003).
For your convenience, preprints of the above articles may be obtained from
http://arXiv.org/abs/cond-mat/9909130, 9805262, and 0204411, respectively.
*/
/* convsets.f -- translated by f2c (version 20010821).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
//#include "f2c.h"
// instead of f2c use only what needed
#define FALSE_ (0)
/* Table of constant values */
static double c_b16 = 0.;
static int c__4 = 4;
static int c__3 = 3;
static int c__2 = 2;
static int c__1 = 1;
static int c__0 = 0;
static int c__6 = 6;
static int c__8 = 8;
static int c__10 = 10;
static int c__12 = 12;
static int c__14 = 14;
static int c__15 = 15;
static int c__13 = 13;
static int c__11 = 11;
static int c__9 = 9;
/* $Id: convsets.cpp,v 1.1.2.5 2003/05/29 18:54:16 ivan Exp $ */
/* $Log: convsets.cpp,v $
/* Revision 1.1.2.5 2003/05/29 18:54:16 ivan
/* Update the PRL reference
/*
/* Revision 1.1.2.4 2003/05/23 01:49:00 sahak
/* added the new header
/*
/* Revision 1.1.2.3 2003/04/03 01:50:49 ivan
/* minor cosmetic fixes to compile on balboa (Alpha) and fiber (SGI)
/* with the native compilers
/*
/* Revision 1.1.2.2 2003/04/02 23:35:14 ivan
/* Removing C-dependencies from the code. Now it needs only C++ compiler
*/
/* Revision 1.1.2.1 2003/03/13 18:56:54 muchomas */
/* Move convsets.f to fortran-source text file convsetsf for use with f2c. */
/* Put f2c result in convsets.c, so as to eliminate need for fortran compilers. */
/* Revision 1.1.2.2 2002/10/28 22:07:29 matthew */
/* Non-orthogonal Laplaciancvs add coeff_d.c! Works for orthorhombic cells, still in need of debugging fo
r other unit cells. */
/* Revision 1.1.1.1 2001/04/18 17:15:14 daykov */
/* starting project */
/* Revision 1.1 1999/09/14 19:40:30 torkel */
/* Initial revision */
/* Revision 1.12 1999/07/04 20:48:04 torkel */
/* Full O implemented. */
/* Revision 1.11 1999/07/02 07:04:05 torkel */
/* the convolution in Osame seems to work for this version. */
/* Revision 1.10 1999/07/02 03:11:48 torkel */
/* In the middle of making Osame work */
/* Revision 1.9 1999/06/23 22:20:05 torkel */
/* revised version where Obelow is dagger of Oabove for grid without */
/* >> siblings. */
/* Revision 1.8 1999/06/21 18:56:14 torkel */
/* updated version for multilevel operators */
/* Revision 1.7 1999/06/15 23:24:45 torkel */
/* debugging, trying to make Obelow work */
/* Revision 1.6 1999/06/11 12:53:50 torkel */
/* including convolutions for obelow */
/* Revision 1.4 1999/04/25 22:32:56 muchomas */
/* *** empty log message *** */
/* Revision 1.3 1999/04/25 18:35:39 muchomas */
/* Contains new routines for same-level operator convolutions */
/* Working toward more arbitrary length convolution routines */
/* Revision 1.2 1999/04/20 15:59:37 muchomas */
/* Fully documented code */
/* Revision 1.1 1999/04/19 18:12:11 muchomas */
/* Initial revision */
/* Revision 1.6 1999/04/17 22:00:09 muchomas */
/* Working dagger operators, still slow... */
/* Revision 1.5 1999/04/14 19:48:37 muchomas */
/* Working transpose transform, but doesn't take input/output arrays -- yet! */
/* Revision 1.4 1999/04/13 15:11:36 muchomas */
/* Does ghosts, making compact version, but still in larger array */
/* Revision 1.3 1999/04/13 14:32:53 muchomas */
/* Idag that computes only only ghosts, leaving sparse data distribution */
/* Revision 1.2 1999/04/11 19:53:46 muchomas */
/* Working forward/inverse xform package */
/* Revision 1.1 1999/04/10 23:04:29 muchomas */
/* Initial revision */
/* Revision 1.3 1999/04/10 22:40:41 muchomas */
/* Now using log functionality */
/* Back and forth, linear convolution of sets of points using padding */
/* for filters of extent +/- 3 (IN PLACE) */
/* out: input/output array for in place computation */
/* outm3,outm2,outm1,outp1,outp2,outp3: pointers in output array to */
/* start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, */
/* 2 forward, 3 forward */
/* area: total size of each set */
/* stride: separation between data sets (stride=&out-&outm1) */
/* n: actual number of data sets being convolved */
/* npad: size of padded region needed for "leakage" of foward filters */
/* tr: reverse-time filter coefficients */
/* tf: forward-time filter coefficients */
/* Subroutine */ int convsets3_(double *out, double *outm3,
double *outm2, double *outm1, double *outp1, double *
outp2, double *outp3, int *area, int *stride, int *n,
int *npad, double *tr, double *tf)
{
/* System generated locals */
int i__1, i__2;
/* Local variables */
int i__, j;
double t0, t1, t2, t3;
int out2;
/* First, make sure pad is zero */
/* Parameter adjustments */
--tf;
--tr;
--outp3;
--outp2;
--outp1;
--outm1;
--outm2;
--outm3;
--out;
/* Function Body */
i__1 = *npad - 1;
for (i__ = *n; i__ <= i__1; ++i__) {
out2 = i__ * *stride;
i__2 = *area + out2;
for (j = out2 + 1; j <= i__2; ++j) {
out[j] = (double)0.;
}
}
/* Reverse phase */
t0 = tr[1];
t1 = tr[2];
t2 = tr[3];
t3 = tr[4];
for (i__ = *npad - 1; i__ >= 3; --i__) {
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j] + t1 * outm1[j] + t2 * outm2[j] + t3 * outm3[
j];
}
}
out2 = *stride << 1;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j] + t1 * outm1[j] + t2 * outm2[j];
}
out2 = *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j] + t1 * outm1[j];
}
out2 = 0;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j];
}
/* Forward phase (padding now ensures room!) */
t0 = tf[1];
t1 = tf[2];
t2 = tf[3];
t3 = tf[4];
i__1 = *n - 1;
for (i__ = 0; i__ <= i__1; ++i__) {
out2 = i__ * *stride;
i__2 = *area + out2;
for (j = out2 + 1; j <= i__2; ++j) {
out[j] = t0 * out[j] + t1 * outp1[j] + t2 * outp2[j] + t3 * outp3[
j];
}
}
return 0;
} /* convsets3_ */
/* Back and forth, linear convolution of sets of points using padding */
/* for filters of extent +/- 5 (IN PLACE) */
/* out: input/output array for in place computation */
/* outm5-outp5,out: pointers in in/output array to */
/* start of data sets */
/* area: total size of each set */
/* stride: separation between data sets (stride=&out-&outm1) */
/* n: actual number of data sets being convolved */
/* npad: size of padded region needed for "leakage" of foward filters */
/* tr: reverse-time filter coefficients */
/* tf: forward-time filter coefficients */
/* Subroutine */ int convsets5_(double *out, double *outm5,
double *outm4, double *outm3, double *outm2, double *
outm1, double *outp1, double *outp2, double *outp3,
double *outp4, double *outp5, int *area, int *stride,
int *n, int *npad, double *tr, double *tf)
{
/* System generated locals */
int i__1, i__2;
/* Local variables */
int i__, j;
double t0, t1, t2, t3, t4, t5;
int out2;
/* First, make sure pad is zero */
/* Parameter adjustments */
--tf;
--tr;
--outp5;
--outp4;
--outp3;
--outp2;
--outp1;
--outm1;
--outm2;
--outm3;
--outm4;
--outm5;
--out;
/* Function Body */
i__1 = *npad - 1;
for (i__ = *n; i__ <= i__1; ++i__) {
out2 = i__ * *stride;
i__2 = *area + out2;
for (j = out2 + 1; j <= i__2; ++j) {
out[j] = (double)0.;
}
}
/* Reverse phase */
t0 = tr[1];
t1 = tr[2];
t2 = tr[3];
t3 = tr[4];
t4 = tr[5];
t5 = tr[6];
for (i__ = *npad - 1; i__ >= 5; --i__) {
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j] + t1 * outm1[j] + t2 * outm2[j] + t3 * outm3[
j] + t4 * outm4[j] + t5 * outm5[j];
}
}
out2 = *stride << 2;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j] + t1 * outm1[j] + t2 * outm2[j] + t3 * outm3[j]
+ t4 * outm4[j];
}
out2 = *stride * 3;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j] + t1 * outm1[j] + t2 * outm2[j] + t3 * outm3[j];
}
out2 = *stride << 1;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j] + t1 * outm1[j] + t2 * outm2[j];
}
out2 = *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j] + t1 * outm1[j];
}
out2 = 0;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j];
}
/* Forward phase (padding now ensures room!) */
t0 = tf[1];
t1 = tf[2];
t2 = tf[3];
t3 = tf[4];
t4 = tf[5];
t5 = tf[6];
i__1 = *n - 1;
for (i__ = 0; i__ <= i__1; ++i__) {
out2 = i__ * *stride;
i__2 = *area + out2;
for (j = out2 + 1; j <= i__2; ++j) {
out[j] = t0 * out[j] + t1 * outp1[j] + t2 * outp2[j] + t3 * outp3[
j] + t4 * outp4[j] + t5 * outp5[j];
}
}
return 0;
} /* convsets5_ */
/* Back and forth, linear convolution of sets of points using padding */
/* for filters of extent +/- 5 (IN PLACE) */
/* out: input/output array for in place computation */
/* outm5-outp5,out: pointers in in/output array to */
/* start of data sets */
/* area: total size of each set */
/* stride: separation between data sets (stride=&out-&outm1) */
/* n: actual number of data sets being convolved */
/* npad: size of padded region needed for "leakage" of foward filters */
/* tr: reverse-time filter coefficients */
/* tf: forward-time filter coefficients */
/* Subroutine */ int convsetsper5_(double *out, double *outm5,
double *outm4, double *outm3, double *outm2, double *
outm1, double *outp1, double *outp2, double *outp3,
double *outp4, double *outp5, int *area, int *stride,
int *n, int *npad, double *tr, double *tf)
{
/* System generated locals */
int i__1, i__2;
/* Local variables */
int i__, j;
double t0, t1, t2, t3, t4, t5;
int out2;
/* First, make sure pad is zero */
/* Parameter adjustments */
--tf;
--tr;
--outp5;
--outp4;
--outp3;
--outp2;
--outp1;
--outm1;
--outm2;
--outm3;
--outm4;
--outm5;
--out;
/* Function Body */
i__1 = *npad - 1;
for (i__ = *n; i__ <= i__1; ++i__) {
out2 = i__ * *stride;
i__2 = *area + out2;
for (j = out2 + 1; j <= i__2; ++j) {
out[j] = (double)0.;
}
}
/* Reverse phase */
t0 = tr[1];
t1 = tr[2];
t2 = tr[3];
t3 = tr[4];
t4 = tr[5];
t5 = tr[6];
i__ = *npad - 1;
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j - (i__ - 4) * *stride] + t1 * out[j - (i__ - 3) *
*stride] + t2 * out[j - (i__ - 2) * *stride] + t3 * out[j - (
i__ - 1) * *stride] + t4 * out[j - i__ * *stride] + t5 *
outm5[j];
}
i__ = *npad - 2;
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j - (i__ - 3) * *stride] + t1 * out[j - (i__ - 2) *
*stride] + t2 * out[j - (i__ - 1) * *stride] + t3 * out[j -
i__ * *stride] + t4 * outm4[j] + t5 * outm5[j];
}
i__ = *npad - 3;
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j - (i__ - 2) * *stride] + t1 * out[j - (i__ - 1) *
*stride] + t2 * out[j - i__ * *stride] + t3 * outm3[j] + t4 *
outm4[j] + t5 * outm5[j];
}
i__ = *npad - 4;
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j - (i__ - 1) * *stride] + t1 * out[j - i__ * *
stride] + t2 * outm2[j] + t3 * outm3[j] + t4 * outm4[j] + t5 *
outm5[j];
}
i__ = *npad - 5;
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j - i__ * *stride] + t1 * outm1[j] + t2 * outm2[j]
+ t3 * outm3[j] + t4 * outm4[j] + t5 * outm5[j];
}
for (i__ = *npad - 6; i__ >= 5; --i__) {
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j] + t1 * outm1[j] + t2 * outm2[j] + t3 * outm3[
j] + t4 * outm4[j] + t5 * outm5[j];
}
}
out2 = *stride << 2;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = out[j + (*npad - 5) * *stride];
/* & t0*out(j)+t1*outm1(j)+t2*outm2(j)+t3*outm3(j)+t4*outm4(j) */
}
out2 = *stride * 3;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = out[j + (*npad - 5) * *stride];
/* & t0*out(j)+t1*outm1(j)+t2*outm2(j)+t3*outm3(j) */
}
out2 = *stride << 1;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = out[j + (*npad - 5) * *stride];
/* & t0*out(j)+t1*outm1(j)+t2*outm2(j) */
}
out2 = *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = out[j + (*npad - 5) * *stride];
/* & t0*out(j)+t1*outm1(j) */
}
out2 = 0;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = out[j + (*npad - 5) * *stride];
/* & t0*out(j) */
}
/* Forward phase (padding now ensures room!) */
t0 = tf[1];
t1 = tf[2];
t2 = tf[3];
t3 = tf[4];
t4 = tf[5];
t5 = tf[6];
i__1 = *n - 1;
for (i__ = 0; i__ <= i__1; ++i__) {
out2 = i__ * *stride;
i__2 = *area + out2;
for (j = out2 + 1; j <= i__2; ++j) {
out[j] = t0 * out[j] + t1 * outp1[j] + t2 * outp2[j] + t3 * outp3[
j] + t4 * outp4[j] + t5 * outp5[j];
}
}
return 0;
} /* convsetsper5_ */
/* Back and forth, linear convolution of sets of points using padding */
/* for filters of extent +/- 5 (IN PLACE) */
/* out: input/output array for in place computation */
/* outm5-outp5,out: pointers in in/output array to */
/* start of data sets */
/* area: total size of each set */
/* stride: separation between data sets (stride=&out-&outm1) */
/* n: actual number of data sets being convolved */
/* npad: size of padded region needed for "leakage" of foward filters */
/* tr: reverse-time filter coefficients */
/* tf: forward-time filter coefficients */
/* Subroutine */ int convsetsio5_(double *datout, double *work,
double *workm5, double *workm4, double *workm3,
double *workm2, double *workm1, double *workp1,
double *workp2, double *workp3, double *workp4,
double *workp5, int *area, int *stride, int *n,
int *npad, int *stridel, double *tr, double *tf)
{
/* System generated locals */
int i__1, i__2;
/* Local variables */
int work2, i__, j, out2s;
double t0, t1, t2, t3, t4, t5;
int jl;
/* First, make sure pad is zero */
/* Parameter adjustments */
--tf;
--tr;
--workp5;
--workp4;
--workp3;
--workp2;
--workp1;
--workm1;
--workm2;
--workm3;
--workm4;
--workm5;
--work;
--datout;
/* Function Body */
i__1 = *npad - 1;
for (i__ = *n; i__ <= i__1; ++i__) {
work2 = i__ * *stride;
i__2 = *area + work2;
for (j = work2 + 1; j <= i__2; ++j) {
work[j] = (double)0.;
}
}
/* Reverse phase */
t0 = tr[1];
t1 = tr[2];
t2 = tr[3];
t3 = tr[4];
t4 = tr[5];
t5 = tr[6];
for (i__ = *npad - 1; i__ >= 5; --i__) {
work2 = i__ * *stride;
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
work[j] = t0 * work[j] + t1 * workm1[j] + t2 * workm2[j] + t3 *
workm3[j] + t4 * workm4[j] + t5 * workm5[j];
}
}
work2 = *stride << 2;
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
work[j] = t0 * work[j] + t1 * workm1[j] + t2 * workm2[j] + t3 *
workm3[j] + t4 * workm4[j];
}
work2 = *stride * 3;
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
work[j] = t0 * work[j] + t1 * workm1[j] + t2 * workm2[j] + t3 *
workm3[j];
}
work2 = *stride << 1;
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
work[j] = t0 * work[j] + t1 * workm1[j] + t2 * workm2[j];
}
work2 = *stride;
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
work[j] = t0 * work[j] + t1 * workm1[j];
}
work2 = 0;
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
work[j] = t0 * work[j];
}
/* Forward phase (padding now ensures room!) */
t0 = tf[1];
t1 = tf[2];
t2 = tf[3];
t3 = tf[4];
t4 = tf[5];
t5 = tf[6];
i__1 = *n - 1;
for (i__ = 0; i__ <= i__1; ++i__) {
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__2 = *area + work2;
for (j = work2 + 1; j <= i__2; ++j) {
jl = j + out2s;
datout[jl] = t0 * work[j] + t1 * workp1[j] + t2 * workp2[j] + t3 *
workp3[j] + t4 * workp4[j] + t5 * workp5[j];
}
}
return 0;
} /* convsetsio5_ */
/* Back and forth, linear convolution of sets of points using padding */
/* for filters of extent +/- 5 (IN PLACE) */
/* out: input/output array for in place computation */
/* outm5-outp5,out: pointers in in/output array to */
/* start of data sets */
/* area: total size of each set */
/* stride: separation between data sets (stride=&out-&outm1) */
/* n: actual number of data sets being convolved */
/* npad: size of padded region needed for "leakage" of foward filters */
/* tr: reverse-time filter coefficients */
/* tf: forward-time filter coefficients */
/* Subroutine */ int convsetsio5full_(double *datout, double *work,
double *workm5, double *workm4, double *workm3,
double *workm2, double *workm1, double *workp1,
double *workp2, double *workp3, double *workp4,
double *workp5, int *area, int *stride, int *n,
int *npad, int *stridel, double *filter)
{
/* System generated locals */
int i__1, i__2;
/* Local variables */
int work2, i__, j, out2s;
double t0, t1, t2, t3, t4, t5;
int jl;
/* First, make sure pad is zero */
/* Parameter adjustments */
--filter;
--workp5;
--workp4;
--workp3;
--workp2;
--workp1;
--workm1;
--workm2;
--workm3;
--workm4;
--workm5;
--work;
--datout;
/* Function Body */
i__1 = *npad - 1;
for (i__ = *n; i__ <= i__1; ++i__) {
work2 = i__ * *stride;
i__2 = *area + work2;
for (j = work2 + 1; j <= i__2; ++j) {
work[j] = (double)0.;
}
}
/* Reverse phase */
t0 = filter[1];
t1 = filter[2];
t2 = filter[3];
t3 = filter[4];
t4 = filter[5];
t5 = filter[6];
if (t0 != (double)0.) {
/* symmetric filter */
i__ = 0;
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
jl = j + out2s;
datout[jl] = t0 * work[j] + t1 * workp1[j] + t2 * workp2[j] + t3 *
workp3[j] + t4 * workp4[j] + t5 * workp5[j];
}
i__ = 1;
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
jl = j + out2s;
datout[jl] = t1 * workm1[j] + t0 * work[j] + t1 * workp1[j] + t2 *
workp2[j] + t3 * workp3[j] + t4 * workp4[j] + t5 *
workp5[j];
}
i__ = 2;
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
jl = j + out2s;
datout[jl] = t2 * workm2[j] + t1 * workm1[j] + t0 * work[j] + t1 *
workp1[j] + t2 * workp2[j] + t3 * workp3[j] + t4 *
workp4[j] + t5 * workp5[j];
}
i__ = 3;
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
jl = j + out2s;
datout[jl] = t3 * workm3[j] + t2 * workm2[j] + t1 * workm1[j] +
t0 * work[j] + t1 * workp1[j] + t2 * workp2[j] + t3 *
workp3[j] + t4 * workp4[j] + t5 * workp5[j];
}
i__ = 4;
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
jl = j + out2s;
datout[jl] = t4 * workm4[j] + t3 * workm3[j] + t2 * workm2[j] +
t1 * workm1[j] + t0 * work[j] + t1 * workp1[j] + t2 *
workp2[j] + t3 * workp3[j] + t4 * workp4[j] + t5 * workp5[
j];
}
i__1 = *n - 6;
for (i__ = 5; i__ <= i__1; ++i__) {
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__2 = *area + work2;
for (j = work2 + 1; j <= i__2; ++j) {
jl = j + out2s;
datout[jl] = t5 * workm5[j] + t4 * workm4[j] + t3 * workm3[j]
+ t2 * workm2[j] + t1 * workm1[j] + t0 * work[j] + t1
* workp1[j] + t2 * workp2[j] + t3 * workp3[j] + t4 *
workp4[j] + t5 * workp5[j];
}
}
i__ = *n - 5;
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
jl = j + out2s;
datout[jl] = t5 * workm5[j] + t4 * workm4[j] + t3 * workm3[j] +
t2 * workm2[j] + t1 * workm1[j] + t0 * work[j] + t1 *
workp1[j] + t2 * workp2[j] + t3 * workp3[j] + t4 * workp4[
j];
}
i__ = *n - 4;
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
jl = j + out2s;
datout[jl] = t5 * workm5[j] + t4 * workm4[j] + t3 * workm3[j] +
t2 * workm2[j] + t1 * workm1[j] + t0 * work[j] + t1 *
workp1[j] + t2 * workp2[j] + t3 * workp3[j];
}
i__ = *n - 3;
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
jl = j + out2s;
datout[jl] = t5 * workm5[j] + t4 * workm4[j] + t3 * workm3[j] +
t2 * workm2[j] + t1 * workm1[j] + t0 * work[j] + t1 *
workp1[j] + t2 * workp2[j];
}
i__ = *n - 2;
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
jl = j + out2s;
datout[jl] = t5 * workm5[j] + t4 * workm4[j] + t3 * workm3[j] +
t2 * workm2[j] + t1 * workm1[j] + t0 * work[j] + t1 *
workp1[j];
}
i__ = *n - 1;
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
jl = j + out2s;
datout[jl] = t5 * workm5[j] + t4 * workm4[j] + t3 * workm3[j] +
t2 * workm2[j] + t1 * workm1[j] + t0 * work[j];
}
/* end of symmetric part */
} else {
/* antisymmetric filter */
i__ = 0;
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
jl = j + out2s;
datout[jl] = t0 * work[j] + t1 * workp1[j] + t2 * workp2[j] + t3 *
workp3[j] + t4 * workp4[j] + t5 * workp5[j];
}
i__ = 1;
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
jl = j + out2s;
datout[jl] = t1 * (double)-1. * workm1[j] + t0 * work[j] + t1 *
workp1[j] + t2 * workp2[j] + t3 * workp3[j] + t4 * workp4[
j] + t5 * workp5[j];
}
i__ = 2;
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
jl = j + out2s;
datout[jl] = t2 * (double)-1. * workm2[j] - t1 * workm1[j] + t0 *
work[j] + t1 * workp1[j] + t2 * workp2[j] + t3 * workp3[j]
+ t4 * workp4[j] + t5 * workp5[j];
}
i__ = 3;
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
jl = j + out2s;
datout[jl] = t3 * (double)-1. * workm3[j] - t2 * workm2[j] - t1 *
workm1[j] + t0 * work[j] + t1 * workp1[j] + t2 * workp2[j]
+ t3 * workp3[j] + t4 * workp4[j] + t5 * workp5[j];
}
i__ = 4;
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
jl = j + out2s;
datout[jl] = t4 * (double)-1. * workm4[j] - t3 * workm3[j] - t2 *
workm2[j] - t1 * workm1[j] + t0 * work[j] + t1 * workp1[j]
+ t2 * workp2[j] + t3 * workp3[j] + t4 * workp4[j] + t5 *
workp5[j];
}
i__1 = *n - 6;
for (i__ = 5; i__ <= i__1; ++i__) {
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__2 = *area + work2;
for (j = work2 + 1; j <= i__2; ++j) {
jl = j + out2s;
datout[jl] = t5 * (double)-1. * workm5[j] - t4 * workm4[j] -
t3 * workm3[j] - t2 * workm2[j] - t1 * workm1[j] + t0
* work[j] + t1 * workp1[j] + t2 * workp2[j] + t3 *
workp3[j] + t4 * workp4[j] + t5 * workp5[j];
}
}
i__ = *n - 5;
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
jl = j + out2s;
datout[jl] = t5 * (double)-1. * workm5[j] - t4 * workm4[j] - t3 *
workm3[j] - t2 * workm2[j] - t1 * workm1[j] + t0 * work[j]
+ t1 * workp1[j] + t2 * workp2[j] + t3 * workp3[j] + t4 *
workp4[j];
}
i__ = *n - 4;
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
jl = j + out2s;
datout[jl] = t5 * (double)-1. * workm5[j] - t4 * workm4[j] - t3 *
workm3[j] - t2 * workm2[j] - t1 * workm1[j] + t0 * work[j]
+ t1 * workp1[j] + t2 * workp2[j] + t3 * workp3[j];
}
i__ = *n - 3;
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
jl = j + out2s;
datout[jl] = t5 * (double)-1. * workm5[j] - t4 * workm4[j] - t3 *
workm3[j] - t2 * workm2[j] - t1 * workm1[j] + t0 * work[j]
+ t1 * workp1[j] + t2 * workp2[j];
}
i__ = *n - 2;
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
jl = j + out2s;
datout[jl] = t5 * (double)-1. * workm5[j] - t4 * workm4[j] - t3 *
workm3[j] - t2 * workm2[j] - t1 * workm1[j] + t0 * work[j]
+ t1 * workp1[j];
}
i__ = *n - 1;
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
jl = j + out2s;
datout[jl] = t5 * (double)-1. * workm5[j] - t4 * workm4[j] - t3 *
workm3[j] - t2 * workm2[j] - t1 * workm1[j] + t0 * work[j]
;
}
/* end of antisymmetric part */
}
return 0;
} /* convsetsio5full_ */
/* Back and forth, linear convolution of sets of points using padding */
/* for filters of extent +/- 5 (IN PLACE) */
/* out: input/output array for in place computation */
/* outm5-outp5,out: pointers in in/output array to */
/* start of data sets */
/* area: total size of each set */
/* stride: separation between data sets (stride=&out-&outm1) */
/* n: actual number of data sets being convolved */
/* npad: size of padded region needed for "leakage" of foward filters */
/* tr: reverse-time filter coefficients */
/* tf: forward-time filter coefficients */
/* Subroutine */ int convsetsioppercum5_(double *datout, double *work,
double *workm5, double *workm4, double *workm3,
double *workm2, double *workm1, double *workp1,
double *workp2, double *workp3, double *workp4,
double *workp5, int *area, int *stride, int *n,
int *npad, int *stridel, double *tr, double *tf)
{
/* System generated locals */
int i__1, i__2;
/* Local variables */
int work2, i__, j, out2s;
double t0, t1, t2, t3, t4, t5;
int jl;
/* First, make sure pad is zero */
/* Parameter adjustments */
--tf;
--tr;
--workp5;
--workp4;
--workp3;
--workp2;
--workp1;
--workm1;
--workm2;
--workm3;
--workm4;
--workm5;
--work;
--datout;
/* Function Body */
i__1 = *npad - 1;
for (i__ = *n; i__ <= i__1; ++i__) {
work2 = i__ * *stride;
i__2 = *area + work2;
for (j = work2 + 1; j <= i__2; ++j) {
work[j] = (double)0.;
}
}
/* Reverse phase */
t0 = tr[1];
t1 = tr[2];
t2 = tr[3];
t3 = tr[4];
t4 = tr[5];
t5 = tr[6];
i__ = *npad - 1;
work2 = i__ * *stride;
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
work[j] = t0 * work[j - (i__ - 4) * *stride] + t1 * work[j - (i__ - 3)
* *stride] + t2 * work[j - (i__ - 2) * *stride] + t3 * work[
j - (i__ - 1) * *stride] + t4 * work[j - i__ * *stride] + t5 *
workm5[j];
}
i__ = *npad - 2;
work2 = i__ * *stride;
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
work[j] = t0 * work[j - (i__ - 3) * *stride] + t1 * work[j - (i__ - 2)
* *stride] + t2 * work[j - (i__ - 1) * *stride] + t3 * work[
j - i__ * *stride] + t4 * workm4[j] + t5 * workm5[j];
}
i__ = *npad - 3;
work2 = i__ * *stride;
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
work[j] = t0 * work[j - (i__ - 2) * *stride] + t1 * work[j - (i__ - 1)
* *stride] + t2 * work[j - i__ * *stride] + t3 * workm3[j] +
t4 * workm4[j] + t5 * workm5[j];
}
i__ = *npad - 4;
work2 = i__ * *stride;
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
work[j] = t0 * work[j - (i__ - 1) * *stride] + t1 * work[j - i__ * *
stride] + t2 * workm2[j] + t3 * workm3[j] + t4 * workm4[j] +
t5 * workm5[j];
}
i__ = *npad - 5;
work2 = i__ * *stride;
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
work[j] = t0 * work[j - i__ * *stride] + t1 * workm1[j] + t2 * workm2[
j] + t3 * workm3[j] + t4 * workm4[j] + t5 * workm5[j];
}
for (i__ = *npad - 6; i__ >= 5; --i__) {
work2 = i__ * *stride;
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
work[j] = t0 * work[j] + t1 * workm1[j] + t2 * workm2[j] + t3 *
workm3[j] + t4 * workm4[j] + t5 * workm5[j];
}
}
work2 = *stride << 2;
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
work[j] = work[j + (*npad - 5) * *stride];
/* & t0*work(j)+t1*workm1(j)+t2*workm2(j)+t3*workm3(j)+t4*workm4(j) */
}
work2 = *stride * 3;
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
work[j] = work[j + (*npad - 5) * *stride];
/* & t0*work(j)+t1*workm1(j)+t2*workm2(j)+t3*workm3(j) */
}
work2 = *stride << 1;
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
work[j] = work[j + (*npad - 5) * *stride];
/* & t0*work(j)+t1*workm1(j)+t2*workm2(j) */
}
work2 = *stride;
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
work[j] = work[j + (*npad - 5) * *stride];
/* & t0*work(j)+t1*workm1(j) */
}
work2 = 0;
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
work[j] = work[j + (*npad - 5) * *stride];
/* & t0*work(j) */
}
/* Forward phase (padding now ensures room!) */
t0 = tf[1];
t1 = tf[2];
t2 = tf[3];
t3 = tf[4];
t4 = tf[5];
t5 = tf[6];
i__1 = *n - 1;
for (i__ = 0; i__ <= i__1; ++i__) {
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__2 = *area + work2;
for (j = work2 + 1; j <= i__2; ++j) {
jl = j + out2s;
datout[jl] = t0 * work[j] + t1 * workp1[j] + t2 * workp2[j] + t3 *
workp3[j] + t4 * workp4[j] + t5 * workp5[j];
}
}
return 0;
} /* convsetsioppercum5_ */
/* Back and forth, linear convolution of sets of points using padding */
/* for filters of extent +/- 5 (IN PLACE) */
/* out: input/output array for in place computation */
/* outm5-outp5,out: pointers in in/output array to */
/* start of data sets */
/* area: total size of each set */
/* stride: separation between data sets (stride=&out-&outm1) */
/* n: actual number of data sets being convolved */
/* npad: size of padded region needed for "leakage" of foward filters */
/* tr: reverse-time filter coefficients */
/* tf: forward-time filter coefficients */
/* Subroutine */ int convsetsiop5_(double *datout, double *work,
double *workm5, double *workm4, double *workm3,
double *workm2, double *workm1, double *workp1,
double *workp2, double *workp3, double *workp4,
double *workp5, int *area, int *stride, int *n,
int *npad, int *stridel, double *tr, double *tf)
{
/* System generated locals */
int i__1, i__2;
/* Local variables */
int work2, i__, j, out2s;
double t0, t1, t2, t3, t4, t5;
int jl;
/* First, make sure pad is zero */
/* Parameter adjustments */
--tf;
--tr;
--workp5;
--workp4;
--workp3;
--workp2;
--workp1;
--workm1;
--workm2;
--workm3;
--workm4;
--workm5;
--work;
--datout;
/* Function Body */
i__1 = *npad - 1;
for (i__ = *n; i__ <= i__1; ++i__) {
work2 = i__ * *stride;
i__2 = *area + work2;
for (j = work2 + 1; j <= i__2; ++j) {
work[j] = (double)0.;
}
}
/* Reverse phase */
t0 = tr[1];
t1 = tr[2];
t2 = tr[3];
t3 = tr[4];
t4 = tr[5];
t5 = tr[6];
for (i__ = *npad - 1; i__ >= 5; --i__) {
work2 = i__ * *stride;
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
work[j] = t0 * work[j] + t1 * workm1[j] + t2 * workm2[j] + t3 *
workm3[j] + t4 * workm4[j] + t5 * workm5[j];
}
}
work2 = *stride << 2;
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
work[j] = t0 * work[j] + t1 * workm1[j] + t2 * workm2[j] + t3 *
workm3[j] + t4 * workm4[j];
}
work2 = *stride * 3;
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
work[j] = t0 * work[j] + t1 * workm1[j] + t2 * workm2[j] + t3 *
workm3[j];
}
work2 = *stride << 1;
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
work[j] = t0 * work[j] + t1 * workm1[j] + t2 * workm2[j];
}
work2 = *stride;
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
work[j] = t0 * work[j] + t1 * workm1[j];
}
work2 = 0;
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
work[j] = t0 * work[j];
}
/* Forward phase (padding now ensures room!) */
t0 = tf[1];
t1 = tf[2];
t2 = tf[3];
t3 = tf[4];
t4 = tf[5];
t5 = tf[6];
i__1 = *n - 1;
for (i__ = 0; i__ <= i__1; ++i__) {
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__2 = *area + work2;
for (j = work2 + 1; j <= i__2; ++j) {
jl = j + out2s;
datout[jl] = datout[jl] + t0 * work[j] + t1 * workp1[j] + t2 *
workp2[j] + t3 * workp3[j] + t4 * workp4[j] + t5 * workp5[
j];
}
}
return 0;
} /* convsetsiop5_ */
/* Back and forth, linear convolution of sets of points using padding */
/* for filters of extent +/- 5 (IN PLACE) */
/* out: input/output array for in place computation */
/* outm5-outp5,out: pointers in in/output array to */
/* start of data sets */
/* area: total size of each set */
/* stride: separation between data sets (stride=&out-&outm1) */
/* n: actual number of data sets being convolved */
/* npad: size of padded region needed for "leakage" of foward filters */
/* tr: reverse-time filter coefficients */
/* tf: forward-time filter coefficients */
/* Subroutine */ int convsetsiop5full_(double *datout, double *work,
double *workm5, double *workm4, double *workm3,
double *workm2, double *workm1, double *workp1,
double *workp2, double *workp3, double *workp4,
double *workp5, int *area, int *stride, int *n,
int *npad, int *stridel, double *filter)
{
/* System generated locals */
int i__1, i__2;
/* Local variables */
int work2, i__, j, out2s;
double t0, t1, t2, t3, t4, t5;
int jl;
/* First, make sure pad is zero */
/* Parameter adjustments */
--filter;
--workp5;
--workp4;
--workp3;
--workp2;
--workp1;
--workm1;
--workm2;
--workm3;
--workm4;
--workm5;
--work;
--datout;
/* Function Body */
i__1 = *npad - 1;
for (i__ = *n; i__ <= i__1; ++i__) {
work2 = i__ * *stride;
i__2 = *area + work2;
for (j = work2 + 1; j <= i__2; ++j) {
work[j] = (double)0.;
}
}
/* Reverse phase */
t0 = filter[1];
t1 = filter[2];
t2 = filter[3];
t3 = filter[4];
t4 = filter[5];
t5 = filter[6];
if (t0 != (double)0.) {
/* symmetric part */
i__ = 0;
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
jl = j + out2s;
datout[jl] = datout[jl] + t0 * work[j] + t1 * workp1[j] + t2 *
workp2[j] + t3 * workp3[j] + t4 * workp4[j] + t5 * workp5[
j];
}
i__ = 1;
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
jl = j + out2s;
datout[jl] = datout[jl] + t1 * workm1[j] + t0 * work[j] + t1 *
workp1[j] + t2 * workp2[j] + t3 * workp3[j] + t4 * workp4[
j] + t5 * workp5[j];
}
i__ = 2;
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
jl = j + out2s;
datout[jl] = datout[jl] + t2 * workm2[j] + t1 * workm1[j] + t0 *
work[j] + t1 * workp1[j] + t2 * workp2[j] + t3 * workp3[j]
+ t4 * workp4[j] + t5 * workp5[j];
}
i__ = 3;
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
jl = j + out2s;
datout[jl] = datout[jl] + t3 * workm3[j] + t2 * workm2[j] + t1 *
workm1[j] + t0 * work[j] + t1 * workp1[j] + t2 * workp2[j]
+ t3 * workp3[j] + t4 * workp4[j] + t5 * workp5[j];
}
i__ = 4;
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
jl = j + out2s;
datout[jl] = datout[jl] + t4 * workm4[j] + t3 * workm3[j] + t2 *
workm2[j] + t1 * workm1[j] + t0 * work[j] + t1 * workp1[j]
+ t2 * workp2[j] + t3 * workp3[j] + t4 * workp4[j] + t5 *
workp5[j];
}
i__1 = *n - 6;
for (i__ = 5; i__ <= i__1; ++i__) {
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__2 = *area + work2;
for (j = work2 + 1; j <= i__2; ++j) {
jl = j + out2s;
datout[jl] = datout[jl] + t5 * workm5[j] + t4 * workm4[j] +
t3 * workm3[j] + t2 * workm2[j] + t1 * workm1[j] + t0
* work[j] + t1 * workp1[j] + t2 * workp2[j] + t3 *
workp3[j] + t4 * workp4[j] + t5 * workp5[j];
}
}
i__ = *n - 5;
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
jl = j + out2s;
datout[jl] = datout[jl] + t5 * workm5[j] + t4 * workm4[j] + t3 *
workm3[j] + t2 * workm2[j] + t1 * workm1[j] + t0 * work[j]
+ t1 * workp1[j] + t2 * workp2[j] + t3 * workp3[j] + t4 *
workp4[j];
}
i__ = *n - 4;
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
jl = j + out2s;
datout[jl] = datout[jl] + t5 * workm5[j] + t4 * workm4[j] + t3 *
workm3[j] + t2 * workm2[j] + t1 * workm1[j] + t0 * work[j]
+ t1 * workp1[j] + t2 * workp2[j] + t3 * workp3[j];
}
i__ = *n - 3;
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
jl = j + out2s;
datout[jl] = datout[jl] + t5 * workm5[j] + t4 * workm4[j] + t3 *
workm3[j] + t2 * workm2[j] + t1 * workm1[j] + t0 * work[j]
+ t1 * workp1[j] + t2 * workp2[j];
}
i__ = *n - 2;
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
jl = j + out2s;
datout[jl] = datout[jl] + t5 * workm5[j] + t4 * workm4[j] + t3 *
workm3[j] + t2 * workm2[j] + t1 * workm1[j] + t0 * work[j]
+ t1 * workp1[j];
}
i__ = *n - 1;
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
jl = j + out2s;
datout[jl] = datout[jl] + t5 * workm5[j] + t4 * workm4[j] + t3 *
workm3[j] + t2 * workm2[j] + t1 * workm1[j] + t0 * work[j]
;
}
/* end of symmetric part */
} else {
/* antisymmetric part */
i__ = 0;
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
jl = j + out2s;
datout[jl] = datout[jl] + t0 * work[j] + t1 * workp1[j] + t2 *
workp2[j] + t3 * workp3[j] + t4 * workp4[j] + t5 * workp5[
j];
}
i__ = 1;
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
jl = j + out2s;
datout[jl] = datout[jl] - t1 * workm1[j] + t0 * work[j] + t1 *
workp1[j] + t2 * workp2[j] + t3 * workp3[j] + t4 * workp4[
j] + t5 * workp5[j];
}
i__ = 2;
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
jl = j + out2s;
datout[jl] = datout[jl] - t2 * workm2[j] - t1 * workm1[j] + t0 *
work[j] + t1 * workp1[j] + t2 * workp2[j] + t3 * workp3[j]
+ t4 * workp4[j] + t5 * workp5[j];
}
i__ = 3;
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
jl = j + out2s;
datout[jl] = datout[jl] - t3 * workm3[j] - t2 * workm2[j] - t1 *
workm1[j] + t0 * work[j] + t1 * workp1[j] + t2 * workp2[j]
+ t3 * workp3[j] + t4 * workp4[j] + t5 * workp5[j];
}
i__ = 4;
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
jl = j + out2s;
datout[jl] = datout[jl] - t4 * workm4[j] - t3 * workm3[j] - t2 *
workm2[j] - t1 * workm1[j] + t0 * work[j] + t1 * workp1[j]
+ t2 * workp2[j] + t3 * workp3[j] + t4 * workp4[j] + t5 *
workp5[j];
}
i__1 = *n - 6;
for (i__ = 5; i__ <= i__1; ++i__) {
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__2 = *area + work2;
for (j = work2 + 1; j <= i__2; ++j) {
jl = j + out2s;
datout[jl] = datout[jl] - t5 * workm5[j] - t4 * workm4[j] -
t3 * workm3[j] - t2 * workm2[j] - t1 * workm1[j] + t0
* work[j] + t1 * workp1[j] + t2 * workp2[j] + t3 *
workp3[j] + t4 * workp4[j] + t5 * workp5[j];
}
}
i__ = *n - 5;
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
jl = j + out2s;
datout[jl] = datout[jl] - t5 * workm5[j] - t4 * workm4[j] - t3 *
workm3[j] - t2 * workm2[j] - t1 * workm1[j] + t0 * work[j]
+ t1 * workp1[j] + t2 * workp2[j] + t3 * workp3[j] + t4 *
workp4[j];
}
i__ = *n - 4;
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
jl = j + out2s;
datout[jl] = datout[jl] - t5 * workm5[j] - t4 * workm4[j] - t3 *
workm3[j] - t2 * workm2[j] - t1 * workm1[j] + t0 * work[j]
+ t1 * workp1[j] + t2 * workp2[j] + t3 * workp3[j];
}
i__ = *n - 3;
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
jl = j + out2s;
datout[jl] = datout[jl] - t5 * workm5[j] - t4 * workm4[j] - t3 *
workm3[j] - t2 * workm2[j] - t1 * workm1[j] + t0 * work[j]
+ t1 * workp1[j] + t2 * workp2[j];
}
i__ = *n - 2;
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
jl = j + out2s;
datout[jl] = datout[jl] - t5 * workm5[j] - t4 * workm4[j] - t3 *
workm3[j] - t2 * workm2[j] - t1 * workm1[j] + t0 * work[j]
+ t1 * workp1[j];
}
i__ = *n - 1;
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
jl = j + out2s;
datout[jl] = datout[jl] - t5 * workm5[j] - t4 * workm4[j] - t3 *
workm3[j] - t2 * workm2[j] - t1 * workm1[j] + t0 * work[j]
;
}
/* end of antisymmetric part */
}
return 0;
} /* convsetsiop5full_ */
/* Back and forth, linear convolution of sets of points using padding */
/* for filters of extent +/- 5 (IN PLACE) */
/* out: input/output array for in place computation */
/* outm5-outp5,out: pointers in in/output array to */
/* start of data sets */
/* area: total size of each set */
/* stride: separation between data sets (stride=&out-&outm1) */
/* n: actual number of data sets being convolved */
/* npad: size of padded region needed for "leakage" of foward filters */
/* tr: reverse-time filter coefficients */
/* tf: forward-time filter coefficients */
/* Subroutine */ int convsetsiopper5_(double *datout, double *work,
double *workm5, double *workm4, double *workm3,
double *workm2, double *workm1, double *workp1,
double *workp2, double *workp3, double *workp4,
double *workp5, int *area, int *stride, int *n,
int *npad, int *stridel, double *tr, double *tf)
{
/* System generated locals */
int i__1, i__2;
/* Local variables */
int work2, i__, j, out2s;
double t0, t1, t2, t3, t4, t5;
int jl;
/* First, make sure pad is zero */
/* Parameter adjustments */
--tf;
--tr;
--workp5;
--workp4;
--workp3;
--workp2;
--workp1;
--workm1;
--workm2;
--workm3;
--workm4;
--workm5;
--work;
--datout;
/* Function Body */
i__1 = *npad - 1;
for (i__ = *n; i__ <= i__1; ++i__) {
work2 = i__ * *stride;
i__2 = *area + work2;
for (j = work2 + 1; j <= i__2; ++j) {
work[j] = (double)0.;
}
}
/* Reverse phase */
t0 = tr[1];
t1 = tr[2];
t2 = tr[3];
t3 = tr[4];
t4 = tr[5];
t5 = tr[6];
i__ = *npad - 1;
work2 = i__ * *stride;
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
work[j] = t0 * work[j - (i__ - 4) * *stride] + t1 * work[j - (i__ - 3)
* *stride] + t2 * work[j - (i__ - 2) * *stride] + t3 * work[
j - (i__ - 1) * *stride] + t4 * work[j - i__ * *stride] + t5 *
workm5[j];
}
i__ = *npad - 2;
work2 = i__ * *stride;
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
work[j] = t0 * work[j - (i__ - 3) * *stride] + t1 * work[j - (i__ - 2)
* *stride] + t2 * work[j - (i__ - 1) * *stride] + t3 * work[
j - i__ * *stride] + t4 * workm4[j] + t5 * workm5[j];
}
i__ = *npad - 3;
work2 = i__ * *stride;
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
work[j] = t0 * work[j - (i__ - 2) * *stride] + t1 * work[j - (i__ - 1)
* *stride] + t2 * work[j - i__ * *stride] + t3 * workm3[j] +
t4 * workm4[j] + t5 * workm5[j];
}
i__ = *npad - 4;
work2 = i__ * *stride;
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
work[j] = t0 * work[j - (i__ - 1) * *stride] + t1 * work[j - i__ * *
stride] + t2 * workm2[j] + t3 * workm3[j] + t4 * workm4[j] +
t5 * workm5[j];
}
i__ = *npad - 5;
work2 = i__ * *stride;
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
work[j] = t0 * work[j - i__ * *stride] + t1 * workm1[j] + t2 * workm2[
j] + t3 * workm3[j] + t4 * workm4[j] + t5 * workm5[j];
}
for (i__ = *npad - 6; i__ >= 5; --i__) {
work2 = i__ * *stride;
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
work[j] = t0 * work[j] + t1 * workm1[j] + t2 * workm2[j] + t3 *
workm3[j] + t4 * workm4[j] + t5 * workm5[j];
}
}
work2 = *stride << 2;
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
work[j] = work[j + (*npad - 5) * *stride];
/* & t0*work(j)+t1*workm1(j)+t2*workm2(j)+t3*workm3(j)+t4*workm4(j) */
}
work2 = *stride * 3;
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
work[j] = work[j + (*npad - 5) * *stride];
/* & t0*work(j)+t1*workm1(j)+t2*workm2(j)+t3*workm3(j) */
}
work2 = *stride << 1;
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
work[j] = work[j + (*npad - 5) * *stride];
/* & t0*work(j)+t1*workm1(j)+t2*workm2(j) */
}
work2 = *stride;
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
work[j] = work[j + (*npad - 5) * *stride];
/* & t0*work(j)+t1*workm1(j) */
}
work2 = 0;
i__1 = *area + work2;
for (j = work2 + 1; j <= i__1; ++j) {
work[j] = work[j + (*npad - 5) * *stride];
/* & t0*work(j) */
}
/* Forward phase (padding now ensures room!) */
t0 = tf[1];
t1 = tf[2];
t2 = tf[3];
t3 = tf[4];
t4 = tf[5];
t5 = tf[6];
i__1 = *n - 1;
for (i__ = 0; i__ <= i__1; ++i__) {
work2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__2 = *area + work2;
for (j = work2 + 1; j <= i__2; ++j) {
jl = j + out2s;
datout[jl] = datout[jl] + t0 * work[j] + t1 * workp1[j] + t2 *
workp2[j] + t3 * workp3[j] + t4 * workp4[j] + t5 * workp5[
j];
}
}
return 0;
} /* convsetsiopper5_ */
/* Back and forth, linear convolution of sets of points using padding, */
/* computing data ONLY on ghostpoints (IN PLACE) */
/* out: input/output array for in place computation */
/* outm3,outm2,outm1,outp1,outp2,outp3: pointers in output array to */
/* start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, */
/* 2 forward, 3 forward */
/* area: total size of each set */
/* stride: separation between data sets (stride=&out-&outm1) */
/* n: actual number of data sets being convolved */
/* npad: size of padded region needed for "leakage" of foward filters */
/* tr: reverse-time filter coefficients */
/* tf: forward-time filter coefficients */
/* Subroutine */ int pconvsets3_(double *out, double *outm3,
double *outm2, double *outm1, double *outp1, double *
outp2, double *outp3, int *area, int *stride, int *n,
int *npad, double *tr, double *tf)
{
/* System generated locals */
int i__1, i__2;
/* Local variables */
int i__, j;
double t0, t1, t2, t3;
int out2;
/* First, make sure pad is zero for forward phase. */
/* For pretty code the last number would be npad-1, but technically */
/* npad-2 (for extra speed) is all that is needed. */
/* Parameter adjustments */
--tf;
--tr;
--outp3;
--outp2;
--outp1;
--outm1;
--outm2;
--outm3;
--out;
/* Function Body */
i__1 = *npad - 2;
for (i__ = *n; i__ <= i__1; ++i__) {
out2 = i__ * *stride;
i__2 = *area + out2;
for (j = out2 + 1; j <= i__2; ++j) {
out[j] = (double)0.;
}
}
/* Reverse phase */
t0 = tr[1];
t1 = tr[2];
t2 = tr[3];
t3 = tr[4];
for (i__ = *n - 1; i__ >= 3; --i__) {
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j] + t1 * outm1[j] + t2 * outm2[j] + t3 * outm3[
j];
}
}
out2 = *stride << 1;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j] + t1 * outm1[j] + t2 * outm2[j];
}
out2 = *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j] + t1 * outm1[j];
}
out2 = 0;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j];
}
/* Forward phase (padding now ensures room!) */
t0 = tf[1];
t1 = tf[2];
t2 = tf[3];
t3 = tf[4];
i__1 = *n - 1;
for (i__ = 0; i__ <= i__1; i__ += 2) {
out2 = i__ * *stride;
i__2 = *area + out2;
for (j = out2 + 1; j <= i__2; ++j) {
out[j] = t0 * out[j] + t1 * outp1[j] + t2 * outp2[j] + t3 * outp3[
j];
}
}
return 0;
} /* pconvsets3_ */
/* Back and forth, linear convolution of sets of points using padding, */
/* computing data ONLY on ghostpoints (IN PLACE) */
/* out: input/output array for in place computation */
/* outm3,outm2,outm1,outp1,outp2,outp3: pointers in output array to */
/* start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, */
/* 2 forward, 3 forward */
/* area: total size of each set */
/* stride: separation between data sets (stride=&out-&outm1) */
/* n: actual number of data sets being convolved */
/* npad: size of padded region needed for "leakage" of foward filters */
/* tr: reverse-time filter coefficients */
/* tf: forward-time filter coefficients */
/* Subroutine */ int pconvsetsio3full_(double *out, double *in,
double *inm3, double *inm2, double *inm1, double *
inp1, double *inp2, double *inp3, int *area, int *
stride, int *n, int *npad, double *filt)
{
/* System generated locals */
int i__1, i__2;
/* Local variables */
int i__, j;
double t0, t1, t3;
int out2;
/* Forward phase (padding now ensures room!) */
/* Parameter adjustments */
--filt;
--inp3;
--inp2;
--inp1;
--inm1;
--inm2;
--inm3;
--in;
--out;
/* Function Body */
t0 = filt[1];
t1 = filt[2];
t3 = filt[4];
i__ = 0;
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t3 * inp3[j];
}
i__ = 2;
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = in[j] + t1 * inp1[j] + t3 * inp3[j];
}
i__1 = *n - 6;
for (i__ = 4; i__ <= i__1; i__ += 2) {
out2 = i__ * *stride;
i__2 = *area + out2;
for (j = out2 + 1; j <= i__2; ++j) {
out[j] = t3 * inm3[j] + t1 * inm1[j] + in[j] + t1 * inp1[j] + t3 *
inp3[j];
}
}
i__ = *n - 4;
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t3 * inm3[j] + t1 * inm1[j] + in[j];
}
i__ = *n - 2;
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t3 * inm3[j];
}
return 0;
} /* pconvsetsio3full_ */
/* Back and forth, linear convolution of sets of points using padding, */
/* computing data ONLY on ghostpoints (IN PLACE) */
/* out: input/output array for in place computation */
/* outm3,outm2,outm1,outp1,outp2,outp3: pointers in output array to */
/* start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, */
/* 2 forward, 3 forward */
/* area: total size of each set */
/* stride: separation between data sets (stride=&out-&outm1) */
/* n: actual number of data sets being convolved */
/* npad: size of padded region needed for "leakage" of foward filters */
/* tr: reverse-time filter coefficients */
/* tf: forward-time filter coefficients */
/* Subroutine */ int pconvsetsio8full_(double *out, double *in,
double *inm8, double *inm7, double *inm6, double *
inm5, double *inm4, double *inm3, double *inm2,
double *inm1, double *inp1, double *inp2, double *
inp3, double *inp4, double *inp5, double *inp6,
double *inp7, double *inp8, int *area, int *stride,
int *n, int *npad, double *filt)
{
/* System generated locals */
int i__1, i__2;
/* Local variables */
int i__, j;
double t0, t1, t2, t3, t4, t5, t6, t7, t8;
int out2;
/* Forward phase (padding now ensures room!) */
/* Parameter adjustments */
--filt;
--inp8;
--inp7;
--inp6;
--inp5;
--inp4;
--inp3;
--inp2;
--inp1;
--inm1;
--inm2;
--inm3;
--inm4;
--inm5;
--inm6;
--inm7;
--inm8;
--in;
--out;
/* Function Body */
t0 = filt[1];
t1 = filt[2];
t2 = filt[3];
t3 = filt[4];
t4 = filt[5];
t5 = filt[6];
t6 = filt[7];
t7 = filt[8];
t8 = filt[9];
if (t0 != (double)0.) {
/* symmetric part */
i__ = 0;
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t8 * inp8[j];
}
i__ = 2;
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t6 * inp6[j] + t7 * inp7[j] + t8 * inp8[j];
}
i__ = 4;
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t4 * inp4[j] + t5 * inp5[j] + t6 * inp6[j] + t7 * inp7[j]
+ t8 * inp8[j];
}
i__ = 6;
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t2 * inp2[j] + t3 * inp3[j] + t4 * inp4[j] + t5 * inp5[j]
+ t6 * inp6[j] + t7 * inp7[j] + t8 * inp8[j];
}
i__ = 8;
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * in[j] + t1 * inp1[j] + t2 * inp2[j] + t3 * inp3[j]
+ t4 * inp4[j] + t5 * inp5[j] + t6 * inp6[j] + t7 * inp7[
j] + t8 * inp8[j];
}
i__ = 10;
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t2 * inm2[j] + t1 * inm1[j] + t0 * in[j] + t1 * inp1[j]
+ t2 * inp2[j] + t3 * inp3[j] + t4 * inp4[j] + t5 * inp5[
j] + t6 * inp6[j] + t7 * inp7[j] + t8 * inp8[j];
}
i__ = 12;
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t4 * inm4[j] + t3 * inm3[j] + t2 * inm2[j] + t1 * inm1[j]
+ t0 * in[j] + t1 * inp1[j] + t2 * inp2[j] + t3 * inp3[j]
+ t4 * inp4[j] + t5 * inp5[j] + t6 * inp6[j] + t7 * inp7[
j] + t8 * inp8[j];
}
i__ = 14;
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t6 * inm6[j] + t5 * inm5[j] + t4 * inm4[j] + t3 * inm3[j]
+ t2 * inm2[j] + t1 * inm1[j] + t0 * in[j] + t1 * inp1[j]
+ t2 * inp2[j] + t3 * inp3[j] + t4 * inp4[j] + t5 * inp5[
j] + t6 * inp6[j] + t7 * inp7[j] + t8 * inp8[j];
}
i__1 = *n - 18;
for (i__ = 16; i__ <= i__1; i__ += 2) {
out2 = i__ * *stride;
i__2 = *area + out2;
for (j = out2 + 1; j <= i__2; ++j) {
out[j] = t8 * inm8[j] + t7 * inm7[j] + t6 * inm6[j] + t5 *
inm5[j] + t4 * inm4[j] + t3 * inm3[j] + t2 * inm2[j]
+ t1 * inm1[j] + t0 * in[j] + t1 * inp1[j] + t2 *
inp2[j] + t3 * inp3[j] + t4 * inp4[j] + t5 * inp5[j]
+ t6 * inp6[j] + t7 * inp7[j] + t8 * inp8[j];
}
}
i__ = *n - 16;
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t7 * inp7[j] + t6 * inp6[j] + t5 * inp5[j] + t4 * inp4[j]
+ t3 * inp3[j] + t2 * inp2[j] + t1 * inp1[j] + t0 * in[j]
+ t1 * inm1[j] + t2 * inm2[j] + t3 * inm3[j] + t4 * inm4[
j] + t5 * inm5[j] + t6 * inm6[j] + t7 * inm7[j] + t8 *
inm8[j];
}
i__ = *n - 14;
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t5 * inp5[j] + t4 * inp4[j] + t3 * inp3[j] + t2 * inp2[j]
+ t1 * inp1[j] + t0 * in[j] + t1 * inm1[j] + t2 * inm2[j]
+ t3 * inm3[j] + t4 * inm4[j] + t5 * inm5[j] + t6 * inm6[
j] + t7 * inm7[j] + t8 * inm8[j];
}
i__ = *n - 12;
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t3 * inp3[j] + t2 * inp2[j] + t1 * inp1[j] + t0 * in[j]
+ t1 * inm1[j] + t2 * inm2[j] + t3 * inm3[j] + t4 * inm4[
j] + t5 * inm5[j] + t6 * inm6[j] + t7 * inm7[j] + t8 *
inm8[j];
}
i__ = *n - 10;
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t1 * inp1[j] + t0 * in[j] + t1 * inm1[j] + t2 * inm2[j]
+ t3 * inm3[j] + t4 * inm4[j] + t5 * inm5[j] + t6 * inm6[
j] + t7 * inm7[j] + t8 * inm8[j];
}
i__ = *n - 8;
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t1 * inm1[j] + t2 * inm2[j] + t3 * inm3[j] + t4 * inm4[j]
+ t5 * inm5[j] + t6 * inm6[j] + t7 * inm7[j] + t8 * inm8[
j];
}
i__ = *n - 6;
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t3 * inm3[j] + t4 * inm4[j] + t5 * inm5[j] + t6 * inm6[j]
+ t7 * inm7[j] + t8 * inm8[j];
}
i__ = *n - 4;
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t5 * inm5[j] + t6 * inm6[j] + t7 * inm7[j] + t8 * inm8[j]
;
}
i__ = *n - 2;
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t7 * inm7[j] + t8 * inm8[j];
}
/* end of symmetric part */
} else {
/* antisymmetric part */
i__ = 0;
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t8 * inp8[j];
}
i__ = 2;
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t6 * inp6[j] + t7 * inp7[j] + t8 * inp8[j];
}
i__ = 4;
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t4 * inp4[j] + t5 * inp5[j] + t6 * inp6[j] + t7 * inp7[j]
+ t8 * inp8[j];
}
i__ = 6;
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t2 * inp2[j] + t3 * inp3[j] + t4 * inp4[j] + t5 * inp5[j]
+ t6 * inp6[j] + t7 * inp7[j] + t8 * inp8[j];
}
i__ = 8;
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * in[j] + t1 * inp1[j] + t2 * inp2[j] + t3 * inp3[j]
+ t4 * inp4[j] + t5 * inp5[j] + t6 * inp6[j] + t7 * inp7[
j] + t8 * inp8[j];
}
i__ = 10;
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t2 * (double)-1. * inm2[j] - t1 * inm1[j] + t0 * in[j] +
t1 * inp1[j] + t2 * inp2[j] + t3 * inp3[j] + t4 * inp4[j]
+ t5 * inp5[j] + t6 * inp6[j] + t7 * inp7[j] + t8 * inp8[
j];
}
i__ = 12;
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t4 * (double)-1. * inm4[j] - t3 * inm3[j] - t2 * inm2[j]
- t1 * inm1[j] + t0 * in[j] + t1 * inp1[j] + t2 * inp2[j]
+ t3 * inp3[j] + t4 * inp4[j] + t5 * inp5[j] + t6 * inp6[
j] + t7 * inp7[j] + t8 * inp8[j];
}
i__ = 14;
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t6 * (double)-1. * inm6[j] - t5 * inm5[j] - t4 * inm4[j]
- t3 * inm3[j] - t2 * inm2[j] - t1 * inm1[j] + t0 * in[j]
+ t1 * inp1[j] + t2 * inp2[j] + t3 * inp3[j] + t4 * inp4[
j] + t5 * inp5[j] + t6 * inp6[j] + t7 * inp7[j] + t8 *
inp8[j];
}
i__1 = *n - 18;
for (i__ = 16; i__ <= i__1; i__ += 2) {
out2 = i__ * *stride;
i__2 = *area + out2;
for (j = out2 + 1; j <= i__2; ++j) {
out[j] = t8 * (double)-1. * inm8[j] - t7 * inm7[j] - t6 * inm6[
j] - t5 * inm5[j] - t4 * inm4[j] - t3 * inm3[j] - t2 *
inm2[j] - t1 * inm1[j] + t0 * in[j] + t1 * inp1[j] +
t2 * inp2[j] + t3 * inp3[j] + t4 * inp4[j] + t5 *
inp5[j] + t6 * inp6[j] + t7 * inp7[j] + t8 * inp8[j];
}
}
i__ = *n - 16;
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t7 * inp7[j] + t6 * inp6[j] + t5 * inp5[j] + t4 * inp4[j]
+ t3 * inp3[j] + t2 * inp2[j] + t1 * inp1[j] + t0 * in[j]
- t1 * inm1[j] - t2 * inm2[j] - t3 * inm3[j] - t4 * inm4[
j] - t5 * inm5[j] - t6 * inm6[j] - t7 * inm7[j] - t8 *
inm8[j];
}
i__ = *n - 14;
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t5 * inp5[j] + t4 * inp4[j] + t3 * inp3[j] + t2 * inp2[j]
+ t1 * inp1[j] + t0 * in[j] - t1 * inm1[j] - t2 * inm2[j]
- t3 * inm3[j] - t4 * inm4[j] - t5 * inm5[j] - t6 * inm6[
j] - t7 * inm7[j] - t8 * inm8[j];
}
i__ = *n - 12;
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t3 * inp3[j] + t2 * inp2[j] + t1 * inp1[j] + t0 * in[j]
- t1 * inm1[j] - t2 * inm2[j] - t3 * inm3[j] - t4 * inm4[
j] - t5 * inm5[j] - t6 * inm6[j] - t7 * inm7[j] - t8 *
inm8[j];
}
i__ = *n - 10;
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t1 * inp1[j] + t0 * in[j] - t1 * inm1[j] - t2 * inm2[j]
- t3 * inm3[j] - t4 * inm4[j] - t5 * inm5[j] - t6 * inm6[
j] - t7 * inm7[j] - t8 * inm8[j];
}
i__ = *n - 8;
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t1 * (double)-1. * inm1[j] - t2 * inm2[j] - t3 * inm3[j]
- t4 * inm4[j] - t5 * inm5[j] - t6 * inm6[j] - t7 * inm7[
j] - t8 * inm8[j];
}
i__ = *n - 6;
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t3 * (double)-1. * inm3[j] - t4 * inm4[j] - t5 * inm5[j]
- t6 * inm6[j] - t7 * inm7[j] - t8 * inm8[j];
}
i__ = *n - 4;
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t5 * (double)-1. * inm5[j] - t6 * inm6[j] - t7 * inm7[j]
- t8 * inm8[j];
}
i__ = *n - 2;
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t7 * (double)-1. * inm7[j] - t8 * inm8[j];
}
/* end of antisymmetric part */
}
return 0;
} /* pconvsetsio8full_ */
/* Back and forth, linear convolution of sets of points using padding, */
/* computing data ONLY on ghostpoints (IN PLACE) */
/* out: input/output array for in place computation */
/* outm3,outm2,outm1,outp1,outp2,outp3: pointers in output array to */
/* start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, */
/* 2 forward, 3 forward */
/* area: total size of each set */
/* stride: separation between data sets (stride=&out-&outm1) */
/* n: actual number of data sets being convolved */
/* npad: size of padded region needed for "leakage" of foward filters */
/* tr: reverse-time filter coefficients */
/* tf: forward-time filter coefficients */
/* Subroutine */ int pconvsets3new_(double *outx, double *out,
double *outm3, double *outm2, double *outm1, double *
outp1, double *outp2, double *outp3, int *area, int *
stride, int *n, int *npad, double *tr, double *tf)
{
/* System generated locals */
int i__1, i__2;
/* Local variables */
int i__, j;
double t0, t1, t2, t3;
int out2;
/* First, make sure pad is zero for forward phase. */
/* For pretty code the last number would be npad-1, but technically */
/* npad-2 (for extra speed) is all that is needed. */
/* Parameter adjustments */
--tf;
--tr;
--outp3;
--outp2;
--outp1;
--outm1;
--outm2;
--outm3;
--out;
--outx;
/* Function Body */
i__1 = *npad - 2;
for (i__ = *n; i__ <= i__1; ++i__) {
out2 = i__ * *stride;
i__2 = *area + out2;
for (j = out2 + 1; j <= i__2; ++j) {
out[j] = (double)0.;
}
}
/* Reverse phase */
t0 = tr[1];
t1 = tr[2];
t2 = tr[3];
t3 = tr[4];
for (i__ = *n - 1; i__ >= 3; --i__) {
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j] + t1 * outm1[j] + t2 * outm2[j] + t3 * outm3[
j];
}
}
out2 = *stride << 1;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j] + t1 * outm1[j] + t2 * outm2[j];
}
out2 = *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j] + t1 * outm1[j];
}
out2 = 0;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j];
}
/* Forward phase (padding now ensures room!) */
t0 = tf[1];
t1 = tf[2];
t2 = tf[3];
t3 = tf[4];
i__1 = *n - 1;
for (i__ = 0; i__ <= i__1; i__ += 2) {
out2 = i__ * *stride;
i__2 = *area + out2;
for (j = out2 + 1; j <= i__2; ++j) {
outx[j] = t0 * out[j] + t1 * outp1[j] + t2 * outp2[j] + t3 *
outp3[j];
}
}
return 0;
} /* pconvsets3new_ */
/* Back and forth, linear convolution of sets of points using padding, */
/* computing data ONLY on ghostpoints (IN PLACE) */
/* out: input/output array for in place computation */
/* outm3,outm2,outm1,outp1,outp2,outp3: pointers in output array to */
/* start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, */
/* 2 forward, 3 forward */
/* area: total size of each set */
/* stride: separation between data sets (stride=&out-&outm1) */
/* n: actual number of data sets being convolved */
/* npad: size of padded region needed for "leakage" of foward filters */
/* tr: reverse-time filter coefficients */
/* tf: forward-time filter coefficients */
/* Subroutine */ int pconvsets8_(double *out, double *outm8,
double *outm7, double *outm6, double *outm5, double *
outm4, double *outm3, double *outm2, double *outm1,
double *outp1, double *outp2, double *outp3, double *
outp4, double *outp5, double *outp6, double *outp7,
double *outp8, int *area, int *stride, int *n,
int *npad, double *tr, double *tf)
{
/* System generated locals */
int i__1, i__2;
/* Local variables */
int i__, j;
double t0, t1, t2, t3, t4, t5, t6, t7, t8;
int out2;
/* First, make sure pad is zero for forward phase. */
/* For pretty code the last number would be npad-1, but technically */
/* npad-2 (for extra speed) is all that is needed. */
/* Parameter adjustments */
--tf;
--tr;
--outp8;
--outp7;
--outp6;
--outp5;
--outp4;
--outp3;
--outp2;
--outp1;
--outm1;
--outm2;
--outm3;
--outm4;
--outm5;
--outm6;
--outm7;
--outm8;
--out;
/* Function Body */
i__1 = *npad - 2;
for (i__ = *n; i__ <= i__1; ++i__) {
out2 = i__ * *stride;
i__2 = *area + out2;
for (j = out2 + 1; j <= i__2; ++j) {
out[j] = (double)0.;
}
}
/* Reverse phase */
t0 = tr[1];
t1 = tr[2];
t2 = tr[3];
t3 = tr[4];
t4 = tr[5];
t5 = tr[6];
t6 = tr[7];
t7 = tr[8];
t8 = tr[9];
for (i__ = *n - 1; i__ >= 8; --i__) {
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j] + t1 * outm1[j] + t2 * outm2[j] + t3 * outm3[
j] + t4 * outm4[j] + t5 * outm5[j] + t6 * outm6[j] + t7 *
outm7[j] + t8 * outm8[j];
}
}
out2 = *stride * 7;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j] + t1 * outm1[j] + t2 * outm2[j] + t3 * outm3[j]
+ t4 * outm4[j] + t5 * outm5[j] + t6 * outm6[j] + t7 * outm7[
j];
}
out2 = *stride * 6;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j] + t1 * outm1[j] + t2 * outm2[j] + t3 * outm3[j]
+ t4 * outm4[j] + t5 * outm5[j] + t6 * outm6[j];
}
out2 = *stride * 5;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j] + t1 * outm1[j] + t2 * outm2[j] + t3 * outm3[j]
+ t4 * outm4[j] + t5 * outm5[j];
}
out2 = *stride << 2;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j] + t1 * outm1[j] + t2 * outm2[j] + t3 * outm3[j]
+ t4 * outm4[j];
}
out2 = *stride * 3;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j] + t1 * outm1[j] + t2 * outm2[j] + t3 * outm3[j];
}
out2 = *stride << 1;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j] + t1 * outm1[j] + t2 * outm2[j];
}
out2 = *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j] + t1 * outm1[j];
}
out2 = 0;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j];
}
/* Forward phase (padding now ensures room!) */
t0 = tf[1];
t1 = tf[2];
t2 = tf[3];
t3 = tf[4];
t4 = tf[5];
t5 = tf[6];
t6 = tf[7];
t7 = tf[8];
t8 = tf[9];
i__1 = *n - 1;
for (i__ = 0; i__ <= i__1; i__ += 2) {
out2 = i__ * *stride;
i__2 = *area + out2;
for (j = out2 + 1; j <= i__2; ++j) {
out[j] = t0 * out[j] + t1 * outp1[j] + t2 * outp2[j] + t3 * outp3[
j] + t4 * outp4[j] + t5 * outp5[j] + t6 * outp6[j] + t7 *
outp7[j] + t8 * outp8[j];
}
}
return 0;
} /* pconvsets8_ */
/* Back and forth, linear convolution of sets of points using padding */
/* computing data ONLY on ghostpoints, AND folding in input of a */
/* different data set (IN PLACE) */
/* out: input/output array for in place computation */
/* outm3,outm2,outm1,outp1,outp2,outp3: pointers in output array to */
/* start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, */
/* 2 forward, 3 forward */
/* in: input/output array for in place computation */
/* inm3,inm2,inm1,inp1,inp2,inp3: pointers in input array to */
/* start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, */
/* 2 forward, 3 forward */
/* area: total size of out sets */
/* stride: separation between out data sets (stride=&out-&outm1) */
/* n: actual number of data sets being convolved */
/* npad: size of padded region needed for "leakage" of foward filters */
/* stridel: separation between in data sets (stride=&in-&inm1) */
/* fptl, lptl: first and last data components actually present in input */
/* tr: reverse-time filter coefficients */
/* tf: forward-time filter coefficients */
/* eflag: mod(,2) of present line (tells if special ghost handling needed) */
/* zcflag: =1 means we can count on zero ghosts and need special handling */
/* Subroutine */ int pconvsetsi3old_(double *out, double *outm3,
double *outm2, double *outm1, double *outp1, double *
outp2, double *outp3, double *in, double *inm3,
double *inm2, double *inm1, double *inp1, double *
inp2, double *inp3, int *area, int *stride, int *n,
int *npad, int *stridel, int *fptl, int *lptl,
double *tr, double *tf, int *eflag, int *zcflag)
{
/* System generated locals */
int i__1, i__2;
/* Local variables */
extern /* Subroutine */ int lineioo3_(double *, double *,
double *, double *, double *, int *, int *,
int *, int *, int *, int *, double *,
double *, double *, double *);
int i__, j;
extern /* Subroutine */ int zeroline3_(double *, double *,
double *, double *, double *, int *, int *,
int *, int *, int *, int *, double *,
double *, double *, double *);
double t0, t1, t2, t3;
extern /* Subroutine */ int lineio3_(double *, double *,
double *, double *, double *, int *, int *,
int *, int *, int *, int *, double *,
double *, double *, double *), lineioe3_(double *,
double *, double *, double *, double *, int *
, int *, int *, int *, int *, int *,
double *, double *, double *, double *);
int out2;
/* Reverse phase */
/* Parameter adjustments */
--tf;
--tr;
--inp3;
--inp2;
--inp1;
--inm1;
--inm2;
--inm3;
--in;
--outp3;
--outp2;
--outp1;
--outm1;
--outm2;
--outm3;
--out;
/* Function Body */
t0 = tr[1];
t1 = tr[2];
t2 = tr[3];
t3 = tr[4];
/* First, make sure pad is zero for forward phase. */
/* For pretty code the last number would be npad-1, but technically */
/* npad-2 (for extra speed) is all that is needed. */
i__1 = *npad - 2;
for (i__ = *n; i__ <= i__1; ++i__) {
zeroline3_(&out[1], &in[1], &inm1[1], &inm2[1], &inm3[1], &i__,
stridel, stride, fptl, lptl, area, &t0, &t1, &t2, &t3);
}
i__1 = *n - 1;
zeroline3_(&out[1], &in[1], &inm1[1], &inm2[1], &inm3[1], &i__1, stridel,
stride, fptl, lptl, area, &t0, &t1, &t2, &t3);
if (*eflag == 0 && *zcflag == 1) {
/* For Idag (zcflag=1) we are coming from real space, and so can't */
/* depend on the ghostpoints being zero. Thus, on the even planes (eflag), */
/* we must skip over the data in the appropriate lines for the even */
/* points. This is what the line??[eo] routines do. */
/* For the UNUSED lines below, we point to "inm3" to avoid seg faults */
/* because we know "inm3" will always be in range (unless the box is */
/* ridiculously small). Although multiplied by zero, the code still */
/* accesses the dummy data! */
i__1 = *n - 2;
lineioe3_(&out[1], &inm3[1], &inm3[1], &inm3[1], &inm3[1], &i__1,
stridel, stride, fptl, lptl, area, &c_b16, &c_b16, &c_b16, &
t3);
i__1 = *n - 3;
lineioo3_(&out[1], &inm3[1], &inm3[1], &inm2[1], &inm3[1], &i__1,
stridel, stride, fptl, lptl, area, &c_b16, &c_b16, &t2, &t3);
i__1 = *n - 4;
lineioe3_(&out[1], &inm3[1], &inm1[1], &inm2[1], &inm3[1], &i__1,
stridel, stride, fptl, lptl, area, &c_b16, &t1, &t2, &t3);
for (i__ = *n - 5; i__ >= 5; --i__) {
if (i__ % 2 == 0) {
lineioe3_(&out[1], &in[1], &inm1[1], &inm2[1], &inm3[1], &i__,
stridel, stride, fptl, lptl, area, &t0, &t1, &t2, &
t3);
} else {
lineioo3_(&out[1], &in[1], &inm1[1], &inm2[1], &inm3[1], &i__,
stridel, stride, fptl, lptl, area, &t0, &t1, &t2, &
t3);
}
}
lineioe3_(&out[1], &in[1], &inm1[1], &inm2[1], &in[1], &c__4, stridel,
stride, fptl, lptl, area, &t0, &t1, &t2, &c_b16);
lineioo3_(&out[1], &in[1], &inm1[1], &in[1], &in[1], &c__3, stridel,
stride, fptl, lptl, area, &t0, &t1, &c_b16, &c_b16);
lineioe3_(&out[1], &in[1], &in[1], &in[1], &in[1], &c__2, stridel,
stride, fptl, lptl, area, &t0, &c_b16, &c_b16, &c_b16);
} else {
/* For the UNUSED lines below, we point to "inm3" to avoid seg faults */
/* because we know "inm3" will always be in range (unless the box is */
/* ridiculously small). Although multiplied by zero, the code still */
/* accesses the dummy data! */
i__1 = *n - 2;
lineio3_(&out[1], &inm3[1], &inm3[1], &inm3[1], &inm3[1], &i__1,
stridel, stride, fptl, lptl, area, &c_b16, &c_b16, &c_b16, &
t3);
i__1 = *n - 3;
lineio3_(&out[1], &inm3[1], &inm3[1], &inm2[1], &inm3[1], &i__1,
stridel, stride, fptl, lptl, area, &c_b16, &c_b16, &t2, &t3);
i__1 = *n - 4;
lineio3_(&out[1], &inm3[1], &inm1[1], &inm2[1], &inm3[1], &i__1,
stridel, stride, fptl, lptl, area, &c_b16, &t1, &t2, &t3);
for (i__ = *n - 5; i__ >= 5; --i__) {
lineio3_(&out[1], &in[1], &inm1[1], &inm2[1], &inm3[1], &i__,
stridel, stride, fptl, lptl, area, &t0, &t1, &t2, &t3);
}
/* For the UNUSED lines below, we must point to "in" to avoid seg */
/* faults. On this end of the convolution, we know that "in" will */
/* always be in range. */
lineio3_(&out[1], &in[1], &inm1[1], &inm2[1], &in[1], &c__4, stridel,
stride, fptl, lptl, area, &t0, &t1, &t2, &c_b16);
lineio3_(&out[1], &in[1], &inm1[1], &in[1], &in[1], &c__3, stridel,
stride, fptl, lptl, area, &t0, &t1, &c_b16, &c_b16);
lineio3_(&out[1], &in[1], &in[1], &in[1], &in[1], &c__2, stridel,
stride, fptl, lptl, area, &t0, &c_b16, &c_b16, &c_b16);
}
zeroline3_(&out[1], &in[1], &inm1[1], &inm2[1], &inm3[1], &c__1, stridel,
stride, fptl, lptl, area, &t0, &t1, &t2, &t3);
zeroline3_(&out[1], &in[1], &inm1[1], &inm2[1], &inm3[1], &c__0, stridel,
stride, fptl, lptl, area, &t0, &t1, &t2, &t3);
/* Forward phase (padding now ensures room!) */
t0 = tf[1];
t1 = tf[2];
t2 = tf[3];
t3 = tf[4];
i__1 = *n - 1;
for (i__ = 0; i__ <= i__1; i__ += 2) {
out2 = i__ * *stride;
i__2 = *area + out2;
for (j = out2 + 1; j <= i__2; ++j) {
out[j] = t0 * out[j] + t1 * outp1[j] + t2 * outp2[j] + t3 * outp3[
j];
}
}
return 0;
} /* pconvsetsi3old_ */
/* Back and forth, linear convolution of sets of points using padding */
/* computing data ONLY on ghostpoints, AND folding in input of a */
/* different data set (IN PLACE) */
/* out: input/output array for in place computation */
/* outm3,outm2,outm1,outp1,outp2,outp3: pointers in output array to */
/* start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, */
/* 2 forward, 3 forward */
/* in: input/output array for in place computation */
/* inm3,inm2,inm1,inp1,inp2,inp3: pointers in input array to */
/* start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, */
/* 2 forward, 3 forward */
/* area: total size of out sets */
/* stride: separation between out data sets (stride=&out-&outm1) */
/* n: actual number of data sets being convolved */
/* npad: size of padded region needed for "leakage" of foward filters */
/* stridel: separation between in data sets (stride=&in-&inm1) */
/* fptl, lptl: first and last data components actually present in input */
/* tr: reverse-time filter coefficients */
/* tf: forward-time filter coefficients */
/* eflag: mod(,2) of present line (tells if special ghost handling needed) */
/* zcflag: =1 means we can count on zero ghosts and need special handling */
/* Subroutine */ int pconvsetsi3new_(double *out, double *in,
double *inm3, double *inm2, double *inm1, double *
inp1, double *inp2, double *inp3, int *area, int *
stride, int *n, int *npad, int *stridel, int *fptl,
int *lptl, double *filt, int *eflag, int *zcflag)
{
/* System generated locals */
int i__1;
/* Local variables */
int i__;
extern /* Subroutine */ int zeroline3_(double *, double *,
double *, double *, double *, int *, int *,
int *, int *, int *, int *, double *,
double *, double *, double *);
double t0, t1, t2, t3;
extern /* Subroutine */ int lineio3fullacc_(double *, double *,
double *, double *, double *, double *,
double *, double *, int *, int *, int *,
int *, int *, int *, double *, double *,
double *, double *, double *, double *,
double *), lineio3full_(double *, double *,
double *, double *, double *, double *,
double *, double *, int *, int *, int *,
int *, int *, int *, double *, double *,
double *, double *, double *, double *,
double *);
/* Filter coefficients */
/* Parameter adjustments */
--filt;
--inp3;
--inp2;
--inp1;
--inm1;
--inm2;
--inm3;
--in;
--out;
/* Function Body */
t0 = filt[1];
t1 = filt[2];
t2 = filt[3];
t3 = filt[4];
i__1 = *npad - 2;
for (i__ = 0; i__ <= i__1; ++i__) {
zeroline3_(&out[1], &in[1], &inm1[1], &inm2[1], &inm3[1], &i__,
stridel, stride, fptl, lptl, area, &t0, &t1, &t2, &t3);
}
/* For the UNUSED lines below, we point to "inm3" or "inp3" to avoid */
/* seg faults because we know "inm3" or "inp3" will always be in range */
/* (unless the box is ridiculously small). Although multiplied */
/* by zero, the code still accesses the dummy data! */
if (*eflag == 0 && *zcflag == 1) {
lineio3full_(&out[1], &inp3[1], &inp3[1], &inp3[1], &inp3[1], &inp3[1]
, &inp3[1], &inp3[1], &c__0, stridel, stride, fptl, lptl,
area, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &t3);
lineio3full_(&out[1], &inp3[1], &inp3[1], &inp3[1], &in[1], &inp1[1],
&inp2[1], &inp3[1], &c__2, stridel, stride, fptl, lptl, area,
&c_b16, &c_b16, &c_b16, &t0, &t1, &t2, &t3);
lineio3full_(&out[1], &inp3[1], &inp3[1], &inm1[1], &in[1], &inp1[1],
&inp2[1], &inp3[1], &c__4, stridel, stride, fptl, lptl, area,
&c_b16, &t2, &t1, &t0, &t1, &t2, &t3);
i__1 = *n - 8;
for (i__ = 6; i__ <= i__1; i__ += 2) {
lineio3full_(&out[1], &inm3[1], &inm2[1], &inm1[1], &in[1], &inp1[
1], &inp2[1], &inp3[1], &i__, stridel, stride, fptl, lptl,
area, &t3, &t2, &t1, &t0, &t1, &t2, &t3);
}
i__1 = *n - 6;
lineio3full_(&out[1], &inm3[1], &inm2[1], &inm1[1], &in[1], &inp1[1],
&inm3[1], &inm3[1], &i__1, stridel, stride, fptl, lptl, area,
&t3, &t2, &t1, &t0, &t1, &c_b16, &c_b16);
i__1 = *n - 4;
lineio3full_(&out[1], &inm3[1], &inm2[1], &inm1[1], &inm3[1], &inm3[1]
, &inm3[1], &inm3[1], &i__1, stridel, stride, fptl, lptl,
area, &t3, &t2, &t1, &c_b16, &c_b16, &c_b16, &c_b16);
i__1 = *n - 2;
lineio3full_(&out[1], &inm3[1], &inm3[1], &inm3[1], &inm3[1], &inm3[1]
, &inm3[1], &inm3[1], &i__1, stridel, stride, fptl, lptl,
area, &t3, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16);
} else {
lineio3fullacc_(&out[1], &inp3[1], &inp3[1], &inp3[1], &inp3[1], &
inp3[1], &inp3[1], &inp3[1], &c__0, stridel, stride, fptl,
lptl, area, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
t3);
lineio3fullacc_(&out[1], &inp3[1], &inp3[1], &inp3[1], &in[1], &inp1[
1], &inp2[1], &inp3[1], &c__2, stridel, stride, fptl, lptl,
area, &c_b16, &c_b16, &c_b16, &t0, &t1, &t2, &t3);
lineio3fullacc_(&out[1], &inp3[1], &inp3[1], &inm1[1], &in[1], &inp1[
1], &inp2[1], &inp3[1], &c__4, stridel, stride, fptl, lptl,
area, &c_b16, &t2, &t1, &t0, &t1, &t2, &t3);
i__1 = *n - 8;
for (i__ = 6; i__ <= i__1; i__ += 2) {
lineio3fullacc_(&out[1], &inm3[1], &inm2[1], &inm1[1], &in[1], &
inp1[1], &inp2[1], &inp3[1], &i__, stridel, stride, fptl,
lptl, area, &t3, &t2, &t1, &t0, &t1, &t2, &t3);
}
i__1 = *n - 6;
lineio3fullacc_(&out[1], &inm3[1], &inm2[1], &inm1[1], &in[1], &inp1[
1], &inm3[1], &inm3[1], &i__1, stridel, stride, fptl, lptl,
area, &t3, &t2, &t1, &t0, &t1, &c_b16, &c_b16);
i__1 = *n - 4;
lineio3fullacc_(&out[1], &inm3[1], &inm2[1], &inm1[1], &inm3[1], &
inm3[1], &inm3[1], &inm3[1], &i__1, stridel, stride, fptl,
lptl, area, &t3, &t2, &t1, &c_b16, &c_b16, &c_b16, &c_b16);
i__1 = *n - 2;
lineio3fullacc_(&out[1], &inm3[1], &inm3[1], &inm3[1], &inm3[1], &
inm3[1], &inm3[1], &inm3[1], &i__1, stridel, stride, fptl,
lptl, area, &t3, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16);
}
return 0;
} /* pconvsetsi3new_ */
/* Back and forth, linear convolution of sets of points using padding */
/* computing data ONLY on ghostpoints, AND folding in input of a */
/* different data set (IN PLACE) */
/* out: input/output array for in place computation */
/* outm3,outm2,outm1,outp1,outp2,outp3: pointers in output array to */
/* start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, */
/* 2 forward, 3 forward */
/* in: input/output array for in place computation */
/* inm3,inm2,inm1,inp1,inp2,inp3: pointers in input array to */
/* start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, */
/* 2 forward, 3 forward */
/* area: total size of out sets */
/* stride: separation between out data sets (stride=&out-&outm1) */
/* n: actual number of data sets being convolved */
/* npad: size of padded region needed for "leakage" of foward filters */
/* stridel: separation between in data sets (stride=&in-&inm1) */
/* fptl, lptl: first and last data components actually present in input */
/* tr: reverse-time filter coefficients */
/* tf: forward-time filter coefficients */
/* eflag: mod(,2) of present line (tells if special ghost handling needed) */
/* zcflag: =1 means we can count on zero ghosts and need special handling */
/* Subroutine */ int pconvsetsi8new_(double *out, double *in,
double *inm8, double *inm7, double *inm6, double *
inm5, double *inm4, double *inm3, double *inm2,
double *inm1, double *inp1, double *inp2, double *
inp3, double *inp4, double *inp5, double *inp6,
double *inp7, double *inp8, int *area, int *stride,
int *n, int *npad, int *stridel, int *fptl, int *
lptl, double *filt, int *eflag, int *zcflag)
{
/* System generated locals */
int i__1;
double d__1, d__2, d__3, d__4, d__5, d__6, d__7, d__8;
/* Local variables */
int i__;
double t0, t1, t2, t3, t4, t5, t6, t7, t8;
extern /* Subroutine */ int lineio8fullacc_(double *, double *,
double *, double *, double *, double *,
double *, double *, double *, double *,
double *, double *, double *, double *,
double *, double *, double *, double *, int *,
int *, int *, int *, int *, int *,
double *, double *, double *, double *,
double *, double *, double *, double *,
double *, double *, double *, double *,
double *, double *, double *, double *,
double *), lineio8full_(double *, double *,
double *, double *, double *, double *,
double *, double *, double *, double *,
double *, double *, double *, double *,
double *, double *, double *, double *, int *,
int *, int *, int *, int *, int *,
double *, double *, double *, double *,
double *, double *, double *, double *,
double *, double *, double *, double *,
double *, double *, double *, double *,
double *);
/* Filter coefficients */
/* Parameter adjustments */
--filt;
--inp8;
--inp7;
--inp6;
--inp5;
--inp4;
--inp3;
--inp2;
--inp1;
--inm1;
--inm2;
--inm3;
--inm4;
--inm5;
--inm6;
--inm7;
--inm8;
--in;
--out;
/* Function Body */
t0 = filt[1];
t1 = filt[2];
t2 = filt[3];
t3 = filt[4];
t4 = filt[5];
t5 = filt[6];
t6 = filt[7];
t7 = filt[8];
t8 = filt[9];
/* For the UNUSED lines below, we point to "inm3" or "inp3" to avoid */
/* seg faults because we know "inm3" or "inp3" will always be in range */
/* (unless the box is ridiculously small). Although multiplied */
/* by zero, the code still accesses the dummy data! */
if (t0 != (double)0.) {
/* symmetric filter */
if (FALSE_) {
/* if (eflag.eq.0 .and. zcflag.eq.1) then */
lineio8full_(&out[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &
inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1],
&inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp8[
1], &inp8[1], &c__0, stridel, stride, fptl, lptl, area, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16, &c_b16, &t3);
lineio8full_(&out[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &
inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1],
&inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp6[1], &inp7[
1], &inp8[1], &c__2, stridel, stride, fptl, lptl, area, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
t6, &t7, &t8);
lineio8full_(&out[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &
inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1],
&inp8[1], &inp8[1], &inp4[1], &inp5[1], &inp6[1], &inp7[
1], &inp8[1], &c__4, stridel, stride, fptl, lptl, area, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &t4, &t5, &t6, &t7,
&t8);
lineio8full_(&out[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &
inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1],
&inp2[1], &inp3[1], &inp4[1], &inp5[1], &inp6[1], &inp7[
1], &inp8[1], &c__6, stridel, stride, fptl, lptl, area, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16, &c_b16, &c_b16, &t2, &t3, &t4, &t5, &t6, &t7, &t8);
lineio8full_(&out[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &
inp8[1], &inp8[1], &inp8[1], &inp8[1], &in[1], &inp1[1], &
inp2[1], &inp3[1], &inp4[1], &inp5[1], &inp6[1], &inp7[1],
&inp8[1], &c__8, stridel, stride, fptl, lptl, area, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16, &t0, &t1, &t2, &t3, &t4, &t5, &t6, &t7, &t8);
lineio8full_(&out[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &
inp8[1], &inp8[1], &inm2[1], &inm1[1], &in[1], &inp1[1], &
inp2[1], &inp3[1], &inp4[1], &inp5[1], &inp6[1], &inp7[1],
&inp8[1], &c__10, stridel, stride, fptl, lptl, area, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &t2, &t1, &
t0, &t1, &t2, &t3, &t4, &t5, &t6, &t7, &t8);
lineio8full_(&out[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &
inm4[1], &inm3[1], &inm2[1], &inm1[1], &in[1], &inp1[1], &
inp2[1], &inp3[1], &inp4[1], &inp5[1], &inp6[1], &inp7[1],
&inp8[1], &c__12, stridel, stride, fptl, lptl, area, &
c_b16, &c_b16, &c_b16, &c_b16, &t4, &t3, &t2, &t1, &t0, &
t1, &t2, &t3, &t4, &t5, &t6, &t7, &t8);
lineio8full_(&out[1], &inp8[1], &inp8[1], &inm6[1], &inm5[1], &
inm4[1], &inm3[1], &inm2[1], &inm1[1], &in[1], &inp1[1], &
inp2[1], &inp3[1], &inp4[1], &inp5[1], &inp6[1], &inp7[1],
&inp8[1], &c__14, stridel, stride, fptl, lptl, area, &
c_b16, &c_b16, &t6, &t5, &t4, &t3, &t2, &t1, &t0, &t1, &
t2, &t3, &t4, &t5, &t6, &t7, &t8);
i__1 = *n - 18;
for (i__ = 16; i__ <= i__1; i__ += 2) {
lineio8full_(&out[1], &inm8[1], &inm7[1], &inm6[1], &inm5[1],
&inm4[1], &inm3[1], &inm2[1], &inm1[1], &in[1], &inp1[
1], &inp2[1], &inp3[1], &inp4[1], &inp5[1], &inp6[1],
&inp7[1], &inp8[1], &i__, stridel, stride, fptl, lptl,
area, &t8, &t7, &t6, &t5, &t4, &t3, &t2, &t1, &t0, &
t1, &t2, &t3, &t4, &t5, &t6, &t7, &t8);
}
i__1 = *n - 16;
lineio8full_(&out[1], &inm8[1], &inm7[1], &inm6[1], &inm5[1], &
inm4[1], &inm3[1], &inm2[1], &inm1[1], &in[1], &inp1[1], &
inp2[1], &inp3[1], &inp4[1], &inp5[1], &inp6[1], &inp7[1],
&inm8[1], &i__1, stridel, stride, fptl, lptl, area, &t8,
&t7, &t6, &t5, &t4, &t3, &t2, &t1, &t0, &t1, &t2, &t3, &
t4, &t5, &t6, &t7, &c_b16);
i__1 = *n - 14;
lineio8full_(&out[1], &inm8[1], &inm7[1], &inm6[1], &inm5[1], &
inm4[1], &inm3[1], &inm2[1], &inm1[1], &in[1], &inp1[1], &
inp2[1], &inp3[1], &inp4[1], &inp5[1], &inm8[1], &inm8[1],
&inm8[1], &i__1, stridel, stride, fptl, lptl, area, &t8,
&t7, &t6, &t5, &t4, &t3, &t2, &t1, &t0, &t1, &t2, &t3, &
t4, &t5, &c_b16, &c_b16, &c_b16);
i__1 = *n - 12;
lineio8full_(&out[1], &inm8[1], &inm7[1], &inm6[1], &inm5[1], &
inm4[1], &inm3[1], &inm2[1], &inm1[1], &in[1], &inp1[1], &
inp2[1], &inp3[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1],
&inm8[1], &i__1, stridel, stride, fptl, lptl, area, &t8,
&t7, &t6, &t5, &t4, &t3, &t2, &t1, &t0, &t1, &t2, &t3, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16);
i__1 = *n - 10;
lineio8full_(&out[1], &inm8[1], &inm7[1], &inm6[1], &inm5[1], &
inm4[1], &inm3[1], &inm2[1], &inm1[1], &in[1], &inp1[1], &
inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1],
&inm8[1], &i__1, stridel, stride, fptl, lptl, area, &t8,
&t7, &t6, &t5, &t4, &t3, &t2, &t1, &t0, &t1, &c_b16, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16);
i__1 = *n - 8;
lineio8full_(&out[1], &inm8[1], &inm7[1], &inm6[1], &inm5[1], &
inm4[1], &inm3[1], &inm2[1], &inm1[1], &inm8[1], &inm8[1],
&inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[
1], &inm8[1], &i__1, stridel, stride, fptl, lptl, area, &
t8, &t7, &t6, &t5, &t4, &t3, &t2, &t1, &c_b16, &c_b16, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16);
i__1 = *n - 6;
lineio8full_(&out[1], &inm8[1], &inm7[1], &inm6[1], &inm5[1], &
inm4[1], &inm3[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1],
&inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[
1], &inm8[1], &i__1, stridel, stride, fptl, lptl, area, &
t8, &t7, &t6, &t5, &t4, &t3, &c_b16, &c_b16, &c_b16, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16);
i__1 = *n - 4;
lineio8full_(&out[1], &inm8[1], &inm7[1], &inm6[1], &inm5[1], &
inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1],
&inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[
1], &inm8[1], &i__1, stridel, stride, fptl, lptl, area, &
t8, &t7, &t6, &t5, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16,
&c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16);
i__1 = *n - 2;
lineio8full_(&out[1], &inm8[1], &inm7[1], &inm8[1], &inm8[1], &
inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1],
&inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[
1], &inm8[1], &i__1, stridel, stride, fptl, lptl, area, &
t8, &t7, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16, &c_b16);
} else {
lineio8fullacc_(&out[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &
inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1],
&inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp8[
1], &inp8[1], &c__0, stridel, stride, fptl, lptl, area, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16, &c_b16, &t8);
lineio8fullacc_(&out[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &
inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1],
&inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp6[1], &inp7[
1], &inp8[1], &c__2, stridel, stride, fptl, lptl, area, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
t6, &t7, &t8);
lineio8fullacc_(&out[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &
inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1],
&inp8[1], &inp8[1], &inp4[1], &inp5[1], &inp6[1], &inp7[
1], &inp8[1], &c__4, stridel, stride, fptl, lptl, area, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &t4, &t5, &t6, &t7,
&t8);
lineio8fullacc_(&out[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &
inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1],
&inp2[1], &inp3[1], &inp4[1], &inp5[1], &inp6[1], &inp7[
1], &inp8[1], &c__6, stridel, stride, fptl, lptl, area, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16, &c_b16, &c_b16, &t2, &t3, &t4, &t5, &t6, &t7, &t8);
lineio8fullacc_(&out[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &
inp8[1], &inp8[1], &inp8[1], &inp8[1], &in[1], &inp1[1], &
inp2[1], &inp3[1], &inp4[1], &inp5[1], &inp6[1], &inp7[1],
&inp8[1], &c__8, stridel, stride, fptl, lptl, area, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16, &t0, &t1, &t2, &t3, &t4, &t5, &t6, &t7, &t8);
lineio8fullacc_(&out[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &
inp8[1], &inp8[1], &inm2[1], &inm1[1], &in[1], &inp1[1], &
inp2[1], &inp3[1], &inp4[1], &inp5[1], &inp6[1], &inp7[1],
&inp8[1], &c__10, stridel, stride, fptl, lptl, area, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &t2, &t1, &
t0, &t1, &t2, &t3, &t4, &t5, &t6, &t7, &t8);
lineio8fullacc_(&out[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &
inm4[1], &inm3[1], &inm2[1], &inm1[1], &in[1], &inp1[1], &
inp2[1], &inp3[1], &inp4[1], &inp5[1], &inp6[1], &inp7[1],
&inp8[1], &c__12, stridel, stride, fptl, lptl, area, &
c_b16, &c_b16, &c_b16, &c_b16, &t4, &t3, &t2, &t1, &t0, &
t1, &t2, &t3, &t4, &t5, &t6, &t7, &t8);
lineio8fullacc_(&out[1], &inp8[1], &inp8[1], &inm6[1], &inm5[1], &
inm4[1], &inm3[1], &inm2[1], &inm1[1], &in[1], &inp1[1], &
inp2[1], &inp3[1], &inp4[1], &inp5[1], &inp6[1], &inp7[1],
&inp8[1], &c__14, stridel, stride, fptl, lptl, area, &
c_b16, &c_b16, &t6, &t5, &t4, &t3, &t2, &t1, &t0, &t1, &
t2, &t3, &t4, &t5, &t6, &t7, &t8);
i__1 = *n - 18;
for (i__ = 16; i__ <= i__1; i__ += 2) {
lineio8fullacc_(&out[1], &inm8[1], &inm7[1], &inm6[1], &inm5[
1], &inm4[1], &inm3[1], &inm2[1], &inm1[1], &in[1], &
inp1[1], &inp2[1], &inp3[1], &inp4[1], &inp5[1], &
inp6[1], &inp7[1], &inp8[1], &i__, stridel, stride,
fptl, lptl, area, &t8, &t7, &t6, &t5, &t4, &t3, &t2, &
t1, &t0, &t1, &t2, &t3, &t4, &t5, &t6, &t7, &t8);
}
i__1 = *n - 16;
lineio8fullacc_(&out[1], &inm8[1], &inm7[1], &inm6[1], &inm5[1], &
inm4[1], &inm3[1], &inm2[1], &inm1[1], &in[1], &inp1[1], &
inp2[1], &inp3[1], &inp4[1], &inp5[1], &inp6[1], &inp7[1],
&inm8[1], &i__1, stridel, stride, fptl, lptl, area, &t8,
&t7, &t6, &t5, &t4, &t3, &t2, &t1, &t0, &t1, &t2, &t3, &
t4, &t5, &t6, &t7, &c_b16);
i__1 = *n - 14;
lineio8fullacc_(&out[1], &inm8[1], &inm7[1], &inm6[1], &inm5[1], &
inm4[1], &inm3[1], &inm2[1], &inm1[1], &in[1], &inp1[1], &
inp2[1], &inp3[1], &inp4[1], &inp5[1], &inm8[1], &inm8[1],
&inm8[1], &i__1, stridel, stride, fptl, lptl, area, &t8,
&t7, &t6, &t5, &t4, &t3, &t2, &t1, &t0, &t1, &t2, &t3, &
t4, &t5, &c_b16, &c_b16, &c_b16);
i__1 = *n - 12;
lineio8fullacc_(&out[1], &inm8[1], &inm7[1], &inm6[1], &inm5[1], &
inm4[1], &inm3[1], &inm2[1], &inm1[1], &in[1], &inp1[1], &
inp2[1], &inp3[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1],
&inm8[1], &i__1, stridel, stride, fptl, lptl, area, &t8,
&t7, &t6, &t5, &t4, &t3, &t2, &t1, &t0, &t1, &t2, &t3, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16);
i__1 = *n - 10;
lineio8fullacc_(&out[1], &inm8[1], &inm7[1], &inm6[1], &inm5[1], &
inm4[1], &inm3[1], &inm2[1], &inm1[1], &in[1], &inp1[1], &
inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1],
&inm8[1], &i__1, stridel, stride, fptl, lptl, area, &t8,
&t7, &t6, &t5, &t4, &t3, &t2, &t1, &t0, &t1, &c_b16, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16);
i__1 = *n - 8;
lineio8fullacc_(&out[1], &inm8[1], &inm7[1], &inm6[1], &inm5[1], &
inm4[1], &inm3[1], &inm2[1], &inm1[1], &inm8[1], &inm8[1],
&inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[
1], &inm8[1], &i__1, stridel, stride, fptl, lptl, area, &
t8, &t7, &t6, &t5, &t4, &t3, &t2, &t1, &c_b16, &c_b16, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16);
i__1 = *n - 6;
lineio8fullacc_(&out[1], &inm8[1], &inm7[1], &inm6[1], &inm5[1], &
inm4[1], &inm3[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1],
&inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[
1], &inm8[1], &i__1, stridel, stride, fptl, lptl, area, &
t8, &t7, &t6, &t5, &t4, &t3, &c_b16, &c_b16, &c_b16, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16);
i__1 = *n - 4;
lineio8fullacc_(&out[1], &inm8[1], &inm7[1], &inm6[1], &inm5[1], &
inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1],
&inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[
1], &inm8[1], &i__1, stridel, stride, fptl, lptl, area, &
t8, &t7, &t6, &t5, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16,
&c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16);
i__1 = *n - 2;
lineio8fullacc_(&out[1], &inm8[1], &inm7[1], &inm8[1], &inm8[1], &
inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1],
&inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[
1], &inm8[1], &i__1, stridel, stride, fptl, lptl, area, &
t8, &t7, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16, &c_b16);
}
/* end of symmetric part */
} else {
/* antisymmetric filter */
if (FALSE_) {
/* if (eflag.eq.0 .and. zcflag.eq.1) then */
lineio8full_(&out[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &
inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1],
&inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp8[
1], &inp8[1], &c__0, stridel, stride, fptl, lptl, area, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16, &c_b16, &t3);
lineio8full_(&out[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &
inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1],
&inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp6[1], &inp7[
1], &inp8[1], &c__2, stridel, stride, fptl, lptl, area, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
t6, &t7, &t8);
lineio8full_(&out[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &
inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1],
&inp8[1], &inp8[1], &inp4[1], &inp5[1], &inp6[1], &inp7[
1], &inp8[1], &c__4, stridel, stride, fptl, lptl, area, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &t4, &t5, &t6, &t7,
&t8);
lineio8full_(&out[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &
inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1],
&inp2[1], &inp3[1], &inp4[1], &inp5[1], &inp6[1], &inp7[
1], &inp8[1], &c__6, stridel, stride, fptl, lptl, area, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16, &c_b16, &c_b16, &t2, &t3, &t4, &t5, &t6, &t7, &t8);
lineio8full_(&out[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &
inp8[1], &inp8[1], &inp8[1], &inp8[1], &in[1], &inp1[1], &
inp2[1], &inp3[1], &inp4[1], &inp5[1], &inp6[1], &inp7[1],
&inp8[1], &c__8, stridel, stride, fptl, lptl, area, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16, &t0, &t1, &t2, &t3, &t4, &t5, &t6, &t7, &t8);
lineio8full_(&out[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &
inp8[1], &inp8[1], &inm2[1], &inm1[1], &in[1], &inp1[1], &
inp2[1], &inp3[1], &inp4[1], &inp5[1], &inp6[1], &inp7[1],
&inp8[1], &c__10, stridel, stride, fptl, lptl, area, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &t2, &t1, &
t0, &t1, &t2, &t3, &t4, &t5, &t6, &t7, &t8);
lineio8full_(&out[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &
inm4[1], &inm3[1], &inm2[1], &inm1[1], &in[1], &inp1[1], &
inp2[1], &inp3[1], &inp4[1], &inp5[1], &inp6[1], &inp7[1],
&inp8[1], &c__12, stridel, stride, fptl, lptl, area, &
c_b16, &c_b16, &c_b16, &c_b16, &t4, &t3, &t2, &t1, &t0, &
t1, &t2, &t3, &t4, &t5, &t6, &t7, &t8);
lineio8full_(&out[1], &inp8[1], &inp8[1], &inm6[1], &inm5[1], &
inm4[1], &inm3[1], &inm2[1], &inm1[1], &in[1], &inp1[1], &
inp2[1], &inp3[1], &inp4[1], &inp5[1], &inp6[1], &inp7[1],
&inp8[1], &c__14, stridel, stride, fptl, lptl, area, &
c_b16, &c_b16, &t6, &t5, &t4, &t3, &t2, &t1, &t0, &t1, &
t2, &t3, &t4, &t5, &t6, &t7, &t8);
i__1 = *n - 18;
for (i__ = 16; i__ <= i__1; i__ += 2) {
lineio8full_(&out[1], &inm8[1], &inm7[1], &inm6[1], &inm5[1],
&inm4[1], &inm3[1], &inm2[1], &inm1[1], &in[1], &inp1[
1], &inp2[1], &inp3[1], &inp4[1], &inp5[1], &inp6[1],
&inp7[1], &inp8[1], &i__, stridel, stride, fptl, lptl,
area, &t8, &t7, &t6, &t5, &t4, &t3, &t2, &t1, &t0, &
t1, &t2, &t3, &t4, &t5, &t6, &t7, &t8);
}
i__1 = *n - 16;
lineio8full_(&out[1], &inm8[1], &inm7[1], &inm6[1], &inm5[1], &
inm4[1], &inm3[1], &inm2[1], &inm1[1], &in[1], &inp1[1], &
inp2[1], &inp3[1], &inp4[1], &inp5[1], &inp6[1], &inp7[1],
&inm8[1], &i__1, stridel, stride, fptl, lptl, area, &t8,
&t7, &t6, &t5, &t4, &t3, &t2, &t1, &t0, &t1, &t2, &t3, &
t4, &t5, &t6, &t7, &c_b16);
i__1 = *n - 14;
lineio8full_(&out[1], &inm8[1], &inm7[1], &inm6[1], &inm5[1], &
inm4[1], &inm3[1], &inm2[1], &inm1[1], &in[1], &inp1[1], &
inp2[1], &inp3[1], &inp4[1], &inp5[1], &inm8[1], &inm8[1],
&inm8[1], &i__1, stridel, stride, fptl, lptl, area, &t8,
&t7, &t6, &t5, &t4, &t3, &t2, &t1, &t0, &t1, &t2, &t3, &
t4, &t5, &c_b16, &c_b16, &c_b16);
i__1 = *n - 12;
lineio8full_(&out[1], &inm8[1], &inm7[1], &inm6[1], &inm5[1], &
inm4[1], &inm3[1], &inm2[1], &inm1[1], &in[1], &inp1[1], &
inp2[1], &inp3[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1],
&inm8[1], &i__1, stridel, stride, fptl, lptl, area, &t8,
&t7, &t6, &t5, &t4, &t3, &t2, &t1, &t0, &t1, &t2, &t3, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16);
i__1 = *n - 10;
lineio8full_(&out[1], &inm8[1], &inm7[1], &inm6[1], &inm5[1], &
inm4[1], &inm3[1], &inm2[1], &inm1[1], &in[1], &inp1[1], &
inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1],
&inm8[1], &i__1, stridel, stride, fptl, lptl, area, &t8,
&t7, &t6, &t5, &t4, &t3, &t2, &t1, &t0, &t1, &c_b16, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16);
i__1 = *n - 8;
lineio8full_(&out[1], &inm8[1], &inm7[1], &inm6[1], &inm5[1], &
inm4[1], &inm3[1], &inm2[1], &inm1[1], &inm8[1], &inm8[1],
&inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[
1], &inm8[1], &i__1, stridel, stride, fptl, lptl, area, &
t8, &t7, &t6, &t5, &t4, &t3, &t2, &t1, &c_b16, &c_b16, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16);
i__1 = *n - 6;
lineio8full_(&out[1], &inm8[1], &inm7[1], &inm6[1], &inm5[1], &
inm4[1], &inm3[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1],
&inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[
1], &inm8[1], &i__1, stridel, stride, fptl, lptl, area, &
t8, &t7, &t6, &t5, &t4, &t3, &c_b16, &c_b16, &c_b16, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16);
i__1 = *n - 4;
lineio8full_(&out[1], &inm8[1], &inm7[1], &inm6[1], &inm5[1], &
inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1],
&inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[
1], &inm8[1], &i__1, stridel, stride, fptl, lptl, area, &
t8, &t7, &t6, &t5, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16,
&c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16);
i__1 = *n - 2;
lineio8full_(&out[1], &inm8[1], &inm7[1], &inm8[1], &inm8[1], &
inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1],
&inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[
1], &inm8[1], &i__1, stridel, stride, fptl, lptl, area, &
t8, &t7, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16, &c_b16);
/* Is any of the above code ever called? Is 2 ever greater than 3? */
} else {
lineio8fullacc_(&out[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &
inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1],
&inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp8[
1], &inp8[1], &c__0, stridel, stride, fptl, lptl, area, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16, &c_b16, &t8);
lineio8fullacc_(&out[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &
inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1],
&inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp6[1], &inp7[
1], &inp8[1], &c__2, stridel, stride, fptl, lptl, area, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
t6, &t7, &t8);
lineio8fullacc_(&out[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &
inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1],
&inp8[1], &inp8[1], &inp4[1], &inp5[1], &inp6[1], &inp7[
1], &inp8[1], &c__4, stridel, stride, fptl, lptl, area, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &t4, &t5, &t6, &t7,
&t8);
lineio8fullacc_(&out[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &
inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1],
&inp2[1], &inp3[1], &inp4[1], &inp5[1], &inp6[1], &inp7[
1], &inp8[1], &c__6, stridel, stride, fptl, lptl, area, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16, &c_b16, &c_b16, &t2, &t3, &t4, &t5, &t6, &t7, &t8);
lineio8fullacc_(&out[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &
inp8[1], &inp8[1], &inp8[1], &inp8[1], &in[1], &inp1[1], &
inp2[1], &inp3[1], &inp4[1], &inp5[1], &inp6[1], &inp7[1],
&inp8[1], &c__8, stridel, stride, fptl, lptl, area, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16, &t0, &t1, &t2, &t3, &t4, &t5, &t6, &t7, &t8);
d__1 = t2 * (double)-1.;
d__2 = t1 * (double)-1.;
lineio8fullacc_(&out[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &
inp8[1], &inp8[1], &inm2[1], &inm1[1], &in[1], &inp1[1], &
inp2[1], &inp3[1], &inp4[1], &inp5[1], &inp6[1], &inp7[1],
&inp8[1], &c__10, stridel, stride, fptl, lptl, area, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &d__1, &
d__2, &t0, &t1, &t2, &t3, &t4, &t5, &t6, &t7, &t8);
d__1 = t4 * (double)-1.;
d__2 = t3 * (double)-1.;
d__3 = t2 * (double)-1.;
d__4 = t1 * (double)-1.;
lineio8fullacc_(&out[1], &inp8[1], &inp8[1], &inp8[1], &inp8[1], &
inm4[1], &inm3[1], &inm2[1], &inm1[1], &in[1], &inp1[1], &
inp2[1], &inp3[1], &inp4[1], &inp5[1], &inp6[1], &inp7[1],
&inp8[1], &c__12, stridel, stride, fptl, lptl, area, &
c_b16, &c_b16, &c_b16, &c_b16, &d__1, &d__2, &d__3, &d__4,
&t0, &t1, &t2, &t3, &t4, &t5, &t6, &t7, &t8);
d__1 = t6 * (double)-1.;
d__2 = t5 * (double)-1.;
d__3 = t4 * (double)-1.;
d__4 = t3 * (double)-1.;
d__5 = t2 * (double)-1.;
d__6 = t1 * (double)-1.;
lineio8fullacc_(&out[1], &inp8[1], &inp8[1], &inm6[1], &inm5[1], &
inm4[1], &inm3[1], &inm2[1], &inm1[1], &in[1], &inp1[1], &
inp2[1], &inp3[1], &inp4[1], &inp5[1], &inp6[1], &inp7[1],
&inp8[1], &c__14, stridel, stride, fptl, lptl, area, &
c_b16, &c_b16, &d__1, &d__2, &d__3, &d__4, &d__5, &d__6, &
t0, &t1, &t2, &t3, &t4, &t5, &t6, &t7, &t8);
i__1 = *n - 18;
for (i__ = 16; i__ <= i__1; i__ += 2) {
d__1 = t8 * (double)-1.;
d__2 = t7 * (double)-1.;
d__3 = t6 * (double)-1.;
d__4 = t5 * (double)-1.;
d__5 = t4 * (double)-1.;
d__6 = t3 * (double)-1.;
d__7 = t2 * (double)-1.;
d__8 = t1 * (double)-1.;
lineio8fullacc_(&out[1], &inm8[1], &inm7[1], &inm6[1], &inm5[
1], &inm4[1], &inm3[1], &inm2[1], &inm1[1], &in[1], &
inp1[1], &inp2[1], &inp3[1], &inp4[1], &inp5[1], &
inp6[1], &inp7[1], &inp8[1], &i__, stridel, stride,
fptl, lptl, area, &d__1, &d__2, &d__3, &d__4, &d__5, &
d__6, &d__7, &d__8, &t0, &t1, &t2, &t3, &t4, &t5, &t6,
&t7, &t8);
}
i__1 = *n - 16;
d__1 = t8 * (double)-1.;
d__2 = t7 * (double)-1.;
d__3 = t6 * (double)-1.;
d__4 = t5 * (double)-1.;
d__5 = t4 * (double)-1.;
d__6 = t3 * (double)-1.;
d__7 = t2 * (double)-1.;
d__8 = t1 * (double)-1.;
lineio8fullacc_(&out[1], &inm8[1], &inm7[1], &inm6[1], &inm5[1], &
inm4[1], &inm3[1], &inm2[1], &inm1[1], &in[1], &inp1[1], &
inp2[1], &inp3[1], &inp4[1], &inp5[1], &inp6[1], &inp7[1],
&inm8[1], &i__1, stridel, stride, fptl, lptl, area, &
d__1, &d__2, &d__3, &d__4, &d__5, &d__6, &d__7, &d__8, &
t0, &t1, &t2, &t3, &t4, &t5, &t6, &t7, &c_b16);
i__1 = *n - 14;
d__1 = t8 * (double)-1.;
d__2 = t7 * (double)-1.;
d__3 = t6 * (double)-1.;
d__4 = t5 * (double)-1.;
d__5 = t4 * (double)-1.;
d__6 = t3 * (double)-1.;
d__7 = t2 * (double)-1.;
d__8 = t1 * (double)-1.;
lineio8fullacc_(&out[1], &inm8[1], &inm7[1], &inm6[1], &inm5[1], &
inm4[1], &inm3[1], &inm2[1], &inm1[1], &in[1], &inp1[1], &
inp2[1], &inp3[1], &inp4[1], &inp5[1], &inm8[1], &inm8[1],
&inm8[1], &i__1, stridel, stride, fptl, lptl, area, &
d__1, &d__2, &d__3, &d__4, &d__5, &d__6, &d__7, &d__8, &
t0, &t1, &t2, &t3, &t4, &t5, &c_b16, &c_b16, &c_b16);
i__1 = *n - 12;
d__1 = t8 * (double)-1.;
d__2 = t7 * (double)-1.;
d__3 = t6 * (double)-1.;
d__4 = t5 * (double)-1.;
d__5 = t4 * (double)-1.;
d__6 = t3 * (double)-1.;
d__7 = t2 * (double)-1.;
d__8 = t1 * (double)-1.;
lineio8fullacc_(&out[1], &inm8[1], &inm7[1], &inm6[1], &inm5[1], &
inm4[1], &inm3[1], &inm2[1], &inm1[1], &in[1], &inp1[1], &
inp2[1], &inp3[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1],
&inm8[1], &i__1, stridel, stride, fptl, lptl, area, &
d__1, &d__2, &d__3, &d__4, &d__5, &d__6, &d__7, &d__8, &
t0, &t1, &t2, &t3, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16)
;
i__1 = *n - 10;
d__1 = t8 * (double)-1.;
d__2 = t7 * (double)-1.;
d__3 = t6 * (double)-1.;
d__4 = t5 * (double)-1.;
d__5 = t4 * (double)-1.;
d__6 = t3 * (double)-1.;
d__7 = t2 * (double)-1.;
d__8 = t1 * (double)-1.;
lineio8fullacc_(&out[1], &inm8[1], &inm7[1], &inm6[1], &inm5[1], &
inm4[1], &inm3[1], &inm2[1], &inm1[1], &in[1], &inp1[1], &
inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1],
&inm8[1], &i__1, stridel, stride, fptl, lptl, area, &
d__1, &d__2, &d__3, &d__4, &d__5, &d__6, &d__7, &d__8, &
t0, &t1, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16);
i__1 = *n - 8;
d__1 = t8 * (double)-1.;
d__2 = t7 * (double)-1.;
d__3 = t6 * (double)-1.;
d__4 = t5 * (double)-1.;
d__5 = t4 * (double)-1.;
d__6 = t3 * (double)-1.;
d__7 = t2 * (double)-1.;
d__8 = t1 * (double)-1.;
lineio8fullacc_(&out[1], &inm8[1], &inm7[1], &inm6[1], &inm5[1], &
inm4[1], &inm3[1], &inm2[1], &inm1[1], &inm8[1], &inm8[1],
&inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[
1], &inm8[1], &i__1, stridel, stride, fptl, lptl, area, &
d__1, &d__2, &d__3, &d__4, &d__5, &d__6, &d__7, &d__8, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16, &c_b16);
i__1 = *n - 6;
d__1 = t8 * (double)-1.;
d__2 = t7 * (double)-1.;
d__3 = t6 * (double)-1.;
d__4 = t5 * (double)-1.;
d__5 = t4 * (double)-1.;
d__6 = t3 * (double)-1.;
lineio8fullacc_(&out[1], &inm8[1], &inm7[1], &inm6[1], &inm5[1], &
inm4[1], &inm3[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1],
&inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[
1], &inm8[1], &i__1, stridel, stride, fptl, lptl, area, &
d__1, &d__2, &d__3, &d__4, &d__5, &d__6, &c_b16, &c_b16, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16, &c_b16);
i__1 = *n - 4;
d__1 = t8 * (double)-1.;
d__2 = t7 * (double)-1.;
d__3 = t6 * (double)-1.;
d__4 = t5 * (double)-1.;
lineio8fullacc_(&out[1], &inm8[1], &inm7[1], &inm6[1], &inm5[1], &
inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1],
&inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[
1], &inm8[1], &i__1, stridel, stride, fptl, lptl, area, &
d__1, &d__2, &d__3, &d__4, &c_b16, &c_b16, &c_b16, &c_b16,
&c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16, &c_b16);
i__1 = *n - 2;
d__1 = t8 * (double)-1.;
d__2 = t7 * (double)-1.;
lineio8fullacc_(&out[1], &inm8[1], &inm7[1], &inm8[1], &inm8[1], &
inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1],
&inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[
1], &inm8[1], &i__1, stridel, stride, fptl, lptl, area, &
d__1, &d__2, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16, &c_b16, &c_b16);
}
/* end of antisymmetric part */
}
return 0;
} /* pconvsetsi8new_ */
/* Back and forth, linear convolution of sets of points using padding */
/* computing data ONLY on ghostpoints, AND folding in input of a */
/* different data set (IN PLACE) */
/* out: input/output array for in place computation */
/* outm3,outm2,outm1,outp1,outp2,outp3: pointers in output array to */
/* start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, */
/* 2 forward, 3 forward */
/* in: input/output array for in place computation */
/* inm3,inm2,inm1,inp1,inp2,inp3: pointers in input array to */
/* start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, */
/* 2 forward, 3 forward */
/* area: total size of out sets */
/* stride: separation between out data sets (stride=&out-&outm1) */
/* n: actual number of data sets being convolved */
/* npad: size of padded region needed for "leakage" of foward filters */
/* stridel: separation between in data sets (stride=&in-&inm1) */
/* fptl, lptl: first and last data components actually present in input */
/* tr: reverse-time filter coefficients */
/* tf: forward-time filter coefficients */
/* eflag: mod(,2) of present line (tells if special ghost handling needed) */
/* zcflag: =1 means we can count on zero ghosts and need special handling */
/* Subroutine */ int pconvsetsi8_(double *out, double *outm8,
double *outm7, double *outm6, double *outm5, double *
outm4, double *outm3, double *outm2, double *outm1,
double *outp1, double *outp2, double *outp3, double *
outp4, double *outp5, double *outp6, double *outp7,
double *outp8, double *in, double *inm8, double *inm7,
double *inm6, double *inm5, double *inm4, double *
inm3, double *inm2, double *inm1, double *inp1,
double *inp2, double *inp3, double *inp4, double *
inp5, double *inp6, double *inp7, double *inp8, int *
area, int *stride, int *n, int *npad, int *stridel,
int *fptl, int *lptl, double *tr, double *tf, int
*eflag, int *zcflag)
{
/* System generated locals */
int i__1, i__2;
/* Local variables */
int i__, j;
extern /* Subroutine */ int zeroline8_(double *, double *,
double *, double *, double *, double *,
double *, double *, double *, double *, int *,
int *, int *, int *, int *, int *,
double *, double *, double *, double *,
double *, double *, double *, double *,
double *);
double t0, t1, t2, t3, t4, t5, t6, t7, t8;
extern /* Subroutine */ int lineio8_(double *, double *,
double *, double *, double *, double *,
double *, double *, double *, double *, int *,
int *, int *, int *, int *, int *,
double *, double *, double *, double *,
double *, double *, double *, double *,
double *);
int out2;
/* Reverse phase */
/* Parameter adjustments */
--tf;
--tr;
--inp8;
--inp7;
--inp6;
--inp5;
--inp4;
--inp3;
--inp2;
--inp1;
--inm1;
--inm2;
--inm3;
--inm4;
--inm5;
--inm6;
--inm7;
--inm8;
--in;
--outp8;
--outp7;
--outp6;
--outp5;
--outp4;
--outp3;
--outp2;
--outp1;
--outm1;
--outm2;
--outm3;
--outm4;
--outm5;
--outm6;
--outm7;
--outm8;
--out;
/* Function Body */
t0 = tr[1];
t1 = tr[2];
t2 = tr[3];
t3 = tr[4];
t4 = tr[5];
t5 = tr[6];
t6 = tr[7];
t7 = tr[8];
t8 = tr[9];
/* First, make sure pad is zero for forward phase. */
/* For pretty code the last number would be npad-1, but technically */
/* npad-2 (for extra speed) is all that is needed. */
i__1 = *npad - 2;
for (i__ = *n; i__ <= i__1; ++i__) {
zeroline8_(&out[1], &in[1], &inm1[1], &inm2[1], &inm3[1], &inm4[1], &
inm5[1], &inm6[1], &inm7[1], &inm8[1], &i__, stridel, stride,
fptl, lptl, area, &t0, &t1, &t2, &t3, &t4, &t5, &t6, &t7, &t8)
;
}
/* call zeroline8(out,in,inm1,inm2,inm3, */
/* $ inm4,inm5,inm6,inm7,inm8,n-1, */
/* $ stridel,stride,fptl,lptl,area, */
/* $ t0,t1,t2,t3,t4,t5,t6,t7,t8) */
/* CALL flush(6) */
/* For the UNUSED lines below, we point to "inm8" to avoid seg faults */
/* because we know "inm8" will always be in range (unless the box is */
/* ridiculously small). Although multiplied by zero, the code still */
/* accesses the dummy data! */
i__1 = *n - 1;
lineio8_(&out[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[
1], &inm8[1], &inm8[1], &inm8[1], &i__1, stridel, stride, fptl,
lptl, area, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16, &c_b16, &t8);
i__1 = *n - 2;
lineio8_(&out[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[
1], &inm8[1], &inm7[1], &inm8[1], &i__1, stridel, stride, fptl,
lptl, area, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &
c_b16, &t7, &t8);
i__1 = *n - 3;
lineio8_(&out[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[
1], &inm6[1], &inm7[1], &inm8[1], &i__1, stridel, stride, fptl,
lptl, area, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &t6, &
t7, &t8);
i__1 = *n - 4;
lineio8_(&out[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm5[
1], &inm6[1], &inm7[1], &inm8[1], &i__1, stridel, stride, fptl,
lptl, area, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &t5, &t6, &t7,
&t8);
i__1 = *n - 5;
lineio8_(&out[1], &inm8[1], &inm8[1], &inm8[1], &inm8[1], &inm4[1], &inm5[
1], &inm6[1], &inm7[1], &inm8[1], &i__1, stridel, stride, fptl,
lptl, area, &c_b16, &c_b16, &c_b16, &c_b16, &t4, &t5, &t6, &t7, &
t8);
i__1 = *n - 6;
lineio8_(&out[1], &inm8[1], &inm8[1], &inm8[1], &inm3[1], &inm4[1], &inm5[
1], &inm6[1], &inm7[1], &inm8[1], &i__1, stridel, stride, fptl,
lptl, area, &c_b16, &c_b16, &c_b16, &t3, &t4, &t5, &t6, &t7, &t8);
i__1 = *n - 7;
lineio8_(&out[1], &inm8[1], &inm8[1], &inm2[1], &inm3[1], &inm4[1], &inm5[
1], &inm6[1], &inm7[1], &inm8[1], &i__1, stridel, stride, fptl,
lptl, area, &c_b16, &c_b16, &t2, &t3, &t4, &t5, &t6, &t7, &t8);
i__1 = *n - 8;
lineio8_(&out[1], &inm8[1], &inm1[1], &inm2[1], &inm3[1], &inm4[1], &inm5[
1], &inm6[1], &inm7[1], &inm8[1], &i__1, stridel, stride, fptl,
lptl, area, &c_b16, &t1, &t2, &t3, &t4, &t5, &t6, &t7, &t8);
for (i__ = *n - 9; i__ >= 16; --i__) {
lineio8_(&out[1], &in[1], &inm1[1], &inm2[1], &inm3[1], &inm4[1], &
inm5[1], &inm6[1], &inm7[1], &inm8[1], &i__, stridel, stride,
fptl, lptl, area, &t0, &t1, &t2, &t3, &t4, &t5, &t6, &t7, &t8)
;
}
/* For the UNUSED lines below, we must point to "in" to avoid seg */
/* faults. On this end of the convolution, we know that "in" will */
/* always be in range. */
lineio8_(&out[1], &in[1], &inm1[1], &inm2[1], &inm3[1], &inm4[1], &inm5[1]
, &inm6[1], &inm7[1], &in[1], &c__15, stridel, stride, fptl, lptl,
area, &t0, &t1, &t2, &t3, &t4, &t5, &t6, &t7, &c_b16);
lineio8_(&out[1], &in[1], &inm1[1], &inm2[1], &inm3[1], &inm4[1], &inm5[1]
, &inm6[1], &in[1], &in[1], &c__14, stridel, stride, fptl, lptl,
area, &t0, &t1, &t2, &t3, &t4, &t5, &t6, &c_b16, &c_b16);
lineio8_(&out[1], &in[1], &inm1[1], &inm2[1], &inm3[1], &inm4[1], &inm5[1]
, &in[1], &in[1], &in[1], &c__13, stridel, stride, fptl, lptl,
area, &t0, &t1, &t2, &t3, &t4, &t5, &c_b16, &c_b16, &c_b16);
lineio8_(&out[1], &in[1], &inm1[1], &inm2[1], &inm3[1], &inm4[1], &in[1],
&in[1], &in[1], &in[1], &c__12, stridel, stride, fptl, lptl, area,
&t0, &t1, &t2, &t3, &t4, &c_b16, &c_b16, &c_b16, &c_b16);
lineio8_(&out[1], &in[1], &inm1[1], &inm2[1], &inm3[1], &in[1], &in[1], &
in[1], &in[1], &in[1], &c__11, stridel, stride, fptl, lptl, area,
&t0, &t1, &t2, &t3, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16);
lineio8_(&out[1], &in[1], &inm1[1], &inm2[1], &in[1], &in[1], &in[1], &in[
1], &in[1], &in[1], &c__10, stridel, stride, fptl, lptl, area, &
t0, &t1, &t2, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16);
lineio8_(&out[1], &in[1], &inm1[1], &in[1], &in[1], &in[1], &in[1], &in[1]
, &in[1], &in[1], &c__9, stridel, stride, fptl, lptl, area, &t0, &
t1, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16);
lineio8_(&out[1], &in[1], &in[1], &in[1], &in[1], &in[1], &in[1], &in[1],
&in[1], &in[1], &c__8, stridel, stride, fptl, lptl, area, &t0, &
c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16, &c_b16);
for (i__ = 7; i__ >= 0; --i__) {
zeroline8_(&out[1], &in[1], &inm1[1], &inm2[1], &inm3[1], &inm4[1], &
inm5[1], &inm6[1], &inm7[1], &inm8[1], &i__, stridel, stride,
fptl, lptl, area, &t0, &t1, &t2, &t3, &t4, &t5, &t6, &t7, &t8)
;
}
/* Forward phase (padding now ensures room!) */
t0 = tf[1];
t1 = tf[2];
t2 = tf[3];
t3 = tf[4];
t4 = tf[5];
t5 = tf[6];
t6 = tf[7];
t7 = tf[8];
t8 = tf[9];
i__1 = *n - 1;
for (i__ = 0; i__ <= i__1; i__ += 2) {
out2 = i__ * *stride;
i__2 = *area + out2;
for (j = out2 + 1; j <= i__2; ++j) {
out[j] = t0 * out[j] + t1 * outp1[j] + t2 * outp2[j] + t3 * outp3[
j] + t4 * outp4[j] + t5 * outp5[j] + t6 * outp6[j] + t7 *
outp7[j] + t8 * outp8[j];
}
}
return 0;
} /* pconvsetsi8_ */
/* Computes convolution for a single data set ("line") during reverse */
/* phases (INPUT/OUTPUT) */
/* out: output set */
/* in,inm1,inm2,inm3: input set */
/* i: index of line to be done */
/* stridel: stride of input */
/* stride: stride of output */
/* fptl: first point of input actually used */
/* lptl: last point of input actually used */
/* area: full area to be generated */
/* t0,t1,t2,t3: filter coefficients */
/* Subroutine */ int lineio3_(double *out, double *in, double *
inm1, double *inm2, double *inm3, int *i__, int *
stridel, int *stride, int *fptl, int *lptl, int *area,
double *t0, double *t1, double *t2, double *t3)
{
/* System generated locals */
int i__1;
/* Local variables */
int dout2, j, out2;
/* Parameter adjustments */
--inm3;
--inm2;
--inm1;
--in;
--out;
/* Function Body */
out2 = *i__ * *stridel;
dout2 = *i__ * (*stride - *stridel);
/* Zero out data before fptl */
i__1 = *fptl + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j + dout2] = (double)0.;
}
/* Do the line for "active" points */
i__1 = *lptl + 1 + out2;
for (j = *fptl + 1 + out2; j <= i__1; ++j) {
out[j + dout2] = *t0 * in[j] + *t1 * inm1[j] + *t2 * inm2[j] + *t3 *
inm3[j];
}
/* Zero out data after lptl */
i__1 = *area + out2;
for (j = *lptl + 1 + out2 + 1; j <= i__1; ++j) {
out[j + dout2] = (double)0.;
}
return 0;
} /* lineio3_ */
/* Computes convolution for a single data set ("line") during reverse */
/* phases (INPUT/OUTPUT) */
/* out: output set */
/* in,inm1,inm2,inm3: input set */
/* i: index of line to be done */
/* stridel: stride of input */
/* stride: stride of output */
/* fptl: first point of input actually used */
/* lptl: last point of input actually used */
/* area: full area to be generated */
/* t0,t1,t2,t3: filter coefficients */
/* Subroutine */ int lineio3full_(double *out, double *inm3,
double *inm2, double *inm1, double *in, double *inp1,
double *inp2, double *inp3, int *i__, int *stridel,
int *stride, int *fptl, int *lptl, int *area,
double *tl3, double *tl2, double *tl1, double *t0,
double *tr1, double *tr2, double *tr3)
{
/* System generated locals */
int i__1;
/* Local variables */
int dout2, j, out2;
/* Parameter adjustments */
--inp3;
--inp2;
--inp1;
--in;
--inm1;
--inm2;
--inm3;
--out;
/* Function Body */
out2 = *i__ * *stridel;
dout2 = *i__ * (*stride - *stridel);
i__1 = *lptl + out2 + 1;
for (j = out2 + 3; j <= i__1; j += 2) {
/* do j=1+out2+1,lptl+out2+1 */
out[j + dout2] = *tr3 * inp3[j] + *tr1 * inp1[j] + *tl1 * inm1[j] + *
tl3 * inm3[j];
out[j + dout2 + 1] = *tr3 * inp3[j + 1] + *tr1 * inp1[j + 1] + *t0 *
in[j + 1] + *tl1 * inm1[j + 1] + *tl3 * inm3[j + 1];
}
/* Zero out data after lptl */
i__1 = *area + out2;
for (j = *lptl + 1 + out2 + 1; j <= i__1; ++j) {
out[j + dout2] = (double)0.;
}
return 0;
} /* lineio3full_ */
/* Computes convolution for a single data set ("line") during reverse */
/* phases (INPUT/OUTPUT) */
/* out: output set */
/* in,inm1,inm2,inm3: input set */
/* i: index of line to be done */
/* stridel: stride of input */
/* stride: stride of output */
/* fptl: first point of input actually used */
/* lptl: last point of input actually used */
/* area: full area to be generated */
/* t0,t1,t2,t3: filter coefficients */
/* Subroutine */ int lineio8full_(double *out, double *inm8,
double *inm7, double *inm6, double *inm5, double *
inm4, double *inm3, double *inm2, double *inm1,
double *in, double *inp1, double *inp2, double *inp3,
double *inp4, double *inp5, double *inp6, double *
inp7, double *inp8, int *i__, int *stridel, int *
stride, int *fptl, int *lptl, int *area, double *tl8,
double *tl7, double *tl6, double *tl5, double *tl4,
double *tl3, double *tl2, double *tl1, double *t0,
double *tr1, double *tr2, double *tr3, double *tr4,
double *tr5, double *tr6, double *tr7, double *tr8)
{
/* System generated locals */
int i__1;
/* Local variables */
int dout2, j, out2;
/* Parameter adjustments */
--inp8;
--inp7;
--inp6;
--inp5;
--inp4;
--inp3;
--inp2;
--inp1;
--in;
--inm1;
--inm2;
--inm3;
--inm4;
--inm5;
--inm6;
--inm7;
--inm8;
--out;
/* Function Body */
out2 = *i__ * *stridel;
dout2 = *i__ * (*stride - *stridel);
i__1 = *lptl + out2 + 1;
for (j = out2 + 3; j <= i__1; j += 2) {
out[j + dout2] = *tr8 * inp8[j] + *tr7 * inp7[j] + *tr6 * inp6[j] + *
tr5 * inp5[j] + *tr4 * inp4[j] + *tr3 * inp3[j] + *tr2 * inp2[
j] + *tr1 * inp1[j] + *tl1 * inm1[j] + *tl2 * inm2[j] + *tl3 *
inm3[j] + *tl3 * inm3[j] + *tl4 * inm4[j] + *tl5 * inm4[j] +
*tl6 * inm6[j] + *tl7 * inm7[j] + *tl8 * inm8[j];
out[j + dout2 + 1] = *tr8 * inp8[j + 1] + *tr7 * inp7[j + 1] + *tr6 *
inp6[j + 1] + *tr5 * inp5[j + 1] + *tr4 * inp4[j + 1] + *tr3 *
inp3[j + 1] + *tr2 * inp2[j + 1] + *tr1 * inp1[j + 1] + *tl1
* inm1[j + 1] + *tl2 * inm2[j + 1] + *tl3 * inm3[j + 1] + *
tl3 * inm3[j + 1] + *tl4 * inm4[j + 1] + *tl5 * inm4[j + 1] +
*tl6 * inm6[j + 1] + *tl7 * inm7[j + 1] + *tl8 * inm8[j + 1];
}
/* Zero out data after lptl */
i__1 = *area + out2;
for (j = *lptl + 1 + out2 + 1; j <= i__1; ++j) {
out[j + dout2] = (double)0.;
}
return 0;
} /* lineio8full_ */
/* Computes convolution for a single data set ("line") during reverse */
/* phases (INPUT/OUTPUT) */
/* out: output set */
/* in,inm1,inm2,inm3: input set */
/* i: index of line to be done */
/* stridel: stride of input */
/* stride: stride of output */
/* fptl: first point of input actually used */
/* lptl: last point of input actually used */
/* area: full area to be generated */
/* t0,t1,t2,t3: filter coefficients */
/* Subroutine */ int lineio3fullacc_(double *out, double *inm3,
double *inm2, double *inm1, double *in, double *inp1,
double *inp2, double *inp3, int *i__, int *stridel,
int *stride, int *fptl, int *lptl, int *area,
double *tl3, double *tl2, double *tl1, double *t0,
double *tr1, double *tr2, double *tr3)
{
/* System generated locals */
int i__1;
/* Local variables */
int dout2, j, out2;
/* Parameter adjustments */
--inp3;
--inp2;
--inp1;
--in;
--inm1;
--inm2;
--inm3;
--out;
/* Function Body */
out2 = *i__ * *stridel;
dout2 = *i__ * (*stride - *stridel);
i__1 = *lptl + out2 + 1;
for (j = out2 + 1 + *fptl; j <= i__1; ++j) {
/* do j=1+out2+1,lptl+out2+1 */
out[j + dout2] = *tr3 * inp3[j] + *tr1 * inp1[j] + *t0 * in[j] + *tl1
* inm1[j] + *tl3 * inm3[j];
}
/* Zero out data after lptl */
i__1 = *area + out2;
for (j = *lptl + 1 + out2 + 1; j <= i__1; ++j) {
out[j + dout2] = (double)0.;
}
return 0;
} /* lineio3fullacc_ */
/* Computes convolution for a single data set ("line") during reverse */
/* phases (INPUT/OUTPUT) */
/* out: output set */
/* in,inm1,inm2,inm3: input set */
/* i: index of line to be done */
/* stridel: stride of input */
/* stride: stride of output */
/* fptl: first point of input actually used */
/* lptl: last point of input actually used */
/* area: full area to be generated */
/* t0,t1,t2,t3: filter coefficients */
/* Subroutine */ int lineio8fullacc_(double *out, double *inm8,
double *inm7, double *inm6, double *inm5, double *
inm4, double *inm3, double *inm2, double *inm1,
double *in, double *inp1, double *inp2, double *inp3,
double *inp4, double *inp5, double *inp6, double *
inp7, double *inp8, int *i__, int *stridel, int *
stride, int *fptl, int *lptl, int *area, double *tl8,
double *tl7, double *tl6, double *tl5, double *tl4,
double *tl3, double *tl2, double *tl1, double *t0,
double *tr1, double *tr2, double *tr3, double *tr4,
double *tr5, double *tr6, double *tr7, double *tr8)
{
/* System generated locals */
int i__1;
/* Local variables */
int dout2, j, out2;
/* Parameter adjustments */
--inp8;
--inp7;
--inp6;
--inp5;
--inp4;
--inp3;
--inp2;
--inp1;
--in;
--inm1;
--inm2;
--inm3;
--inm4;
--inm5;
--inm6;
--inm7;
--inm8;
--out;
/* Function Body */
out2 = *i__ * *stridel;
dout2 = *i__ * (*stride - *stridel);
/* Zero out data before fptl */
i__1 = *fptl + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j + dout2] = (double)0.;
}
i__1 = *lptl + out2 + 1;
for (j = out2 + 1 + *fptl; j <= i__1; ++j) {
/* do j=1+out2+2,lptl+out2+1 */
/* do j=1+out2+1,lptl+out2+1 */
out[j + dout2] = *tr8 * inp8[j] + *tr7 * inp7[j] + *tr6 * inp6[j] + *
tr5 * inp5[j] + *tr4 * inp4[j] + *tr3 * inp3[j] + *tr2 * inp2[
j] + *tr1 * inp1[j] + *t0 * in[j] + *tl1 * inm1[j] + *tl2 *
inm2[j] + *tl3 * inm3[j] + *tl4 * inm4[j] + *tl5 * inm5[j] + *
tl6 * inm6[j] + *tl7 * inm7[j] + *tl8 * inm8[j];
}
/* Zero out data after lptl */
i__1 = *area + out2;
for (j = *lptl + 1 + out2 + 1; j <= i__1; ++j) {
out[j + dout2] = (double)0.;
}
return 0;
} /* lineio8fullacc_ */
/* Computes convolution for a single data set ("line") during reverse */
/* phases (INPUT/OUTPUT) */
/* out: output set */
/* in,inm1,inm2,inm3,inm4,inm5,inm6,inm7,inm8: input set */
/* i: index of line to be done */
/* stridel: stride of input */
/* stride: stride of output */
/* fptl: first point of input actually used */
/* lptl: last point of input actually used */
/* area: full area to be generated */
/* t0,t1,t2,t3,t4,t5,t6,t7,t8: filter coefficients */
/* Subroutine */ int lineio8_(double *out, double *in, double *
inm1, double *inm2, double *inm3, double *inm4,
double *inm5, double *inm6, double *inm7, double *
inm8, int *i__, int *stridel, int *stride, int *fptl,
int *lptl, int *area, double *t0, double *t1,
double *t2, double *t3, double *t4, double *t5,
double *t6, double *t7, double *t8)
{
/* System generated locals */
int i__1;
/* Local variables */
int dout2, j, out2;
/* Parameter adjustments */
--inm8;
--inm7;
--inm6;
--inm5;
--inm4;
--inm3;
--inm2;
--inm1;
--in;
--out;
/* Function Body */
out2 = *i__ * *stridel;
dout2 = *i__ * (*stride - *stridel);
/* Zero out data before fptl */
i__1 = *fptl + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j + dout2] = (double)0.;
}
/* Do the line for "active" points */
i__1 = *lptl + 1 + out2;
for (j = *fptl + 1 + out2; j <= i__1; ++j) {
out[j + dout2] = *t0 * in[j] + *t1 * inm1[j] + *t2 * inm2[j] + *t3 *
inm3[j] + *t4 * inm4[j] + *t5 * inm5[j] + *t6 * inm6[j] + *t7
* inm7[j] + *t8 * inm8[j];
}
/* Zero out data after lptl */
i__1 = *area + out2;
for (j = *lptl + 1 + out2 + 1; j <= i__1; ++j) {
out[j + dout2] = (double)0.;
}
return 0;
} /* lineio8_ */
/* Zeros an output line (OUTPUT) */
/* consistent calling sequence with lineio */
/* out: output set */
/* in,inm1,inm2,inm3: input set */
/* i: index of line to be done */
/* stridel: stride of input */
/* stride: stride of output */
/* fptl: first point of input actually used */
/* lptl: last point of input actually used */
/* area: full area to be generated */
/* t0,t1,t2,t3: filter coefficients */
/* Subroutine */ int zeroline3_(double *out, double *in, double *
inm1, double *inm2, double *inm3, int *i__, int *
stridel, int *stride, int *fptl, int *lptl, int *area,
double *t0, double *t1, double *t2, double *t3)
{
/* System generated locals */
int i__1;
/* Local variables */
int j, out2;
/* Parameter adjustments */
--inm3;
--inm2;
--inm1;
--in;
--out;
/* Function Body */
out2 = *i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = (double)0.;
}
return 0;
} /* zeroline3_ */
/* Zeros an output line (OUTPUT) */
/* consistent calling sequence with lineio */
/* out: output set */
/* i: index of line to be done */
/* stride: stride of output */
/* area: full area to be generated */
/* Subroutine */ int zeroline8_(double *out, double *in, double *
inm1, double *inm2, double *inm3, double *inm4,
double *inm5, double *inm6, double *inm7, double *
inm8, int *i__, int *stridel, int *stride, int *fptl,
int *lptl, int *area, double *t0, double *t1,
double *t2, double *t3, double *t4, double *t5,
double *t6, double *t7, double *t8)
{
/* System generated locals */
int i__1;
/* Local variables */
int j, out2;
/* Parameter adjustments */
--inm8;
--inm7;
--inm6;
--inm5;
--inm4;
--inm3;
--inm2;
--inm1;
--in;
--out;
/* Function Body */
out2 = *i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = (double)0.;
}
return 0;
} /* zeroline8_ */
/* Computes convolution for a single "even" data set ("line") during */
/* reverse phases (INPUT/OUTPUT) */
/* out: output set */
/* in,inm1,inm2,inm3: input set */
/* i: index of line to be done */
/* stridel: stride of input */
/* stride: stride of output */
/* fptl: first point of input actually used */
/* lptl: last point of input actually used */
/* area: full area to be generated */
/* t0,t1,t2,t3: filter coefficients */
/* Subroutine */ int lineioe3_(double *out, double *in, double *
inm1, double *inm2, double *inm3, int *i__, int *
stridel, int *stride, int *fptl, int *lptl, int *area,
double *t0, double *t1, double *t2, double *t3)
{
/* System generated locals */
int i__1;
/* Local variables */
int dout2, j, out2;
/* Parameter adjustments */
--inm3;
--inm2;
--inm1;
--in;
--out;
/* Function Body */
out2 = *i__ * *stridel;
dout2 = *i__ * (*stride - *stridel);
/* Zero out data before fptl */
i__1 = *fptl + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j + dout2] = (double)0.;
}
/* Do the line for "active" points: */
/* Ignore any even input data on even input lines! */
i__1 = *lptl + 1 + out2;
for (j = *fptl + 1 + out2; j <= i__1; j += 2) {
out[j + dout2] = *t1 * inm1[j] + *t3 * inm3[j];
out[j + 1 + dout2] = *t0 * in[j + 1] + *t1 * inm1[j + 1] + *t2 * inm2[
j + 1] + *t3 * inm3[j + 1];
}
/* Zero out data after lptl */
i__1 = *area + out2;
for (j = *lptl + 1 + out2 + 1; j <= i__1; ++j) {
out[j + dout2] = (double)0.;
}
return 0;
} /* lineioe3_ */
/* Computes convolution for a single "even" data set ("line") during */
/* reverse phases (INPUT/OUTPUT) */
/* out: output set */
/* in,inm1,inm2,inm3,inm4,inm5,inm6,inm7,inm8: input set */
/* i: index of line to be done */
/* stridel: stride of input */
/* stride: stride of output */
/* fptl: first point of input actually used */
/* lptl: last point of input actually used */
/* area: full area to be generated */
/* t0,t1,t2,t3,t4,t5,t6,t7,t8: filter coefficients */
/* Subroutine */ int lineioe8_(double *out, double *in, double *
inm1, double *inm2, double *inm3, double *inm4,
double *inm5, double *inm6, double *inm7, double *
inm8, int *i__, int *stridel, int *stride, int *fptl,
int *lptl, int *area, double *t0, double *t1,
double *t2, double *t3, double *t4, double *t5,
double *t6, double *t7, double *t8)
{
/* System generated locals */
int i__1;
/* Local variables */
int dout2, j, out2;
/* Parameter adjustments */
--inm8;
--inm7;
--inm6;
--inm5;
--inm4;
--inm3;
--inm2;
--inm1;
--in;
--out;
/* Function Body */
out2 = *i__ * *stridel;
dout2 = *i__ * (*stride - *stridel);
/* Zero out data before fptl */
i__1 = *fptl + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j + dout2] = (double)0.;
}
/* Do the line for "active" points: */
/* Ignore any even input data on even input lines! */
i__1 = *lptl + 1 + out2;
for (j = *fptl + 1 + out2; j <= i__1; j += 2) {
out[j + dout2] = *t1 * inm1[j] + *t3 * inm3[j];
out[j + 1 + dout2] = *t0 * in[j + 1] + *t1 * inm1[j + 1] + *t2 * inm2[
j + 1] + *t3 * inm3[j + 1] + *t4 * inm4[j + 1] + *t5 * inm5[j
+ 1] + *t6 * inm6[j + 1] + *t7 * inm7[j + 1] + *t8 * inm8[j +
1];
}
/* Zero out data after lptl */
i__1 = *area + out2;
for (j = *lptl + 1 + out2 + 1; j <= i__1; ++j) {
out[j + dout2] = (double)0.;
}
return 0;
} /* lineioe8_ */
/* Computes convolution for a single "odd" data set ("line") during */
/* reverse phases (INPUT/OUTPUT) */
/* out: output set */
/* in,inm1,inm2,inm3: input set */
/* i: index of line to be done */
/* stridel: stride of input */
/* stride: stride of output */
/* fptl: first point of input actually used */
/* lptl: last point of input actually used */
/* area: full area to be generated */
/* t0,t1,t2,t3: filter coefficients */
/* Subroutine */ int lineioo3_(double *out, double *in, double *
inm1, double *inm2, double *inm3, int *i__, int *
stridel, int *stride, int *fptl, int *lptl, int *area,
double *t0, double *t1, double *t2, double *t3)
{
/* System generated locals */
int i__1;
/* Local variables */
int dout2, j, out2;
/* Parameter adjustments */
--inm3;
--inm2;
--inm1;
--in;
--out;
/* Function Body */
out2 = *i__ * *stridel;
dout2 = *i__ * (*stride - *stridel);
/* Zero out data before fptl */
i__1 = *fptl + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j + dout2] = (double)0.;
}
/* Do the line for "active" points: */
/* Ignore any even input data on even input lines! */
i__1 = *lptl + 1 + out2;
for (j = *fptl + 1 + out2; j <= i__1; j += 2) {
out[j + dout2] = *t0 * in[j] + *t2 * inm2[j];
out[j + 1 + dout2] = *t0 * in[j + 1] + *t1 * inm1[j + 1] + *t2 * inm2[
j + 1] + *t3 * inm3[j + 1];
}
/* Zero out data after lptl */
i__1 = *area + out2;
for (j = *lptl + 1 + out2 + 1; j <= i__1; ++j) {
out[j + dout2] = (double)0.;
}
return 0;
} /* lineioo3_ */
/* Computes convolution for a single "odd" data set ("line") during */
/* reverse phases (INPUT/OUTPUT) */
/* out: output set */
/* in,inm1,inm2,inm3,inm4,inm5,inm6,inm7,inm8: input set */
/* i: index of line to be done */
/* stridel: stride of input */
/* stride: stride of output */
/* fptl: first point of input actually used */
/* lptl: last point of input actually used */
/* area: full area to be generated */
/* t0,t1,t2,t3,t4,t5,t6,t7,t8: filter coefficients */
/* Subroutine */ int lineioo8_(double *out, double *in, double *
inm1, double *inm2, double *inm3, double *inm4,
double *inm5, double *inm6, double *inm7, double *
inm8, int *i__, int *stridel, int *stride, int *fptl,
int *lptl, int *area, double *t0, double *t1,
double *t2, double *t3, double *t4, double *t5,
double *t6, double *t7, double *t8)
{
/* System generated locals */
int i__1;
/* Local variables */
int dout2, j, out2;
/* Parameter adjustments */
--inm8;
--inm7;
--inm6;
--inm5;
--inm4;
--inm3;
--inm2;
--inm1;
--in;
--out;
/* Function Body */
out2 = *i__ * *stridel;
dout2 = *i__ * (*stride - *stridel);
/* Zero out data before fptl */
i__1 = *fptl + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j + dout2] = (double)0.;
}
/* Do the line for "active" points: */
/* Ignore any even input data on even input lines! */
i__1 = *lptl + 1 + out2;
for (j = *fptl + 1 + out2; j <= i__1; j += 2) {
out[j + dout2] = *t0 * in[j] + *t2 * inm2[j];
out[j + 1 + dout2] = *t0 * in[j + 1] + *t1 * inm1[j + 1] + *t2 * inm2[
j + 1] + *t3 * inm3[j + 1] + *t4 * inm4[j + 1] + *t5 * inm5[j
+ 1] + *t6 * inm6[j + 1] + *t8 * inm8[j + 1];
}
/* Zero out data after lptl */
i__1 = *area + out2;
for (j = *lptl + 1 + out2 + 1; j <= i__1; ++j) {
out[j + dout2] = (double)0.;
}
return 0;
} /* lineioo8_ */
/* Back and forth, linear convolution of sets of points for +/- 3 */
/* onvolutions using padding, and ignoring input data on non-ghost */
/* points (IN PLACE) */
/* out: input/output array for in place computation */
/* outm3,outm2,outm1,outp1,outp2,outp3: pointers in output array to */
/* start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, */
/* 2 forward, 3 forward */
/* area: total size of each set */
/* stride: separation between data sets (stride=&out-&outm1) */
/* n: actual number of data sets being convolved */
/* npad: size of padded region needed for "leakage" of foward filters */
/* tr: reverse-time filter coefficients */
/* tf: forward-time filter coefficients */
/* Subroutine */ int convsetsp3_(double *out, double *in, double *
inm3, double *inm2, double *inm1, double *inp1,
double *inp2, double *inp3, int *area, int *stride,
int *n, int *npad, double *tr, double *tf)
{
/* System generated locals */
int i__1, i__2;
/* Local variables */
int i__, j;
double t0, t1, t3;
int out2;
/* Forward phase (padding now ensures room!) */
/* No change here relative to convsets */
/* Parameter adjustments */
--tf;
--tr;
--inp3;
--inp2;
--inp1;
--inm1;
--inm2;
--inm3;
--in;
--out;
/* Function Body */
t0 = tf[1];
t1 = tf[2];
t3 = tf[4];
i__1 = *n - 5;
for (i__ = 3; i__ <= i__1; i__ += 2) {
out2 = i__ * *stride;
i__2 = *area + out2;
for (j = out2 + 1; j <= i__2; ++j) {
out[j] = t3 * inm3[j] + t1 * inm1[j] + t1 * inp1[j] + t3 * inp3[j]
;
}
}
i__1 = *n - 5;
for (i__ = 2; i__ <= i__1; i__ += 2) {
out2 = i__ * *stride;
i__2 = *area + out2;
for (j = out2 + 1; j <= i__2; ++j) {
out[j] = in[j];
}
}
return 0;
} /* convsetsp3_ */
/* Back and forth, linear convolution of sets of points for +/- 8 */
/* onvolutions using padding, and ignoring input data on non-ghost */
/* points (IN PLACE) */
/* out: input/output array for in place computation */
/* outm3,outm2,outm1,outp1,outp2,outp3: pointers in output array to */
/* start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, */
/* 2 forward, 3 forward */
/* area: total size of each set */
/* stride: separation between data sets (stride=&out-&outm1) */
/* n: actual number of data sets being convolved */
/* npad: size of padded region needed for "leakage" of foward filters */
/* tr: reverse-time filter coefficients */
/* tf: forward-time filter coefficients */
/* Subroutine */ int convsetsp8_(double *out, double *outm8,
double *outm7, double *outm6, double *outm5, double *
outm4, double *outm3, double *outm2, double *outm1,
double *outp1, double *outp2, double *outp3, double *
outp4, double *outp5, double *outp6, double *outp7,
double *outp8, int *area, int *stride, int *n,
int *npad, double *tr, double *tf)
{
/* System generated locals */
int i__1, i__2;
/* Local variables */
int i__, j;
double t0, t1, t2, t3, t4, t5, t6, t7, t8;
int out2;
/* Reverse direction */
/* Relative to convsets, we wiped out the non-ghost references */
/* Parameter adjustments */
--tf;
--tr;
--outp8;
--outp7;
--outp6;
--outp5;
--outp4;
--outp3;
--outp2;
--outp1;
--outm1;
--outm2;
--outm3;
--outm4;
--outm5;
--outm6;
--outm7;
--outm8;
--out;
/* Function Body */
t0 = tr[1];
t1 = tr[2];
t2 = tr[3];
t3 = tr[4];
t4 = tr[5];
t5 = tr[6];
t6 = tr[7];
t7 = tr[8];
t8 = tr[9];
/* Pad region */
/* out2=(npad-1)*stride */
/* do j=1+out2,area+out2 */
/* out(j)=0. */
/* enddo */
/* out2=(npad-2)*stride */
/* do j=1+out2,area+out2 */
/* out(j)=t8*outm8(j) */
/* enddo */
/* out2=(npad-3)*stride */
/* do j=1+out2,area+out2 */
/* out(j)=t7*outm7(j) */
/* enddo */
/* out2=(npad-4)*stride */
/* do j=1+out2,area+out2 */
/* out(j)=t6*outm6(j)+t8*outm8(j) */
/* enddo */
/* out2=(npad-5)*stride */
/* do j=1+out2,area+out2 */
/* out(j)=t5*outm5(j)+t7*outm7(j) */
/* enddo */
/* out2=(npad-6)*stride */
/* do j=1+out2,area+out2 */
/* out(j)=t4*outm4(j)+t6*outm6(j)+t8*outm8(j) */
/* enddo */
/* out2=(npad-7)*stride */
/* do j=1+out2,area+out2 */
/* out(j)=t3*outm3(j)+t5*outm5(j)+t7*outm7(j) */
/* enddo */
/* out2=(npad-8)*stride */
/* do j=1+out2,area+out2 */
/* out(j)=t2*outm2(j)+t4*outm4(j)+t6*outm6(j)+t8*outm8(j) */
/* enddo */
/* Main reverse phase exploiting zeros! */
for (i__ = *n - 1; i__ >= 8; i__ += -2) {
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t1 * outm1[j] + t3 * outm3[j] + t5 * outm5[j] + t7 *
outm7[j];
}
out2 = (i__ - 1) * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j] + t2 * outm2[j] + t4 * outm4[j] + t6 * outm6[
j] + t8 * outm8[j];
}
}
/* Fence posting */
/* out2=7*stride */
/* do j=1+out2,area+out2 */
/* out(j)=t1*outm1(j)+t3*outm3(j)+t5*outm5(j)+t7*outm7(j) */
/* enddo */
/* out2=6*stride */
/* do j=1+out2,area+out2 */
/* out(j)=t0*out(j)+t2*outm2(j)+t4*outm4(j)+t6*outm6(j) */
/* enddo */
/* out2=5*stride */
/* do j=1+out2,area+out2 */
/* out(j)=t1*outm1(j)+t3*outm3(j)+t5*outm5(j) */
/* enddo */
/* out2=4*stride */
/* do j=1+out2,area+out2 */
/* out(j)=t0*out(j)+t2*outm2(j)+t4*outm4(j) */
/* enddo */
/* out2=3*stride */
/* do j=1+out2,area+out2 */
/* out(j)=t1*outm1(j)+t3*outm3(j) */
/* enddo */
/* out2=2*stride */
/* do j=1+out2,area+out2 */
/* out(j)=t0*out(j)+t2*outm2(j) */
/* enddo */
/* out2=1*stride */
/* do j=1+out2,area+out2 */
/* out(j)=t1*outm1(j) */
/* enddo */
/* out2=0*stride */
/* do j=1+out2,area+out2 */
/* out(j)=t0*out(j) */
/* enddo */
/* Forward phase (padding now ensures room!) */
/* No change here relative to convsets */
t0 = tf[1];
t1 = tf[2];
t2 = tf[3];
t3 = tf[4];
t4 = tf[5];
t5 = tf[6];
t6 = tf[7];
t7 = tf[8];
t8 = tf[9];
/* Need to verify that this is OK..... */
/* do i=0,n-1 */
i__1 = *n - 8;
for (i__ = 8; i__ <= i__1; ++i__) {
out2 = i__ * *stride;
i__2 = *area + out2;
for (j = out2 + 1; j <= i__2; ++j) {
out[j] = t0 * out[j] + t1 * outp1[j] + t2 * outp2[j] + t3 * outp3[
j] + t4 * outp4[j] + t5 * outp5[j] + t6 * outp6[j] + t7 *
outp7[j] + t8 * outp8[j];
}
}
return 0;
} /* convsetsp8_ */
/* Back and forth, linear convolution of sets of points for +/- 3 */
/* onvolutions using padding, and ignoring input data on non-ghost */
/* points (IN PLACE) */
/* out: input/output array for in place computation */
/* outm3,outm2,outm1,outp1,outp2,outp3: pointers in output array to */
/* start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, */
/* 2 forward, 3 forward */
/* area: total size of each set */
/* stride: separation between data sets (stride=&out-&outm1) */
/* n: actual number of data sets being convolved */
/* npad: size of padded region needed for "leakage" of foward filters */
/* tr: reverse-time filter coefficients */
/* tf: forward-time filter coefficients */
/* Subroutine */ int convsetsp8full_(double *out, double *in,
double *inm8, double *inm7, double *inm6, double *
inm5, double *inm4, double *inm3, double *inm2,
double *inm1, double *inp1, double *inp2, double *
inp3, double *inp4, double *inp5, double *inp6,
double *inp7, double *inp8, int *area, int *stride,
int *n, int *npad, double *filter)
{
/* System generated locals */
int i__1, i__2;
/* Local variables */
int i__, j;
double t0, t1, t2, t3, t4, t5, t6, t7, t8;
int out2;
/* Forward phase (padding now ensures room!) */
/* No change here relative to convsets */
/* Parameter adjustments */
--filter;
--inp8;
--inp7;
--inp6;
--inp5;
--inp4;
--inp3;
--inp2;
--inp1;
--inm1;
--inm2;
--inm3;
--inm4;
--inm5;
--inm6;
--inm7;
--inm8;
--in;
--out;
/* Function Body */
t0 = filter[1];
t1 = filter[2];
t2 = filter[3];
t3 = filter[4];
t4 = filter[5];
t5 = filter[6];
t6 = filter[7];
t7 = filter[8];
t8 = filter[9];
if (t0 != (double)0.) {
/* symmetric filter */
/* thin func */
i__1 = *n - 8;
for (i__ = 9; i__ <= i__1; i__ += 2) {
out2 = i__ * *stride;
i__2 = *area + out2;
for (j = out2 + 1; j <= i__2; ++j) {
out[j] = t7 * inm7[j] + t5 * inm5[j] + t3 * inm3[j] + t1 *
inm1[j] + t1 * inp1[j] + t3 * inp3[j] + t5 * inp5[j]
+ t7 * inp7[j];
}
}
/* fat func */
i__1 = *n - 8;
for (i__ = 8; i__ <= i__1; i__ += 2) {
out2 = i__ * *stride;
i__2 = *area + out2;
for (j = out2 + 1; j <= i__2; ++j) {
out[j] = t8 * inm8[j] + t6 * inm6[j] + t4 * inm4[j] + t2 *
inm2[j] + t0 * in[j] + t2 * inp2[j] + t4 * inp4[j] +
t6 * inp6[j] + t8 * inp8[j];
}
}
/* end of symmetric part */
} else {
/* antisymmetric filter */
/* thin func */
i__1 = *n - 8;
for (i__ = 9; i__ <= i__1; i__ += 2) {
out2 = i__ * *stride;
i__2 = *area + out2;
for (j = out2 + 1; j <= i__2; ++j) {
out[j] = t7 * (double)-1. * inm7[j] - t5 * inm5[j] - t3 * inm3[
j] - t1 * inm1[j] + t1 * inp1[j] + t3 * inp3[j] + t5 *
inp5[j] + t7 * inp7[j];
}
}
/* fat func */
i__1 = *n - 8;
for (i__ = 8; i__ <= i__1; i__ += 2) {
out2 = i__ * *stride;
i__2 = *area + out2;
for (j = out2 + 1; j <= i__2; ++j) {
out[j] = t8 * (double)-1. * inm8[j] - t6 * inm6[j] - t4 * inm4[
j] - t2 * inm2[j] + t0 * in[j] + t2 * inp2[j] + t4 *
inp4[j] + t6 * inp6[j] + t8 * inp8[j];
}
}
/* end of antisymmetric part */
}
return 0;
} /* convsetsp8full_ */
/* Back and forth, linear convolution of sets of points using padding, */
/* ignoring input data on non-ghost points (IN/OUT) */
/* lower: output array */
/* outm3,outm2,outm1,outp1,outp2,outp3: pointers to input work array to */
/* start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, */
/* 2 forward, 3 forward */
/* area: total size of each set */
/* stride: separation between data sets (stride=&out-&outm1) */
/* n: actual number of data sets being convolved */
/* npad: size of padded region needed for "leakage" of foward filters */
/* stridel: stride on output */
/* fptl,lptl: first,last actual points in output set */
/* fsetl,lsetl: first,last actual sets in output */
/* tr: reverse-time filter coefficients */
/* tf: forward-time filter coefficients */
/* Subroutine */ int convsetspo3_(double *lower, double *out,
double *outm3, double *outm2, double *outm1, double *
outp1, double *outp2, double *outp3, int *area, int *
stride, int *n, int *npad, int *stridel, int *fptl,
int *lptl, int *fsetl, int *lsetl, double *tr,
double *tf)
{
/* System generated locals */
int i__1, i__2;
/* Local variables */
int i__, j, out2s;
double t0, t1, t2, t3;
int jl, out2;
/* Reverse direction */
/* No change here relative to convsetsp */
/* Parameter adjustments */
--tf;
--tr;
--outp3;
--outp2;
--outp1;
--outm1;
--outm2;
--outm3;
--out;
--lower;
/* Function Body */
t0 = tr[1];
t1 = tr[2];
t2 = tr[3];
t3 = tr[4];
/* Pad region */
out2 = (*npad - 1) * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = (double)0.;
}
out2 = (*npad - 2) * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t3 * outm3[j];
}
out2 = (*npad - 3) * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t2 * outm2[j];
}
/* Main reverse phase exploiting zeros! */
for (i__ = *n - 1; i__ >= 3; i__ += -2) {
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t1 * outm1[j] + t3 * outm3[j];
}
out2 = (i__ - 1) * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j] + t2 * outm2[j];
}
}
/* Fence posting */
out2 = *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t1 * outm1[j];
}
out2 = 0;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j];
}
/* Forward phase (padding now ensures room!) */
/* Relative to convsetsp, output now places into lower instead */
t0 = tf[1];
t1 = tf[2];
t2 = tf[3];
t3 = tf[4];
i__1 = *lsetl;
for (i__ = *fsetl; i__ <= i__1; ++i__) {
out2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__2 = out2 + 1 + *lptl;
for (j = out2 + 1 + *fptl; j <= i__2; ++j) {
jl = j + out2s;
/* $$$ out(j)= */
/* $$$ $ t0*out(j)+t1*outp1(j)+t2*outp2(j)+t3*outp3(j) */
lower[jl] = t0 * out[j] + t1 * outp1[j] + t2 * outp2[j] + t3 *
outp3[j];
}
}
return 0;
} /* convsetspo3_ */
/* lower: output array */
/* outm3,outm2,outm1,outp1,outp2,outp3: pointers to input work array to */
/* start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, */
/* 2 forward, 3 forward */
/* area: total size of each set */
/* stride: separation between data sets (stride=&out-&outm1) */
/* n: actual number of data sets being convolved */
/* npad: size of padded region needed for "leakage" of foward filters */
/* stridel: stride on output */
/* fptl,lptl: first,last actual points in output set */
/* fsetl,lsetl: first,last actual sets in output */
/* tr: reverse-time filter coefficients */
/* tf: forward-time filter coefficients */
/* Subroutine */ int convsetspoadd3_(double *lower, double *addon,
double *out, double *outm3, double *outm2, double *
outm1, double *outp1, double *outp2, double *outp3,
int *area, int *stride, int *n, int *npad, int *
stridel, int *fptl, int *lptl, int *fsetl, int *lsetl,
double *tr, double *tf)
{
/* System generated locals */
int i__1, i__2;
/* Local variables */
int i__, j, out2s;
double t0, t1, t2, t3;
int jl, out2;
/* Reverse direction */
/* No change here relative to convsetsp */
/* Parameter adjustments */
--tf;
--tr;
--outp3;
--outp2;
--outp1;
--outm1;
--outm2;
--outm3;
--out;
--addon;
--lower;
/* Function Body */
t0 = tr[1];
t1 = tr[2];
t2 = tr[3];
t3 = tr[4];
/* Pad region */
out2 = (*npad - 1) * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = (double)0.;
}
out2 = (*npad - 2) * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t3 * outm3[j];
}
out2 = (*npad - 3) * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t2 * outm2[j];
}
/* Main reverse phase exploiting zeros! */
for (i__ = *n - 1; i__ >= 3; i__ += -2) {
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t1 * outm1[j] + t3 * outm3[j];
}
out2 = (i__ - 1) * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j] + t2 * outm2[j];
}
}
/* Fence posting */
out2 = *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t1 * outm1[j];
}
out2 = 0;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j];
}
/* Forward phase (padding now ensures room!) */
/* Relative to convsetsp, output now places into lower instead */
t0 = tf[1];
t1 = tf[2];
t2 = tf[3];
t3 = tf[4];
i__1 = *lsetl;
for (i__ = *fsetl; i__ <= i__1; ++i__) {
out2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__2 = out2 + 1 + *lptl;
for (j = out2 + 1 + *fptl; j <= i__2; ++j) {
jl = j + out2s;
/* $$$ out(j)= */
/* $$$ $ t0*out(j)+t1*outp1(j)+t2*outp2(j)+t3*outp3(j) */
lower[jl] = addon[jl] + t0 * out[j] + t1 * outp1[j] + t2 * outp2[
j] + t3 * outp3[j];
}
}
return 0;
} /* convsetspoadd3_ */
/* Back and forth, linear convolution of sets of points for +/- 3 */
/* filters using padding, ignoring input data on non-ghost points */
/* (IN/+=OUT) */
/* lower: output array for accumulation */
/* outm3,outm2,outm1,outp1,outp2,outp3: pointers to input work array to */
/* start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, */
/* 2 forward, 3 forward */
/* area: total size of each set */
/* stride: separation between data sets (stride=&out-&outm1) */
/* n: actual number of data sets being convolved */
/* npad: size of padded region needed for "leakage" of foward filters */
/* stridel: stride on output */
/* fptl,lptl: first,last actual points in output set */
/* fsetl,lsetl: first,last actual sets in output */
/* tr: reverse-time filter coefficients */
/* tf: forward-time filter coefficients */
/* Subroutine */ int convsetspop3_(double *lower, double *out,
double *outm3, double *outm2, double *outm1, double *
outp1, double *outp2, double *outp3, int *area, int *
stride, int *n, int *npad, int *stridel, int *fptl,
int *lptl, int *fsetl, int *lsetl, double *filt)
{
/* System generated locals */
int i__1, i__2;
/* Local variables */
int i__, j, out2s;
double t0, t1, t3;
int jl, out2;
/* Forward phase (padding now ensures room!) */
/* Relative to convsetspo, only change is accumulation into lower */
/* Parameter adjustments */
--filt;
--outp3;
--outp2;
--outp1;
--outm1;
--outm2;
--outm3;
--out;
--lower;
/* Function Body */
t0 = filt[1];
t1 = filt[2];
t3 = filt[4];
i__1 = *lsetl;
for (i__ = *fsetl + 1; i__ <= i__1; i__ += 2) {
out2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__2 = out2 + 1 + *lptl;
for (j = out2 + 1 + *fptl; j <= i__2; ++j) {
jl = j + out2s;
lower[jl] = lower[jl] + t3 * outm3[j] + t1 * outm1[j] + t1 *
outp1[j] + t3 * outp3[j];
}
}
i__1 = *lsetl;
for (i__ = *fsetl; i__ <= i__1; i__ += 2) {
out2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__2 = out2 + 1 + *lptl;
for (j = out2 + 1 + *fptl; j <= i__2; ++j) {
jl = j + out2s;
lower[jl] += out[j];
}
}
return 0;
} /* convsetspop3_ */
/* Back and forth, linear convolution of sets of points for +/- 3 */
/* filters using padding, ignoring input data on non-ghost points */
/* (IN/+=OUT) */
/* lower: output array for accumulation */
/* outm3,outm2,outm1,outp1,outp2,outp3: pointers to input work array to */
/* start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, */
/* 2 forward, 3 forward */
/* area: total size of each set */
/* stride: separation between data sets (stride=&out-&outm1) */
/* n: actual number of data sets being convolved */
/* npad: size of padded region needed for "leakage" of foward filters */
/* stridel: stride on output */
/* fptl,lptl: first,last actual points in output set */
/* fsetl,lsetl: first,last actual sets in output */
/* tr: reverse-time filter coefficients */
/* tf: forward-time filter coefficients */
/* Subroutine */ int convsetspom3_(double *lower, double *out,
double *outm3, double *outm2, double *outm1, double *
outp1, double *outp2, double *outp3, int *area, int *
stride, int *n, int *npad, int *stridel, int *fptl,
int *lptl, int *fsetl, int *lsetl, double *filt)
{
/* System generated locals */
int i__1, i__2;
/* Local variables */
int i__, j, out2s;
double t0, t1, t3;
int jl, out2;
/* Forward phase (padding now ensures room!) */
/* Relative to convsetspo, only change is accumulation into lower */
/* Parameter adjustments */
--filt;
--outp3;
--outp2;
--outp1;
--outm1;
--outm2;
--outm3;
--out;
--lower;
/* Function Body */
t0 = filt[1];
t1 = filt[2];
t3 = filt[4];
i__1 = *lsetl;
for (i__ = *fsetl + 1; i__ <= i__1; i__ += 2) {
out2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__2 = out2 + 1 + *lptl;
for (j = out2 + 1 + *fptl; j <= i__2; ++j) {
jl = j + out2s;
lower[jl] = lower[jl] - t3 * outm3[j] - t1 * outm1[j] - t1 *
outp1[j] - t3 * outp3[j];
}
}
i__1 = *lsetl;
for (i__ = *fsetl; i__ <= i__1; i__ += 2) {
out2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__2 = out2 + 1 + *lptl;
for (j = out2 + 1 + *fptl; j <= i__2; ++j) {
jl = j + out2s;
lower[jl] -= out[j];
}
}
return 0;
} /* convsetspom3_ */
/* Back and forth, linear convolution of sets of points for +/- 8 */
/* filters using padding, ignoring input data on non-ghost points */
/* (IN/+=OUT) */
/* lower: output array for accumulation */
/* outm3,outm2,outm1,outp1,outp2,outp3: pointers to input work array to */
/* start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, */
/* 2 forward, 3 forward */
/* area: total size of each set */
/* stride: separation between data sets (stride=&out-&outm1) */
/* n: actual number of data sets being convolved */
/* npad: size of padded region needed for "leakage" of foward filters */
/* stridel: stride on output */
/* fptl,lptl: first,last actual points in output set */
/* fsetl,lsetl: first,last actual sets in output */
/* tr: reverse-time filter coefficients */
/* tf: forward-time filter coefficients */
/* Subroutine */ int convsetspop8_(double *lower, double *out,
double *outm8, double *outm7, double *outm6, double *
outm5, double *outm4, double *outm3, double *outm2,
double *outm1, double *outp1, double *outp2, double *
outp3, double *outp4, double *outp5, double *outp6,
double *outp7, double *outp8, int *area, int *stride,
int *n, int *npad, int *stridel, int *fptl, int *
lptl, int *fsetl, int *lsetl, double *tr, double *tf)
{
/* System generated locals */
int i__1, i__2;
/* Local variables */
int i__, j, out2s;
double t0, t1, t2, t3, t4, t5, t6, t7, t8;
int jl, out2;
/* Reverse direction */
/* Relative to convsets, we wiped out the non-ghost references */
/* Parameter adjustments */
--tf;
--tr;
--outp8;
--outp7;
--outp6;
--outp5;
--outp4;
--outp3;
--outp2;
--outp1;
--outm1;
--outm2;
--outm3;
--outm4;
--outm5;
--outm6;
--outm7;
--outm8;
--out;
--lower;
/* Function Body */
t0 = tr[1];
t1 = tr[2];
t2 = tr[3];
t3 = tr[4];
t4 = tr[5];
t5 = tr[6];
t6 = tr[7];
t7 = tr[8];
t8 = tr[9];
/* Pad region */
out2 = (*npad - 1) * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = (double)0.;
}
out2 = (*npad - 2) * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t8 * outm8[j];
}
out2 = (*npad - 3) * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t7 * outm7[j];
}
out2 = (*npad - 4) * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t6 * outm6[j] + t8 * outm8[j];
}
out2 = (*npad - 5) * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t5 * outm5[j] + t7 * outm7[j];
}
out2 = (*npad - 6) * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t4 * outm4[j] + t6 * outm6[j] + t8 * outm8[j];
}
out2 = (*npad - 7) * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t3 * outm3[j] + t5 * outm5[j] + t7 * outm7[j];
}
out2 = (*npad - 8) * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t2 * outm2[j] + t4 * outm4[j] + t6 * outm6[j] + t8 * outm8[j]
;
}
/* Main reverse phase exploiting zeros! */
for (i__ = *n - 1; i__ >= 8; i__ += -2) {
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t7 * outm7[j] + t5 * outm5[j] + t3 * outm3[j] + t1 *
outm1[j];
}
out2 = (i__ - 1) * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j] + t2 * outm2[j] + t4 * outm4[j] + t6 * outm6[
j] + t8 * outm8[j];
}
}
/* Fence posting */
out2 = *stride * 7;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t1 * outm1[j] + t3 * outm3[j] + t5 * outm5[j] + t7 * outm7[j]
;
}
out2 = *stride * 6;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j] + t2 * outm2[j] + t4 * outm4[j] + t6 * outm6[j];
}
out2 = *stride * 5;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t1 * outm1[j] + t3 * outm3[j] + t5 * outm5[j];
}
out2 = *stride << 2;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j] + t2 * outm2[j] + t4 * outm4[j];
}
out2 = *stride * 3;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t1 * outm1[j] + t3 * outm3[j];
}
out2 = *stride << 1;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j] + t2 * outm2[j];
}
out2 = *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t1 * outm1[j];
}
out2 = 0;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j];
}
/* Forward phase (padding now ensures room!) */
/* No change here relative to convsets */
t0 = tf[1];
t1 = tf[2];
t2 = tf[3];
t3 = tf[4];
t4 = tf[5];
t5 = tf[6];
t6 = tf[7];
t7 = tf[8];
t8 = tf[9];
i__1 = *lsetl;
for (i__ = *fsetl; i__ <= i__1; ++i__) {
out2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__2 = out2 + 1 + *lptl;
for (j = out2 + 1 + *fptl; j <= i__2; ++j) {
jl = j + out2s;
/* $$$ out(j)= */
/* $$$ $ t0*out(j)+t1*outp1(j)+t2*outp2(j)+t3*outp3(j) */
lower[jl] = lower[jl] + t0 * out[j] + t1 * outp1[j] + t2 * outp2[
j] + t3 * outp3[j] + t4 * outp4[j] + t5 * outp5[j] + t6 *
outp6[j] + t7 * outp7[j] + t8 * outp8[j];
}
}
return 0;
} /* convsetspop8_ */
/* Back and forth, linear convolution of sets of points for +/- 3 */
/* filters using padding, ignoring input data on non-ghost points */
/* (IN/+=OUT) */
/* lower: output array for accumulation */
/* outm3,outm2,outm1,outp1,outp2,outp3: pointers to input work array to */
/* start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, */
/* 2 forward, 3 forward */
/* area: total size of each set */
/* stride: separation between data sets (stride=&out-&outm1) */
/* n: actual number of data sets being convolved */
/* npad: size of padded region needed for "leakage" of foward filters */
/* stridel: stride on output */
/* fptl,lptl: first,last actual points in output set */
/* fsetl,lsetl: first,last actual sets in output */
/* tr: reverse-time filter coefficients */
/* tf: forward-time filter coefficients */
/* Subroutine */ int convsetspop8full_(double *lower, double *out,
double *outm8, double *outm7, double *outm6, double *
outm5, double *outm4, double *outm3, double *outm2,
double *outm1, double *outp1, double *outp2, double *
outp3, double *outp4, double *outp5, double *outp6,
double *outp7, double *outp8, int *area, int *stride,
int *n, int *npad, int *stridel, int *fptl, int *
lptl, int *fsetl, int *lsetl, double *filt)
{
/* System generated locals */
int i__1, i__2;
/* Local variables */
int i__, j, out2s;
double t0, t1, t2, t3, t4, t5, t6, t7, t8;
int jl, out2;
/* Forward phase (padding now ensures room!) */
/* Relative to convsetspo, only change is accumulation into lower */
/* Parameter adjustments */
--filt;
--outp8;
--outp7;
--outp6;
--outp5;
--outp4;
--outp3;
--outp2;
--outp1;
--outm1;
--outm2;
--outm3;
--outm4;
--outm5;
--outm6;
--outm7;
--outm8;
--out;
--lower;
/* Function Body */
t0 = filt[1];
t1 = filt[2];
t2 = filt[3];
t3 = filt[4];
t4 = filt[5];
t5 = filt[6];
t6 = filt[7];
t7 = filt[8];
t8 = filt[9];
if (t0 != (double)0.) {
/* symmetric filter */
i__1 = *lsetl;
for (i__ = *fsetl + 1; i__ <= i__1; i__ += 2) {
out2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__2 = out2 + 1 + *lptl;
for (j = out2 + 1 + *fptl; j <= i__2; ++j) {
jl = j + out2s;
lower[jl] = lower[jl] + t7 * outm7[j] + t5 * outm5[j] + t3 *
outm3[j] + t1 * outm1[j] + t1 * outp1[j] + t3 * outp3[
j] + t5 * outp5[j] + t7 * outp7[j];
}
}
i__1 = *lsetl;
for (i__ = *fsetl; i__ <= i__1; i__ += 2) {
out2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__2 = out2 + 1 + *lptl;
for (j = out2 + 1 + *fptl; j <= i__2; ++j) {
jl = j + out2s;
lower[jl] = lower[jl] + t8 * outm8[j] + t6 * outm6[j] + t4 *
outm4[j] + t2 * outm2[j] + t0 * out[j] + t2 * outp2[j]
+ t4 * outp4[j] + t6 * outp6[j] + t8 * outp8[j];
}
}
/* end of symmetric part */
} else {
/* antisymmetric filter */
i__1 = *lsetl;
for (i__ = *fsetl + 1; i__ <= i__1; i__ += 2) {
out2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__2 = out2 + 1 + *lptl;
for (j = out2 + 1 + *fptl; j <= i__2; ++j) {
jl = j + out2s;
lower[jl] = lower[jl] - t7 * outm7[j] - t5 * outm5[j] - t3 *
outm3[j] - t1 * outm1[j] + t1 * outp1[j] + t3 * outp3[
j] + t5 * outp5[j] + t7 * outp7[j];
}
}
i__1 = *lsetl;
for (i__ = *fsetl; i__ <= i__1; i__ += 2) {
out2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__2 = out2 + 1 + *lptl;
for (j = out2 + 1 + *fptl; j <= i__2; ++j) {
jl = j + out2s;
lower[jl] = lower[jl] - t8 * outm8[j] - t6 * outm6[j] - t4 *
outm4[j] - t2 * outm2[j] + t0 * out[j] + t2 * outp2[j]
+ t4 * outp4[j] + t6 * outp6[j] + t8 * outp8[j];
}
}
/* end of antisymmetric part */
}
return 0;
} /* convsetspop8full_ */
/* Back and forth, linear convolution of sets of points for +/- 3 */
/* filters using padding, ignoring input data on non-ghost points */
/* (IN/+=OUT) */
/* lower: output array for accumulation */
/* outm3,outm2,outm1,outp1,outp2,outp3: pointers to input work array to */
/* start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, */
/* 2 forward, 3 forward */
/* area: total size of each set */
/* stride: separation between data sets (stride=&out-&outm1) */
/* n: actual number of data sets being convolved */
/* npad: size of padded region needed for "leakage" of foward filters */
/* stridel: stride on output */
/* fptl,lptl: first,last actual points in output set */
/* fsetl,lsetl: first,last actual sets in output */
/* tr: reverse-time filter coefficients */
/* tf: forward-time filter coefficients */
/* Subroutine */ int convsetspo8full_(double *lower, double *out,
double *outm8, double *outm7, double *outm6, double *
outm5, double *outm4, double *outm3, double *outm2,
double *outm1, double *outp1, double *outp2, double *
outp3, double *outp4, double *outp5, double *outp6,
double *outp7, double *outp8, int *area, int *stride,
int *n, int *npad, int *stridel, int *fptl, int *
lptl, int *fsetl, int *lsetl, double *filt)
{
/* System generated locals */
int i__1, i__2;
/* Local variables */
int i__, j, out2s;
double t0, t1, t2, t3, t4, t5, t6, t7, t8;
int jl, out2;
/* Forward phase (padding now ensures room!) */
/* Relative to convsetspo, only change is accumulation into lower */
/* Parameter adjustments */
--filt;
--outp8;
--outp7;
--outp6;
--outp5;
--outp4;
--outp3;
--outp2;
--outp1;
--outm1;
--outm2;
--outm3;
--outm4;
--outm5;
--outm6;
--outm7;
--outm8;
--out;
--lower;
/* Function Body */
t0 = filt[1];
t1 = filt[2];
t2 = filt[3];
t3 = filt[4];
t4 = filt[5];
t5 = filt[6];
t6 = filt[7];
t7 = filt[8];
t8 = filt[9];
if (t0 != (double)0.) {
/* symmetric filter */
i__1 = *lsetl;
for (i__ = *fsetl + 1; i__ <= i__1; i__ += 2) {
out2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__2 = out2 + 1 + *lptl;
for (j = out2 + 1 + *fptl; j <= i__2; ++j) {
jl = j + out2s;
lower[jl] = t7 * outm7[j] + t5 * outm5[j] + t3 * outm3[j] +
t1 * outm1[j] + t1 * outp1[j] + t3 * outp3[j] + t5 *
outp5[j] + t7 * outp7[j];
}
}
i__1 = *lsetl;
for (i__ = *fsetl; i__ <= i__1; i__ += 2) {
out2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__2 = out2 + 1 + *lptl;
for (j = out2 + 1 + *fptl; j <= i__2; ++j) {
jl = j + out2s;
lower[jl] = t8 * outm8[j] + t6 * outm6[j] + t4 * outm4[j] +
t2 * outm2[j] + t0 * out[j] + t2 * outp2[j] + t4 *
outp4[j] + t6 * outp6[j] + t8 * outp8[j];
}
}
/* end of symmetric part */
} else {
/* antisymmetric filter */
i__1 = *lsetl;
for (i__ = *fsetl + 1; i__ <= i__1; i__ += 2) {
out2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__2 = out2 + 1 + *lptl;
for (j = out2 + 1 + *fptl; j <= i__2; ++j) {
jl = j + out2s;
lower[jl] = t7 * (double)-1. * outm7[j] - t5 * outm5[j] - t3 *
outm3[j] - t1 * outm1[j] + t1 * outp1[j] + t3 * outp3[
j] + t5 * outp5[j] + t7 * outp7[j];
}
}
i__1 = *lsetl;
for (i__ = *fsetl; i__ <= i__1; i__ += 2) {
out2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__2 = out2 + 1 + *lptl;
for (j = out2 + 1 + *fptl; j <= i__2; ++j) {
jl = j + out2s;
lower[jl] = t8 * (double)-1. * outm8[j] - t6 * outm6[j] - t4 *
outm4[j] - t2 * outm2[j] + t0 * out[j] + t2 * outp2[j]
+ t4 * outp4[j] + t6 * outp6[j] + t8 * outp8[j];
}
}
/* end of antisymmetric part */
}
return 0;
} /* convsetspo8full_ */
/* Back and forth, linear convolution of sets of points for +/- 8 */
/* filters using padding, ignoring input data on non-ghost points */
/* (IN/+=OUT) */
/* lower: output array for accumulation */
/* outm3,outm2,outm1,outp1,outp2,outp3: pointers to input work array to */
/* start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, */
/* 2 forward, 3 forward */
/* area: total size of each set */
/* stride: separation between data sets (stride=&out-&outm1) */
/* n: actual number of data sets being convolved */
/* npad: size of padded region needed for "leakage" of foward filters */
/* stridel: stride on output */
/* fptl,lptl: first,last actual points in output set */
/* fsetl,lsetl: first,last actual sets in output */
/* tr: reverse-time filter coefficients */
/* tf: forward-time filter coefficients */
/* Subroutine */ int convsetspo8_(double *lower, double *out,
double *outm8, double *outm7, double *outm6, double *
outm5, double *outm4, double *outm3, double *outm2,
double *outm1, double *outp1, double *outp2, double *
outp3, double *outp4, double *outp5, double *outp6,
double *outp7, double *outp8, int *area, int *stride,
int *n, int *npad, int *stridel, int *fptl, int *
lptl, int *fsetl, int *lsetl, double *tr, double *tf)
{
/* System generated locals */
int i__1, i__2;
/* Local variables */
int i__, j, out2s;
double t0, t1, t2, t3, t4, t5, t6, t7, t8;
int jl, out2;
/* Reverse direction */
/* Relative to convsets, we wiped out the non-ghost references */
/* Parameter adjustments */
--tf;
--tr;
--outp8;
--outp7;
--outp6;
--outp5;
--outp4;
--outp3;
--outp2;
--outp1;
--outm1;
--outm2;
--outm3;
--outm4;
--outm5;
--outm6;
--outm7;
--outm8;
--out;
--lower;
/* Function Body */
t0 = tr[1];
t1 = tr[2];
t2 = tr[3];
t3 = tr[4];
t4 = tr[5];
t5 = tr[6];
t6 = tr[7];
t7 = tr[8];
t8 = tr[9];
/* Pad region */
out2 = (*npad - 1) * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = (double)0.;
}
out2 = (*npad - 2) * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t8 * outm8[j];
}
out2 = (*npad - 3) * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t7 * outm7[j];
}
out2 = (*npad - 4) * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t6 * outm6[j] + t8 * outm8[j];
}
out2 = (*npad - 5) * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t5 * outm5[j] + t7 * outm7[j];
}
out2 = (*npad - 6) * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t4 * outm4[j] + t6 * outm6[j] + t8 * outm8[j];
}
out2 = (*npad - 7) * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t3 * outm3[j] + t5 * outm5[j] + t7 * outm7[j];
}
out2 = (*npad - 8) * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t2 * outm2[j] + t4 * outm4[j] + t6 * outm6[j] + t8 * outm8[j]
;
}
/* Main reverse phase exploiting zeros! */
for (i__ = *n - 1; i__ >= 8; i__ += -2) {
out2 = i__ * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t7 * outm7[j] + t5 * outm5[j] + t3 * outm3[j] + t1 *
outm1[j];
}
out2 = (i__ - 1) * *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j] + t2 * outm2[j] + t4 * outm4[j] + t6 * outm6[
j] + t8 * outm8[j];
}
}
/* Fence posting */
out2 = *stride * 7;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t1 * outm1[j] + t3 * outm3[j] + t5 * outm5[j] + t7 * outm7[j]
;
}
out2 = *stride * 6;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j] + t2 * outm2[j] + t4 * outm4[j] + t6 * outm6[j];
}
out2 = *stride * 5;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t1 * outm1[j] + t3 * outm3[j] + t5 * outm5[j];
}
out2 = *stride << 2;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j] + t2 * outm2[j] + t4 * outm4[j];
}
out2 = *stride * 3;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t1 * outm1[j] + t3 * outm3[j];
}
out2 = *stride << 1;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j] + t2 * outm2[j];
}
out2 = *stride;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t1 * outm1[j];
}
out2 = 0;
i__1 = *area + out2;
for (j = out2 + 1; j <= i__1; ++j) {
out[j] = t0 * out[j];
}
/* Forward phase (padding now ensures room!) */
/* No change here relative to convsets */
t0 = tf[1];
t1 = tf[2];
t2 = tf[3];
t3 = tf[4];
t4 = tf[5];
t5 = tf[6];
t6 = tf[7];
t7 = tf[8];
t8 = tf[9];
i__1 = *lsetl;
for (i__ = *fsetl; i__ <= i__1; ++i__) {
out2 = i__ * *stride;
out2s = i__ * (*stridel - *stride);
i__2 = out2 + 1 + *lptl;
for (j = out2 + 1 + *fptl; j <= i__2; ++j) {
jl = j + out2s;
/* $$$ out(j)= */
/* $$$ $ t0*out(j)+t1*outp1(j)+t2*outp2(j)+t3*outp3(j) */
lower[jl] = t0 * out[j] + t1 * outp1[j] + t2 * outp2[j] + t3 *
outp3[j] + t4 * outp4[j] + t5 * outp5[j] + t6 * outp6[j]
+ t7 * outp7[j] + t8 * outp8[j];
}
}
return 0;
} /* convsetspo8_ */
syntax highlighted by Code2HTML, v. 0.9.1