%Tests csd_qc.m

example = 2
switch example
	case (1)
		U = rand_unitary(8, true);
	case (2)
		U = kron(rand_unitary(2,true), rand_unitary(2, true));
	otherwise
		error("example number is out of range");
end

ns=rows(U);
nsh = ns/2;
U00 = U(1:nsh, 1:nsh);
U01 = U(1:nsh, nsh+1:ns);
U10 = U(nsh+1:ns, 1:nsh);
U11 = U(nsh+1:ns, nsh+1:ns);

[L0,L1,cc,ss,R0,R1]=csd_qc(U);

error=(norm(U00 -L0*cc*R0)
+norm(U01 -L0*ss*R1)
+norm(U10 +L1*ss*R0)
+norm(U11 -L1*cc*R1))


fi = fopen ("csd_out.txt", "w", "native");

fprintf(fi, "error = %10.4e\n", error);

fdisp(fi, "input U ="),fdisp(fi, U);

fdisp(fi, "L0 ="),fdisp(fi, L0);

fdisp(fi, "L1 ="),fdisp(fi, L1);

fdisp(fi, "R0 ="),fdisp(fi, R0);

fdisp(fi, "R1 ="),fdisp(fi, R1);

fdisp(fi, "cc ="),fdisp(fi, diag(cc));

fdisp(fi, "ss ="),fdisp(fi, diag(ss));

fclose(fi);