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

public static (double,vector) raygrad(matrix H,vector x){
	vector Hx=H*x;
	double xx=x%x;
	double R=x%Hx/xx;
	vector G=(2/xx)*(Hx-R*x);
	return (R,G);
	}

public static (double,vector) eigen0
(matrix H, vector x, double acc=1e-6){
	matrix B=matrix.id(x.size);
	(double fx,vector gx)=raygrad(H,x);
	int nsteps=0;
	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.0;
		while(true){// backtracking linesearch
			z=x+Dx*lambda;
			(fz,gz)=raygrad(H,z);
			if(fz<fx){
				break; // good step
				}
			if(lambda<0.01){
				B.setid();
Error.WriteLine($"error: lambda={lambda}");
				break; // accept anyway
				}
			lambda/=2;
		}
Error.WriteLine($"lambda={lambda}");
		vector s=z-x;
		vector y=gz-gx;
		vector u=s-B*y;
		double uy=u%y;

// symmetric Broyden update
		double sy=s%y;
		if(Abs(sy)>1e-6){
		double gamma=uy/sy/2;
		vector a=(u-gamma*s)/sy;
		B.update(a,s);
		B.update(s,a);}
/*
// SR1 update
		B.update(u,u,1/uy);
*/

		double normz=z.norm();
		double max=2;
		if(normz>max || normz<1/max){
Error.WriteLine($"nsteps={nsteps} normz={normz}");
			x=z/normz;
			(fx,gx)=raygrad(H,x);
			//B.setid();
		}
		else { x=z; fx=fz; gx=gz; }
	}
	Error.WriteLine($"H.size={H.size1} nstpes={nsteps}");
	return (fx,x);
}//SR1

}//class