/*  Programming template for Final Part 2 Question 3

    Refer to Exercise 4.8 from the text and also the Wikipedia entry
    on the Thomas Algorithm. Suppose A is a n x n tridiagonal matrix
    with the vector v along the diagonal, w on the upper diagonal and
    u on the lower diagonal. Thus,

            [ v[0]   w[0]                        ]
            [ u[1]   v[1]   w[1]                 ]
            [        u[2]   v[2]   w[2]          ]
        A = [           .      .      .          ]
            [             .      .      .        ]
            [               .      .      .      ]
            [               u[n-2] v[n-2] w[n-2] ]
            [                      u[n-1] v[n-1] ]

    Assume |u[i]| + |w[i]| < |v[i]| for every i so that Gaussian
    elimination can be performed without any pivoting.  Write a
    subroutine called tdthomas to solve Ax = b in O(n) time.

    Feel free to use the matrix library developed in class or any other
    code that you find useful as part of your program.
*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <sys/time.h>
#include <sys/resource.h>

static double tic_time;
void tic() {
    struct timeval tp;
    gettimeofday(&tp,0);
    tic_time=tp.tv_sec+tp.tv_usec*1.0e-6;
}
double toc() {
    struct timeval tp;
    gettimeofday(&tp,0);
    return (tp.tv_sec+tp.tv_usec*1.0e-6)-tic_time;
}

// Multiply the tridiagonal system Ax to find y assume n>3.
void tdmult(int n,double u[n],double v[n],double w[n-1],
    double x[n],double y[n]){
    y[0]=v[0]*x[0]+w[0]*x[1];
    for(int i=1;i<n-1;i++){
        y[i]=u[i]*x[i-1]+v[i]*x[i]+w[i]*x[i+1];
    }
    y[n-1]=u[n-1]*x[n-2]+v[n-1]*x[n-1];
}

// Solve the tridiagonal system Ax=b for x.
void tdthomas(int n,double u[n],double v[n],double w[n-1],
    double x[n],double b[n]){
    w[0]/=v[0];
    for(int i=1;i<n;i++) w[i]/=v[i]-u[i]*w[i-1];
    b[0]/=v[0];
    for(int i=1;i<n;i++) b[i]=(b[i]-u[i]*b[i-1])/(v[i]-u[i]*w[i-1]);
    x[n-1]=b[n-1];
    for(int i=n-2;i>=0;i--) x[i]=b[i]-w[i]*x[i+1];
}

#define K 10
int main(){
    printf("Final Part 2 Question 3.\n");
    setrlimit(RLIMIT_STACK,
        &(const struct rlimit)
        {RLIM_INFINITY,RLIM_INFINITY});
    int n=512;
    double T[K];
    // Use different values of n to make sure O(n) time complexity.
    for(int k=0;k<K;k++){
        n*=2;
        double ugold[n],vgold[n],wgold[n],bgold[n],xgold[n];
        double u[n],v[n],w[n],b[n],x[n];
        for(int i=0;i<n;i++) {
            ugold[i]=2.0*random()/RAND_MAX-1.0;
            wgold[i]=2.0*random()/RAND_MAX-1.0;
            vgold[i]=fabs(ugold[i])+fabs(wgold[i])+2.0*random()/RAND_MAX;
            if(random()<RAND_MAX/2) vgold[i]=-vgold[i];
            xgold[i]=2.0*random()/RAND_MAX-1.0;
        }
        tdmult(n,ugold,vgold,wgold,xgold,bgold);
        double elapsed=0;
        int c=0;
        // For accurate timing repeat until 0.5 seconds have elapsed.
        while(elapsed<0.5){
            memcpy(u,ugold,sizeof(double)*n);
            memcpy(v,vgold,sizeof(double)*n);
            memcpy(w,wgold,sizeof(double)*n);
            memcpy(b,bgold,sizeof(double)*n);
            tic();
            tdthomas(n,u,v,w,x,b);
               // At this point x should equal xgold.
            elapsed+=toc();
            c++;
        }
        T[k]=elapsed/c;
        double r=0;
        for(int i=0;i<n;i++){
            double t=x[i]-xgold[i];
            r+=t*t;
        }
        r=sqrt(r);
           if(r<1e-7){
            printf("n=%d |x-xgold|=%g T=%g sec\n",n,r,T[k]);
        } else {
            printf("Incorrect!  For n=%d computation of x is wrong!\n",n);
            exit(-1);
        }
    }
    double order=log(T[K-1]/T[K-2])/log(2.0);
    printf("Running time is O(n^%g)\n",order);
    if(fabs(order-1.0)<.2) {
        printf("Correct!  Computations of x are consistent and O(n)!\n");
    }
    return 0;
}