In[]:=
deploy
Mon 10 May 2021 16:15:39
In[]:=
On[Assert];(*Partitionsmatrixintoblocks{{axa,axb},{bxa,bxb}}*)partitionMatrix[mat_,{a_,b_}]:=Module[{},Assert[a+bLength@mat];Assert;[a+bLength@mat];Internal`PartitionRagged[mat,{{a,b},{a,b}}]];randomRotation2[n_,eps_]:=Module[{z,q,r,d,ph},z=IdentityMatrix[n]+epsRandomVariate[NormalDistribution[0,1],{n,n}];{q,r}=QRDecomposition[z];d=Diagonal[r];ph=d/Abs[d];q*ph];(*GeneratesrealisticsetofA,B,C*)prepareRealistic[d_]:=(diag=Table[1/k,{k,1,d/2}];diag=DiagonalMatrix[diag~Join~diag];rot=randomRotation2[d,.1];sigma=rot.diag.Inverse[rot];{{AA,AB},{BA,BB}}=partitionMatrix[sigma,{d/2,d/2}];b0=BB;a0=AA;c0=BA;d0=AB;e0=RandomReal[{-1,1},{d/2,d/2}];A0=Inverse[c0].b0;B0=d0.Inverse[a0];C0=Inverse[c0].e0.Inverse[a0];{A0,B0,C0});prepareRandom[d_]:=(a=RandomReal[{-3,3},{d,d}];b=RandomReal[{-3,3},{d,d}];c=RandomReal[{-3,3},{d,d}];)(*SolveAX+XB=Cusingdiagonalization*)sylvesterDiag[A_,B_,C0_]:=Module[{da,db,DA,T,DB,U,denom,cutoff,sdiv,Y},{da,db}=Length/@{A,B};{DA,T}=Eigensystem[A+$MachineEpsilonIdentityMatrix[da]];T=Transpose[T];{DB,U}=Eigensystem[B+$MachineEpsilonIdentityMatrix[db]];U=Transpose[U];denom=Outer[Plus,DA,DB];cutoff=Max@Abs[denom]**$MachineEpsilon;sdiv=Map[If[Abs[#]>cutoff,1/#,#]&,denom,{2}];Y=Inverse[T].C0.U*sdiv;T.Y.Inverse[U]];rr[a_,b_]:=RandomInteger[{-5,5},{a,b}];{a,b,c}={rr[3,3],rr[4,4],rr[3,4]};Assert[Norm[sylvesterDiag[a,b,c]-LyapunovSolve[a,b,c]]<](*SolveAX+X'B=CusingdiagonalizationandEddyreduction*)tSylvesterDiagEddy[a_,b_,c_]:=(g=a+b;ig=Inverse[g];h=(c+c)/2;u=sylvesterDiag[a.ig,-ig.b,c-a.ig.h-h.ig.b];u=(u-u)/2;x=ig.(h+u);Print["error is ",Norm[a.x+Transpose[x].b-c]];x);(*SolveAX+X'B=Cusingbuilt-inandEddyreduction*)tSylvesterBuiltinEddy[a_,b_,c_]:=(g=a+b;ig=Inverse[g];h=(c+c)/2;u=LyapunovSolve[a.ig,-ig.b,c-a.ig.h-h.ig.b];u=(u-u)/2;x=ig.(h+u);Print["error is ",Norm[a.x+Transpose[x].b-c]];x);(*SolveAX+X'B=Cusingbuilt-inandDopicoreduction*)tSylvesterBuiltinDopico[a_,b_,c_]:=Module[{},dims=Length[a];X=LyapunovSolve[Inverse[b].a,-Inverse[a].b,Inverse[b].c-Inverse[b].c.Inverse[a].b];Print["error is ",Norm[a.X+Transpose[X].b-c]];X];(*SolveAX+X'B=CusingdiagonalizationandEddyreduction*)tSylvesterDiagDopico[a_,b_,c_]:=Module[{},dims=Length[a];X=sylvesterDiag[Inverse[b].a,-Inverse[a].b,Inverse[b].c-Inverse[b].c.Inverse[a].b];Print["error is ",Norm[a.X+Transpose[X].b-c]];X];d=1000;{a,b,c}=prepareRealistic[d];Print["eddy+spectral"];tSylvesterDiagEddy[a,b,c];//TimingPrint["eddy+builtin"];tSylvesterBuiltinEddy[a,b,c];//TimingPrint["dopico+spectral"];tSylvesterDiagDopico[a,b,c];//TimingPrint["dopico+builtin"];tSylvesterBuiltinDopico[a,b,c];//Timing
6
10
-10
10
eddy+spectral
error is 0.0000721691
Out[]=
{2.28724,Null}
eddy+builtin
error is 0.0000664644
Out[]=
{2.63042,Null}
dopico+spectral
error is 0.00453468
Out[]=
{2.10813,Null}
dopico+builtin
error is 0.0065095
Out[]=
{2.70095,Null}