From owner-octave-sources at bevo dot che dot wisc dot edu Thu Aug 22 13:55:24 1996 Subject: Band and tri-diagonal linear systems solver From: Mario Storti To: octave-sources at bevo dot che dot wisc dot edu Date: Thu, 22 Aug 1996 15:52:25 -0300 Hi all! I had to develop this for solving systems arising from Finite Element discretizations. They are comparatively much slower than the built-in full-matrix versions, since they are, at most, "vectorized" at the scale of the band-width. Then, they are useful only for very large systems, specially if the full version doesn't fit in memory. The test.m script generates an example with a random generated matrix and solves it with the full and banded version, comparing the elapsed times. If someone knows of something better, drop me a line, please... :-) Mario --------Cut here and run with /bin/sh ----------------------------- #!/bin/sh # shar: Shell Archiver (v1.22) # # Run the following text with /bin/sh to create: # bandfac.m # bandsol.m # test.m # trifac.m # trisol.m # sed 's/^X//' << 'SHAR_EOF' > bandfac.m && Xfunction afac=bandfac(a); X% afac=BANDFAC(aband); X% Band matrix factorization. X% aband(n,2*nd+1)=non-null co-diagonals of A X% nd=number of non-null co-diagonals X% afac(n,nd)=diagonals of factorized A (Gauss elimination) X% Once factorized the system is solved by calling BANDSOL X% Band storage must be performed by shifting rows. (Complete description X% in the BANDFAC.M file.) X% X% see also: TRIFAC, TRISOL, BANDSOL X X%$Id: bandfac.m,v 1.3 1996/08/22 18:33:16 mstorti Exp $ X X% Bugs or comments to Mario Storti X X% Let: X% X% A x = b X% X% be your banded linear system A in R^(n x n) and A_ij=0 for abs(i-j)>nd X% then you should store the non-zero coefficients of A in a rectangular matrix of X% Aband size [n,2*nd+1] in the following way: X% X% Aband = [A(1,-nd:1+nd); X% A(2,1+(-nd:1+nd)); X% . X% . X% . X% A(n-1,n-2+(-nd:1+nd)); X% A(n,n-1+(-nd:1+nd));] X% X% To fix ideas, if n=5, nd=2 and: X% X% A= [1 2 4 0 0; Aband= [0 0 1 2 4; X% 3 2 1 5 0; 0 3 2 1 5; X% 2 4 3 1 2; then 2 4 3 1 2; X% 0 1 3 2 3; 1 3 2 3 0; X% 0 0 3 2 1]; 3 2 1 0 0]; X% X% To solve the linear system, you should first factorize the matrix wiht X% bandfac and then solve the system with bandsol: X% X% Afac=bandfac(Aband); X% x=bandsol(Afac,b); X% X X[n,nd]=size(a); X Xncod=round((nd-1)/2); Xif 2*ncod+1~=nd X% disp("bandfac: nd=column dimension of a not odd") Xend X Xafac=a; Xupdiag=ncod+(1:ncod+1)'; Xafac(1,updiag)=a(1,updiag); Xfor j=1:n-1 X for k=1:min([ncod n-j]) X pivote=afac(j+k,ncod+1-k)/afac(j,ncod+1); X afac(j+k-1,ncod+1-k)=pivote; X afac(j+k,updiag-k)=afac(j+k,updiag-k)-pivote*afac(j,updiag); X end Xend X SHAR_EOF chmod 0444 bandfac.m || echo "restore of bandfac.m fails" sed 's/^X//' << 'SHAR_EOF' > bandsol.m && Xfunction x=bandsol(ap,b) X% x=BANDSOL(ap,b) X% Solution of a tridiagonal system previously factorized by BANDFAC. X% ap(n,nd)= codiagonals of the factorized matrix, as obtained from BANDFAC X% nd=number of diagonals (must be an odd number) X% b(n)= right hand side vector X% x(n)= unknown X% X% see also: TRIFAC, TRISOL, BANDFAC X X% Bugs or comments to Mario Storti X% $Id: bandsol.m,v 1.2 1996/08/22 18:33:33 mstorti Exp $ X X[n,nd]=size(ap); X Xncod=round((nd-1)/2); Xif 2*ncod+1~=nd X disp("bandsol: nd=column dimension of a not odd") Xend X Xx=b; Xx(1,:)=b(1,:); Xfor j=1:n-1 X for k=1:min([ncod n-j]) X x(j+k,:)=x(j+k,:)-ap(j+k-1,ncod-k+1)*x(j,:); X end Xend X Xx(n,:)=x(n,:)/ap(n,ncod+1); Xfor j=n-1:-1:1 X ncodj=min([ncod n-j]); X x(j,:)=(x(j,:)-ap(j,ncod+1+(1:ncodj))*x(j+(1:ncodj)',:))/ap(j,ncod+1); Xend Xendfunction SHAR_EOF chmod 0444 bandsol.m || echo "restore of bandsol.m fails" sed 's/^X//' << 'SHAR_EOF' > test.m && X%$Id: test.m,v 1.1 1996/08/21 20:15:45 mstorti Exp mstorti $ X Xn=200; Xm=15; XA=rand(n)+5*eye(n); XAband=zeros(n,2*m+1); X Xfor k=1:n X for j=1:n X if abs(j-k)>m X A(j,k)=0; X else X if k-j+m+1>2*m+1 X disp("que paso?") X break,break X endif X Aband(j,k-j+m+1)=A(j,k); X endif X endfor Xendfor X Xb=rand(n,1); X Xtic;x=A\b; Xdisp("Full system solution (secs):"); toc X Xtic; XAfac=bandfac(Aband); Xxband=bandsol(Afac,b); Xdisp("Banded solution (secs):"); toc X Xdisp("Difference between band solution and full matrix solution:") Xmax(abs(xband-x)) SHAR_EOF chmod 0644 test.m || echo "restore of test.m fails" sed 's/^X//' << 'SHAR_EOF' > trifac.m && Xfunction ap=trifac(a); X% ap=TRIFAC(a); X% Tridiagonal factorization. X% a(n,3)=diagonals of A X% ap(n,3)=diagonals of factorized A X% Band storage must be performed by shifting rows. (Complete description X% in the BANDFAC.M file.) X% X% see also: TRISOL, BANDFAC, BANDSOL X X% Bugs or comments to Mario Storti X%$Id: trifac.m,v 1.2 1996/08/22 18:31:04 mstorti Exp $ X X[n,bid]=size(a); X Xap=a; Xap(1,2)=a(1,2); Xfor j=1:n-1 X j1=j+1; X ap(j,1)=a(j1,1)/ap(j,2); X ap(j1,2)=ap(j1,2)-ap(j,1)*ap(j,3); Xend SHAR_EOF chmod 0444 trifac.m || echo "restore of trifac.m fails" sed 's/^X//' << 'SHAR_EOF' > trisol.m && Xfunction x=trisol(ap,b) X% x=TRISOL(ap,b) X% Solution of a tridiagonal system previously factorized by TRIFAC. X% ap(n,3)= codiagonals of the factorized matrix, as obtained from TRIFAC X% b(n)= right hand side vector X% x(n)= unknown X% Band storage must be performed by shifting rows. (Complete description X% in the BANDFAC.M file.) X% X% see also: TRIFAC, BANDSOL, BANDFAC X X% Bugs or comments to Mario Storti X%$Id: trisol.m,v 1.2 1996/08/22 18:31:12 mstorti Exp $ X X X[n,bid]=size(ap); X Xx=b; Xx(1,:)=b(1,:); Xfor j=1:n-1 X x(j+1,:)=b(j+1,:)-ap(j,1)*x(j,:); Xend X Xx(n,:)=x(n,:)/ap(n,2); Xfor j=n-1:-1:1 X x(j,:)=(x(j,:)-ap(j,3)*x(j+1,:))/ap(j,2); Xend X SHAR_EOF chmod 0444 trisol.m || echo "restore of trisol.m fails" exit 0 -- %%%%%%<>%%%%%%<>%%%%%%<>%%%%%%<>%%%%%%<>%%%%%%<>%%%%%%<>%%%%%%<>% Mario Alberto Storti | Fax: (54)(42) 55.09.44 | Grupo de Tecnologia Mecanica | Tel: (54)(42) 55.91.75 | INTEC, Guemes 3450 - 3000 Santa Fe | http://venus.unl.edu.ar/ | Argentina | Home: Gob. Vera 3161 | Reply: mstorti at galileo dot unl dot edu dot ar| | (54)(42) 55 dot 00 dot 23 | (54)(42) 55.00.23 |