Unit Matrices; {$A+,B-,D+,E+,F-,I+,L+,N+,O+,R+,S+,V+} {$M 65520,0,655360} Interface Uses Definitions, Sorties, { pour so_rt } Environnement; { pour progression } Const dimMat = npmax; { dimMat … aller piocher dans D‚finitions... } { ne peut atteindre 70 en matrices statiques... } Type Vmat = Array[1..dimMat] Of DataMat; colMatInt = Array[1..dimMat] Of Integer; { possibilit‚ Byte } matrice = Array[1..dimMat] Of Vmat; Procedure matmult ( a, b : matrice; n, m, p : integer; var c : matrice ); Procedure mattransp ( a : matrice; n, m : integer; var b : matrice ); procedure matread(var a:matrice;n,m:integer); procedure matwrite(var a:matrice;n,m:integer); (*procedure rechpiv(var a:matrice;n,k:integer;var i,j:integer;var pivot:DataMat);*) procedure solvetriang(var a,b:matrice;n,m:integer;var x:matrice); procedure mattriang(var a,b:matrice;n,m:integer;var x:matrice;var det:Data); procedure matdet(a:matrice;n:integer;var det:Data); procedure systlin(a:matrice;b:vmat;n:integer;var x:vmat;var det:Data); procedure matinv(a:matrice;n:integer;var ainv:matrice;var det:Data); procedure diagsym(a:matrice;n,niteration:integer;precision:Data; var valpropre:vmat;var vectpropre:matrice;var precatt:Data); Implementation { Var poubMat1, poubMat2 : matrice; solution trop gourmande } procedure matmult; { (a,b:matrice;n,m,p:integer;var c:matrice); } var i,j,k:integer; { trans:DataMat; utilit‚ ??? } begin for i:=1 to n do for j:=1 to p do begin c[i,j]:=0; for k:=1 to m do c[i,j]:=c[i,j]+a[i,k]*b[k,j] end; end; procedure mattransp; { (a:matrice;n,m:integer;var b:matrice); } var i,j:integer; begin for i:=1 to m do for j:=1 to n do b[i,j]:=a[j,i] end; procedure matread; { (var a:matrice;n,m:integer); } var i,j:integer; begin for j:=1 to m do begin writeln('Entrer la colonne ',j,' : '); for i:=1 to n do begin write(' terme ',i,',',j,' = '); readln(a[i,j]) end; writeln end; end; procedure matwrite; { (var a:matrice;n,m:integer); } var i,j:integer; begin for i:=1 to n do begin for j:=1 to m do write(so_rt,a[i,j]:10:4); writeln(so_rt) end; end; procedure rechpiv ( var a:matrice;n,k:integer;var i,j:integer;var pivot:DataMat); var l,m:integer; begin i:=k;j:=k; pivot:=a[k,k]; for l:=k to n do for m:=k to n do begin if abs(a[l,m])>abs(pivot) then begin pivot:=a[l,m]; i:=l; j:=m end; end; end; procedure echligne(var a:matrice;n,i,j,k:integer); var l:integer; trans:DataMat; begin for l:=k to n do begin trans:=a[i,l]; a[i,l]:=a[j,l]; a[j,l]:=trans end; end; procedure echcol(var a:matrice;n,i,j:integer); var l:integer; trans:DataMat; begin for l:=1 to n do begin trans:=a[l,j]; a[l,j]:=a[l,i]; a[l,i]:=trans end; end; procedure solvetriang; { (var a,b:matrice;n,m:integer;var x:matrice); } var i,j,k:integer; trans:Extended; begin for i:=1 to m do begin x[n,i]:=b[n,i]/a[n,n]; for j:=1 to n-1 do begin trans:=0; for k:=n-j+1 to n do trans:=trans+a[n-j,k]*x[k,i]; x[n-j,i]:=(b[n-j,i]-trans)/a[n-j,n-j] end; end; end; procedure mattriang; { (var a,b:matrice;n,m:integer;var x:matrice;var det:Data); } var trans,pivot:DataMat; memcol:colMatInt; i,j,k,ik,jk:integer; det2 : Extended; begin det2:=1.0; for k:=1 to n-1 do begin rechpiv(a,n,k,ik,jk,pivot); if pivot=0 then begin writeln('Matrice non-inversible'); writeln('D‚terminant nul'); exit end; memcol[k]:=jk; det2:=det2*pivot; if ik<>k then begin det2:=-det2; echligne(a,n,ik,k,k); echligne(b,m,ik,k,1) end; if jk<>k then begin det2:=-det2; echcol(a,n,jk,k) end; for i:=k+1 to n do begin trans:=a[i,k]/a[k,k]; for j:=k to n do a[i,j]:=a[i,j]-trans*a[k,j]; for j:=1 to m do b[i,j]:=b[i,j]-trans*b[k,j]; end; end; det2:=det2*a[n,n]; if det2=0 then begin writeln('Matrice non-inversible'); writeln('d‚terminant nul'); exit end; solvetriang(a,b,n,m,x); for j:=n-1 downto 1 do begin k:=memcol[j]; if k<>j then begin for i:=1 to m do begin trans:=x[j,i]; x[j,i]:=x[k,i]; x[k,i]:=trans end; end; end; If det2<1.0E+30 Then det:=det2 Else det:=1.0E+30; { trŠs gros... } end; procedure matdet; { (a:matrice;n:integer;var det:Data); } var i:integer; x,b:matrice; begin for i:=1 to n do b[i,1]:=0; { b[] } mattriang(a,b,n,1,x,det) { mattriang(a,b,n,1,x,det) } end; procedure systlin; { (a:matrice;b:vmat;n:integer;var x:vmat;var det:Data); } var c,y:matrice; i:integer; begin for i:=1 to n do c[i,1]:=b; { c devient poubMat1 et y devient poubMat2 } mattriang(a,c,n,1,y,det); for i:=1 to n do x:=y[i,1] end; procedure matinv; { (a:matrice;n:integer;var ainv:matrice;var det:Data); } var b:matrice; i,j:integer; begin for i:=1 to n do { poubMat1 remplace b } for j:=1 to n do if i=j then b[i,j]:=1 else b[i,j]:=0; mattriang(a,b,n,n,ainv,det) end; procedure pivotsym(var a:matrice;n:integer;var lpivot,cpivot:integer); var i,j:integer; pivot:DataMat; begin pivot:=0; for i:=1 to n-1 do for j:=i+1 to n do if abs(a[i,j])>pivot then begin pivot:=abs(a[i,j]); lpivot:=i; cpivot:=j end; end; procedure initdiag(var a,vectpropre:matrice;n:integer;var norme,tracecarre:Data); var i,j:integer; begin tracecarre:=0; for i:=1 to n do begin vectpropre[i,i]:=1; tracecarre:=tracecarre+sqr(a[i,i]) end; norme:=0; for i:=1 to n-1 do for j:=i+1 to n do begin vectpropre[i,j]:=0; vectpropre[j,i]:=0; norme:=norme+sqr(a[i,j]) end; norme:=norme*2+tracecarre end; procedure rotation(var a:matrice;n,lpivot,cpivot:integer;var uelem:matrice); var i,j:integer; p,q,b:DataMat; c,s:DataMat; begin b:=a[lpivot,lpivot]-a[cpivot,cpivot]; if b=0 then begin c:=1/sqrt(2); s:=a[lpivot,cpivot]/sqrt(2)/abs(a[lpivot,cpivot]) end else begin q:=abs(b); p:=2*b*a[lpivot,cpivot]/q; b:=sqrt(p*p+q*q); c:=sqrt((1+q/b)/2); s:=p/(2*b*c) end; for i:=1 to n do uelem[i,i]:=1; for i:=1 to n-1 do for j:=i+1 to n do begin uelem[i,j]:=0; uelem[j,i]:=0; end; uelem[lpivot,lpivot]:=c; uelem[cpivot,cpivot]:=c; uelem[lpivot,cpivot]:=-s; uelem[cpivot,lpivot]:=s end; procedure diagsym; { (a:matrice;n,niteration:integer;precision:Data; var valpropre:vmat;var vectpropre:matrice;var precatt:Data); } var iteration,lpivot,cpivot,i:integer; norme,tracecarre:Data; rotelem:matrice; begin iteration:=0; precatt:=1.0; { pour Progression... } initdiag(a,vectpropre,n,norme,tracecarre); { poubMat1 remplace rotelem } repeat Progression(Trunc(100*precision/precatt),1); iteration:=iteration+1; pivotsym(a,n,lpivot,cpivot); tracecarre:=tracecarre+2*sqr(a[lpivot,cpivot]); rotation(a,n,lpivot,cpivot,rotelem); matmult(vectpropre,rotelem,n,n,n,vectpropre); matmult(a,rotelem,n,n,n,a); mattransp(rotelem,n,n,rotelem); matmult(rotelem,a,n,n,n,a); precatt:=abs(1-(norme/tracecarre)) until (precatt<=precision) or (iteration=niteration); for i:=1 to n do valpropre:=a[i,i]; end; Begin { Init Part } WriteLn('Matrices Statiques'); End. { Unit Matrices }