using System;
using static System.Math;
using static System.Console;
using System.Threading;
public class Hessenberg{

public matrix H=null;
public matrix Q=null;

public static void timesJ(matrix A, int p, int q, double c, double s){
        for(int i=0;i<A.size1;i++){
                double new_aip=c*A[i,p]-s*A[i,q];
                double new_aiq=s*A[i,p]+c*A[i,q];
                A[i,p]=new_aip;
                A[i,q]=new_aiq;
                }

	}//timesJ

public static void Jtimes(matrix A, int p, int q, double c,double s){
        for(int j=0;j<A.size2;j++){
                double new_apj= c*A[p,j]+s*A[q,j];
                double new_aqj=-s*A[p,j]+c*A[q,j];
                A[p,j]=new_apj;
                A[q,j]=new_aqj;
                }
	}//Jtimes

public Hessenberg(matrix A){
	H=A.copy();
	int n=A.size1;
	Q = matrix.id(n);
	for(int p=1;p<n;p++)
	for(int q=p+1;q<n;q++) {
		double theta=Atan2(-H[p-1,q],H[p-1,p]);
		double c=Cos(theta);
		double s=Sin(theta);
		timesJ(H,p,q,c,s);
		Jtimes(H,p,q,c,-s);
		timesJ(Q,p,q,c,s);
		}
	}//ctor
}//class Hessenberg

class main{
static void Main(){
	int n=6;
	var A=new matrix(n,n);
	var rnd=new System.Random();
	for(int i=0;i<A.size1;i++)
	for(int j=0;j<A.size2;j++)
		A[i,j]=rnd.NextDouble();
	var hess=new Hessenberg(A);
	hess.H.print("H=");
	var tmp=hess.Q.T*A*hess.Q;
	tmp.print("tmp=");
	tmp=hess.Q.T*hess.Q;
	tmp.print("=1 ?");
	}
}