2.6 Singular Value Decomposition 69 f=s*rv1[1] rvi[i]=c*rvi[i]; if ((float)(fabs(f)+anorm) =anorm)break; g=w[i]; h=pythag(f,g); w[i]=h; h=1.0/h: c=g*h; s =-fwh; for (j=1;j<=m;j++) y=a[j][nm]; z=a[j][i]; http://www.nr.com or call 1-800-872- Permission is a[j门[nm]=y*c+z*s; a[j门[i]=z*c-y*s; 2 z-w[k]i if (1==k){ Convergence 1f(z<0.0)[ Singular value is made nonnegative w]=-z; for （j=1j<n;j+) v[j][k]=-v[j]k]: break; if (its ==30)nrerror("no convergence in 30 svdcmp iterations"); -7423 (North America x=w[1]; Shift from bottom 2-by-2 minor. nm=k-1; y=w[nm]; granted for internet users to make one paper copy for their g=rv1 [nm]: h=rv1[k]; f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y); g=pythag(f,1.0): f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x; c=g=1.0; Next QR transformation: for (j=l;j<=nm;j++) 1=j+1; grv1[1]; y=w[i]; Sample page from NUMERICAL RECIPES IN C:THE ART OF SCIENTIFIC COMPUTING(ISBN 0-521-43108-5) h=s*g; Copyright (C)1988-1992 by Cambridge University Press.Programs Copyright(C)1988-1992 by Numerical Recipes 8=C*g; z=pythag(f,h); rv1[j]=z; c=f/z: s=h/z; f=x*c+g*s; gg*c-x*s h=y*s; y*=c; for (jj=1;jj<=mijj++){ only),orsend email to directcustserv@cambridge.org (outside North America) x=v[jj][j]; Software. z=v[jj][1; v[jj][j]=x*c+z*s; v[ji][i]=z*c-x*s; ying of machine z=pythag(f,h); w[j]-z; Rotation can be arbitrary if z =0. if (z){ z=1.0/2: c=f*2; s=h*z; f=c*g+s*y; x=c*y-s*gi2.6 Singular Value Decomposition 69 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). f=s*rv1[i]; rv1[i]=c*rv1[i]; if ((float)(fabs(f)+anorm) == anorm) break; g=w[i]; h=pythag(f,g); w[i]=h; h=1.0/h; c=g*h; s = -f*h; for (j=1;j<=m;j++) { y=a[j][nm]; z=a[j][i]; a[j][nm]=y*c+z*s; a[j][i]=z*c-y*s; } } } z=w[k]; if (l == k) { Convergence. if (z < 0.0) { Singular value is made nonnegative. w[k] = -z; for (j=1;j<=n;j++) v[j][k] = -v[j][k]; } break; } if (its == 30) nrerror("no convergence in 30 svdcmp iterations"); x=w[l]; Shift from bottom 2-by-2 minor. 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=pythag(f,1.0); f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x; c=s=1.0; Next QR transformation: for (j=l;j<=nm;j++) { i=j+1; g=rv1[i]; y=w[i]; h=s*g; g=c*g; z=pythag(f,h); rv1[j]=z; c=f/z; s=h/z; f=x*c+g*s; g = g*c-x*s; h=y*s; y *= c; for (jj=1;jj<=n;jj++) { x=v[jj][j]; z=v[jj][i]; v[jj][j]=x*c+z*s; v[jj][i]=z*c-x*s; } z=pythag(f,h); w[j]=z; Rotation can be arbitrary if z = 0. if (z) { z=1.0/z; c=f*z; s=h*z; } f=c*g+s*y; x=c*y-s*g;