(* ::Package:: *)

BeginPackage["Eigenvalues`"]


Maximum::usage = "Maximum[u_]"
PowerMethod::usage = "PowerMethod[{A_,itmax_, toll_,t_]"
InversePowerMethod::usage = "PowerMethod[{A_,itmax_, toll_,t_]"
InversePowerMethodImprove::usage = "InversePowerMethodImprove[A_,mu_,itmax_, toll_,t_]"

HouseHolderMatrix::usage = "HouseHolderMatrix[t_]"
HouseHolderDeflatedMatrix::usage = "HouseHolderDeflatedMatrix[A_,x_]"


Off[General::spelll]


Begin["Private`"]


Maximum[u_]:= Module[{m,pos,i},
 m=u[[1]];
 pos=1;
 For [i=1,i<= Length[u],i++,
     If[u[[i]]>= m,
        m=u[[i]];
        pos=i;
     ];
 ];
 Return[{m,pos}];
];


PowerMethod[A_,itmax_, toll_,x_] := Module[{t=x,k,err,lamv,lam,u,m,mas},
    k=1;
	err=1;
	lamv=0;
    While[k<itmax && err>toll, 
         u=Flatten[A.Transpose[{t}]];
         {mas,m}=Maximum[Abs[u]];
         lam=u[[m]];
         t=u/lam;
         err=Abs[lamv-lam]/Abs[lam];
         lamv=lam;
         k=k+1;
	];
    k=k-1; 
    Return[{lam,t,k}];
     
];
    
InversePowerMethod[A_,itmax_, toll_,x_] := Module[{t=x,n,m,mas,k,err,lamv,lam,u,y,L,U,lu,p,c},
    {n,m}=Dimensions[A]; 
    {lu,p,c}=LUDecomposition[A];
    L=lu SparseArray[{i_,j_}/;j<i->1,{n,n}]+IdentityMatrix[n];
    U=lu SparseArray[{i_,j_}/;j>=i->1,{n,n}];
    k=1;
	err=1;
	lamv=0;
	While[k<itmax && err>toll, 
         y=LinearSolve[L,t[[p]]];
         u=LinearSolve[U,y];
         {mas,m}=Maximum[Abs[u]];
         lam=1/u[[m]];
         t=u/u[[m]];
         err=Abs[lamv-lam]/Abs[lam];
         lamv=lam;
         k=k+1;
	];
    Return[{lam,t,k}];
     
];   

InversePowerMethodImprove[A_,mu_,itmax_, toll_,x_] := Module[{t=x,n,m,mas,k,err,lamv,lam,u,y,L,U,lu,p,c},
    {n,m}=Dimensions[A]; 
    {lu,p,c}=LUDecomposition[A-mu IdentityMatrix[n]];
    L=lu SparseArray[{i_,j_}/;j<i->1,{n,n}]+IdentityMatrix[n];
    U=lu SparseArray[{i_,j_}/;j>=i->1,{n,n}];
    k=1;
	err=1;
	lamv=0;
	While[k<itmax && err>toll, 
         y=LinearSolve[L,t[[p]]];
         u=LinearSolve[U,y];
         {mas,m}=Maximum[Abs[u]];
         lam=mu+1/u[[m]];
         t=u/u[[m]];
         err=Abs[lamv-lam]/Abs[lam];
         lamv=lam;
         k=k+1;
	];
    Return[{lam,t,k}];
     
];   


HouseHolderMatrix[t_]:=Module[{x=t,n,m,eta,sigma,beta},
	n=Length[x];
    {eta,m}=Maximum[Abs[x]];
	sigma=0;
    For [i=1,i<= n,i++,
         If[ Abs[x[[i]]]>= Sqrt[$MachineEpsilon] eta,
            sigma=sigma+(x[[i]]/eta)^2;
         ];
    ];
    sigma=Sign[x[[1]]] Sqrt[sigma] eta;
    x[[1]]=x[[1]]+sigma;
    beta=sigma x[[1]];
    Return[{x,beta,sigma}];
];


HouseHolderDeflatedMatrix[A_,x_]:=Module[{n,m,beta,sigma,B,C,tau,i,j,u},

(* x is the dominant eigenvector di A *)

	{n,m}=Dimensions[A];
	{u,beta,sigma} =HouseHolderMatrix[x];

(* H=IdentityMatrix[n]-u*u^T /beta;*)

(*Compute C=(HA')'*)
	B=Transpose[A];
    For [i=1, i<= n, i++,
         tau=u.Transpose[B][[i]]/beta;
         For [j=1,j<= n,j++,
             B[[j,i]]=B[[j,i]]-tau u[[j]];
         ];        
    ];
   C=Transpose[B];
(* Compute HC=HAH *)
    For [i=1, i<= n,i++,
         tau=u.Transpose[C][[i]]/beta;
         For [j=1,j<= n,j++,
             C[[j,i]]=C[[j,i]]-tau u[[j]];
         ];  
    ];
    A2=Take[C,{2,n},{2,n}];

	Return[A2];

];


End[]


EndPackage[]