code 34138;
procedure LSQREFSOL(A, QR, N, M, N1, AUX, AID, CI, B,LDX,X,RES);
value N,M,N1;integer N,M,N1;integer array CI;real LDX;
array QR,A,AID,AUX,X,RES,B;
begin integer I,J,K,S;
       real C1,NEXVE,NDX,NDR,D,DD,OP,OPL,CORRNORM;
       array F[1:N],G[1:M];
       procedure HOUSEHOLDER(P, Q, R, E);
       value P,Q,R,E;integer P,Q,R,E;
       begin for S:=P step Q until R do 
               ELMVECCOL(S, E, S, F, QR,TAMVEC(S, E, S,QR, F)/(QR[S,S]*
               AID[S]))
       end;
       for J:=1 step 1 until M do 
       begin S:=CI[J];if S^=J"THEN" ICHCOL(1,N,J,S,A) end;
       for J:=1step 1until M"DO"X[J]:=G[J]:=0;
       for I:=1step 1until N"DO"
       begin RES[I]:=0;F[I]:=B[I]end;
       for K:=0,1,K+1
          while (CORRNORM>AUX[2]*NEXVE & K<=AUX[6])
       do 
       begin NDX:=NDR:=0;
          if K^=0then 
          begin for I:=1step 1until N"DO"RES[I]:=RES[I]+F[I];
                 for S:=1step 1until M"DO"
                 begin X[S]:=X[S]+G[S];
                        LNGTAMVEC(1,N,S,A,RES,0,0,D,DD);
                        G[S]:=(-D-TAMVEC(1,S-1,S,QR,G))/AID[S]
                 end;
                 for I:=1step 1until N"DO"
                 begin LNGMATVEC(1, M, I, A, X,
                         if I>N1 then RES[I] else 0, 0, D, DD);
                        LNG SUB(B[I],0,D,DD,OP,OPL);
                        F[I]:=OP
                 end 
          end;
          NEXVE:=SQRT(VECVEC(1,M,0,X,X)+VECVEC(1,N,0,RES,RES));
          HOUSEHOLDER(1, 1, N1, N1);
          for I:=N1+1 step 1 until N do 
          F[I]:=F[I]-MATVEC(1,N1, I, QR, F);
          HOUSEHOLDER(N1+1, 1, M, N);
          for I:=1 step 1 until M do 
          begin C1:=F[I];F[I]:=G[I];
              G[I]:=if I>N1 then C1-G[I] else C1
          end;
          for S:=M"STEP"-1until 1do 
          begin G[S]:=(G[S]-MATVEC(S+1,M,S,QR,G))/AID[S];
                 NDX:=NDX+G[S]**2
          end;
          HOUSEHOLDER(M, -1, N1+1, N);
          for S:=1 step 1 until N1 do 
          F[S]:=F[S]-TAMVEC(N1+1, N, S, QR, F);
          HOUSEHOLDER(N1, -1, 1, N1);
          AUX[7]:=K;
          for I:=1step 1until N"DO"NDR:=NDR+F[I]**2;
          CORRNORM:=SQRT(NDX+NDR)
       end FOR K;
       LDX:=SQRT(NDX);
       for S:=M step -1 until 1 do 
       begin J:=CI[S];if J^=S then 
          begin C1:=X[J];X[J]:=X[S];X[S]:=C1;ICHCOL(1,N,J,S,A)end 
       end 
end LSQREFSOL;
        eop