/* This is the template file for HW4 problem 4.
 *
 * To compile, use UNIX command 
 *      cc -g -Wall -o hw4_4 hw4_4.c -llapack -lm 
 * This will generate an executable hw4_4. The -g option is optional and 
 * is needed only if you need to debug the program in a debugger such as ddd.
 % The -Wall option would enable compiler's warning messages.
 *
 * Run the program 
 *      hw4_4 
 * would create two MATLAB/Octave files "results_hw4_4a.m" and "results_hw4_4b.m".
 *
 * To produce the plot, start octave by issuing UNIX command 
 *      octave
 * Then within octave, issue commands
 *      results_hw4_4a
 *      results_hw4_4b
 *      exit
 * This will produce an EPS file results_hw4_4[ab].eps. Alternatively, 
 * you can also run Octave in batch mode by issuing the UNIX command
 *     octave -q -f --eval "results_hw4_4a; results_hw4_4ab; exit"
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

#define N     100
#define MFNAME1   "results_hw4_4a.m"
#define MFNAME2   "results_hw4_4b.m"
#define EPSFNAME1 "results_hw4_4a.eps"
#define EPSFNAME2 "results_hw4_4b.eps"
#define TITLE1    "Computer Problem 7.6(a)-(c) -- Interpolate to Square Root"
#define TITLE2    "Computer Problem 7.6(d) -- Interpolate to Square Root"

/* Function prototypes */
void polyfit8(const double *t, const double *y, double cs[9]);
void polyval8(const double p[9], int n, const double* ts, double* yp);

void spline(const double *t, const double *y, double pp[8][4]);
void ppval(double pp[8][4], const double* t_input, int n, const double* ts, double* ys);

void  write_plot_script_1(const char *mfname, const char *epsfname, 
                          const char *title, const double *t, const double *y,
                          const double *ts, const double *yr, const double *yp, 
                          const double * ys, int n);

void  write_plot_script_2(const char *mfname, const char *epsfname, 
                          const char *title, const double *t, const double *y,
                          const double *ts, const double *yr, const double *yp, 
                          const double * ys, int n);

int solve9x9(double A[9][9], double x[9]);

int solve32x32(double A[32][32], double x[32]);

int main(int argc, char **argv) {
    
    double t_input[] = {0, 1, 4, 9, 16, 25, 36, 49, 64};
    double y_input[] = {0, 1, 2, 3,  4,  5,  6,  7,  8};

    /* cs is an array containing the coefficients of the polynomial
       in ascending power for monomial basis functions:
       p(t) = cs(0) + cs(1)*t + cs(2)*t^2 + ... + cs(8)*t^8 */
    double cs[9];

    /* PP is an 8x4 matrix containing the eight piecewise polynomial coefficients in ascending powers,
       f(t) = pp[i][0] + pp[i][1]*t + pp[i][2]*t^2 + pp[i][3]*t^3, i is in [0, 7] */
    double pp[8][4];

    double ts[N]; /* x-value data points */
    double yr[N]; /* the y-values by the build-in sqrt function */
    double yp[N]; /* the y-values by the interpolation polynomial of degree eight */
    double ys[N]; /* the y-values by the interpolation cubic spline function */

    int i;
    double h;

    /* Computer Problem 7.6(a)-(c) -- Interpolate to Square Root */
    /* Compute the corresponding values given by the built-in sqrt function over the domain [0, 64] */
    for (i=0, h=64.0/(N-1); i<N; i++)
    {
        ts[i] = i * h;
        yr[i] = sqrt(ts[i]);
    }
	
    /* 1. Polynomial interpolation */
    /* Compute the polynomial of degree eight that interpolates these nine data points. */
    polyfit8(t_input,  y_input, cs);

    /* POLYVAL8 computes the values yp[] of a polynomial p evaluated at ts[].*/
    polyval8(cs, N, ts, yp);


    /* 2. Cubic spline routine */
    /* Compute the piecewise polynomial form of the cubic spline that interpolates these nine data points. */
    spline(t_input,  y_input, pp);

    /* PPVAL returns the values yp[] of a piecewise polynomial pp evaluated at ts[].*/
    ppval(pp, t_input, N, ts, ys);

    /* Write results into an Octave script for plotting */
    write_plot_script_1( MFNAME1, EPSFNAME1, TITLE1, 
                         t_input, y_input, ts, yr, yp, ys, N);

    /* Computer Problem 7.6(d) -- Interpolate to Square Root */
    for (i=0, h=1.0/(N-1); i<N; i++)
    {
        ts[i] = i * h;
        yr[i] = sqrt(ts[i]);
    }

    polyval8(cs, N, ts, yp);
    ppval(pp, t_input, N, ts, ys);

    write_plot_script_2(MFNAME2, EPSFNAME2, TITLE2, 
                        t_input, y_input, ts, yr, yp, ys, N);

    return 0;
}

