PROCEDURE svdcmp(VAR a: glmpbynp; m,n,mp,np: integer;
		 VAR w: glnparray; VAR v: glnpbynp);


{
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % 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
% %
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
}


(* 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,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
      IF (a > b) THEN max := a ELSE max := b
   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;
      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
	       BEGIN
		  ll:=l;
		  GOTO 2;
	       END;
	    IF nm>0 THEN {* Put in by me - see book *}
	       IF ((abs(w[nm])+anorm) = anorm) THEN
		  BEGIN
		     ll:=l;
		     GOTO 1
		  END;
         END;
1:         c := 0.0;
         s := 1.0;
         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 := pythag(f,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 (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
         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); 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 := 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 := pythag(f,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);**}
	    z := pythag(f,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);

         END;
         rv1[l] := 0.0;
         rv1[k] := f;
         w[k] := x
      END;
3:   END
END;