trlinsolve.pas

maths lib with source

                            XJ := 1;
                            S := 0;
                            XMAX := 0;
                        end;
                    end;
                end;
                
                //
                // Scale x if necessary to avoid overflow when adding a
                // multiple of column j of A.
                //
                if XJ>1 then
                begin
                    REC := 1/XJ;
                    if CNORM[J]>(BIGNUM-XMAX)*REC then
                    begin
                        
                        //
                        // Scale x by 1/(2*abs(x(j))).
                        //
                        REC := REC*0.5;
                        APVMul(@X[0], 1, N, REC);
                        S := S*REC;
                    end;
                end
                else
                begin
                    if XJ*CNORM[J]>BIGNUM-XMAX then
                    begin
                        
                        //
                        // Scale x by 1/2.
                        //
                        APVMul(@X[0], 1, N, 0.5);
                        S := S*0.5;
                    end;
                end;
                if UPPER then
                begin
                    if J>1 then
                    begin
                        
                        //
                        // Compute the update
                        // x(1:j-1) := x(1:j-1) - x(j) * A(1:j-1,j)
                        //
                        V := X[J]*TSCAL;
                        JM1 := J-1;
                        for i_ := 1 to JM1 do
                        begin
                            X[i_] := X[i_] - V*A[i_,J];
                        end;
                        I := 1;
                        K:=2;
                        while K<=J-1 do
                        begin
                            if AbsReal(X[K])>AbsReal(X[I]) then
                            begin
                                I := K;
                            end;
                            Inc(K);
                        end;
                        XMAX := ABSReal(X[I]);
                    end;
                end
                else
                begin
                    if J<N then
                    begin
                        
                        //
                        // Compute the update
                        // x(j+1:n) := x(j+1:n) - x(j) * A(j+1:n,j)
                        //
                        JP1 := J+1;
                        V := X[J]*TSCAL;
                        for i_ := JP1 to N do
                        begin
                            X[i_] := X[i_] - V*A[i_,J];
                        end;
                        I := J+1;
                        K:=J+2;
                        while K<=N do
                        begin
                            if AbsReal(X[K])>AbsReal(X[I]) then
                            begin
                                I := K;
                            end;
                            Inc(K);
                        end;
                        XMAX := ABSReal(X[I]);
                    end;
                end;
                J := J+JINC;
            end;
        end
        else
        begin
            
            //
            // Solve A' * x = b
            //
            J := JFIRST;
            while (JINC>0) and (J<=JLAST) or (JINC<0) and (J>=JLAST) do
            begin
                
                //
                // Compute x(j) = b(j) - sum A(k,j)*x(k).
                //   k<>j
                //
                XJ := ABSReal(X[J]);
                USCAL := TSCAL;
                REC := 1/Max(XMAX, 1);
                if CNORM[J]>(BIGNUM-XJ)*REC then
                begin
                    
                    //
                    // If x(j) could overflow, scale x by 1/(2*XMAX).
                    //
                    REC := REC*0.5;
                    if NOUNIT then
                    begin
                        TJJS := A[J,J]*TSCAL;
                    end
                    else
                    begin
                        TJJS := TSCAL;
                    end;
                    TJJ := ABSReal(TJJS);
                    if TJJ>1 then
                    begin
                        
                        //
                        // Divide by A(j,j) when scaling x if A(j,j) > 1.
                        //
                        REC := Min(1, REC*TJJ);
                        USCAL := USCAL/TJJS;
                    end;
                    if REC<1 then
                    begin
                        APVMul(@X[0], 1, N, REC);
                        S := S*REC;
                        XMAX := XMAX*REC;
                    end;
                end;
                SUMJ := 0;
                if USCAL=1 then
                begin
                    
                    //
                    // If the scaling needed for A in the dot product is 1,
                    // call DDOT to perform the dot product.
                    //
                    if UPPER then
                    begin
                        if J>1 then
                        begin
                            JM1 := J-1;
                            SUMJ := 0.0;
                            for i_ := 1 to JM1 do
                            begin
                                SUMJ := SUMJ + A[i_,J]*X[i_];
                            end;
                        end
                        else
                        begin
                            SUMJ := 0;
                        end;
                    end
                    else
                    begin
                        if J<N then
                        begin
                            JP1 := J+1;
                            SUMJ := 0.0;
                            for i_ := JP1 to N do
                            begin
                                SUMJ := SUMJ + A[i_,J]*X[i_];
                            end;
                        end;
                    end;
                end
                else
                begin
                    
                    //
                    // Otherwise, use in-line code for the dot product.
                    //
                    if UPPER then
                    begin
                        I:=1;
                        while I<=J-1 do
                        begin
                            V := A[I,J]*USCAL;
                            SUMJ := SUMJ+V*X[I];
                            Inc(I);
                        end;
                    end
                    else
                    begin
                        if J<N then
                        begin
                            I:=J+1;
                            while I<=N do
                            begin
                                V := A[I,J]*USCAL;
                                SUMJ := SUMJ+V*X[I];
                                Inc(I);
                            end;
                        end;
                    end;
                end;
                if USCAL=TSCAL then
                begin
                    
                    //
                    // Compute x(j) := ( x(j) - sumj ) / A(j,j) if 1/A(j,j)
                    // was not used to scale the dotproduct.
                    //
                    X[J] := X[J]-SUMJ;
                    XJ := ABSReal(X[J]);
                    Flg := 0;
                    if NOUNIT then
                    begin
                        TJJS := A[J,J]*TSCAL;
                    end
                    else
                    begin
                        TJJS := TSCAL;
                        if TSCAL=1 then
                        begin
                            Flg := 150;
                        end;
                    end;
                    
                    //
                    // Compute x(j) = x(j) / A(j,j), scaling if necessary.
                    //
                    if Flg<>150 then
                    begin
                        TJJ := ABSReal(TJJS);
                        if TJJ>SMLNUM then
                        begin
                            
                            //
                            // abs(A(j,j)) > SMLNUM:
                            //
                            if TJJ<1 then
                            begin
                                if XJ>TJJ*BIGNUM then
                                begin
                                    
                                    //
                                    // Scale X by 1/abs(x(j)).
                                    //
                                    REC := 1/XJ;
                                    APVMul(@X[0], 1, N, REC);
                                    S := S*REC;
                                    XMAX := XMAX*REC;
                                end;
                            end;
                            X[J] := X[J]/TJJS;
                        end
                        else
                        begin
                            if TJJ>0 then
                            begin
                                
                                //
                                // 0 < abs(A(j,j)) <= SMLNUM:
                                //
                                if XJ>TJJ*BIGNUM then
                                begin
                                    
                                    //
                                    // Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.
                                    //
                                    REC := TJJ*BIGNUM/XJ;
                                    APVMul(@X[0], 1, N, REC);
                                    S := S*REC;
                                    XMAX := XMAX*REC;
                                end;
                                X[J] := X[J]/TJJS;
                            end
                            else
                            begin
                                
                                //
                                // A(j,j) = 0:  Set x(1:n) = 0, x(j) = 1, and
                                // scale = 0, and compute a solution to A'*x = 0.
                                //
                                I:=1;
                                while I<=N do
                                begin
                                    X[I] := 0;
                                    Inc(I);
                                end;
                                X[J] := 1;
                                S := 0;
                                XMAX := 0;
                            end;
                        end;
                    end;
                end
                else
                begin
                    
                    //
                    // Compute x(j) := x(j) / A(j,j)  - sumj if the dot
                    // product has already been divided by 1/A(j,j).
                    //
                    X[J] := X[J]/TJJS-SUMJ;
                end;
                XMAX := Max(XMAX, ABSReal(X[J]));
                J := J+JINC;
            end;
        end;
        S := S/TSCAL;
    end;
    
    //
    // Scale the column norms by 1/TSCAL for return.
    //
    if TSCAL<>1 then
    begin
        V := 1/TSCAL;
        APVMul(@CNORM[0], 1, N, V);
    end;
end;


end.