Test Matrices

I downloaded three FEA matrices (positive definate and probably symmetric) from Matrix Market.  I added stuff to the diagonals to make them diagonally dominant.

PorkDiagonals[A_]:= Module[
{m=Max[Abs[Diagonal[A]]]},
A+SparseArray[Band[{1,1}]->m, {Length[A],Length[A]}]
]
SetDirectory["/home/campus03/struther/Desktop/Link to Courses/5580/CodeFragments/LinearAlgebra/"];
names=FileNames["*mtx.gz"];
{A1,A2,A3}=Sort[Map[Import,names], Length[#1]<Length[#2]&];
{A1,A2,A3}=Map[PorkDiagonals, {A1,A2,A3}];
Map[MatrixPlot,{A1,A2,A3}]
{n1,n2,n3}=Map[Length, {A1,A2,A3}];
{u1,u2,u3}=Map[ Sin[Range[0.0,1,(1/(#-1))]]&, {n1,n2,n3}];
{b1,b2,b3}={A1.u1,A2.u2,A3.u3};
{Norm[LinearSolve[A1,b1]-u1],Norm[LinearSolve[A2,b2]-u2],Norm[LinearSolve[A3,b3]-u3]}

Apparently the condition number of some of the matrices is relatively large.

Jacobi and Related Iterations

The following defines a Jacobi Iteration with a specified order, an inplace Jacobi like Iterate and Gauss-Seidel iteration.  The standard order is just 1,2,3,4,…

JacobiIterate[xIn_,{{Dia_,LPlusU_},b_},order_]:=Module[
{i,j},
Table[
i=order[[j]];
(b[[i]]-LPlusU[[i]].x)/Dia[[i]],
{j,Length[Dia]}] 
]
SetAttributes[JacobiIterateInPlace, HoldAll]
JacobiIterateInPlace[x_,{{Dia_,LPlusU_},b_},order_]:=Module[
{i,j},
Do[
i=order[[j]];
x[[i]]=(b[[i]]-LPlusU[[i]].x)/Dia[[i]],
{j,Length[Dia]}] 
]
GaussSeidelIterate[xIn_,{{L_,DiaPlusU_},b_},order_]:=Module[
{},
LinearSolve[DiaPlusU,b-L.xIn]
]

All three of these should converge on the test problems.

NIter=23;
{TestMatrix,TestRHS,TestSol}={A1,b1,u1};
order=Range[1,Length[TestMatrix]];
Dia=Diagonal[TestMatrix];
LPlusU=TestMatrix-DiagonalMatrix[Dia];
L=LowerTriangularize[TestMatrix,-1];
DiaPlusU=TestMatrix-L;
x=ConstantArray[0,Length[TestMatrix]];
DataJ=Table[
x=JacobiIterate[x,{{Dia,LPlusU},TestRHS},order];
Norm[x-TestSol],
{NIter}];
x=ConstantArray[0,Length[TestMatrix]];
DataJIP=Table[
JacobiIterateInPlace[x,{{Dia,LPlusU},TestRHS},order];
Norm[x-TestSol],
{NIter}];
x=ConstantArray[0,Length[TestMatrix]];
DataJIPRandom=Table[
JacobiIterateInPlace[x,{{Dia,LPlusU},TestRHS},RandomSample[order,Length[TestMatrix]]];
Norm[x-TestSol],
{NIter}];
x=ConstantArray[0,Length[TestMatrix]];
FixedRandomOrder1=RandomSample[order,Length[TestMatrix]];
DataJIPFixedRandom1=Table[
JacobiIterateInPlace[x,{{Dia,LPlusU},TestRHS},FixedRandomOrder1];
Norm[x-TestSol],
{NIter}];
x=ConstantArray[0,Length[TestMatrix]];
FixedRandomOrder2=RandomSample[order,Length[TestMatrix]];
DataJIPFixedRandom2=Table[
JacobiIterateInPlace[x,{{Dia,LPlusU},TestRHS},FixedRandomOrder2];
Norm[x-TestSol],
{NIter}];
x=ConstantArray[0,Length[TestMatrix]];
DataGS=Table[
x=GaussSeidelIterate[x,{{L,DiaPlusU},TestRHS},order];
Norm[x-TestSol],
{NIter}];
ListLogPlot[{DataJ, DataJIP,DataGS,DataJIPRandom,DataJIPFixedRandom1,DataJIPFixedRandom2 },
PlotStyle->{Red, Blue,Green,Pink, Cyan, Black}]

ListLogPlot[{DataGS,DataJIPFixedRandom1,DataJIPFixedRandom2 },
PlotStyle->{Green,Cyan, Black}]

The order does not seem to matter too much provided it is consistent from run to run. I wonder howmuach it matters how consistent.

Block Ideas

I can implement block versions of these as well. I need to be able to block a sparse array and keeps the pieces sparse.  Unfortunately, Partition makes dense arrays but there is no real problem.

SparsePartition[A_SparseArray,{nSub_,mSub_}]:= Module[
{n,m},
{n,m}=Dimensions[A];
Table[
A[[i;;i-1+nSub,j;;j-1+mSub]],
{i, 1, n+1-nSub, nSub},
{j,1, m+1-mSub, mSub}
]
]

The thing above merely takes what fits.  I need to make my partitioning fit my problem.  Well actually I need to fix that.  This is just a test that the thing is working correctly.   It also shows that there are a lot of zero blocks.

SplitA1=SparsePartition[A1,{13,13}];
{Dimensions[A1],Dimensions[SplitA1]}
Norm[A1-ArrayFlatten[SplitA1]]
MatrixPlot[Map[Norm, SplitA1,{2}]]

To make things work in block land I need to block my right hand side and vectors.  I should also check that it is working correctly

Splitu1=Partition[u1,13];
Splitb1=Partition[b1,13];
test=Table[
Sum[SplitA1[[i,j]].Splitu1[[j]],{j,1,Length[Splitu1]}],
{i,1, Length[SplitA1]}
];
Norm[Flatten[test]-b1]

SetAttributes[JacobiIPBlock, HoldAll]
JacobiIPBlock[Splitx_,{{SplitADia_,SplitLPlusU_},Splitb_},order_]:=Module[
{i,j,n,m,RHS},
{n,m}=Take[Dimensions[SplitLPlusU],2];
Do[
i=order[[k]];
RHS= Sum[SplitLPlusU[[i,j]].Splitx[[j]],{j,1,i-1}]+ Sum[SplitLPlusU[[i,j]].Splitx[[j]],{j,i+1,m}];
RHS=Splitb[[i]]-RHS;
Splitx[[i]]=LinearSolve[SplitADia[[i]],RHS],
{k,n}
] 
]

SplitA1Dia=Diagonal[SplitA1];
SplitA1LPlusU = SplitA1;
Do[SplitA1LPlusU[[i,i]]=SparseArray[{},Dimensions[SplitA1LPlusU[[i,i]]]],
{i,1, Length[SplitA1LPlusU]}]

NIter=15;
Splitx=0*Splitu1;
order=Range[1,36];
DataJacobiIPBlock=Table[
JacobiIPBlock[Splitx,{{SplitA1Dia,SplitA1LPlusU},Splitb1},order];
Norm[Flatten[Splitx]-u1],
{NIter}];
Splitx=0*Splitu1;
DataJacobiIPBlockRO=Table[
JacobiIPBlock[Splitx,{{SplitA1Dia,SplitA1LPlusU},Splitb1},RandomSample[order,Length[order]]];
Norm[Flatten[Splitx]-u1],
{NIter}];
Clear[FixedRO, DataJacobiIPBlockFRO,p]
Do[
Splitx=0*Splitu1;
FixedRO[p]=RandomSample[order,Length[order]];
DataJacobiIPBlockFRO[p]=Table[
JacobiIPBlock[Splitx,{{SplitA1Dia,SplitA1LPlusU},Splitb1},FixedRO[p]];
Norm[Flatten[Splitx]-u1],
{NIter}],
{p,1,100}]
ListLogPlot[{DataJacobiIPBlock, DataJacobiIPBlockRO,DataJacobiIPBlockFRO[1],DataGS},
PlotStyle->{Red, Green,Pink, Blue}]

ListLogPlot[Table[DataJacobiIPBlockFRO[p],{p,1,100}],
PlotStyle->{Red, Green,Pink, Blue}]