/*************************************************************
 * FUNCTION: polyfit8
 *
 * Compute the polynomial of degree eight that interpolates
 * these nine data points.
 * The coefficients of the interpolating polynomial are stored
 * in the array cs in ascending powers, 
 *     f(t) = cs(0) + cs(1)*t + cs(2)*t^2 + ... + cs(8)*t^8
 *************************************************************/
void polyfit8(const double *t, const double *y, double cs[9])
{
    double mat[9][9];
    int i, j;
    double s;
   
    /* Construct 9x9 linear system. 
       We must store matrix in column major to be compatible with LAPACK. */
    for (j=0; j<9; ++j) {
        mat[0][j] = 1;
        cs[j] = y[j];
    }

    for (j=0; j<9; ++j) {
        /* We must scale the matrix to make it stable. */
        s = t[j]/64.;
        for (i=1; i<9; ++i) {
            mat[i][j] = mat[i-1][j] * s;
        }
    }

    /* Solve linear system. */
    solve9x9( mat, cs);

    /* Rescale the coefficients */
    s = 1;
    for (i=1; i<9; ++i) {
        s /= 64.;
        cs[i] *= s;
    }
}

/*************************************************************
 * FUNCTION: polyval8
 *
 * Evaluate polynomial of degree eight at given values
 *************************************************************/
void polyval8(const double cs[9], int n, const double* ts, double* yp)
{
    int i, j;

    /* FIXME: Implement this function to evaluate polynomial interpolation. */
    assert(0);
    for (i=0; i<n; ++i) {
        /* Evaluate using Horner's nested formula */

    }
}

/*************************************************************
 * FUNCTION: spline
 *
 * Compute coefficients of the cubic spline that interpolates 
 * the given data points. The computed coefficients are stored
 * in the 2-dimensional array pp[][4] in ascending powers.
 *************************************************************/
void spline(const double *t, const double *y, double pp[8][4])
{
    double mat[32][32], b[32];
    int i, j;
    double s;

    /* Initialize mat to 0 */
    memset( &mat[0][0], 0, sizeof(double)*32*32);

    /* Construct 32x32 linear system. Note that mat is supposed to be a sparse 
       matrix, but it is small here so we use a dense matrix for simplicity.
       We store matrix in column major to be compatible with LAPACK. */
    for (j=0; j<8; ++j) {
        /* Require interpolation of left vertex  */
        s = 1.;
        for (i=0; i<4; ++i) {
            mat[4*j+i][4*j] = s;
            s *= t[j];
        }
        b[4*j] = y[j];

        /* Require interpolation of right vertex  */
        s = 1.;
        for (i=0; i<4; ++i) {
            mat[4*j+i][4*j+1] = s;
            s *= t[j+1];
        }
        b[4*j+1] = y[j+1];

        if (j>0) {
            /* Require C^1 continuity at left vertex  */
            s = t[j];

            mat[4*(j-1)+1][4*j+2] = 1;
            mat[4*j+1][4*j+2] = -1;

            mat[4*(j-1)+2][4*j+2] = 2*s;
            mat[4*j+2][4*j+2] = -2*s;

            mat[4*(j-1)+3][4*j+2] = 3*s*s;
            mat[4*j+3][4*j+2] = -3*s*s;
            b[4*j+2] = 0;

            /* Require C^2 continuity at left vertex  */
            mat[4*(j-1)+2][4*j+3] = 2;
            mat[4*j+2][4*j+3] = -2;

            mat[4*(j-1)+3][4*j+3] = 6*s;
            mat[4*j+3][4*j+3] = -6*s;
            b[4*j+3] = 0;
        }
        else {
            /* Require second derivative to be 0 at two ends. */
            mat[2][2] = 2;
            mat[3][2] = 0;
            b[2] = 0;

            mat[4*7+2][3] = 2;
            mat[4*7+3][3] = 6*t[8];
            b[3] = 0;
        }
    }

    /* Solve linear system. */
    memcpy( &pp[0][0], b, sizeof(double)*32);
    solve32x32( mat, &pp[0][0]);
}


