code 36401; procedure SYMEIGIMP(N,A,VEC,VAL1,VAL2,LBOUND,UBOUND,AUX); value N;integer N;array A,VEC,VAL1,VAL2,LBOUND,UBOUND,AUX; begin integer K,I,J,I0,I1,ITER,MAXITP1;real S,HEAD,TAIL,MAX,TOL, MATEPS,RELERRA,RELTOLR,NORMA;integer array PERM[1:N]; array R,P,Y[1:N,1:N],RQ,RQT,EPS,Z,VAL3,ETA[1:N]; procedure BOUNDS(I0,I1,N,LBOUND,UBOUND);value I0,I1,N; integer I0,I1,N;array LBOUND,UBOUND; begin integer K,I,J,I01;real EPS2,DL,DR; for I:=I0,I01 while I<=I1 do begin J:=I01:=I; for J:=J+1 while if J>I1 then false else RQ[J]-RQ[J-1]<=EPS[J]+EPS[J-1] do I01:=J; if I = I01 then begin if I<N then begin if I=1 then DL:=DR:=RQ[I+1]-RQ[I]-EPS[I+1] else begin DL:=RQ[I]-RQ[I-1]-EPS[I-1]; DR:=RQ[I+1]-RQ[I]-EPS[I+1] end end else DL:=DR:=RQ[I]-RQ[I-1]-EPS[I-1]; EPS2:=EPS[I]*EPS[I];LBOUND[I]:=EPS2/DR+MATEPS; UBOUND[I]:=EPS2/DL+MATEPS end else begin for K:=I step 1 until I01 do LBOUND[K]:=UBOUND[K]:=EPS[K]+MATEPS end;I01:=I01+1 end end BOUNDS; boolean STOP;STOP:=false;NORMA:=INFNRMMAT(1,N,1,N,J,A); RELERRA:=AUX[0];RELTOLR:=AUX[2];MAXITP1:=AUX[4]+1; MATEPS:=RELERRA*NORMA;TOL:=RELTOLR*NORMA; for ITER:=1 step 1 until MAXITP1 do begin STOP:=true;MAX:=0; for J:=1 step 1 until N do for I:=1 step 1 until N do begin LNGMATVEC(J,J,I,VEC,VAL1,0,0,HEAD,TAIL); LNGMATMAT(1,N,I,J,A,VEC,-HEAD,-TAIL,R[I,J],TAIL); if ABS(R[I,J])>MAX then MAX:=ABS(R[I,J]) end;if MAX > TOL then STOP:=false; if not STOP and ITER<MAXITP1 then begin for I:=1 step 1 until N do LNGTAMMAT(1,N,I,I,VEC,R,VAL1[I],0,RQ[I],RQT[I]); for J:=1 step 1 until N do begin for I:=1 step 1 until N do ETA[I]:=R[I,J]-(RQ[J]-VAL1[J])*VEC[I,J]; Z[J]:=SQRT(VECVEC(1,N,0,ETA,ETA)) end; MERGESORT(RQ,PERM,1,N);VECPERM(PERM,1,N,RQ); for I:=1 step 1 until N do begin EPS[I]:=Z[PERM[I]];VAL3[I]:=VAL1[PERM[I]]; ROWPERM(PERM,1,N,I,VEC);ROWPERM(PERM,1,N,I,R) end; for I:=1 step 1 until N do for J:=I step 1 until N do P[I,J]:=P[J,I]:=TAMMAT(1,N,I,J,VEC,R); end; for I0:=1,I1+1 while I0<=N do begin J:=I1:=I0; for J:=J+1 while if J>N then false else RQ[J]-RQ[J-1]<=SQRT((EPS[J]+EPS[J-1])*NORMA) do I1:=J; if STOP or ITER=MAXITP1 then BOUNDS(I0,I1,N,LBOUND,UBOUND) else begin if I0=I1 then begin for K:=1 step 1 until N do if K=I0 then Y[K,I0]:=1 else R[K,I0]:=P[K,I0]; VAL1[I0]:=RQ[I0];VAL2[I0]:=RQT[PERM[I0]] end else begin integer N1,I0M1,I1P1;real M1;array EM[0:5]; N1:=I1-I0+1;EM[0]:=EM[2]:="-14;EM[4]:=10*N1; begin array PP[1:N1,1:N1],VAL4[1:N1];M1:=0; for K:=I0 step 1 until I1 do M1:=M1+VAL3[K];M1:=M1/N1; for I:=1 step 1 until N1 do for J:=1 step 1 until N1 do begin PP[I,J]:=P[I+I0-1,J+I0-1]; if I=J then PP[I,J]:=PP[I,J]+VAL3[J+I0-1]-M1 end;for I:=I0 step 1 until I1 do begin VAL3[I]:=M1;VAL1[I]:=RQ[I]; VAL2[I]:=RQT[PERM[I]] end; QRISYM(PP,N1,VAL4,EM); MERGESORT(VAL4,PERM,1,N1); for I:=1 step 1 until N1 do for J:=1 step 1 until N1 do P[I+I0-1,J+I0-1]:=PP[I,PERM[J]]; I0M1:=I0-1;I1P1:=I1+1; for J:=I0 step 1 until I1 do begin for I:=1 step 1 until I0M1, I1P1 step 1 until N do begin S:=0; for K:=I0 step 1 until I1 do S:=S+P[I,K]*P[K,J]; R[I,J]:=S end;for I:=I0 step 1 until I1 do Y[I,J]:=P[I,J] end FOR J end INNERBLOCK end I1>I0 end NOT STOP end FOR I0; if not STOP and ITER<MAXITP1 then begin for J:=1 step 1 until N do for I:=1 step 1 until N do if VAL3[I]^=VAL3[J] then Y[I,J]:=R[I,J]/(VAL3[J]-VAL3[I]); for I:=1 step 1 until N do begin for J:=1 step 1 until N do Z[J]:=MATMAT(1,N,I,J,VEC,Y); for J:=1 step 1 until N do VEC[I,J]:=Z[J] end;ORTHOG(N,1,N,VEC) end else begin AUX[5]:=ITER-1;goto EXIT end end FOR ITER; EXIT: AUX[1]:=NORMA;AUX[3]:=MAX end SYMEIGIMP; eop