|
% ***************************************************************************
% J. Dankert: Gauss-Algorithmus mit Spaltenpivotisierung
% fuer unsymmetrische Bandmatrix
%
% Vorsicht: Dies ist ein Algorithmus zum Gehirn-Verrenken! Der Autor
% warnt ausdruecklich vor dem Versuch, diesen Algorithmus
% nachempfinden zu wollen.
%
% Es wird das lineare Gleichungssystem A*X=B nach dem Gaussschen
% Algorithmus mit Spaltenpivotisierung gelöst. Dabei ist die N*N-Matrix A
% als unsymmetrische Bandmatrix entweder als kompakte Rechteckmatrix oder
% kompakt als Vektor gespeichert.
%
% BANDMATRIZEN haben nur innerhalb eines Bandes in der Naehe der
% Hauptdiagonalen von Null verschiedene Elemente. Die "Bandweiten" IBWL
% und IBWR charakterisieren die Breite dieses Bandes links bzw. rechts
% von der Hauptdiagonalen, wobei die "weiteste Entfernung" eines
% Elementes von der Hauptdiagonalen entscheidend ist und die
% Hauptdiagonalelement-Position jeweils mitgezaehlt wird.
%
% Beispiel: Wenn das Element A(7,10) einer Matrix am weitesten nach rechts
% von der Hauptdiagonalen entfernt ist, waere IBWR = 4. Dieser Wert
% gilt dann fuer die gesamte Matrix, auch wenn in anderen Zeilen
% die von Null verschiedenen Elemente naeher an der
% Hauptiagonalen liegen.
%
% In der folgenden Matrix ist die 3 in der zweiten Zeile verantwortlich
% fuer die rechte Bandweite IBWR = 4, die 8 in der vierten Zeile fuer die
% linke Bandweite IBWL = 3 (Diagonalelement jeweils mitgezaehlt). Damit ist
% die Gesamtbandweite IBW = IBWL + IBWR -1 = 6 die Anzahl der Elemente,
% die fuer jede Zeile gespeichert werden. Die ersten Zeilen werden durch
% fuehrende Nullen, die letzten Zeilen durch nachgestellte Nullen
% auf diese Anzahl aufgefuellt, wie die Darstellung rechts zeigt (bei einer
% kleinen Matrix wie in diesem Beispiel ergibt sich so natuerlich keine
% Einsparung an Speicherplatz).
% __ __
% | 4 2 0 0 0 0 | 0 0 4 2 0 0
% | 1 0 2 0 3 0 | 0 1 0 2 0 3
% | 0 0 7 1 4 0 | 0 0 7 1 4 0
% | 0 8 2 0 9 9 | -----> 8 2 0 9 9 0
% | 0 0 0 1 1 2 | 0 1 1 2 0 0
% | 0 0 0 0 4 2 | 0 4 2 0 0 0
% --- ---
%
% Die rechts dargestellten Elemente werden entweder als Rechteckmatrix
% (Speichervariante a) oder als Vektor (Speichervariante b) gespeichert.
%
% Speichervariante a
% ------------------
% Es wird folgende Rechteckmatrix gespeichert (N Zeilen, IBW Spalten, die
% Hauptdiagonalelemente stehen hier in der Spalte IBWL=hdsp=3):
% __ __
% | 0 0 4 2 0 0 |
% | 0 1 0 2 0 3 |
% | 0 0 7 1 4 0 |
% | 8 2 0 9 9 0 |
% | 0 1 1 2 0 0 |
% | 0 4 2 0 0 0 |
% --- ---
%
% Speichervariante b:
% -------------------
% Die Matrixelemente werden in folgender Reihenfolge im Vektor A
% gespeichert (N*IBW Elemente):
%
% 0 0 4 2 0 0 0 1 0 2 0 3 0 0 7 1 4 0 8 2 0 9 9 0 0 1 1 2 0 0 0 4 2 0 0 0
%
% Das Element A(i,j) der Matrix landet im Vektor A auf der Position
% (i-1)*IBW+IBWL+j-i.
%
% Aufrufmoeglichkeiten aus MATLAB fuer Speichervariante a:
% ========================================================
%
% Die Anzahl der Gleichungen (Zeilenzahl der Matrix A) und die
% gesamte Bandweite IBW werden aus dem Format der Matrix A entnommen,
% die Anzahl der "rechten Seiten" wird dem Format (Spaltenanzahl) der
% Matrix B entnommen. Der Aufruf fuer Speichervariante a muss drei oder
% zwei Eingabeparameter haben. Der Aufruf
%
% X = gabamp (A , B , hdsp)
%
% kann eine Matrix A mit zwei unterschiedlichen Bandweiten IBWL und IBWR
% uebergeben. Die Hauptdiagonalelemente von A muessen in der Spalte hdsp
% (entspricht der linken Bandweite IBWL) stehen.
%
% Der dritte Parameter kann entsprechend
%
% X = gabamp (A , B)
%
% weggelassen werden, wenn die beiden Bandweiten IBWL und IBWR gleich sind.
% In diesem Fall muss die Matrix A eine ungerade Spaltenanzahl haben, so
% dass die Hauptdiagonalelemente in der Mittelspalte stehen.
%
% Aufrufmoeglichkeiten aus MATLAB fuer Speichervariante a:
% ========================================================
%
% Der einfachste Aufruf der Funktion gabamp fuer Speichervariante b
% (Matrix A ist als kompakter Vektor gespeichert)
%
% X = gabamp (A , B , N , IBWR)
%
% loest ein Gleichungssystem A*X=B mit der N*N-Matrix A. Fuer die beiden
% Bandweiten gilt IBWL=IBWR, so dass fuer A wie oben beschrieben
% N*(2*IBWR-1) Elemente gespeichert sein muessen, B ist ein Vektor mit N
% Elementen.
%
% Wenn die beiden Bandweiten unterschiedlich sind, muss der Aufruf 5
% Eingabeparameter enthalten:
%
% X = gabamp (A , B , N , IBWR , IBWL)
%
% Das Gleichungssystem kann simultan fuer mehrere "rechte Seiten" geloest
% werden. Der zweite Parameter muss dann eine Matrix mit N Zeilen und NB
% Spalten sein, wobei NB als 6. Eingabeparameter angegeben werden muss
% (das Ergebnis wird dann in einer Matrix gleichen Formats abgeliefert):
%
% X = gabamp (A , B , N , IBWR , IBWL , NB)
%
% Die Matrix A wird als singulaer angenommen, wenn der Absolutbetrag eines
% Pivotelements beim Gauss-Algorithmus kleiner als das EPS-fache erste
% Pivotelement wird. Als Standardwert wird EPS=1.e-10 verwendet. Wenn man
% mit einem anderen EPS-Wert arbeiten moechte, kann dieser als 7.
% Eingabeparameter hinzugefuegt werden, z. B.:
%
% X = gabamp (A , B , N , IBWR , IBWL , NB , 1.e-6)
%
% Wenn die Matrix A singulaer ist, wird die Rechnung mit einer
% Fehlerausschrift abgebrochen. Wenn man das nicht moechte, dann muss man
% gabamp mit zwei Ausgabeparametern aufrufen, z. B.:
%
% [X, IREG] = gabamp (A , B , N , IBWR , IBWL)
%
% Dann erfolgt bei singulaerer Matrix A kein Abbruch, Singularitaet wird
% durch den Wert 0 fuer IREG signalisiert, was gegebenenfalls im aufrufenden
% Script abgeprueft werden muss. Bei regulaerer Matrix A hat IREG entweder
% den Wert 1 oder den Wert -1.
%
% ***************************************************************************
%
function [X , IREG] = gabamp (A , B , N , IBWR , IBWL , NB , EPS)
%
if nargin < 4
[n , ibw] = size (A) ;
[m , NB ] = size (B) ;
if m*NB == n
NB = 1 ;
end
if nargin < 3
IBWR = round (ibw/2 - 0.1) + 1 ;
IBWL = IBWR ;
if ibw ~= IBWR+IBWL-1
IREG = 0 ;
error ('Fehler: Anzahl der Spalten der Koeffizientenmatrix muss eine ungerade Zahl sein!')
end
else
IBWL = N ;
IBWR = ibw - IBWL + 1 ;
end
N = n ;
A = A' ;
else
if nargin < 6, NB = 1 ; end
if nargin < 5, IBWL = IBWR ; end
end
%
if nargin < 7, EPS = 1.e-10 ; end
%
[msa,nsa] = size (A) ;
if (msa*nsa ~= (IBWR+IBWL-1)*N)
IREG = 0 ;
error ('Fehler: Anzahl der Elemente der Koeffizientenmatrix!')
end
%
[msb,nsb] = size (B) ;
if (msb*nsb ~= N*NB)
IREG = 0 ;
error ('Fehler: Anzahl der Elemente der "rechten Seite B"!')
end
%
ZERO=0. ;
EINS=1. ;
IREG=1 ;
%
% Definition von Hilfsgroessen:
%
MBH =N ;
IBW =IBWL+IBWR-1 ;
IBW1=IBW-1 ;
NIBW=N*IBW-IBW1 ;
NMBH=NB*MBH ;
IOG =(IBWL-1)*IBW ;
%
if (IBWL ~= 1)
%
% Vorbereitendes Linksverschieben der ersten (IBWL-1) Zeilen:
%
JANF=IBWL ;
IGR =IBWR-1 ;
%
for I=1:IBW:IOG
JJ =JANF ;
JOG=I+IGR ;
JUG=JOG+1 ;
for J=I:JOG
A(J)=A(JJ) ;
JJ=JJ+1 ;
end % for J=I:JOG
JOG=I+IBW1 ;
for J=JUG:JOG
A(J)= ZERO ;
end % for J=JUG:JOG
IGR=IGR+1 ;
JANF=JANF+IBW1 ;
end % for I=1:IBW:IOG
end % if (IBWL ~= 1)
II=1;
%
% Schleife ueber allen Matrixzeilen (II - Zeilennummer,
% I - Position des ersten Elements der I-ten Zeile):
%
for I=1:IBW:NIBW
IZEIE=I+IBW1 ;
%
% PIVOTSUCHE:
%
JOG=I+IOG ;
if (JOG > NIBW)
JOG=NIBW ;
end % if (JOG > NIBW)
%
AMX=ZERO ;
III=II ;
for J=I:IBW:JOG
if (AMX <= abs(A(J)))
AMX=abs(A(J)) ;
IBI=III ;
IMX=J ;
end % if (AMX <= abs(A(J)))
III=III+1 ;
end % for J=I:IBW:JOG
%
% Definition der Abbruchschranke:
%
if (I == 1)
AMX1=AMX*abs(EPS) ;
end % if (I == 1)
%
% Abbruchtest:
%
if (AMX <= AMX1)
IREG = 0 ;
if (nargout == 1)
error ('Matrix A ist singulaer!')
end
return
end % if (AMX <= AMX1)
%
% Zeilentausch:
%
if (I == NIBW)
break
end
if (IMX ~= I)
%
IREG=-IREG ;
for J=I:IZEIE
S =A(J) ;
A(J) =A(IMX) ;
A(IMX)=S ;
IMX=IMX+1 ;
end % for J=I:IZEIE
%
if (NB ~= 0)
for J=II:MBH:NMBH
S =B(J) ;
B(J) =B(IBI) ;
B(IBI)=S ;
IBI =IBI+MBH ;
end % for J=II:MBH:NMBH
end % if (NB ~= 0)
end % if (IMX ~= I)
%
% Elimination:
%
S =EINS/A(I) ;
JUG=I+IBW ;
III=II ;
KUG=I+1 ;
for J=JUG:IBW:JOG
SM=A(J)*S ;
KK=J ;
for K=KUG:IZEIE
A(KK)=A(KK+1)-A(K)*SM ;
KK =KK+1 ;
end % for K=KUG:IZEIE
A(KK)=ZERO ;
if (NB ~= 0)
III=III+1 ;
IK =III ;
for K=II:MBH:NMBH
B(IK)=B(IK)-B(K)*SM ;
IK =IK+MBH ;
end % for K=II:MBH:NMBH
end % if (NB ~= 0)
end
%
II=II+1 ;
end % for I=1:IBW:NIBW
%
% ***********************************************************
% * Rueckwaertseinsetzen fuer Gauss-Algorithmus *
% * mit unsymmetrischer Bandmatrix *
% ***********************************************************
%
% Definition von Hilfsgroessen:
%
MXH = N ;
I =(N-1)*IBW+1 ;
X = zeros (N,NB) ;
%
% Letzte Zeile:
%
K=N ;
S=EINS/A(I) ;
for J=N:MBH:NMBH
X(K)=B(J)*S ;
K =K+MXH ;
end % for J=N:MBH:NMBH
%
% zeilen N-1 ... 1:
%
II=1 ;
KB=N ;
%
while (1 == 1)
KB=KB-1 ;
if (KB < 1)
return
end
JOG=I-1 ;
I =I-IBW ;
if (II < IBW)
JOG=I+II ;
end
II=II+1 ;
JUG=I+1 ;
%
% JUG,JOG - Grenzindizes der Nichtdiagonalelemente in der
% aktuellen Zeile von A
%
D =EINS/A(I) ;
KK=KB ;
%
% Schleife ueber die rechten Seiten:
%
for K=KB:MBH:NMBH
KKK=KK+1 ;
S =B(K) ;
if (JUG <= JOG)
for J=JUG:JOG
S =S-A(J)*X(KKK) ;
KKK=KKK+1 ;
end
end % if (JUG <= JOG)
X(KK)=S*D ;
KK =KK+MXH ;
end % for K=KB:MBH:NMBH
end % while (1 == 1)
|
