#include <stdio.h>
#include <stdlib.h>
#include <strings.h>
#include <string.h>
#include <math.h>
#include <lapacke.h>
#include "matlib.h"

#define N 4
int main(){
    double A[N][N],B[N][N];
    mkA(N,A,12);
    printf("A=[\n");
    pmatrix(N,A);
    AtA(N,B,A);
    printf("]\nB=[\n");
    pmatrix(N,B);
    printf("]\n");
    double x[N],y[N],lambda;
    for(int i=0;i<N;i++) x[i]=1;
    int ipv[N];
    int info=LAPACKE_dgetrf(LAPACK_ROW_MAJOR,
        N,N,B[0],N,ipv);
    for(int c=0;c<10;c++){
        memcpy(y,x,sizeof(double)*N);
        info=LAPACKE_dgetrs(LAPACK_ROW_MAJOR,
            'N',N,1,B[0],N,ipv,y,1);
        double t=sqrt(dotprod(N,y,y));
        lambda=dotprod(N,x,y)/dotprod(N,x,x);
        printf("lambda=%22.14e\n",lambda);
        for(int i=0;i<N;i++){
            x[i]=y[i]/t;
        }
    }
    printf("|inv-A|=%22.14e\n",sqrt(lambda));
    return 0;
}