Programok

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% inkoga %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{verbatim}
function y=inkogaus(a,b,idd1,idd2,n,maxit,epsi,S)
% Kiszamitja az inkomplett Gauss-eliminaciot, majd iteracio-
% val oldja meg az Sy=b rendszert. Itt S az a matrix, amely
% a spdiags(a,idd1,n,n) utasitassal keszithet"o. Ennek meg-
% felel"oen a idd1 megadja a foglalt atlok helyet.
% Hivasa:
%         y=inkogaus(a,b,idd1,idd2,n,maxit,epsi,S)
% ahol:
% y       a lin. rendszer k"ozelit"o megoldasa,
% a       az atlonkent tarolt matrix,
% b       a jobboldali (oszlop-)vektor,
% idd1    az atlok helyet megado (oszlop-)vektor,
% idd2    az iteracio jobb konvergencia celjabol el"oirhato,
% kiegeszit"o atlok helye, mint oszlopvektor
% n       az a matrix es a b vektor merete
% maxit   a maximalis iteracioszam
% epsi    a kilepesi feltetel abszolut pontossagi korlatja
% S       kihagyhato: S=spdiags(a,idd1,n,n)
%         Ha ez a matrix meg nincs, akkor inkoga kesziti.

[dda dda1]=size(idd1);
[idd,indd]=sort([idd1' idd2']); % a teljes indexvektor
%idd,indd,dda                   % "osszeallitasa

[dda1 dd]=size(idd);idi=0;
for i=1:dd
   if idd(i)==0            % a f"oatlo helyenek
     idi=i;                % meghatarozasa
     break
   end
end
%idi
if idi==0
  inputhiba=idi
else
  idip=idi+1;idim=idi-1;
  ddd=dd-idi;

  for j=1:idim
     k=j+1;llij=idd(j);
     for i=1:ddd
        iae=-1;iai=idd(i+idi)+llij;
        for l=k:dd         % az itt kesz"ul"o ia
           if idd(l)==iai  % indext"omb megadja azon
             k=l+1;        % atlok indexet, amelyen
             iae=1;        % valtozhatnak a matrix
             ia(j,i)=l;    % elemei (a f"oatlo alatti
             break         % atlokrol nezve)
           elseif idd(l)>iai
             break
           end
        end
        if iae==-1
          ia(j,i)=0;
        end
     end
  end
% ia                       % elkesz"ult az indext"omb

id1=1;                     % bef"uzve a kiegeszit"o
for id=1:idi               % nulla oszlopokat, olyan
   if indd(id)>dda         % alakra hozzuk a matrixot,
     aa(:,id)=zeros(n,1);  % hogy a felbontasnal keves
   else                    % i>n-fele teszt kelljen
     aa(:,id)=[zeros(-idd(id),1);a(1:n+idd(id),id1);];
     id1=id1+1;
   end
end
for id=idip:dd
   if indd(id)>dda
     aa(:,id)=zeros(n,1);
   else
     aa(:,id)=[a(1+idd(id):n,id1);zeros(idd(id),1)];
     id1=id1+1;
   end
end
clear indd
%aa                        % ezen matrix helyen
                           % t"ortenik a felbontas
zerus=0;
for i=1:n
   if aa(i,idi)==0
     zerus=i
     break
   else
      d(i)=1/aa(i,idi);
      for j=1:idim
         ij=i-idd(j);
         if ij<=n
           llij=aa(ij,j)*d(i);
           ll(ij,j)=llij;
           for k=1:ddd     % a felbontas legbels"obb
              l=ia(j,k);   % ciklusa
              if l>0
                aa(ij,l)=aa(ij,l)-llij*aa(i,idi+k);
              end
           end
         end
      end % j-ciklus vege
      for id=1:ddd
         uu(i,id)=aa(i,id+idi);
      end
   end
end

%aa,d,uu,ll                % a felbontas matrixai
clear aa ia                % elkesz"ultek

