% ***************************************************************************
% 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)
