/*
    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_ */