/*************************************************************************\
* Copyright (c) 2002 The University of Chicago, as Operator of Argonne
* National Laboratory.
* Copyright (c) 2002 The Regents of the University of California, as
* Operator of Los Alamos National Laboratory.
* This file is distributed subject to a Software License Agreement found
* in the file LICENSE that is included with this distribution. 
\*************************************************************************/

/* cost.f -- translated by f2c (version of 30 January 1990  16:02:04).
   You must link the resulting object file with the libraries:
	-lF77 -lI77 -lm -lc   (in that order)
 $Log: cost.c,v $
 Revision 1.3  2002/08/14 16:41:18  soliday
 Added Open License

 Revision 1.2  1995/09/05 21:12:59  saunders
 First test release of the SDDS1.5 package.

*/

#include "f2c.h"

/* Subroutine */ int cost_(integer *n, doublereal *x, doublereal *wsave)
{
    /* System generated locals */
    integer i_1;

    /* Local variables */
    static integer modn, i, k;
    extern /* Subroutine */ int rfftf_(integer *, doublereal *, doublereal *);

    static doublereal c1, t1, t2;
    static integer kc;
    static doublereal xi;
    static integer nm1, np1;
    static doublereal x1h;
    static integer ns2;
    static doublereal tx2, x1p3, xim2;

    /* Parameter adjustments */
    --x;
    --wsave;

    /* Function Body */
    nm1 = *n - 1;
    np1 = *n + 1;
    ns2 = *n / 2;
    if ((i_1 = *n - 2) < 0) {
	goto L106;
    } else if (i_1 == 0) {
	goto L101;
    } else {
	goto L102;
    }
L101:
    x1h = x[1] + x[2];
    x[2] = x[1] - x[2];
    x[1] = x1h;
    return 0;
L102:
    if (*n > 3) {
	goto L103;
    }
    x1p3 = x[1] + x[3];
    tx2 = x[2] + x[2];
    x[2] = x[1] - x[3];
    x[1] = x1p3 + tx2;
    x[3] = x1p3 - tx2;
    return 0;
L103:
    c1 = x[1] - x[*n];
    x[1] += x[*n];
    i_1 = ns2;
    for (k = 2; k <= i_1; ++k) {
	kc = np1 - k;
	t1 = x[k] + x[kc];
	t2 = x[k] - x[kc];
	c1 += wsave[kc] * t2;
	t2 = wsave[k] * t2;
	x[k] = t1 - t2;
	x[kc] = t1 + t2;
/* L104: */
    }
    modn = *n % 2;
    if (modn != 0) {
	x[ns2 + 1] += x[ns2 + 1];
    }
    rfftf_(&nm1, &x[1], &wsave[*n + 1]);
    xim2 = x[2];
    x[2] = c1;
    i_1 = *n;
    for (i = 4; i <= i_1; i += 2) {
	xi = x[i];
	x[i] = x[i - 2] - x[i - 1];
	x[i - 1] = xim2;
	xim2 = xi;
/* L105: */
    }
    if (modn != 0) {
	x[*n] = xim2;
    }
L106:
    return 0;
} /* cost_ */