if zerus==0
  y=b;
  for i=1:n-1              % el"ore a nulladik
     for j=1:idim          % approximaciohoz
        ij=i-idd(j);
        if ij<=n
          y(ij)=y(ij)-ll(ij,j)*y(i);
        end
     end
  end
  for i=n:-1:1             % vissza
     w=y(i);
       for id=idip:dd
          is=i+idd(id);
          if is<=n
            w=w-uu(i,id-idi)*y(is);
          end
       end
     y(i)=w*d(i);
  end

  if nargin<8
    S=spdiags(a,idd1,n,n);
  end

  for it=1:maxit           % az inkomplett felbontas
                           % utani iteracio kezdete, a
     f=b-S*y;              % maradekvektor kiszamitasa
     z=norm(f,inf);
     if z<epsi
       break
     end
     for i=1:n-1           % el"ore
        for j=1:idim
           ij=i-idd(j);
           if ij<=n
             f(ij)=f(ij)-ll(ij,j)*f(i);
           end
        end
     end
     for i=n:-1:1          % vissza
        w=f(i);
        for id=idip:dd
           is=i+idd(id);
           if is<=n
             w=w-uu(i,id-idi)*f(is);
           end
        end
        f(i)=w*d(i);
     end
     y=y+f;
   end                     % iteracio vege
