Numerik / sources / Diff of /MatLib.mod.m2pp

Diff of /MatLib.mod.m2pp [139f46] .. [28b809]

Switch to unified view

-a/MatLib.mod.m2pp
+b/MatLib.mod.m2pp
 ...
   (* 07.12.17, MRi: Einfuegen von StuerzMN,NeumaierSum & NeumaierProdSum    *)
   (*                aus Testumgebung                                        *)
   (* 30.05.18, MRi: Einfuegen von TriInfNorm                                *)
   (* 10.06.18, MRi: ZufallMat um den nicht-quadratischen Fall ergaenzt      *)
   (* 12.09.18, MRi: CMatVekProd erweitert und an MatVekProd angepasst       *)
+  (* 21.09.18, MRi: CStuerzMN eingefuegt                                    *)
+  (* 28.09.18, MRi: MatNormL1 auf "Rechtecksfall" erweitert                 *)
+  (* 29.09.18, MRi: In StuerzMN & CStuerzMN das Nullsetzen der Restmatrix   *)
+  (*                ueberarbeitet, in CStuerzMN den Parameter IfConj        *)
+  (*                eingefuegt                                              *)
   (*------------------------------------------------------------------------*)
   (* Offene Punkte                                                          *)
   (*                                                                        *)
   (* - Nachsehen ob InitMat nicht "gelitten" hat                            *)
   (* - Komplexe Typen auf ISO M2 umstellen                                  *)
 ...
   (*------------------------------------------------------------------------*)
   (* Implementation : Michael Riedl                                         *)
   (* Licence        : GNU Lesser General Public License (LGPL)              *)
   (*------------------------------------------------------------------------*)
-  (* $Id: MatLib.mod.m2pp,v 1.16 2018/06/09 13:29:13 mriedl Exp mriedl $ *)
+  (* $Id: MatLib.mod.m2pp,v 1.17 2019/02/01 22:24:03 mriedl Exp mriedl $ *)
 FROM SYSTEM    IMPORT TSIZE,ADR;
 FROM Storage   IMPORT ALLOCATE,DEALLOCATE;
 FROM Deklera   IMPORT VEKTOR,MATRIX,SUPERVEKTOR,PMATRIX,MaxDim,
                       CVEKTOR,CMATRIX;
 FROM Errors    IMPORT Fehlerflag,Fehler;
                IMPORT Errors;
 FROM LongMath  IMPORT sqrt;
                IMPORT LongComplexMath;
+FROM LMathLib  IMPORT MaxCard;
 FROM RandomLib IMPORT Zufall;
 FROM LibDBlas  IMPORT idamax,dscal;
                IMPORT LibDBlasM2;
 FROM DynMat    IMPORT AllocMat,DeAllocMat;
 ...
 END Stuerz;
 PROCEDURE StuerzMN(VAR A   : ARRAY OF ARRAY OF LONGREAL;
                        M,N : CARDINAL);
-          VAR i,j : CARDINAL;
+          VAR i,j,mn : CARDINAL;
-              T   : PMATRIX;
+              T      : PMATRIX;
 BEGIN
+      IF (M = N) THEN
+        Stuerz(A,N);
+      ELSE
+        mn := MaxCard(M,N);
-      IF (HIGH(A) < N) OR (HIGH(A[0]) < M) THEN
+        IF (HIGH(A) < mn) OR (HIGH(A[0]) < mn) THEN
-        Errors.FatalError("Matrix A falsch dimensioniert (MatLib.StuerzMN)");
+          Errors.FatalError("Matrix A falsch dimensioniert (MatLib.StuerzMN)");
-      END;
+        END;
-      AllocMat(T,M,N);
+        AllocMat(T,M,N);
-      IF (T = NIL) THEN
+        IF (T = NIL) THEN
-        Errors.FatalError("Kein Freispeicher vorhanden (MatLib.StuerzMN)");
+          Errors.FatalError("Kein Freispeicher vorhanden (MatLib.StuerzMN)");
-      END;
+        END;
-      FOR i:=0 TO M-1 DO
+        FOR i:=0 TO M-1 DO
-        FOR j:=0 TO N-1 DO T^[i+1]^[j+1]:=A[i,j]; END;
+          FOR j:=0 TO N-1 DO T^[i+1]^[j+1]:=A[i,j]; END;
-      END;
+        END;
-      FOR i:=N TO HIGH(A) DO
-        FOR j:=M TO HIGH(A[0]) DO
+        IF (1 = 1) THEN (* This code is not really necessary *)
-          A[i,j]:=0.0;
+          (* Initialisiere den Rest der (quadratischen) Matrix mit 0 *)
+          IF (N < M) THEN
+            FOR i:=N TO M-1 DO
+              FOR j:=0 TO M-1 DO A[i,j]:=0.0; END;
-        END;
+            END;
+          ELSE
+            FOR i:=0 TO N-1 DO
+              FOR j:=M TO N-1 DO A[j,i]:=0.0; END;
-      END;
+            END;
+          END; (* IF M > N *)
+        END; (* IF *)
-      FOR i:=0 TO M-1 DO
+        FOR i:=0 TO M-1 DO
-        FOR j:=0 TO N-1 DO A[j,i]:=T^[i+1]^[j+1]; END;
+          FOR j:=0 TO N-1 DO
+            A[j,i]:=T^[i+1]^[j+1];
+          END;
-      END;
+        END;
-      DeAllocMat(T,M,N);
+        DeAllocMat(T,M,N);
+      END;
 END StuerzMN;
