#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <sys/time.h>
#include <sys/resource.h>
#include <string.h>
#include <strings.h>
#include <cilk/cilk.h>
#include <lapacke.h>
#include "matrixlib.h"
#include "tictoc.h"

#define N 2000
double vecnorm1(int n,double x[n]){
    double r=0;
    for(int i=0;i<n;i++){
        r+=fabs(x[i]);
    }
    return r;
}
double matnorm1(int n,double A[n][n]){
    double r=0;
    for(int j=0;j<n;j++){
        double s=0;
        for(int i=0;i<n;i++){
            s+=fabs(A[i][j]);
        }
        if(s>r) r=s;
    }
    return r;
}
double vecnorm2(int n,double x[n]){
    double r=0;
    for(int i=0;i<n;i++){
        r+=x[i]*x[i];
    }
    return sqrt(r);
}
    
double matnorm2(int n,double A[n][n]){
    double B[n][n];

    bzero(B,sizeof(double)*n*n);
    for(int k=0;k<n;k++){
        for(int i=0;i<n;i++){   // B = A^t * A
            for(int j=0;j<n;j++){
                B[i][j]+=A[k][i]*A[k][j];
            }
        }
    }
    printf("here\n");
    double x[n],y[n],yk[n];
    for(int i=0;i<n;i++){
        x[i]=2.0*random()/RAND_MAX+1.0;
    }
    multAx(n,n,A,x,y);
    double normy,normyk;
    for(int k=2;k<200;k++){
        multAx(n,n,B,y,yk);
        normy=vecnorm2(n,y);
        normyk=vecnorm2(n,yk);
        for(int j=0;j<n;j++){
            y[j]=yk[j]/normyk;
        }
//        memcpy(y,yk,sizeof(double)*n);
        printf("%g/%g = %.15e\n",
            normyk,normy,normyk/normy);
    }
    return sqrt(normyk/normy);
}

double A[N][N];
double B[N],X[N],XX[N],*P[N];
int IP[N];
int main(){
    setrlimit(RLIMIT_STACK,
        &(const struct rlimit)
        {RLIM_INFINITY,RLIM_INFINITY});
//    srandom(1235);
    printf("N=%d\n",N);
//    printf("A=\n"); matprint(N,N,A);
    for(int i=0;i<N;i++){
        for(int j=0;j<N;j++){
            A[i][j]=2.0*random()/RAND_MAX-1.0;
        }
        X[i]=1;
    }
    printf("|A|_1=%g\n",matnorm1(N,A));
    double m2=matnorm2(N,A);
    printf("|A|_2=%g\n",m2);
    multAx(N,N,A,X,B);  // B=AX
//    printf("B=\n"); vecprint(N,B);
    tic();
//  PLUfact(N,A,P);     // now A contains LU fact
//  PLUsolve(N,A,P,XX,B);
    memcpy(XX,B,sizeof(double)*N);
    LAPACKE_dgesv(LAPACK_ROW_MAJOR,N,1,&A[0][0],N,IP,XX,1);
    double t=toc();
    printf("The elapsed time was %g seconds.\n",t);
    double flops=2.0*(N/3.0+1)*N*N/t;
    printf("The speed is %g GFLOPS.\n",flops/1e9);
//    printf("XX=\n"); vecprint(N,XX);
    double r=0;
    for(int i=0;i<N;i++){
        double rd=XX[i]-1;
        r+=rd*rd;
    }
    r=sqrt(r);
    printf("The error is %g.\n",r);
    return 0;
}