/*************************************************************
 * FUNCTION: ppval
 *
 * Evaluate cubic spline at given values
 *************************************************************/
void ppval(double pp[8][4], const double *t_input, 
           int n, const double* ts, double* ys)
{
    int i,j,k=0;

    /* FIXME: Implement this function to evaluate cubic spline. */
    assert(0);
    for (i=0; i<n; ++i) {
        /* Determine proper interval for evaluation */


        /* Evaluate using Horner's nested formula */

    }
}

/*************************************************************
 * FUNCTION: write_plot_script
 *
 * Write an Octave script for plotting the results
 *************************************************************/
void  write_plot_script_1(const char *mfname, const char *epsfname, 
                          const char *title, const double *t, const double *y,
                          const double *ts, const double *yr, const double *yp, 
                          const double * ys, int n) {
    int i;
    FILE *fid;

    fid = fopen(mfname, "w");

    /* Print out variables */
    fprintf(fid, "t=[");
    for (i=0; i<n-1; i++) {
	fprintf(fid, "%.16g; \n", t[i]);
    }
    fprintf(fid, "%.16g];\n\n", t[n-1]);

    fprintf(fid, "y=[");
    for (i=0; i<n-1; i++) {
	fprintf(fid, "%.16g; \n", y[i]);
    }
    fprintf(fid, "%.16g];\n\n", y[n-1]);

    fprintf(fid, "ts=[");
    for (i=0; i<n-1; i++) {
	fprintf(fid, "%.16g; \n", ts[i]);
    }
    fprintf(fid, "%.16g];\n\n", ts[n-1]);

    fprintf(fid, "yr=[");
    for (i=0; i<n-1; i++) {
	fprintf(fid, "%.16g; \n", yr[i]);
    }
    fprintf(fid, "%.16g];\n\n", yr[n-1]);

    fprintf(fid, "yp=[");
    for (i=0; i<n-1; i++) {
	fprintf(fid, "%.16g; \n", yp[i]);
    }
    fprintf(fid, "%.16g];\n\n", yp[n-1]);

    fprintf(fid, "ys=[");
    for (i=0; i<n-1; i++) {
	fprintf(fid, "%.16g; \n", ys[i]);
    }
    fprintf(fid, "%.16g];\n\n", ys[n-1]);

    /* Write out statements for plotting */
    fprintf(fid, "%% Plot results.\n");
    fprintf(fid, "plot(t, y, 'ko', ts, yr, 'k-', ts, yp, 'b--', ts, ys, 'r-.');\n");
    fprintf(fid, "title('%s');\n", title);
    fprintf(fid, "xlabel('t');\n");
    fprintf(fid, "ylabel('y');\n");
    fprintf(fid, "axis([0, 64, 0, 8]);\n");
    fprintf(fid, "legend('data points', 'square root', 'polynomial', 'spline', 0);\n");

    fprintf(fid, "%% Print figure into EPS file \n");
    fprintf(fid, "print -depsc %s\n", epsfname);

    fclose(fid);

    return;
}

/*************************************************************
 * FUNCTION: write_plot_script
 *
 * Write an Octave script for plotting the results
 *************************************************************/
