Re: High-speed sparse matrix operations with Mathematica 4.2 ???
- To: mathgroup at smc.vnet.net
- Subject: [mg42321] Re: High-speed sparse matrix operations with Mathematica 4.2 ???
- From: carlos at colorado.edu (Carlos Felippa)
- Date: Tue, 1 Jul 2003 08:45:32 -0400 (EDT)
- References: <bdh7oa$o6s$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
heimann at ism.tu-berlin.de (Justus Heimann) wrote in message news:<bdh7oa$o6s$1 at smc.vnet.net>... > Hi, > > I have a speed problem running matrix operations like Dot[...] and > LinearSolve[...] on large sparse numerical matrices. Currently I'm > using Linux-Mathematica 4.2. > > I noticed that the new Mathematica 5 promises to be much faster on > this kind of problems than version 4.2 > (http://www.wolfram.com/products/mathematica/newin5/performance/sparselinear.html). > But since I don't have the chance to upgrade to Rel. 5 yet, I wonder > if there still is a smart way to speed up large sparse matrix > operations under Mathematica 4.2 ?! Can anybody help ??? > > Many thanks > Justus Personally I could never find a way. For example, here is a symmetric skyline factor translated from Fortran into Mathematica for a course on finite elements. It uses Dot and Take in the inner loop. SymmSkyMatrixFactor[S_,tol_]:= Module[ {p,a,fail,i,j,k,l,m,n,ii,ij,jj,jk,jmj,d,s,row,v}, row=SymmSkyMatrixRowLengths[S]; s=Max[row]; {p,a}=S; n=Length[p]-1; v=Table[0,{n}]; fail=0; Do [jj=p[[j+1]]; If [jj<0|row[[j]]==0, Continue[]]; d=a[[jj]]; jmj=Abs[p[[j]]]; jk=jj-jmj; Do [i=j-jk+k; v[[k]]=0; ii=p[[i+1]]; If [ii<0, Continue[]]; m=Min[ii-Abs[p[[i]]],k]-1; ij=jmj+k; v[[k]]=a[[ij]]; v[[k]]-=Take[a,{ii-m,ii-1}].Take[v,{k-m,k-1}]; a[[ij]]=v[[k]]*a[[ii]], {k,1,jk-1}]; d-=Take[a,{jmj+1,jmj+jk-1}].Take[v,{1,jk-1}]; If [Abs[d]<tol*row[[j]], fail=j; a[[jj]]=Infinity; Break[] ]; a[[jj]]=1/d, {j,1,n}]; Return[{{p,a},fail}] ]; SymmSkyMatrixRowLengths[S_]:= Module[ {p,a,i,j,n,ii,jj,m,d,row}, {p,a}=S; n=Length[p]-1; row=Table[0,{n}]; Do [ii=p[[i+1]]; If [ii<0, Continue[]]; m=ii-i; row[[i]]=a[[ii]]^2; Do [If [p[[j+1]]>0, d=a[[m+j]]^2; row[[i]]+=d; row[[j]]+=d], {j,Max[1,Abs[p[[i]]]-m+1],Min[n,i]-1}], {i,1,n}]; Return[Sqrt[row]]; ]; The penalty from Fortran was roughly 10^4 under Mathematica 4.0-4.2. Example: order 1000, mean skyline width 60: f77 ~ 0.03 sec, Mathematica ~ 5 min. SymmSkyMatrixFactor in fact runs slower than LinearSolve on a full matrix (for order 1000, LinearSolve takes about 20 sec). So I never bothered to use it in the course. Mathematica 5.0 looks more competitive.