using System;
using static System.Console;
public static partial class matlib{

public static (double,vector) RG(matrix H,vector x,matrix N){
	vector Hx=H*x;
	vector Nx=N*x;
	double xNx=x%Nx;
	double R=x%Hx/xNx;
	vector G=(2/xNx)*(Hx-R*Nx);
/*
	matrix B=(2/xNx)*(H-R*N);
	B.update(G,Nx,-2/xNx);
	B.update(Nx,G,-2/xNx);
*/
	return (R,G);
	}

public static (double,vector) rayleigh
(matrix H, vector x, matrix N, double acc=1e-3){
	Func<vector,(double,vector)> F = v => RG(H,v,N);
	return gnewton(F,v,acc);
}

public static (double,vector) gnewton
(Func<vector,(double,vector)> F, vector x, double acc=1e-3){
	int nsteps=0,n=x.size;
	matrix B=matrix.id(n);
	(double fx,vector gx)=F(x);
	while(nsteps<1000){
		if(gx.norm()<acc){
			Error.Write($"broyden: job done: |gx|<acc\n");
			break;
			}
		nsteps++;
		vector Dx=-B*gx;
		vector z,gz;
		double fz,lambda=1;
		while(true){// backtracking linesearch
			z=x+Dx*lambda;
			(fz,gz)=F(z);
			if(fz<fx){
				break; // good step
				}
			if(lambda<1e-3){
				Error.WriteLine($"error: 1/lambda={1/lambda}");
				B.setid();
				break; // accept anyway
				}
			lambda/=2;
		}
		Error.WriteLine($"1/lambda={1/lambda}");

		vector s=z-x;
		vector y=gz-gx;
		vector u=s-B*y;
		double uy=u%y;
		if(Abs(uTy)>1e-6){
			B.update(u,u,1/uTy); // SR1 update
		}

		double normz;
		if(N!=null) normz=Math.Sqrt(z%(N*z));
		else normz=z.norm();
		double max=2;
		if(normz>max || normz<1/max){
			Error.WriteLine($"nsteps={nsteps} normz={normz}");
			x=z/normz;
			(fx,gx,Bx)=RG(H,x,N);
		}
		else { x=z; fx=fz; gx=gz;}
	}
	WriteLine($"{H.size1} {nsteps}");
	return (fx,x);
}//SR1

}//class