void  write_plot_script_2(const char *mfname, const char *epsfname, 
                          const char *title, const double *t, const double *y,
                          const double *ts, const double *yr, const double *yp, 
                          const double * ys, int n) {
    int i;
    FILE *fid;

    fid = fopen(mfname, "w");

    /* Print out variables */
    fprintf(fid, "t=[");
    for (i=0; i<n-1; i++) {
	fprintf(fid, "%.16g; \n", t[i]);
    }
    fprintf(fid, "%.16g];\n\n", t[n-1]);

    fprintf(fid, "y=[");
    for (i=0; i<n-1; i++) {
	fprintf(fid, "%.16g; \n", y[i]);
    }
    fprintf(fid, "%.16g];\n\n", y[n-1]);

    fprintf(fid, "ts=[");
    for (i=0; i<n-1; i++) {
	fprintf(fid, "%.16g; \n", ts[i]);
    }
    fprintf(fid, "%.16g];\n\n", ts[n-1]);

    fprintf(fid, "yr=[");
    for (i=0; i<n-1; i++) {
	fprintf(fid, "%.16g; \n", yr[i]);
    }
    fprintf(fid, "%.16g];\n\n", yr[n-1]);

    fprintf(fid, "yp=[");
    for (i=0; i<n-1; i++) {
	fprintf(fid, "%.16g; \n", yp[i]);
    }
    fprintf(fid, "%.16g];\n\n", yp[n-1]);

    fprintf(fid, "ys=[");
    for (i=0; i<n-1; i++) {
	fprintf(fid, "%.16g; \n", ys[i]);
    }
    fprintf(fid, "%.16g];\n\n", ys[n-1]);

    /* Write out statements for plotting */
    fprintf(fid, "%% Plot results.\n");
    fprintf(fid, "plot(t(1:2), y(1:2), 'ko', ts, yr, 'k-', ts, yp, 'b--', ts, ys, 'r-.');\n");
    fprintf(fid, "title('%s');\n", title);
    fprintf(fid, "xlabel('t');\n");
    fprintf(fid, "ylabel('y');\n");
    fprintf(fid, "axis([0, 1, 0, 1]);\n");
    fprintf(fid, "legend('data points', 'square root', 'polynomial', 'spline', 0);\n");

    fprintf(fid, "%% Print figure into EPS file \n");
    fprintf(fid, "print -depsc %s\n", epsfname);

    fclose(fid);

    return;
}

/*****************************************************************************
 * The following are helper functions 
 */

/* Declare function prototype */
extern int dgesv_(int *n, int *nrhs, double *a, int *lda,
                  int *ipiv, double *b, int *ldb, int *info);
/*  -- LAPACK driver routine (version 3.1) --   
    Purpose   
    =======   

    DGESV computes the solution to a real system of linear equations   
       A * X = B,   
    where A is an N-by-N matrix and X and B are N-by-NRHS matrices.   

    The LU decomposition with partial pivoting and row interchanges is   
    used to factor A as   
       A = P * L * U,   
    where P is a permutation matrix, L is unit lower triangular, and U is   
    upper triangular.  The factored form of A is then used to solve the   
    system of equations A * X = B.   

    Arguments   
    =========   

    N       (input) INTEGER   
            The number of linear equations, i.e., the order of the   
            matrix A.  N >= 0.   

    NRHS    (input) INTEGER   
            The number of right hand sides, i.e., the number of columns   
            of the matrix B.  NRHS >= 0.   

    A       (input/output) REAL array, dimension (LDA,N)   
            On entry, the N-by-N coefficient matrix A.   
            On exit, the factors L and U from the factorization   
            A = P*L*U; the unit diagonal elements of L are not stored.   

    LDA     (input) INTEGER   
            The leading dimension of the array A.  LDA >= max(1,N).   

    IPIV    (output) INTEGER array, dimension (N)   
            The pivot indices that define the permutation matrix P;   
            row i of the matrix was interchanged with row IPIV(i).   

    B       (input/output) REAL array, dimension (LDB,NRHS)   
            On entry, the N-by-NRHS matrix of right hand side matrix B.   
            On exit, if INFO = 0, the N-by-NRHS solution matrix X.   

    LDB     (input) INTEGER   
            The leading dimension of the array B.  LDB >= max(1,N).   

    INFO    (output) INTEGER   
            = 0:  successful exit   
            < 0:  if INFO = -i, the i-th argument had an illegal value   
            > 0:  if INFO = i, U(i,i) is exactly zero.  The factorization   
                  has been completed, but the factor U is exactly   
                  singular, so the solution could not be computed.   

    ===================================================================== */

/* Solve 9x9 linear system by calling LAPACK routine dgesv_ */
int solve9x9(double A[9][9], double b[9]) {

    int info,n,nrhs,lda,ipiv[9],ldb;
       
    /* Invoke dgesv_ */
    n = lda = ldb = 9; nrhs = 1;
    dgesv_(&n, &nrhs, &A[0][0], &lda, ipiv, b, &ldb, &info);

    if (info) printf("Error: dgesv returned info %d\n", info);

    return info;
}

/* Solve 32x32 linear system by calling LAPACK routine dgesv_ */
int solve32x32(double A[32][32], double b[32]) {

    int info,n,nrhs,lda,ipiv[32],ldb;
       
    /* Invoke dgesv_ */
    n = lda = ldb = 32; nrhs = 1;
    dgesv_(&n, &nrhs, &A[0][0], &lda, ipiv, b, &ldb, &info);

    if (info) printf("Error: dgesv returned info %d\n", info);

    return info;
}