DEFINITION  MODULE LinLib;

  (*------------------------------------------------------------------------*)
  (* Modul zur L"osung linearer Gleichungssystem                            *)
  (* Modul for solving lineare systems of equations.                        *)
  (*------------------------------------------------------------------------*)
  (* Implementation : Michael Riedl                                         *)
  (* Licence        : GNU Lesser General Public License (LGPL)              *)
  (*------------------------------------------------------------------------*)

  (* $Id: LinLib.def,v 1.11 2018/08/31 14:14:46 mriedl Exp mriedl $ *)

FROM Deklera IMPORT PMATRIX;

PROCEDURE CondHager(VAR A     : ARRAY OF ARRAY OF LONGREAL;
                        N     : INTEGER;
                    VAR Cond  : LONGREAL;
                    VAR iFehl : INTEGER);

         (*-----------------------------------------------------------------*)
         (* Estimates the L1 condition number of a matrix                   *)
         (*                                                                 *)
         (* Hager, William; Condition estimates, SIAM journal on scientific *)
         (* and statistical computing, 5,2, 311-316  (1984)                 *)
         (*                                                                 *)
         (* This procedure is based on a Fortran version by John Burkardt   *)
         (*                                                                 *)
         (* iFehl = 1 : Matrix A cannot be factorizised                     *)
         (* iFehl = 2 : not enought free memory                             *)
         (*                                                                 *)
         (* In case iFehl # 0 Cond is MAX(LONGREAL)                         *)
         (*-----------------------------------------------------------------*)

PROCEDURE Zerlege(VAR A        : ARRAY OF ARRAY OF LONGREAL;
                  VAR IndexI   : ARRAY OF CARDINAL;
                  VAR IndexJ   : ARRAY OF CARDINAL;
                  VAR RPivot   : ARRAY OF LONGREAL;
                  VAR Det      : LONGREAL;   (* Determinante der Matrix     *)
                      dim      : CARDINAL;   (* Dimension der qadr. Matrix  *)
                      MaxPivot : BOOLEAN);   (* w=Totalpivot, f=Zeilenpivot *)

          (*----------------------------------------------------------------*)
          (* Zerlegt die Matrix A nach dem Gaussschen Eliminationsverfahren *)
          (* In den beiden Feldern IndexI und IndexJ werden dabei die       *)
          (* Informationen \"uber das Zeilen- bzw. Spaltenvertauschen f"ur  *)
          (* die Routine BerechneLoesung gespeichert.                       *)
          (* Wenn A singul\"ar oder unl\"osbar wird die globale Variable    *) 
          (* Fehler auf TRUE gesetzt.                                       *)
          (*----------------------------------------------------------------*)

PROCEDURE BerechneLoesung(VAR A      : ARRAY OF ARRAY OF LONGREAL;
                          VAR X      : ARRAY OF LONGREAL;
                              C      : ARRAY OF LONGREAL;
                          VAR RPivot : ARRAY OF LONGREAL;
                          VAR IndexI : ARRAY OF CARDINAL;
                          VAR IndexJ : ARRAY OF CARDINAL;
                              dim    : CARDINAL);

          (*----------------------------------------------------------------*)
          (* Berechnet den L\"osungsvektor X der Gleichung A X = C mit der  *)
          (* in der Zerlege nach Gauss zerlegten Matrix A                   *)
          (*                                                                *)
          (*   A      : Koeffizientenmatrix                                 *)
          (*   X      : Gesuchte Loesung des Gleichungssystems              *)
          (*   C      : Konstantenvektor, rechte Seite d. Gleichungssystems *)
          (*   RPivot : Reziproke Pivotelemente wie von Zerlege berechnet   *)
          (*   IndexI : Zeilenvertauschungsinformationen.                   *)
          (*   IndexJ : Spaltenvertauschungsinformationen. IndexJ_k = k     *)
          (*            wenn nur Zeilenpivot.                               *)
          (*   dim    : Dimension der qadr. Matrix A                        *)
          (*----------------------------------------------------------------*)

PROCEDURE Det(VAR A   : ARRAY OF ARRAY OF LONGREAL;
                  dim : CARDINAL) : LONGREAL;

          (*----------------------------------------------------------------*)
          (* Berechnet die Determinante von A, wobei die Koeffizienten-     *)
          (* matrix A nicht zerst\"ort wird.                                *)
          (* Wenn Fehler auftreten, wird MAX(LONGREAL) zur\"uckgegeben und  *)
          (* die globale Variable Fehler auf TRUE gesetzt.                  *)
          (*----------------------------------------------------------------*)

PROCEDURE Gauss(VAR A        : ARRAY OF ARRAY OF LONGREAL;
                VAR X        : ARRAY OF LONGREAL;
                    C        : ARRAY OF LONGREAL;
                    dim      : CARDINAL;
                VAR Det      : LONGREAL;  (* Determinante von A      *)
                    MaxIter  : CARDINAL;
                    PivStrat : CARDINAL;  (* Pivotstrategie          *)
                VAR iFehl    : INTEGER); (* Fehlernummer *)

          (*----------------------------------------------------------------*)
          (* Berechnet die L\"osung des Gleichngssystems Ax = c durch       *)
          (* Gausselemination. Die Pivotstrategie ist dabei entweder eine   *)
          (* Totalpivotierung f\"ur PivStrat = 1, ansonsten wird (zeit -    *)
          (* sparend) mit Zeilenpivotierung gearbeitet.                     *)
          (* Die Matrix A wird (au\3er bei iFehl = 1) nicht ver\"andert.    *)
          (* Wenn MaxIter > 1 wird der L\"osungsvektor mit dem Residuen-    *) 
          (* vektor r = (Ax - c) nach A x2 = r und x = x-x2 MaxIter-mal     *) 
          (* nachiteriert.                                                  *)
          (* Wenn die Matrix singul\"ar oder sonstwie nicht l\"osbar ist,   *)
          (* wird die globale Variable Fehler auf TRUE gesetzt.             *)
          (*                                                                *)
          (*   A        : Koeffizientenmatrix                               *)
          (*   X        : L\"osungsvektor des LG                            *)
          (*   C        : Konstantenvektor des LG                           *)
          (*   dim      : Dimension von A                                   *)
          (*   Det      : Determinante von A                                *)
          (*   MaxIter  : Maximalzahl der Nachiterationen                   *)
          (*   PivStrat : Pivotstrategie                                    *)
          (*   iFehl    : Fehlernummer                                      *)
          (*                0 : kein Fehler aufgetreten.                    *)
          (*                1 : Koeffizentenmatrix zerst\"ort               *)
          (*                    (kein Nachiterieren moeglich)               *)
          (*                2 : Iteration nicht konvergent.                 *)
          (*                3 : Koeffizentenmatix singul\"ar.               *)
          (*                4 : Dimension der Matrix < 1                    *)
          (*----------------------------------------------------------------*)

PROCEDURE Householder(VAR A  : ARRAY OF ARRAY OF LONGREAL;
                      VAR X  : ARRAY OF LONGREAL;
                          C  : ARRAY OF LONGREAL;
                          n  : CARDINAL;  (* Anzahl der Gleichungen, n >= m *)
                          m  : CARDINAL); (* Anzahl der Unbekannten         *)


          (*----------------------------------------------------------------*)
          (* Lit.: J.Stoer, 'Einf\"uhrung in die numerische Mathematik 1',  *)
          (*       4.Auflage, Berlin, Springer (1983).                      *)
          (* L\"o\3t lineares Gleichungssystem A*X = C.                     *)
          (* Wenn die Matrix singul\"ar oder sonstwie nicht l\"osbar ist,   *)
          (* wird die globale Variable Fehler auf TRUE gesetzt. Das Ver-    *)
          (* fahren kann insbesondere dann angewandt werden, wenn mehr      *)
          (* Gleichungen als Unbekannte vorhanden sind und nur diese zu-    *)
          (* s\"atzlichen Gleichungen eine eindeutige Bestimmung des        *)
          (* L\"osungsvektors zulassen.                                     *)
          (*----------------------------------------------------------------*)

PROCEDURE PHouseholder(VAR A  : PMATRIX;
                       VAR X  : ARRAY OF LONGREAL;
                       VAR C  : ARRAY OF LONGREAL;
                           m  : CARDINAL;  (* Anzahl d. Gleichungen, m >= n *)
                           n  : CARDINAL); (* Anzahl d. Unbekannten         *)

          (*----------------------------------------------------------------*)
          (* Diese Routine ist mit der obenstehenden Householderroutine bis *)
          (* auf die Form der Darstellung der Matrix A identisch.           *)
          (*----------------------------------------------------------------*)

PROCEDURE HhLin(VAR A : ARRAY OF LONGREAL;
                VAR X : ARRAY OF LONGREAL;
                VAR C : ARRAY OF LONGREAL;
                    n : CARDINAL;  (* Zahl d. Gleichungen, n >= m *)
                    m : CARDINAL); (* Zahl d. Unbekannten         *)

          (*----------------------------------------------------------------*)
          (* Die PROCEDURE HhLin ist, bis auf die Form der \"Ubergabe       *)
          (* der Parameter, identisch mit der PROCEDURE Householder.        *)
          (*                                                                *)
          (* Die Koeffizientenmatrix A mu\3 als eindimensionales, ge -      *)
          (* packtes Feld der physikalischen Dimension n*m \"ubergeben      *)
          (* werden.                                                        *)
          (*----------------------------------------------------------------*)

PROCEDURE Invert(VAR Ainv  : ARRAY OF ARRAY OF LONGREAL; (* Invertiertes A  *)
                 VAR A     : ARRAY OF ARRAY OF LONGREAL; (* Zu invertierend *)
                 VAR Det   : LONGREAL;  (* Determinante der Matrix *)
                     dim   : CARDINAL;  (* Dimension der Matrix    *)
                 VAR iFehl : INTEGER);  (* Fehlercode              *)

          (*----------------------------------------------------------------*)
          (* Berechnet die Inverse zu der Matrix A, wobei A durch die       *)
          (* Inverse \"uberschrieben wird, wenn der Parameter Ainv mit A    *)
          (* identisch ist.                                                 *)
          (* Fehlermeldung, wenn A singul\"ar, die globale Variable Fehler  *)
          (* wird dann auf TRUE gesetzt.                                    *)
          (*                                                                *)
          (*   iFehl :  0 - kein Fehler                                     *)
          (*           -1 - Matrix A ist (numerisch) singulart              *)
          (*            1 - Matrix A ist (numerisch) schlekt konditioniert  *)
          (*            2 - Speicherproblem - Abbruch                       *)
          (*                                                                *)
          (* Wenn iFehl eingangsseitig auf den Wert 99 gesetzt wird werden  *)
          (* eventuell Fehlermeldungen auf dem Fehlerkanal ausgegeben.      *)
          (*                                                                *)
          (* Alternative Routine: InvertMat                                 *)
          (*----------------------------------------------------------------*)

PROCEDURE InvertMat(VAR Ainv : ARRAY OF ARRAY OF LONGREAL;
                    VAR A    : ARRAY OF ARRAY OF LONGREAL;
                    VAR det  : LONGREAL;
                        dim  : CARDINAL);

          (*----------------------------------------------------------------*)
          (* Berechnet das Inverse der Matrix A und gibt dieses in Ainv     *)
          (* zurueck. Routine stuetzt sich auf LibDBlas und LinPack ab      *)
          (*----------------------------------------------------------------*)

PROCEDURE InvertSV(    N     : CARDINAL;
                   VAR A     : ARRAY OF LONGREAL; (* SUPERVEKTOR *)
                   VAR iFehl : INTEGER);
  
          (*----------------------------------------------------------------*)
          (* Gauss-Jordan algorithm for the in situ inversion of a          *)
          (* positive definite matrix. Only the lower triangle is           *)
          (* transformed which is stored row by row in the vector           *)
          (* A[0:n*(n+1)/2-1].                                              *)
          (*                                                                *)
          (* iFehl : -1 wenn A singulaer oder nicht positiv definit ist     *)
          (*                                                                *)
          (* Bauer/Reinsch, pp. 47 im "Handbuch", Alogol60 routine gidef2   *)
          (*----------------------------------------------------------------*)

PROCEDURE Jacobi(VAR A       : ARRAY OF ARRAY OF LONGREAL;   (* <==  *)
                 VAR X       : ARRAY OF LONGREAL;            (* <==> *)
                 VAR C       : ARRAY OF LONGREAL;            (* <==  *)
                     dim     : CARDINAL; 
                     Tol     : LONGREAL; 
                     Neu     : CARDINAL;
                     MaxIter : CARDINAL; 
                 VAR iFehl   : INTEGER);

          (*----------------------------------------------------------------*)
          (* Loest den Satz linearer Gleichungssystem A X = C nach der      *)
          (* Jacobi-Methode.                                                *)
          (*                                                                *)
          (* Tol    : Genauigkeitsangabe für                                *)
          (*            $\delta = \sum_{i=1}^dim \| X_i^k - X_i^{k-1} \|$   *)
          (*          der gefundenen L"osungsvektoren zweier aufeinander    *)
          (*          folgenden Iterationen k,k-1. Wenn $\delta < Tol$      *)
          (*          gilt die L"osung als gefunden                         *)
          (* Neu:                                                           *) 
          (*   1    : Der Lösungsvektor wird neu mit Nullen initialisiert   *)
          (*   0    : Der Lösungsvektor wird von der rufenden Prouzedur     *)
          (*          "ubernommen.                                          *)
          (* MaxIter:                                                       *)
          (*   0    : Die Routine setzt MaxIter auf 8*(dim**3)              *)
          (*  >0    : Die Maximalzahl der Iterationen wird durch den Auf-   *)
          (*          rufenden Prozess bestimmt.                            *)
          (* iFehl:                                                         *)
          (*   0    : Es ist kein Fehler aufgetreten.                       *)
          (*   1    : Matrix ist nicht diagonal dominant, aber es wurde     *)
          (*          eine Lösung gefunden (h\"oherer Zeitaufwand).         *)
          (*   2    : Die Matrix ist diagonal dominat aber MaxIter          *)
          (*          wurde "uberschritten. Die L"osung kann ungenau sein.  *)
          (*   3    : Matrix ist nicht diagonal dominant und MaxIter        *)
          (*          wurde "uberschritten.                                 *)
          (*   4    : Der L"oesung divergierte ueber die letzten zwei       *)
          (*          Iterationen.                                          *)
          (*   5    : Bei der Berechnung des L"oesungsvektor droht ein      *)
          (*          Ueberlauf - die L"oeung ist nicht brauchbar.          *)
          (*   6    : Ein Diagonalelement von A ist Null - es kann keine    *)
          (*          L"osung  berechnet werden.                            *)
          (* F"ur iFehl > 0 wird die globale Variable Errors.Fehlerflag     *)
          (* mit einer entsprechenden Fehlermeldung belegt.                 *)
          (* F"ur alle Fehlercodes gr"o\3er 3 wurde kein L"osung gefunden.  *)
          (*----------------------------------------------------------------*)

PROCEDURE GaussSeidel(VAR A        : ARRAY OF ARRAY OF LONGREAL;
                      VAR X        : ARRAY OF LONGREAL;
                      VAR C        : ARRAY OF LONGREAL;
                          dim      : CARDINAL;
                          genau    : LONGREAL;
                          MaxIter  : CARDINAL;
                          Neu      : CARDINAL;
                          FehlMeld : BOOLEAN;
                      VAR Iter     : CARDINAL;
                      VAR iFehl    : INTEGER);

          (*----------------------------------------------------------------*)
          (* Solves the (diagonal dominant) linear system of equations      *)
          (* A x X = C by the (iterativ) Gauss - Seidel Algorithm.          *)
          (* (E - L)x^{k+1} = C + Ux^k mit A = E-U-L.                       *)
          (*                                                                *)
          (* L\"o\3t das (diagonal dominante) lineare Gleichungssystem      *)
          (* A x X = C nach dem (iterativen) Gauss - Seidel Algorithmus.    *)
          (* (E - L)x^{k+1} = C + Ux^k mit A = E-U-L.                       *)
          (*                                                                *)
          (*    A        : Koeffizinten - Matrix                            *)
          (*    X        : L\"osungsvektor                                  *)
          (*    C        : Konstantenvektor                                 *)
          (*    dim      : Dimension der Koeffizientenmatrix                *)
          (*    genau    : Abbruchkriterium                                 *)
          (*    MaxIter  : Maximalzahl der Iterationen                      *)
          (*    Neu      : s.u.                                             *)
          (*    FehlMeld : Fehlermeldungen j/n 1/0                          *)
          (*    Iter     : Anzahl der benoetigten Iterationen               *)
          (*    iFehl    : s.u.                                             *)
          (*                                                                *)
          (* Neu:                                                           *)
          (*                                                                *)
          (*   0 : Der Startvektor f\"ur X wird vorgegeben.                 *)
          (*   1 : Der Startvektor f\"ur X wird neu berechnet.              *)
          (*                                                                *)
          (* iFehl:                                                         *)
          (*                                                                *)
          (*   0 : Es ist kein Fehler aufgetreten.                          *)
          (*   1 : Es ist kein Fehler aufgetreten, die Matrix ist aber      *)
          (*       nicht diagonal dominant. (h\"oherer Zeitaufwand)         *)
          (*   2 : MaxIter ist \"uberschritten                              *)
          (*   3 : Problem ist mit dieser Routine nicht l\"o\3bar da ein    *)
          (*       Diagonalelement im Zuge der Berechnung zu klein geworden *)
          (*       ist                                                      *)
          (*   4 : Die Matrix ist nicht diagonal dominant und die           *)
          (*       Berechnung der Loesung divergiert                        *)
          (*   6 : Ein Diagonalelement von A exakt Null - es kann keine     *)
          (*       L"osung  berechnet werden.                               *)
          (*                                                                *)
          (* In allen F\"allen, in denen iFehl > 0 werden auch die globalen *)
          (* Variablen Fehler und Fehlerflag entsprechend gesetzt. Wenn     *)
          (* iFehl=1 kann die Routine mit (Neu = 0) wieder gerufen werden.  *)
          (*                                                                *)
          (* FehlMeld:                                                      *)
          (*                                                                *)
          (*   FALSE : Fehler werden nicht auf dem Bildschirm angezeigt.    *)
          (*   TRUE  : Alle Fehler werden auf dem Bildschirm angezeigt.     *)
          (*                                                                *)
          (* Hinweis: Fuer ein ueberdeterminiertes Gleichungssystem         *)
          (*          konvergiert der Algorithmus zu eine nicht eindeutigen *)
          (*          Loesung. Diese Situation kann nicht erkannt werden.   *)
          (*                                                                *)
          (* Lit.: A.Jenninger, 'Matrix Computations for Engineers and      *)
          (*       Sientists', J.Wiley & Sons, London 1977.                 *)
          (*       S.J.Leon, 'Linear Algebra', Maxmillian Pub. Co.,         *)
          (*       New York 1980.                                           *)
          (*----------------------------------------------------------------*)

PROCEDURE GaussJordan(    N,M   : CARDINAL;
                      VAR A     : ARRAY OF ARRAY OF LONGREAL;
                      VAR B     : ARRAY OF ARRAY OF LONGREAL;
                          inv   : BOOLEAN;
                      VAR Det   : LONGREAL;
                      VAR iFehl : INTEGER);

          (*----------------------------------------------------------------*)
          (* Solving a linear matrix system Ax = B with Gauss-Jordan method *)
          (* using full pivoting at each step. During the process, original *)
          (* matrices A and B are destroyed                                 *)
          (*                                                                *)
          (* On Input:                                                      *)
          (*                                                                *)
          (*   N   : size of matrix A                                       *)
          (*   M   : number of simultaneous equations to be solved          *)
          (*   A   : matrix of size NxN                                     *)
          (*   B   : matrix of size NxM                                     *)
          (*   inv : if true only the inverse of A is calculated            *)
          (*                                                                *)
          (* On Output:                                                     *)
          (*                                                                *)
          (*   A   : inverse of A                                           *)
          (*   Det : determinant of A                                       *)
          (*   B   : solution matrix of size NxM                            *)
          (*                                                                *)
          (* Based on a Pascal version provided in tpmath library           *)
          (*----------------------------------------------------------------*)

PROCEDURE Cholesky(VAR A     : ARRAY OF LONGREAL; (* SUPERVEKTOR *)
                   VAR X     : ARRAY OF LONGREAL; (* L\"osungsvektor *)
                       C     : ARRAY OF LONGREAL;
                       dim   : CARDINAL;
                       neu   : BOOLEAN;
                   VAR Det   : LONGREAL;   (* Deteriminate von A *)
                   VAR iFehl : INTEGER);

          (*----------------------------------------------------------------*)
          (* Calculates the solution of a linear system of equations        *)
          (* A x X = C with symmetric, positiv definite coefficien matrix   *)
          (* A stored in super linear form.                                 *)
          (*                                                                *)
          (* Berechnet die L\"osung X eines linearen Gleichungsystem        *)
          (* A X = C mit symmetrischer, positiv definiten Koeffizienten -   *)
          (* matrix A in Superlinearform.                                   *)
          (* Ist A nicht positiv definit, so werden die globale Variabel    *)
          (* Fehler und Fehlerflag entsprechend gesetzt und auf Det wird    *)
          (* MAX(LONGREAL) zur\"uckgegeben.                                 *)
          (*                                                                *)
          (*   A     : Koeffizientenmatrix, ausgangsseitig die Cholesky-    *)
          (*           zerlegung von A                                      *)
          (*   C     : Koeffizientenvektor (rechte Seite des Gl.Sy.)        *)
          (*   X     : Loesungsvektor                                       *)
          (*   dim   : Dimension des Problems                               *)
          (*   neu   : Wenn wahr wird die Zerlegung durchgefuehrt und der   *)
          (*           der Loesungsvektor X berechnet                       *)
          (*         : Wenn falsch wird die Zerlegung nicht durchgefuehrt   *)
          (*           und nur ein neuer Loesungsvektor berechnet. Dazu     *)
          (*           muss zumindest ein Aufruf mit neu=TRUE vorausge-     *)
          (*           gangen sein und A nicht mehr veraendert worden sein  *)
          (*   Det   : Determinate der Matrix A                             *)
          (*   iFehl : iFehl =  0  - keine Fehler                           *)
          (*           iFehl >  0  - A ist (eventuell aus numerischen       *)
          (*                         Gruenden, nicht positiv definit. iFehl *)
          (*                         ist dann gleich dem Index der Gleich-  *)
          (*                         Gleichung bei der der Fehler aufge-    *)
          (*                         treten ist                             *)
          (*           iFehl = -1  - Dimensionierungsfehler                 *)
          (*----------------------------------------------------------------*)

PROCEDURE CholDet1(    N     : CARDINAL;
                   VAR A     : ARRAY OF ARRAY OF LONGREAL;
                   VAR p     : ARRAY OF LONGREAL;
                   VAR d1    : LONGREAL;
                   VAR d2    : INTEGER;
                   VAR iFehl : INTEGER);

          (*----------------------------------------------------------------*)
          (* The upper triangle of a positive definite symmetrie matrix,    *)
          (* A, is stored in the upper triangle of an NxN array             *)
          (* A[i,j], i= 1 (1) n,j= 1 (1) n.                                 *)
          (* The Cholesky deeomposition A =LU, where U is the transpose of  *)
          (* L, is performed and stored in the remainder of the array A     *)
          (* except for the reciproeals of the diagonal elements which are  *)
          (* stored in p[i], i = 1 (1) n, instead of the elements           *)
          (* themselves. A is retained so that the solution obtained can    *)
          (* be subsequently improved. The determinant,                     *)
          (* d1 * 2 ^ d2, of A is also computed. The proeedure will fail if *)
          (* A, modified by the rounding errors, is not positive definite.  *)
          (* In that case iFehl will be equal to the index of the equation  *)
          (* where the error occured, 0 otherwise                           *)
          (*                                                                *)
          (* "Handbuch" routine choldet1                                    *)
          (*----------------------------------------------------------------*)

PROCEDURE CholSol1(    N : CARDINAL;
                       r : CARDINAL;
                   VAR A : ARRAY OF ARRAY OF LONGREAL;
                   VAR p : ARRAY OF LONGREAL;
                   VAR B : ARRAY OF ARRAY OF LONGREAL;
                   VAR X : ARRAY OF ARRAY OF LONGREAL);

          (*----------------------------------------------------------------*)
          (* Solves A X = B, where A is a positive definite symmetrie       *)
          (* matrix and B is an n x r matrix of r right-hand sides. The     *)
          (* proeedure cholsol1 must be preeeded by choldetl in which L is  *)
          (* produced in A[i,j] and p[i], from A. AX=B is solved in two     *)
          (* steps, Ly=B and UX=y. The matrix B is retained in order to     *)
          (* facilitate the refinement of X, but X is overwritten on y.     *)
          (* However, X and B can be identified in the call of the          *)
          (* procedure                                                      *)
          (*                                                                *)
          (* "Handbuch" routine cholsol1                                    *)
          (*----------------------------------------------------------------*)

PROCEDURE CZerlege(VAR A      : ARRAY OF ARRAY OF LONGCOMPLEX;
                   VAR Index  : ARRAY OF CARDINAL;
                   VAR Det    : LONGCOMPLEX;
                   VAR EDet   : INTEGER;
                       dim    : CARDINAL);

          (*----------------------------------------------------------------*)
          (* Zerlegt die komplexe Matrix A. (LU = A) mit Zeilenpivotierung  *)
          (* Vertauschungen werden in Index angezeigt                       *)
          (* Fehlermeldungen wenn A singul\"ar oder unl\"osbar.             *)
          (*                                                                *)
          (*   A      : Koeffizientenmatrix                                 *)
          (*   Index  : Information ueber Vertauschungen durch Pivotierung  *)
          (*   Det    : Basis der Determinante von A                        *)
          (*   EDet   : Det(A) = Det*2^EDet                                 *)
          (*   dim    : Dimension der qadr. Matrix A                        *)
          (*----------------------------------------------------------------*)

PROCEDURE CBerechneLoesung(VAR A     : ARRAY OF ARRAY OF LONGCOMPLEX;
                           VAR X     : ARRAY OF LONGCOMPLEX;
                               C     : ARRAY OF LONGCOMPLEX;
                           VAR Index : ARRAY OF CARDINAL;
                               dim   : CARDINAL);

          (*----------------------------------------------------------------*)
          (* Berechnet den L\"osungsvektor X der Gleichung A x = c mit      *)
          (* der in der CZerlege LU-zerlegten Matrix A.                     *)
          (* Die Routine ist das komplexe Analogon zu BerechneLoesung.      *)
          (*----------------------------------------------------------------*)

PROCEDURE CDet(    dim : CARDINAL; 
               VAR A   : ARRAY OF ARRAY OF LONGCOMPLEX) : LONGCOMPLEX;

          (*----------------------------------------------------------------*)
          (* Berechnet die komplexe Determinante von A.                     *)
          (* Die Koeffizientenmatrix A wird nach M\"oglichkeit nicht zer-   *)
          (* st\"ort, andernfalls wird die globale Variable Fehler auf TRUE *)
          (* gesetzt und eine Fehlermeldung ausgeschrieben.                 *)
          (*----------------------------------------------------------------*)
 
PROCEDURE CLoeseGlSy(VAR A     : ARRAY OF ARRAY OF LONGCOMPLEX;
                     VAR X     : ARRAY OF LONGCOMPLEX; (* L\"osungsvektor.  *)
                         C     : ARRAY OF LONGCOMPLEX;
                     VAR Det   : LONGCOMPLEX; (* Determinte von A *)
                         dim   : CARDINAL;
                     VAR iFehl : INTEGER);

          (*----------------------------------------------------------------*)
          (* Berechnet die Loesung des komplexen Gleichngssystems A x X = C *)
          (* Wenn die Matrix singul\"ar oder sonstwie nicht l\"osbar        *)
          (* ist, wird dir globale Variable Fehler auf TRUE gesetzt.        *)
          (*                                                                *)
          (* iFehl : 0 wenn alles in Ordnung                                *)
          (*       : 1 wenn dim < 1                                         *)
          (*       : 2 wenn die Koeffizientenmatrix (ev. aus numerischen    *)
          (*         Gruenden) singulaer ist                                *)
          (*----------------------------------------------------------------*)

END LinLib.