+PROCEDURE CStuerzMN(VAR A      : ARRAY OF ARRAY OF LONGCOMPLEX;
+                        M,N    : CARDINAL;
+                        IfConj : BOOLEAN);
+          VAR i,j,mn  : CARDINAL;
+              Tre,Tim : PMATRIX;
+BEGIN
+      IF (M = N) THEN
+        CStuerz(A,N,IfConj);
+      ELSE
+        mn := MaxCard(M,N);
+        IF (HIGH(A) < mn) OR (HIGH(A[0]) < mn) THEN
+          Errors.FatalError("Matrix A falsch dimensioniert (MatLib.StuerzMN)");
+        END;
+        AllocMat(Tre,M,N);
+        AllocMat(Tim,M,N);
+        IF (Tre = NIL) OR (Tim = NIL) THEN
+          Errors.FatalError("Kein Freispeicher vorhanden (MatLib.CStuerzMN)");
+        END;
+        FOR i:=0 TO M-1 DO
+          FOR j:=0 TO N-1 DO
+            Tre^[i+1]^[j+1] := RE(A[i,j]);
+            Tim^[i+1]^[j+1] := IM(A[i,j]);
+          END;
+        END;
+        IF (1 = 1) THEN (* This code is not really necessary *)
+          (* Initialisiere den Rest der (quadratischen) Matrix mit 0 *)
+          IF (N < M) THEN
+            FOR i:=N TO M-1 DO
+              FOR j:=0 TO M-1 DO A[i,j]:=zero; END;
+            END;
+          ELSE
+            FOR i:=0 TO N-1 DO
+              FOR j:=M TO N-1 DO A[i,j]:=zero; END;
+            END;
+          END; (* IF M > N *)
+        END; (* IF *)
+        FOR i:=0 TO M-1 DO
+          IF IfConj THEN
+            FOR j:=0 TO N-1 DO
+              A[j,i]:=CMPLX(Tre^[i+1]^[j+1], - Tim^[i+1]^[j+1]);
+            END;
+          ELSE
+            FOR j:=0 TO N-1 DO
+              A[j,i]:=CMPLX(Tre^[i+1]^[j+1],   Tim^[i+1]^[j+1]);
+            END;
+          END; (* IF *)
+        END;
+        DeAllocMat(Tre,M,N);
+        DeAllocMat(Tim,M,N);
+      END;
+END CStuerzMN;
 PROCEDURE Transpose(VAR A : ARRAY OF ARRAY OF LONGREAL;
                         n : CARDINAL);
           (*----------------------------------------------------------------*)
 ...
       END;
       RETURN norm;
 END TriInfNorm;
 PROCEDURE MatNormL1(VAR A    : ARRAY OF ARRAY OF LONGREAL;
-                        N    : INTEGER;
+                        M,N  : CARDINAL;
                     VAR Norm : LONGREAL);
-          VAR i,j : INTEGER;
+          VAR i,j : CARDINAL;
               Sum : LONGREAL;
 BEGIN
       Norm:=0.0;
-      FOR i:=0 TO N-1 DO
+      FOR i:=0 TO M-1 DO
         Sum:=0.0;
         FOR j:=0 TO N-1 DO Sum:=Sum + ABS(A[i,j]); END;
         IF (Sum > Norm) THEN Norm:=Sum; END;
       END;
 END MatNormL1;