Index: mttroot/mtt/bin/trans/p/svdcmp.p
==================================================================
--- mttroot/mtt/bin/trans/p/svdcmp.p
+++ mttroot/mtt/bin/trans/p/svdcmp.p
@@ -6,10 +6,14 @@
 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % % Version control history
 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % % $Id$
 % % $Log$
+% % Revision 1.2  1998/08/12 11:05:33  peterg
+% % Taken from NR share library nrpas13 as SVDCMP.PAS
+% % and renamed svdcmp.p
+% %
 % % Revision 1.1  1998/08/12 11:03:57  peterg
 % % Initial revision
 % %
 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 }
@@ -23,14 +27,31 @@
 in the main routine. *)
 LABEL 1,2,3;
 CONST
    nmax=100;
 VAR
-   nm,l,k,j,jj,its,i: integer;
-   z,y,x,scale,s,h,g,f,c,anorm: real;
-   rv1: ARRAY [1..nmax] OF real;
-FUNCTION sign(a,b: real): real;
+   nm,l,k,j,jj,its,i,ll	       : integer;
+   z,y,x,scale,s,h,g,f,c,anorm : real;
+   rv1			       : ARRAY [1..nmax] OF real;
+			  
+FUNCTION pythag(a,b : real): real;
+VAR p,at,bt : REAL;
+BEGIN
+   at:=abs(a);
+   bt:=abs(b);
+   IF at>bt THEN
+      p:= at*sqrt(1.0+sqr(bt/at))
+   ELSE
+      IF bt=0.0 THEN
+	 p := 0.0
+      ELSE
+	 p := bt*sqrt(1.0+sqr(at/bt));
+   pythag := p;
+END{pythag};
+			   
+   
+FUNCTION sign(a,b : real): real;
    BEGIN
       IF (b >= 0.0) THEN sign := abs(a) ELSE sign := -abs(a)
    END;
 FUNCTION max(a,b: real): real;
    BEGIN
@@ -118,12 +139,12 @@
    END;
    FOR i := n DOWNTO 1 DO BEGIN
       IF (i < n) THEN BEGIN
          IF (g <> 0.0) THEN BEGIN
             FOR j := l TO n DO BEGIN
-               v[j,i] := (a[i,j]/a[i,l])/g
-            END;
+	       v[j,i] := (a[i,j]/a[i,l])/g;
+	    END;
             FOR j := l TO n DO BEGIN
                s := 0.0;
                FOR k := l TO n DO BEGIN
                   s := s+a[i,k]*v[k,j]
                END;
@@ -175,21 +196,30 @@
    END;
    FOR k := n DOWNTO 1 DO BEGIN
       FOR its := 1 TO 30 DO BEGIN
          FOR l := k DOWNTO 1 DO BEGIN
             nm := l-1;
-            IF ((abs(rv1[l])+anorm) = anorm) THEN GOTO 2;
+            IF ((abs(rv1[l])+anorm) = anorm) THEN
+	       BEGIN
+		  ll:=l;
+		  GOTO 2;
+	       END;
 	    IF nm>0 THEN {* Put in by me - see book *}
-	       IF ((abs(w[nm])+anorm) = anorm) THEN GOTO 1
+	       IF ((abs(w[nm])+anorm) = anorm) THEN
+		  BEGIN
+		     ll:=l;
+		     GOTO 1
+		  END;
          END;
 1:         c := 0.0;
          s := 1.0;
-         FOR i := l TO k DO BEGIN
+         FOR i := ll TO k DO BEGIN
             f := s*rv1[i];
             IF ((abs(f)+anorm) <> anorm) THEN BEGIN
 	       g := w[i];
-               h := sqrt(f*f+g*g);
+	       {**h := sqrt(f*f+g*g);**}
+	       h := pythag(f,g);
                w[i] := h;
                h := 1.0/h;
                c := (g*h);
                s := -(f*h);
                FOR j := 1 TO m DO BEGIN
@@ -199,11 +229,11 @@
                   a[j,i] := -(y*s)+(z*c)
                END
             END
          END;
 2:         z := w[k];
-         IF (l = k) THEN BEGIN
+         IF (ll = k) THEN BEGIN
             IF (z < 0.0) THEN BEGIN
                w[k] := -z;
                FOR j := 1 TO n DO BEGIN
                v[j,k] := -v[j,k]
             END
@@ -210,28 +240,30 @@
          END;
          GOTO 3
          END;
          IF (its = 30) THEN BEGIN
             writeln ('no convergence in 30 SVDCMP iterations'); readln
-         END;
+	 END;
          x := w[l];
          nm := k-1;
          y := w[nm];
          g := rv1[nm];
          h := rv1[k];
          f := ((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
-         g := sqrt(f*f+1.0);
+	 {***g := sqrt(f*f+1.0); writeln(g);***}
+	 g := pythag(f,1.0);
          f := ((x-z)*(x+z)+h*((y/(f+sign(g,f)))-h))/x;
          c := 1.0;
          s := 1.0;
-         FOR j := l TO nm DO BEGIN
+         FOR j := ll TO nm DO BEGIN
             i := j+1;
             g := rv1[i];
             y := w[i];
             h := s*g;
             g := c*g;
-            z := sqrt(f*f+h*h);
+            {**z := sqrt(f*f+h*h);**}
+	    z := pythag(f,h);
             rv1[j] := z;
             c := f/z;
             s := h/z;
             f := (x*c)+(g*s);
             g := -(x*s)+(g*c);
@@ -241,11 +273,12 @@
                x := v[jj,j];
                z := v[jj,i];
                v[jj,j] := (x*c)+(z*s);
                v[jj,i] := -(x*s)+(z*c)
             END;
-            z := sqrt(f*f+h*h);
+            {**z := sqrt(f*f+h*h);**}
+	    z := pythag(f,h);
             w[j] := z;
             IF (z <> 0.0) THEN BEGIN
                z := 1.0/z;
                c := f*z;
                s := h*z
@@ -263,5 +296,6 @@
          rv1[k] := f;
          w[k] := x
       END;
 3:   END
 END;
+