/* This is the template file for HW2 problem 4. 
 *
 * To compile, use UNIX command 
 *      cc -g -Wall -o hw2_4 hw2_4.c -lm
 * This will generate an executable hw2_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 
 *      ./hw2_4 
 * would print out the errors on the screen.
 */

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

/*************************************************************
* FUNCTION: norm
*
* Compute norm of a 2-vector
*************************************************************/
double norm( const double s[2]) {
    return sqrt( s[0]*s[0] + s[1]*s[1]);
}

/*************************************************************
* FUNCTION: mul_matrix_vec
*
* Compute b=A*x for 2x2 case.
*************************************************************/
void mul_matrix_vec( double A[2][2], const double x[2], double b[2]) {
    b[0] = A[0][0]*x[0] + A[0][1]*x[1];
    b[1] = A[1][0]*x[0] + A[1][1]*x[1];
}

/*************************************************************
* FUNCTION: solve2x2
*
* Solve Ax=b for 2x2 case.
*************************************************************/
void solve2x2( double A[2][2], const double b[2], double x[2]) {
    double d = A[0][0]*A[1][1]-A[0][1]*A[1][0];
    
    x[0] = (b[0]*A[1][1]-b[1]*A[0][1])/d;
    x[1] = (A[0][0]*b[1]-A[1][0]*b[0])/d;
}

/*************************************************************
* FUNCTION: f
*
* Evaluate the function at x.
*************************************************************/
void f(const double x[2], double fx[2]) {  
    /* FIXME: Implement this function. */
}

/*************************************************************
* FUNCTION: Df
*
* Evaluate the Jacobian matrix of function f at x.
*************************************************************/
void Df(const double x[2], double J[2][2]) {  
    /* FIXME: Implement this function. */
}

/*************************************************************
* FUNCTION: Newton
*
* Newton's method. Note that x_true in general is not available, 
*      but we use it here to compute errors.
*************************************************************/
void Newton( const double x0[2], const double x_true[2], 
             double tol, int maxiter, double x[2]) {
    int k=0;
    double err[2];

    /* Initialization */
    x[0] = x0[0]; x[1] = x0[1];
    err[0]=x[0]-x_true[0]; err[1] = x[1]-x_true[1];
    printf("%3d %17.10e %17.10e %17.10e\n", k, x[0], x[1], norm(err)); 

    for (k=1; k<=maxiter; ++k) {
        double J[2][2], fx[2], s[2];
        f(x, fx); 
        Df(x, J);

        /* FIXME: Implement this function. */
        /* Update x */

        err[0]=x[0]-x_true[0]; err[1] = x[1]-x_true[1];
        printf("%3d %17.10e %17.10e %17.10e\n", k, x[0], x[1], norm(err)); 

        if (norm(s) <= tol) break;
    }
}


/*************************************************************
* FUNCTION: Broyden
*
* Broyden's method. Note that x_true in general is not available, 
*      but we use it here to compute errors.
*************************************************************/
void Broyden( const double x0[2], const double x_true[2], 
              double tol, int maxiter, double x[2]) {
    int k=0;
    double err[2], B[2][2], fxnext[2];

    /* Initialization */
    x[0] = x0[0]; x[1] = x0[1];
    err[0]=x[0]-x_true[0]; err[1] = x[1]-x_true[1];
    printf("%3d %17.10e %17.10e %17.10e\n", k, x[0], x[1], norm(err)); 

    Df(x0, B); f(x0, fxnext);
    for (k=1; k<=maxiter; ++k) {
        double fx[2], s[2], y[2], t[2], s_sqr;

        /* FIXME: Implement this function. */
        /* Update x */

        /* Update B */

        err[0]=x[0]-x_true[0]; err[1] = x[1]-x_true[1];
        printf("%3d %17.10e %17.10e %17.10e\n", k, x[0], x[1], norm(err)); 

        if (norm(s) <= tol) break;
    }
}

int main(int argc, char **argv) {
    
    double x0[2] = {-0.5, 1.4}, x_true[2] = {0,1}, x[2];
    double tol = 2.220446049250313e-16; /* Machine precision */
    int maxiter = 20;

    /* First, run Newton's method */
    printf("Newton's method\n");
    Newton( x0, x_true, tol, maxiter, x);

    /* Second, run Broyden's method */
    printf("Broyden's method\n");
    Broyden( x0, x_true, tol, maxiter, x);

    return 0;
}