More on: Newbie question: big matrix calculations

*To*: mathgroup at smc.vnet.net*Subject*: [mg9298] More on: [mg9253] Newbie question: big matrix calculations*From*: Luci Ellis <elisha at dot.net.au>*Date*: Mon, 27 Oct 1997 02:47:08 -0500*Sender*: owner-wri-mathgroup at wolfram.com

Dennis I have looked at your code and have realised there is a lot quicker way to do it. See below: Your original code: In[1]:= size =5; In[2]:= matrxD=Table[Switch[i-j,0,5,_,0],{i,size},{j,size}]; In[3]:= matrxU=Table[Switch[i-j,-1,-3,_,0],{i,size},{j,size}]; In[4]:= matrxL=Table[Switch[i-j,1,-2,_,0],{i,size},{j,size}]; In[5]:= b=Table[1,{i,size}]; In[6]:= x1=Table[1,{i,size}]; In[9]:= x2 = -Inverse[matrxD].(matrxL+matrxU).x1+Inverse[matrxD].b;//Timing Out[9]= {0.716667 Second,Null} My original suggestion about not inverting D twice doesn't buy you much for small sizes, if anything you'll go backwards (-; In[11]:= x2a = With[{dinv=Inverse[matrxD]},-dinv.(matrxL+matrxU).x1+dinv.b];//Timing Out[11]= {0.95 Second,Null} Let's look at what's going on here: In[10]:= x2 Out[10]= {4/5, 6/5, 6/5, 6/5, 3/5} This is the first issue -- you are working in integers, with infinite precision, but you are implementing what is in essence a numerical procedure. So you let's try it in Real numbers not rationals. In[19]:= size = 5; In[22]:= matrxD=Table[Switch[i-j,0,5.,_,0.],{i,size},{j,size}]; In[23]:= matrxU=Table[Switch[i-j,-1,-3.,_,0.],{i,size},{j,size}]; In[24]:= matrxL=Table[Switch[i-j,1.,-2.,_,0.],{i,size},{j,size}]; In[16]:= b=Table[1.,{i,size}]; In[17]:= x1=Table[1.,{i,size}]; In[25]:= x2 = -Inverse[matrxD].(matrxL+matrxU).x1+Inverse[matrxD].b;//Timing Out[25]= {0.05 Second,Null} Much quicker! (But I got a different answer, I don't know why...) In[26]:= x2 Out[26]= {0.8,0.8,0.8,0.8,0.2} But there's more to it than that. Looking at your code I see that what you are trying to do is solve for the matrix X that satisfies the equation: X = -inv(D).(L+U).x +inv(D).b or if you like: inv(D)-inv(D).(L+U).x +inv(D).b-X=0. So why don't you Solve for the equation explicitly instead of using repeated evaluation to find the solution In[32]:= X= Table[x[i],{i,size}] Out[32]= {x[1],x[2],x[3],x[4],x[5]} In[36]:= eqns=Thread[ -Inverse[matrxD].(matrxL+matrxU).X+Inverse[matrxD].b-X==0] Out[36]= {0.2 -1. x[1]+0.6 x[2]+0. x[3]+0. x[4]+0. x[5]==0, 0.2 +0. x[1]-1. x[2]+0.6 x[3]+0. x[4]+0. x[5]==0, 0.2 +0. x[1]+0. x[2]-1. x[3]+0.6 x[4]+0. x[5]==0, 0.2 +0. x[1]+0. x[2]+0. x[3]-1. x[4]+0.6 x[5]==0, 0.2 +0. x[1]+0. x[2]+0. x[3]+0. x[4]-1. x[5]==0} In[40]:= solution=Solve[eqns,X]//Timing Out[40]= {0.0666667 Second,{{x[1]->0.46112,x[2]->0.4352,x[3]->0.392, x[4]->0.32,x[5]->0.2}}} Pretty fast considering I'm doing this on a three-year-old Mac Quadra 630! So you could write a function: In[45]:= solvesparsething[size_Integer]:= Module[{matrxD,matrxU,matrxL,b,X}, b=Table[1.,{i,size}]; matrxD=Table[Switch[i-j,0,5.,_,0.],{i,size},{j,size}]; matrxU=Table[Switch[i-j,-1,-3.,_,0.],{i,size},{j,size}]; matrxL=Table[Switch[i-j,1.,-2.,_,0.],{i,size},{j,size}]; X= Table[x[i],{i,size}]; eqns=Thread[ -Inverse[matrxD].(matrxL+matrxU).X+Inverse[matrxD].b-X==0]; Flatten[X/.Solve[eqns,X]]] In[46]:= solvesparsething[6]//Timing Out[46]= {0.766667 Second,{0.476672,0.46112,0.4352,0.392,0.32,0.2}} In[48]:= solvesparsething[10]//Timing Out[48]= {1. Second,{0.496977,0.494961,0.491602,0.486003,0.476672,0.46112,0.4352,0.392, 0.32,0.2}} In[49]:= solvesparsething[50]//Timing Out[49]= {31.2833 Second,{0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5, 0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.499999,0.499999,0.499999,0.499998, 0.499996,0.499993,0.499989,0.499982,0.49997,0.499949,0.499915,0.499859, 0.499765,0.499608,0.499347,0.498912,0.498186,0.496977,0.494961,0.491602, 0.486003,0.476672,0.46112,0.4352,0.392,0.32,0.2}} In[50]:= solvesparsething[100]//Timing Hope this helps, Regards, Luci -------------- Luci Ellis: elisha at dot.net.au http://www.dot.net.au/~elisha