ADDED   mttroot/mtt/bin/trans/p/svdcmp.p
Index: mttroot/mtt/bin/trans/p/svdcmp.p
==================================================================
--- /dev/null
+++ mttroot/mtt/bin/trans/p/svdcmp.p
@@ -0,0 +1,263 @@
+PROCEDURE svdcmp(VAR a: glmpbynp; m,n,mp,np: integer;
+		 VAR w: glnparray; VAR v: glnpbynp);
+
+{
+% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% % Version control history
+% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% % $Id$
+% % $Log$
+% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+}
+
+
+(* Programs using routine SVDCMP must define the types
+TYPE
+   glnparray = ARRAY [1..np] OF real;
+   glmpbynp = ARRAY [1..mp,1..np] OF real;
+   glnpbynp = ARRAY [1..np,1..np] OF real;
+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;
+   BEGIN
+      IF (b >= 0.0) THEN sign := abs(a) ELSE sign := -abs(a)
+   END;
+FUNCTION max(a,b: real): real;
+   BEGIN
+      IF (a > b) THEN max := a ELSE max := b
+   END;
+BEGIN
+   g := 0.0;
+   scale := 0.0;
+   anorm := 0.0;
+   FOR i := 1 TO n DO BEGIN
+      l := i+1;
+      rv1[i] := scale*g;
+      g := 0.0;
+      s := 0.0;
+      scale := 0.0;
+      IF (i <= m) THEN BEGIN
+         FOR k := i TO m DO BEGIN
+            scale := scale+abs(a[k,i])
+         END;
+         IF (scale <> 0.0) THEN BEGIN
+            FOR k := i TO m DO BEGIN
+               a[k,i] := a[k,i]/scale;
+               s := s+a[k,i]*a[k,i]
+            END;
+            f := a[i,i];
+            g := -sign(sqrt(s),f);
+            h := f*g-s;
+            a[i,i] := f-g;
+            IF (i <> n) THEN BEGIN
+               FOR j := l TO n DO BEGIN
+                  s := 0.0;
+                  FOR k := i TO m DO BEGIN
+                     s := s+a[k,i]*a[k,j]
+                  END;
+                  f := s/h;
+                  FOR k := i TO m DO BEGIN
+                     a[k,j] := a[k,j]+
+                        f*a[k,i]
+                  END
+               END
+            END;
+            FOR k := i TO m DO BEGIN
+               a[k,i] := scale*a[k,i]
+            END
+         END
+      END;
+      w[i] := scale*g;
+      g := 0.0;
+      s := 0.0;
+      scale := 0.0;
+      IF ((i <= m) AND (i <> n)) THEN BEGIN
+         FOR k := l TO n DO BEGIN
+            scale := scale+abs(a[i,k])
+         END;
+         IF (scale <> 0.0) THEN BEGIN
+            FOR k := l TO n DO BEGIN
+               a[i,k] := a[i,k]/scale;
+               s := s+a[i,k]*a[i,k]
+            END;
+            f := a[i,l];
+            g := -sign(sqrt(s),f);
+            h := f*g-s;
+            a[i,l] := f-g;
+            FOR k := l TO n DO BEGIN
+               rv1[k] := a[i,k]/h
+            END;
+            IF (i <> m) THEN BEGIN
+               FOR j := l TO m DO BEGIN
+                  s := 0.0;
+                  FOR k := l TO n DO BEGIN
+                     s := s+a[j,k]*a[i,k]
+                  END;
+                  FOR k := l TO n DO BEGIN
+                     a[j,k] := a[j,k]
+                        +s*rv1[k]
+                  END
+               END
+            END;
+            FOR k := l TO n DO BEGIN
+               a[i,k] := scale*a[i,k]
+            END
+         END
+      END;
+      anorm := max(anorm,(abs(w[i])+abs(rv1[i])));
+   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;
+            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;
+               FOR k := l TO n DO BEGIN
+                  v[k,j] := v[k,j]+s*v[k,i]
+               END
+            END
+         END;
+         FOR j := l TO n DO BEGIN
+            v[i,j] := 0.0;
+            v[j,i] := 0.0
+         END
+      END;
+      v[i,i] := 1.0;
+      g := rv1[i];
+      l := i
+   END;
+   FOR i := n DOWNTO 1 DO BEGIN
+      l := i+1;
+      g := w[i];
+      IF (i < n) THEN BEGIN
+         FOR j := l TO n DO BEGIN
+            a[i,j] := 0.0
+         END
+      END;
+      IF (g <> 0.0) THEN BEGIN
+         g := 1.0/g;
+         IF (i <> n) THEN BEGIN
+            FOR j := l TO n DO BEGIN
+               s := 0.0;
+               FOR k := l TO m DO BEGIN
+                  s := s+a[k,i]*a[k,j]
+               END;
+               f := (s/a[i,i])*g;
+               FOR k := i TO m DO BEGIN
+                  a[k,j] := a[k,j]+f*a[k,i]
+               END
+            END
+         END;
+         FOR j := i TO m DO BEGIN
+            a[j,i] := a[j,i]*g
+         END
+      END ELSE BEGIN
+         FOR j := i TO m DO BEGIN
+            a[j,i] := 0.0
+         END
+      END;
+      a[i,i] := a[i,i]+1.0
+   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 nm>0 THEN {* Put in by me - see book *}
+	       IF ((abs(w[nm])+anorm) = anorm) THEN GOTO 1
+         END;
+1:         c := 0.0;
+         s := 1.0;
+         FOR i := l TO k DO BEGIN
+            f := s*rv1[i];
+            IF ((abs(f)+anorm) <> anorm) THEN BEGIN
+	       g := w[i];
+               h := sqrt(f*f+g*g);
+               w[i] := h;
+               h := 1.0/h;
+               c := (g*h);
+               s := -(f*h);
+               FOR j := 1 TO m DO BEGIN
+                  y := a[j,nm];
+                  z := a[j,i];
+                  a[j,nm] := (y*c)+(z*s);
+                  a[j,i] := -(y*s)+(z*c)
+               END
+            END
+         END;
+2:         z := w[k];
+         IF (l = 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
+         END;
+         GOTO 3
+         END;
+         IF (its = 30) THEN BEGIN
+            writeln ('no convergence in 30 SVDCMP iterations'); readln
+         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);
+         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
+            i := j+1;
+            g := rv1[i];
+            y := w[i];
+            h := s*g;
+            g := c*g;
+            z := sqrt(f*f+h*h);
+            rv1[j] := z;
+            c := f/z;
+            s := h/z;
+            f := (x*c)+(g*s);
+            g := -(x*s)+(g*c);
+            h := y*s;
+            y := y*c;
+            FOR jj := 1 TO n DO BEGIN
+               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);
+            w[j] := z;
+            IF (z <> 0.0) THEN BEGIN
+               z := 1.0/z;
+               c := f*z;
+               s := h*z
+            END;
+            f := (c*g)+(s*y);
+            x := -(s*g)+(c*y);
+            FOR jj := 1 TO m DO BEGIN
+               y := a[jj,j];
+               z := a[jj,i];
+               a[jj,j] := (y*c)+(z*s);
+               a[jj,i] := -(y*s)+(z*c)
+            END
+         END;
+         rv1[l] := 0.0;
+         rv1[k] := f;
+         w[k] := x
+      END;
+3:   END
+END;