norm_y=z,it
end                        % nullaoszto if vege
end                        % inputhiba  if vege
\end{verbatim}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% qfigset %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{verbatim}
 function qfigset(width,height,xlabels,subplab,linew)
 %
 %ke'szitette: Kolla'r Istva'n (BME) e's Gergo' Lajos (ELTE)
 %
 %
 %a qfigset a rajzot (width x height) mereture allitja es
 %          beallitja a vonalvastagsagot linew ertekre
 %   width = az abra szelessege, alapertelmezes: az
 %           aranymetszes szerint a magassag ertekehez
 %           illesztve
 %   height = az abra magassaga, az aranymetszes szabalya
 %            szerint a szelesseg alapjan szamitva vagy
 %            3.5 inch, ha a szelesseg nincs megadva
 %   xlabels = a keret es az x-tengely kozotti teruletet
 %            szabalyozza
 %           'xtick':    megengedi, hogy az x-tengely
 %                       beosztasait megcimkezzuk
 %           'nolabels': nem enged cimkezest
 %           'plot' (alapertelmezes): megengedi az x-tengely
 %                       es a beosztas megcimkezeset is
 %           'subplot': megengedi az x-tengely es a beosztas
 %                      megcimkezeset
 %    subplab = 'felirat' ha elozoleg a 'subplot' parameter
 %                      volt megadva, akkor az x-tengely
 %                      kozepere helyezi a feliratot.
 %  A keret meretenek az alapertelmezese  3.5 x 5.66 inch.
 %
 if nargin<5, linew=1.7; end
 if nargin<3, xlabels='plot'; end
 if strcmp(xlabels,'nolabels'), fcdist=0.1;
    elseif strcmp(xlabels,'xtick'), fcdist=0.25;
    elseif strcmp(xlabels,'plot'), fcdist=0.48;
    elseif strcmp(xlabels,'subplot'), fcdist=0.70;
    else error(['xlabels = ''',xlabels,''' is not valid'])
 end
 if nargin<4, subplab=''; end
 if nargin<1, width = NaN; end
 if nargin<2, height = NaN; end
 gold=(sqrt(5)-1)/2;
 if isempty(height), height=NaN; end
 if isempty(width), width=NaN; end
 %
 %
 % ha nem adtak meg szelesseget vagy magassagot
 % (vagy egyiket sem),akkor beallitjuk az
 % aranymetszes szerint (3.5 inch magassaggal)
 %
 %
 if isnan(height) & isnan(width), height=3.5; end
 if isnan(height), height=gold*width; end
 if isnan(width), width=height/gold; end
 %
 Hp = 11;        Wp = 8.4;
 eval('tp=get(0,''TerminalProtocol'');','tp=''PCWIN'';')
 if strcmp(tp,'PCWIN')
    zunits=get(0,'Units');
    set(0,'Units', 'inches')
    scrdim=get(0,'ScreenSize');
    if ~all(finite(scrdim))
       scrdim
       figure(1), scrdim=get(0,'ScreenSize')
    end
    set(0,'Units',zunits)
 elseif ~strcmp(tp,'none') %ha csak rajz keszulhet
    zunits=get(0,'Units');
    set(0,'Units', 'inches')
    scrdim=get(0,'ScreenSize');
    set(0,'Units',zunits)
 else
    scrdim=get(0,'ScreenSize');
 end
 %
 %
 % itt szamitjuk ki a kepernyo meretet es a papir
 % meretet a papir legfeljebb 8.4x11 inch meretu
 % lehet (kb A4-es meret)
 %
 Hm=scrdim(4);  Wm=scrdim(3);
 Hp=min(Hp,Hm); Wp=min(Wm,Wp);

 TopMarg = 0.1;  BotMarg = fcdist;
 LefMarg = 0.25; RiMarg = 0.25;
 %
 % ha vannak mar grafikus ablakok a kepernyon, es
 % az aktualis ablak 'NextPlot' parametere 'new',
 % akkor uj ablakot nyitunk
 %
 hc=get(0,'Children');
 if ~isempty(hc)
    if strcmp(get(gcf,'NextPlot'),'new'),
    figure(max(hc)+1);
    end
 end
 funits=get(gcf,'Units');
 set( gcf, 'Units', 'inches' )
 set( gcf, 'PaperUnits', 'inches' )

 pwidth = width+RiMarg+LefMarg+0.5;
 pheight = height+TopMarg+BotMarg+0.35;

 % a grafikus ablak meretenek a beallitasa a
 % kepernyon, az ablakot igyekszunk a jobb
 % felso sarokba tenni
 %
 fleft = Wm - pwidth-0.05;
 fbottom = Hm - pheight-0.7;
 if fleft<0, pwidth=pwidth-fleft; fleft=0; end
 if pwidth>Wp, pwidth=Wp; end
 if pheight>Hp, pheight=Hp; end
 if fbottom<0,pheight=pheight-fbottom;fbottom=0;end

 set( gcf, 'Position',...
                  [ fleft fbottom pwidth pheight ])

 % az abra meretenek a beallitasa a papiron

 left = (Wp-pwidth);
 bottom = Hp-pheight;
 set( gcf, 'PaperPosition',...
                 [left-1.0 bottom pwidth pheight])
 set(gcf,'DefaultAxesLineWidth',0.7)
 set(gca,'LineWidth',0.7)
 set(gcf,'DefaultLineLineWidth',linew)
 set(gca,'DefaultLineLineWidth',linew)

 %
 %  a tengelyek es a keret mereteinek a beallitasa
 %
 aunits=get(gca,'Units');
 set( gca, 'Units', 'inches')
 set( gca, 'Position', [ 0.75 fcdist width height ])
 set(gcf,'DefaultAxesUnits','inches')
 set( gcf, 'Units',funits)
 set(gcf,'DefaultAxesPosition',[ 0.75/pwidth ...
     fcdist/pheight width/pwidth height/pheight ])
 set( gca, 'Units',aunits)

 figure(gcf)
 %
 %
 % a szoveg elhelyezese az x-tengely ala
 %
 %
 if strcmp(xlabels,'subplot')
    if ~isempty(subplab)
        ax=axis;
        xc=mean(ax(1:2));
        yc=ax(3)-fcdist/pheight;
        delete(findobj(gca,'Type','text','string',subplab))
        text(xc,yc,subplab,...
          'Horizontalalignment',...
            'center','Verticalalignment','top');
    end
 end
 %vege a qfigset programnak
\end{verbatim}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% gaufile %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{verbatim}
function [x,res]=gaufile(AF,AM,b,me,ep,imax)
%Nagymeretu, n ismeretlenes Ax=b egyenletrendszer megoldasa
%LU-modszerrel, iterativ javitassal
%  Parameterek:
%  AF : az A egyutthatomatrix elemeit oszlopfolytonos
%   sorrendben tartalmazo file-nev (string)
%  AM : segedfile-nev, a permutalt A matrix LU-felbontasat
%   fogja tartalmazni a vegen (string)
%   b : az egyenletrendszer jobboldala
%  me : hany egyutthato fer be egyszerre a memoriaba,
%   me>=n (integer)
%  ep : a relativ modositas tureshatara (double)
%imax : az iteraciok maximalis szama (integer)
%   x : az egyenletrendszer megoldasa
% res : a rezidualis hiba (res=Ax-b)
  n=length(b);
% A memoriaban egyszerre tarthato oszlopok szama: kc
  kc=floor(me/n);
  if kc>n, kc=n;
  elseif kc<1,
    disp('HIBA! Kotelezoen: me>=egyenletek szama'), return
  end
% Az A egyutthatomatrix particioinak szama: mc
  mc=floor(n/kc);
  if mc*kc>n, mc=mc-1;
  elseif kc*mc<n, mc=mc+1; end
  x=zeros(n,1); y=zeros(n,1);
% Az egyenletek permutaciojat mutato pointer vektor: p
  p=1:n;
  am=fopen(AM,'wb'); fclose(am);
  [af,hibajel]=fopen(AF,'rb');
  if af<0 disp('adatfile nyitasi HIBA az LU-felbontas elott!')
     uzenet=hibajel, return
  end
  l=0;
  n1=kc;
  jel=0;
  jels=0;
%Az LU-felbontas mc reszletben:
for k1=1:mc
   if k1==mc, n1=n-(mc-1)*kc; end
%Az egyutthatomatrix kovetkezo n1 db oszlopanak beolvasasa
%(n1=kc, ha nem az utolso particio, egyebkent a maradek):
  jel=jel+8*n*n1;
  A=fread(af,[n,n1],'double');
  if jel~=ftell(af)
    particio=k1, disp('adatolvasasi HIBA!')
    uzenet=ferror(af), return
  end
%%Elmaradt eliminaciok potlasa
%(szukseg van az elozo, mar kesz oszlopokra):
  [am,hibajel]=fopen(AM,'rb');
  if am<0
    particio=k1, disp('segedfile nyitasi HIBA!')
    uzenet=hibajel, return
  end
  sjel=0;
  for k=1:l
    sjel=sjel+8*n;
    col=fread(am,n,'double');
    if sjel~=ftell(am)
      particio=k1, disp('segedolvasasi HIBA!')
      uzenet=ferror(am), return
    end
    for i=k+1:n
      r=col(p(i));
      A(p(i),1:n1)=A(p(i),1:n1)-r*A(p(k),1:n1);
    end;
  end
  fclose(am);
%A memoriaban levo oszlopok tovabbi eliminacioja:
  for k=1:n1
    [s,t]=max(abs(A(p(l+k:n),k)));
    if abs(s)<1e-50
      disp('Figyelem! Az egyutthatomatrix szingularis,...
        vagy majdnem szingularis!')
    t=l+t+k-1; if k<t, p([l+k,t])=p([t,l+k]); end;
    for i=l+k+1:n
      r=A(p(i),k)/A(p(l+k),k); A(p(i),k)=r;
      if k<n1
    A(p(i),k+1:n1)=A(p(i),k+1:n1)-r*A(p(l+k),k+1:n1);
      end;
    end;
  end;
%Elkeszult oszlopok felvitele az AM file-ba:
  [am,hibajel]=fopen(AM,'ab');
  if am<0
    particio=k1, disp('segedfile nyitasi HIBA!')
    uzenet=hibajel, return
  end
  jels=jels+8*n*n1;
  as=fwrite(am,A,'double');
    if jels~=ftell(am)|as~=n*n1
      particio=k1, disp('segedirasi hiba HIBA!')
      uzenet=ferror(am), return
    end
  fclose(am);
  l=l+kc;
end
fclose(af); clear A;
%Befejezodott a permutalt A matrix LU-felbontasa.
%Mind az L, mind az U az AM file-ban van.
%A p mutatja a sorok permutaciojat, b1 jelolje a permutalt b-t.
%Az LUx=b1 megoldasa az eredeti Ax=b megoldasaval egyezik.
%Oldjuk meg elobb az Ly=b1 egyenletet, majd az Ux=y egyenletet!
  y=solvtril(AM,b,p,me);
  x=solvtriu(AM,y,p,me);
%Az x javitasa iteracioval mindaddig amig a relativ modositas
%nagyobb epszilon-nal, de legfeljebb imax lepesben:
  for k=1:imax
%Eloszor az res=Ax-b rezidualis hibat szamitjuk:
    res=-b;
    n1=kc;
    l=0;
    [af,hibajel]=fopen(AF,'rb');
    if af<0 disp('adatfile nyitasi HIBA az iteracio elott!')
      uzenet=hibajel, return
    end
    jel=0;
    for k1=1:mc
      if k1==mc, n1=n-(mc-1)*kc; end
      jel=jel+8*n*n1;
      A=fread(af,[n,n1],'double');
      if jel~=ftell(af)
    disp('adatolvasasi hiba HIBA az iteracioban!')
    particio=k1, uzenet=ferror(af), return
      end
      res=res+A*x(l+1:l+n1);
      l=l+kc;
    end
    fclose(af); clear A;
%Most az Az=res egyenletrendszert oldjuk meg,
%itt is felhasznalva a permutalt A matrix LU felbontasat:
    y=solvtril(AM,res,p,me);
    z=solvtriu(AM,y,p,me);
    d=norm(z)/norm(x);
%A kozelito x modositasa:
    x=x-z;
    if d<=ep return, end
  end
return
\end{verbatim}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% solvtril %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{verbatim}
function [x]=solvtril(AM,b,p,me)
%Lx=b1 megoldasa, ahol L az AM-ben oszlopfolytonosan tarolt
%matrix sorok szerinti permutaltjanak also haromszog resze,
%a foatloban 1-ek feltetelezve, tovabba
%b1 a b parameter megfelelo permutacioja.
%Parameterek:
%  AM : file-nev (string)
%   b : az Lx=b1 egyenletrendszer
%   jobboldalanak inverz permutaltja
%   p : a permutacio mutatoja (integer vektor)
%  me : az U hany eleme lehet egyszerre a memoriaban
  n=length(b);
  x=zeros(n,1);
% A memoriaban egyszerre tarthato oszlopok szama: kc
  kc=floor(me/n);
  if kc>n, kc=n;
  elseif kc<1,
    disp('HIBA! Kotelezoen: me>=egyenletek szama'), return
  end
% Az L egyutthatomatrix particioinak szama: mc
  mc=floor(n/kc);
  if mc*kc>n, mc=mc-1;
  elseif kc*mc<n, mc=mc+1; end
  l=0;
  n1=kc;
  [am,hibajel]=fopen(AM,'rb');
  if am<0 disp('file nyitasi HIBA az also haromszogmatrixban!')
     uzenet=hibajel, return
  end
  jel=0;
%Az egyenletrendszer megoldasa mc reszletben:
  for k1=1:mc
    if k1==mc
      n1=n-(mc-1)*kc;
    end
    jel=jel+8*n*n1;
%Az L egyutthatomatrix kovetkezo n1 db oszlopanak beolvasasa
%(n1=kc, ha nem az utolso particio, egyebkent a maradek):
    L=fread(am,[n,n1],'double'); %AM(1:n,l+1:l+n1);
    if jel~=ftell(am)
      particio=k1,
      disp('adatolvasasi HIBA az also haromszogmatrixban!')
      uzenet=ferror(am), return
    end
    for i=1:n1
      x(l+i)=b(p(l+i));
      b=b-x(l+i)*L(1:n,i);
   end
   clear L;
    l=l+kc;
  end
  fclose(am);
return
\end{verbatim}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% solvtriu %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{verbatim}
function [x]=solvtriu(AM,b,p,me)
%Ux=b megoldasa, ahol U az AM-ben oszlopfolytonosan tarolt
%matrix sorok szerinti permutaltjanak felso haromszog resze
%Parameterek:
%  AM : file-nev (string)
%   b : az egyenletrendszer jobboldala
%   p : a permutacio mutatoja (integer vektor)
%  me : az U hany eleme lehet egyszerre a memoriaban (integer)
  n=length(b);
  x=zeros(n,1);
% A memoriaban egyszerre tarthato oszlopok szama: kc
  kc=floor(me/n);
  if kc>n, kc=n;
  elseif kc<1,
    disp('HIBA! Kotelezoen: me>=egyenletek szama'), return
  end
% Az U egyutthatomatrix particioinak szama: mc
  mc=floor(n/kc);
  if mc*kc>n, mc=mc-1;
  elseif kc*mc<n, mc=mc+1; end
  l=n;
  n1=kc;
  [am,hibajel]=fopen(AM,'rb');
  if am<0 disp('file nyitasi HIBA a felso haromszogmatrixban!')
     uzenet=hibajel, return
  end
  jel=8*n*n;
%Az egyenletrendszer megoldasa mc reszletben:
  for k1=1:mc
    if k1==mc
      n1=n-(mc-1)*kc;
    end
      fseek(am,8*n*(l-n1),-1);
%Az U egyutthatomatrix kovetkezo n1 db oszlopanak beolvasasa
%(n1=kc, ha nem az utolso particio, egyebkent a maradek):
      U=fread(am,[n,n1],'double'); %AM(1:n,l-n1+1:l);
      if jel~=ftell(am)
    particio=k1,
    disp('adatolvasasi HIBA a felso haromszog matrixban!')
    uzenet=ferror(am), return
      end
    for i=1:n1
      s=U(p(l-i+1),n1-i+1)
      if abs(s)<1e-50
    disp('Figyelem! Az egyutthatomatrix szingularis,...
          vagy majdnem szingularis!')
      x(l-i+1)=b(l-i+1)/s;
      b=b-x(l-i+1)*U(p,n1-i+1);
    end
    clear U;
    l=l-kc; jel=jel-8*n*n1;
  end
  fclose(am);
return
\end{verbatim}