%PROCESS MARGINS(1,100); /* SINGULAR VALUE DECOMPOSITION AND LEAST SQUARES SOLUTIONS */ /* From J. H. Wilkinson & C. Reinsch, "Handbook for Automatic Computation", Vol. II, "Linear Algebra", Springer-Verlag, New York, 1971, pp. 141-147. */ /* p. 141 */ /* L1 added to declarations of variables, 22/3/2007. */ SVD:PROC(M,N,WITHU,WITHV,EPS,TOL,A,Q,U,V) OPTIONS (REORDER); /*PAGE 141*/ DCL (M,N,I,J,K,L,L1, LL,LBN,UBN,LBM,UBM) BIN FIXED (31,0), /*00000020*/ (EPS,TOL,C,F,G,H,X,Y,Z) FLOAT BIN(21), /*00000030*/ S FLOAT BIN (53), /*00000040*/ ((A,U)(*,*), /*00000050*/ V(*,*), /*00000060*/ Q(*),E(LBN:UBN) CTL) BIN FLOAT(21), /*00000070*/ (WITHU,WITHV)BIT(8); /*00000080*/ LBN=LBOUND(A,2); UBN=LBN+N-1; /*00000090*/ LBM=LBOUND(A,1); UBM=LBM+M-1; /*00000100*/ ALLOC E; /*00000110*/ DO I=LBN TO UBN; /*00000120*/ DO J=LBM TO UBM; /*00000130*/ U(J,I)=A(J,I); /*00000140*/ END; /*00000150*/ END; /*00000160*/ G,X=0E0; /*00000170*/ DO I=LBN TO UBN; /*00000180*/ LL=LBM+I-LBN; /*00000190*/ E(I)=G; S=0E0; L=I+1; /*00000200*/ DO J=LL TO UBM; /*00000210*/ S=S+U(J,I)**2; /*00000220*/ END; /*00000230*/ IF S=0E0 THEN G=-G; /*00000270*/ H=F*G-S; U(LL,I)=F-G; /*00000280*/ DO J=L TO UBN; /*00000290*/ S=0E0; /*00000300*/ DO K=LL TO UBM; /*00000310*/ S=S+U(K,I)*U(K,J); /*00000320*/ END; /*00000330*/ F=S/H; /*00000340*/ DO K=LL TO UBM; /*00000350*/ U(K,J)=U(K,J)+F*U(K,I); /*00000360*/ END; /*00000370*/ END; /*00000380*/ END; /*00000390*/ Q(I)=G; S=0E0; /*00000400*/ DO J=L TO UBN; /*00000410*/ S=S+U(LL,J)**2; /*00000420*/ END; /*00000430*/ IF S=0E0 THEN G=-G; /*00000470*/ H=F*G-S; U(LL,I+1)=F-G; /*00000480*/ DO J=L TO UBN; /*00000490*/ E(J)=U(I,J)/H; /*00000500*/ END; /*00000510*/ DO J=L TO UBM; /*00000520*/ S=0E0; /*00000530*/ DO K=L TO UBN; /*00000540*/ S=S+U(J,K)*U(I,K); /*00000550*/ END; /*00000560*/ DO K=L TO UBN; /*00000570*/ U(J,K)=U(J,K)+S*E(K); /*00000580*/ END; /*00000590*/ END; /*00000600*/ END; /*00000610*/ Y=ABS(Q(I))+ABS(E(I)); IF Y>X THEN X=Y; /*00000620*/ END; /*00000630*/ /*ACCUMULATION OF RIGHT-HAND TRANSFORMS */ /*00000640*/ IF WITHV^='0'B THEN DO I=UBN TO LBN BY -1; /*00000650*/ LL=LBM+I-LBN; /*00000660*/ IF G^=0E0 THEN DO; /*00000670*/ H=U(LL,I+1)*G; /*00000680*/ DO J=L TO UBN; /*00000690*/ V(J,I)=U(LL,J)/H; /*00000700*/ END; /*00000710*/ DO J=L TO UBN; /*00000720*/ S=0E0; /*00000730*/ DO K=L TO UBN; /*00000740*/ S=S+U(I,K)*V(K,J); /*00000750*/ END; /*00000760*/ DO K=L TO UBN; /*00000770*/ V(K,J)=V(K,J)+S*V(K,I); /*00000780*/ END; /*00000790*/ END; /*00000800*/ END; /*00000810*/ DO J=L TO UBN; /*00000820*/ V(I,J),V(J,I)=0E0; /*00000830*/ END; /*00000840*/ V(I,I)=1E0; G=E(I); L=I; /*00000850*/ END; /*00000860*/ /*ACCUMULATION OF LEFT-HAND TRANSFORMATIONS */ /*00000870*/ IF WITHU^='0'B THEN DO I=UBN TO LBN BY -1; /*00000880*/ LL=LBM+I-LBN; /*00000890*/ L=I+1; G=Q(I); /*00000900*/ DO J=L TO UBN; /*00000910*/ U(I,J)=0E0; /*00000920*/ END; /*00000930*/ IF G^=0E0 THEN DO; /*00000940*/ H=U(LL,I)*G; /*00000950*/ DO J=L TO UBN; /*00000960*/ S=0E0; /*00000970*/ DO K=L TO UBM; /*00000980*/ S=S+U(K,I)*U(K,J); /*00000990*/ END; /*00001000*/ F=S/H; /*00001010*/ DO K=LL TO UBM; /*00001020*/ U(K,J)=U(K,J)+F*U(K,I); /*00001030*/ END; /*00001040*/ END; /*00001050*/ DO J=LL TO UBM; /*00001060*/ U(J,I)=U(J,I)/G; /*00001070*/ END; /*00001080*/ END; /*00001090*/ ELSE DO J=LL TO UBM; /*00001100*/ U(J,I)=0E0; /*00001110*/ END; /*00001120*/ U(LL,I)=U(LL,I)+1; /*00001130*/ END; /*00001140*/ /*DIAGONALIZATION OF THE BIDIAGONAL FORM */ /*00001150*/ EPS=EPS*X; /*00001160*/ DO K=UBN TO LBN BY -1; /*00001170*/ /*TEST F SPLITTING*/ /*00001180*/ TEST_F_SPLITTING: /*00001190*/ DO L=K TO LBN BY -1; /*00001200*/ IF ABS(E(L))<=EPS THEN GO TO TEST_F_CONVERGENCE; /*00001210*/ IF ABS(Q(L-1))<=EPS THEN GO TO CANCELLATION; /*00001220*/ END; /*00001230*/ /*CANCELLATION OF E(L) IF L>1*/ /*00001240*/ CANCELLATION: /*00001250*/ C=0E0; S=1E0; L1=L-1; /*00001260*/ DO I=L TO K; /*00001270*/ F=S*E(I); E(I)=C*E(I); /*00001280*/ IF ABS(F)<=EPS THEN GO TO TEST_F_CONVERGENCE; /*00001290*/ G=Q(I); H,Q(I)=SQRT(F*F+G*G); C=G/H; S=-F/H; /*00001300*/ IF WITHU^='0'B THEN DO J=LBM TO UBM; /*00001310*/ Y=U(J,L1); Z=U(J,I);U(J,L1)=Y*C+Z*S; /*00001320*/ U(J,I)=-Y*S+Z*C; /*00001330*/ END; /*00001340*/ END; /*00001350*/ TEST_F_CONVERGENCE: /*00001360*/ Z=Q(K); IF L=K THEN GO TO CONVERGENCE; /*00001370*/ /*SHIFT FROM BOTTOM 2X2 MINOR */ /*00001380*/ X=Q(L); Y=Q(K-1); G=E(K-1); H=E(K); /*00001390*/ F=((Y-Z)*(Y+Z)+(G-H)*(G+H))/(2E0*H*Y); /*00001400*/ G=SQRT(F*F+1E0); /*00001410*/ IF F<0E0 THEN S=F-G; ELSE S=F+G; /*00001420*/ F=((X+Z)*(X-Z)+H*(Y/S-H))/X; /*00001430*/ /*NEXT QR TRANSFORMATION*/ /*00001440*/ C,S=1E0; /*00001450*/ DO I=L+1 TO K; /*00001460*/ G=E(I); Y=Q(I); H=S*G; G=C*G; /*00001470*/ E(I-1),Z=SQRT(F*F+H*H); C=F/Z; S=H/Z; /*00001480*/ F=X*C+G*S; G=-X*S+G*C; H=Y*S; Y=Y*C; /*00001490*/ IF WITHV^='0'B THEN DO J=LBN TO UBN; /*00001500*/ X=V(J,I-1); Z=V(J,I); /*00001510*/ V(J,I-1)=X*C+Z*S; V(J,I)=-X*S+Z*C; /*00001520*/ END; /*00001530*/ Q(I-1),Z=SQRT(F*F+H*H); C=F/Z; S=H/Z; /*00001540*/ F=C*G+S*Y; X=-S*G+C*Y; /*00001550*/ IF WITHU^='0'B THEN DO J=LBM TO UBM; /*00001560*/ Y=U(J,I-1); Z=U(J,I); /*00001570*/ U(J,I-1)=Y*C+Z*S; U(J,I)=-Y*S+Z*C; /*00001580*/ END; /*00001590*/ END; /*00001600*/ E(L)=0E0;E(K)=F;Q(K)=X; GO TO TEST_F_SPLITTING; /*00001610*/ CONVERGENCE: /*00001620*/ IF Z<0E0 THEN DO; /*00001630*/ /*Q(K) IS MADE NON-NEGATIVE */ /*00001640*/ Q(K)=-Z; /*00001650*/ IF WITHV^='0'B THEN V(*,K)=-V(*,K); /*00001660*/ END; /*00001670*/ END; /*00001680*/ END SVD; /*00001690*/ /* p. 144 */ MINFIT:PROC(M,N,P,EPS,TOL,AB,Q) OPTIONS (REORDER); /*PAGE 144*/ /*00000010*/ DCL (M,N,P,I,J,K,L,L1,N1,NP,LBP,UBP,LBN,UBN,II) BIN FIXED(31,0), /*00000020*/ (EPS,TOL,C,F,G,H,X,Y,Z) BIN FLOAT (21), /*00000030*/ S BIN FLOAT (53), /*00000040*/ (AB(*,*), /*00000050*/ Q(*),E(LBP:LBP+N-1) CTL) BIN FLOAT(21); /*00000060*/ /*HOUSEHOLDERS REDUCTION TO BIDIAGONAL FORM*/ /*00000070*/ LBN=LBOUND(AB,1); UBN=LBN+MAX(N,M)-1; /*00000080*/ LBP=LBOUND(AB,2); UBP=LBP+(N+P)-1; /*00000090*/ ALLOC E; /*00000100*/ G,X=0E0; NP=N+P; /*00000110*/ DO I=LBP TO (LBP+N-1); /*00000120*/ II=(LBN+N-1)+(I-(LBP+N-1)); /*00000130*/ E(I)=G;S=0E0;L=I+1; /*00000140*/ DO J=LBN+(I-LBP) TO (LBN+M-1); /*00000150*/ S=S+AB(J,I)**2; /*00000160*/ END; /*00000170*/ IF S=0E0 THEN G=-G; /*00000210*/ H=F*G-S; AB(II,I)=F-G; /*00000220*/ DO J=L TO UBP; /*00000230*/ S=0E0; /*00000240*/ DO K=LBN+(I-LBP) TO (LBN+M-1); /*00000250*/ S=S+AB(K,I)*AB(K,J); /*00000260*/ END; /*00000270*/ F=S/H; /*00000280*/ DO K=LBN+(I-LBP) TO (LBN+M-1); /*00000290*/ AB(K,J)=AB(K,J)+F*AB(K,I); /*00000300*/ END; /*00000310*/ END; /*00000320*/ END; /*00000330*/ Q(I)=G; S=0E0; /*00000340*/ IF ((I-LBP)+1)<=M THEN DO J=L TO (LBP+N-1); /*00000350*/ S=S+AB(I,J)**2; /*00000360*/ END; /*00000370*/ IF S=0E0 THEN G=-G; /*00000410*/ H=F*G-S; AB(II,I+1)=F-G; /*00000420*/ DO J=L TO (LBP+N-1); /*00000430*/ E(J)=AB(II,J)/H; /*00000440*/ END; /*00000450*/ DO J=LBN+(L-LBP) TO (LBN+M-1); /*00000460*/ S=0E0; /*00000470*/ DO K=L TO (LBP+N-1); /*00000480*/ S=S+AB(J,K)*AB(II,K); /*00000490*/ END; /*00000500*/ DO K=L TO (LBP+N-1); /*00000510*/ AB(J,K)=AB(J,K)+S*E(K); /*00000520*/ END; /*00000530*/ END; /*00000540*/ END; /*00000550*/ Y=ABS(Q(I))+ABS(E(I)); IF Y>X THEN X=Y; /*00000560*/ END; /*00000570*/ /*ACCUMULATION OF RIGHT-HAND TRANSFORMATIONS*/ /*00000580*/ DO I=(LBP+N-1) TO LBP BY -1; /*00000590*/ II=(LBN+N-1)+(I-(LBP+N-1)); /*00000600*/ IF G^=0E0 THEN DO; /*00000610*/ H=AB(II,I+1)*G; /*00000620*/ DO J=L TO (LBP+N-1); /*00000630*/ AB(LBN+(J-LBP),I)=AB(II,J)/H; /*00000640*/ END; /*00000650*/ DO J=L TO (LBP+N-1); /*00000660*/ S=0E0; /*00000670*/ DO K=L TO (LBP+N-1); /*00000680*/ S=S+AB(II,K)*AB(LBN+(K-LBP),J); /*00000690*/ END; /*00000700*/ DO K=LBN+(L-LBP) TO LBN+N-1; /*00000710*/ AB(K,J)=AB(K,J)+S*AB(K,I); /*00000720*/ END; /*00000730*/ END; /*00000740*/ END; /*00000750*/ DO J=L TO (LBP+N-1); /*00000760*/ AB(II,J),AB(LBN+(J-LBP),I)=0E0; /*00000770*/ END; /*00000780*/ AB(II,I)=1E0; G=E(I); L=I; /*00000790*/ END; /*00000800*/ EPS=EPS*X; /*00000810*/ DO I=LBN+M TO LBN+N-1; /*00000820*/ DO J=LBP+N TO LBP+NP-1; /*00000830*/ AB(I,J)=0E0; /*00000840*/ END; /*00000850*/ END; /*00000860*/ /* DIAGONALIZATION OF THE BIDIAGONAL FORM*/ /*00000870*/ DO K=LBP+N-1 TO LBP BY -1; /*00000880*/ TEST_F_SPLITTING: /*00000890*/ DO L=K TO LBP BY -1; /*00000900*/ IF ABS(E(L))<=EPS THEN GO TO TEST_F_CONVERGENCE; /*00000910*/ IF ABS(Q(L-1))<=EPS THEN GO TO CANCELLATION; /*00000920*/ END; /*00000930*/ /*CANCELLATION OF E(L) IF L>1*/ /*00000940*/ CANCELLATION: /*00000950*/ C=0E0; S=1E0; L1=L-1; /*00000960*/ DO I=L TO K; /*00000970*/ F=S*E(I); E(I)=C*E(I); /*00000980*/ IF ABS(F)<=EPS THEN GO TO TEST_F_CONVERGENCE; /*00000990*/ G=Q(I); Q(I),H=SQRT(F*F+G*G); C=G/H; S=-F/H; /*00001000*/ DO J=LBP+N TO LBP+NP-1; /*00001010*/ Y=AB(L1,J); Z=AB(I,J); /*00001020*/ AB(L1,J)=C*Y+S*Z; AB(I,J)=-S*Y+C*Z; /*00001030*/ END; /*00001040*/ END; /*00001050*/ TEST_F_CONVERGENCE: /*00001060*/ Z=Q(K); IF L=K THEN GO TO CONVERGENCE; /*00001070*/ /*SHIFT FROM BOTTOM 2X2 MINOR*/ /*00001080*/ X=Q(L); Y=Q(K-1); G=E(K-1); H=E(K); /*00001090*/ F=((Y-Z)*(Y+Z)+(G-H)*(G+H))/(2E0*H*Y); G=SQRT(F*F+1E0); /*00001100*/ IF F<0E0 THEN S=F-G; ELSE S=F+G; /*00001110*/ F=((X-Z)*(X+Z)+H*(Y/S-H))/X; /*00001120*/ /*NEXT QR TRANSFORMATION*/ /*00001130*/ C,S=1E0; /*00001140*/ DO I=L+1 TO K; /*00001150*/ G=E(I); Y=Q(I); H=S*G; G=C*G; /*00001160*/ E(I-1),Z=SQRT(F*F+H*H); C=F/Z; S=H/Z; /*00001170*/ F=X*C+G*S; G=-X*S+G*C; H=Y*S; Y=Y*C; /*00001180*/ DO J=LBN TO LBN+N-1; /*00001190*/ X=AB(J,I-1); Z=AB(J,I); /*00001200*/ AB(J,I-1)=X*C+Z*S; AB(J,I)=-X*S+Z*C; /*00001210*/ END; /*00001220*/ Q(I-1),Z=SQRT(F*F+H*H); /*00001230*/ C=F/Z; S=H/Z; /*00001240*/ F=C*G+S*Y; X=-S*G+C*Y; /*00001250*/ DO J=LBP+N TO LBP+NP-1; /*00001260*/ Y=AB(I-1,J); Z=AB(I,J); /*00001270*/ AB(I-1,J)=C*Y+S*Z; /*00001280*/ AB(I,J)=-S*Y+C*Z; /*00001290*/ END; /*00001300*/ END; /*00001310*/ E(L)=0E0; E(K)=F; Q(K)=X; GO TO TEST_F_SPLITTING; /*00001320*/ CONVERGENCE: /*00001330*/ IF Z<0E0 THEN DO; /*00001340*/ /* Q(K) IS MADE NON-NEGATIVE */ /*00001350*/ Q(K)=-Z; /*00001360*/ AB(*,K)=-AB(*,K); /*00001370*/ END; /*00001380*/ END; /*00001390*/ FREE E; /*00001400*/ END MINFIT; /*00001410*/