/* SOLUTION OF REAL AND COMPLEX SYSTEMS OF LINEAR EQUATIONS */ /* From J. H. Wilkinson & C. Reinsch, "Handbook for Automatic Computation", Vol. II, "Linear Algebra", Springer-Verlag, New York, 1971, pp. 98-106. */ /* p. 98. */ INNERPROD: PROCEDURE (K, L, S, U, C, A, B, D); DECLARE (K, L, S, U) FIXED BINARY (31); DECLARE C FLOAT (18); DECLARE (A(*), B(*)) FLOAT (18); DECLARE D FLOAT (18); D = C; DO K = L TO U BY S; D = D + A(K) * B(K); END; END INNERPROD; /* P. 99 */ UNSDET:PROC(N,EPS,A,D1,D2,INT,FAIL); /*PAGE 99*/ /*00000010*/ DCL (D2,N,I,J,K,L,LBN,UBN) BIN FIXED(31, 0), /*00000020*/ (X,Y,YY,D1,EPS) BIN FLOAT(21), /*00000030*/ SUMBOX BIN FLOAT(53), /*00000040*/ FAIL LABEL, /*00000050*/ A(*,*) BIN FLOAT(21), /*00000060*/ INT(*) BIN FIXED(31), /*00000070*/ (ATTEMPT(LBN:UBN) CTL) BIN FLOAT(21); /*00000080*/ LBN=LBOUND(A,1); UBN=LBN+N-1; /*00000090*/ ALLOC ATTEMPT; /*00000100*/ DO I=LBN TO UBN; /*00000110*/ SUMBOX=0E0; /*00000120*/ DO J=LBN TO UBN; /*00000130*/ SUMBOX=SUMBOX+A(I,J)**2; /*00000140*/ END; /*00000150*/ ATTEMPT(I)=1E0/SQRT(SUMBOX); /*00000160*/ END; /*00000170*/ D1=1E0; D2=0; /*00000180*/ DO K=LBN TO UBN; /*00000190*/ L=K; X=0E0; /*00000200*/ DO I=K TO UBN; /*00000210*/ SUMBOX= A(I,K); /*00000220*/ DO J=LBN TO K-1; /*00000230*/ SUMBOX=SUMBOX-A(I,J)*A(J,K); /*00000240*/ END; /*00000250*/ A(I,K)= SUMBOX; Y=ABS(SUMBOX*ATTEMPT(I)); /*00000260*/ IF Y>X /*00000270*/ THEN DO; /*00000280*/ X=Y; L=I; /*00000290*/ END; /*00000300*/ END; /*00000310*/ IF L^=K /*00000320*/ THEN DO; /*00000330*/ D1=-D1; /*00000340*/ DO J=LBN TO UBN; /*00000350*/ Y=A(K,J); A(K,J)=A(L,J); A(L,J)=Y; /*00000360*/ END; /*00000370*/ ATTEMPT(L)=ATTEMPT(K); /*00000380*/ END; /*00000390*/ ATTEMPT(K)=L; D1=D1*A(K,K); /*00000400*/ IF X<8E0*EPS /*00000410*/ THEN GO TO FAIL; /*00000420*/ DO WHILE(ABS(D1)>=1E0); /*00000430*/ D1=D1*.0625E0; D2=D2+4; /*00000440*/ END; /*00000450*/ DO WHILE(ABS(D1)<=.0625); /*00000460*/ D1=D1*16E0; D2=D2-4; /*00000470*/ END; /*00000480*/ X=-1E0/A(K,K); /*00000490*/ DO J=K+1 TO UBN; /*00000500*/ SUMBOX=-A(K,J); /*00000510*/ DO I=LBN TO K-1; /*00000520*/ SUMBOX =SUMBOX+A(K,I)*A(I,J); /*00000530*/ END; /*00000540*/ A(K,J)=X*SUMBOX; /*00000550*/ END; /*00000560*/ END; /*00000570*/ INT=FLOOR(ATTEMPT+.5); /*00000580*/ FREE ATTEMPT; /*00000590*/ END UNSDET; /*00000600*/ /* P. 100 */ UNSYM_ACC_SOLVE: PROC(N,R,A,AA,P,B,EPS,X,BB,L,ILL) OPTIONS (REORDER); /*PAGE 100*/ /*00000010*/ DCL (N,R,L,I,J,K,D2,LBN,UBN,LBR,UBR) BIN FIXED(31,0), /*00000020*/ ILL LABEL, /*00000030*/ (EPS,D0,D1,XMAX,BBMAX) BIN FLOAT(21), /*00000040*/ UNSSOL ENTRY, /*00000050*/ C BIN FLOAT(53), /*00000060*/ (A,AA)(*,*) BIN FLOAT(21), /*00000070*/ P(*) BIN FIXED(31,0), /*00000080*/ (X,B,BB)(*,*) BIN FLOAT(21); /*00000090*/ LBN=LBOUND(A,1); UBN=LBN+N-1; /*00000100*/ LBR=LBOUND(X,2); UBR=LBR+R-1; /*00000110*/ X=0E0; BB=B; L=0; D0=0E0; /*00000120*/ L3:CALL UNSYMSOL(N,R,AA,P,BB); /*00000130*/ L=L+1; D2=0; D1=0E0; /*00000140*/ X=X+BB; /*00000150*/ DO J=LBR TO UBR; /*00000160*/ XMAX,BBMAX=0E0; /*00000170*/ DO I=LBN TO UBN; /*00000180*/ IF ABS(X(I,J))>XMAX THEN XMAX=ABS(X(I,J)); /*00000190*/ IF ABS(BB(I,J))>BBMAX THEN BBMAX=ABS(BB(I,J)); /*00000200*/ C= B(I,J); /*00000210*/ DO K=LBN TO UBN; /*00000220*/ C=C-A(I,K)*X(K,J); /*00000230*/ END; /*00000240*/ BB(I,J)= C; /*00000250*/ END; /*00000260*/ IF BBMAX>D1*XMAX THEN D1=BBMAX/XMAX; /*00000270*/ IF BBMAX>2E0*EPS*XMAX THEN D2=1; /*00000280*/ END; /*00000290*/ IF D1>D0*.5&L^=1 THEN GO TO ILL; /*00000300*/ D0=D1; /*00000310*/ IF D2=1 THEN GO TO L3; /*00000320*/ END UNSYM_ACC_SOLVE; /*00000330*/ UNSYMSOL:PROC(N,R,A,INT,B) OPTIONS (REORDER); /*PAGE 100*/ /*00000010*/ DCL (N,R,I,J,K,LBN,UBN,LBR,UBR) BIN FIXED(31,0), /*00000020*/ SUMBOX BIN FLOAT(53), /*00000030*/ X BIN FLOAT(21), /*00000040*/ A(*,*) BIN FLOAT(21), /*00000050*/ B(*,*) BIN FLOAT(21), /*00000060*/ INT(*) BIN FIXED(31,0); /*00000070*/ LBN=LBOUND(A,1); UBN=LBN+N-1; /*00000080*/ LBR=LBOUND(B,2); UBR=LBR+R-1; /*00000090*/ /*INTERCHANGING OF ELEMENTS OF B */ /*00000100*/ DO I=LBN TO UBN; /*00000110*/ IF INT(I)^=I /*00000120*/ THEN DO K=LBR TO UBR; /*00000130*/ X=B(I,K); B(I,K)=B(INT(I),K); /*00000140*/ B(INT(I),K)=X; /*00000150*/ END; /*00000160*/ END; /*00000170*/ DO K=LBR TO UBR; /*00000180*/ /*SOLUTION OF LY=B */ /*00000190*/ DO I=LBN TO UBN; /*00000200*/ SUMBOX=B(I,K); /*00000210*/ DO J=LBN TO I-1; /*00000220*/ SUMBOX=SUMBOX+A(I,J)*B(J,K); /*00000230*/ END; /*00000240*/ B(I,K)=-SUMBOX/A(I,I); /*00000250*/ END; /*00000260*/ /*SOLUTION OF UX=Y */ /*00000270*/ DO I=UBN TO LBN BY -1; /*00000280*/ SUMBOX=B(I,K); /*00000290*/ DO J=I+1 TO UBN; /*00000300*/ SUMBOX=SUMBOX+A(I,J)*B(J,K); /*00000310*/ END; /*00000320*/ B(I,K)=-SUMBOX; /*00000330*/ END; /*00000340*/ END; /*00000350*/ END UNSYMSOL; /*00000360*/