trlinsolve.pas

maths lib with source

                end;
            end
            else
            begin
                
                //
                // A is unit triangular.
                //
                // Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
                //
                GROW := Min(1, 1/Max(XBND, SMLNUM));
                J := JFIRST;
                while (JINC>0) and (J<=JLAST) or (JINC<0) and (J>=JLAST) do
                begin
                    
                    //
                    // Exit the loop if the growth factor is too small.
                    //
                    if GROW<=SMLNUM then
                    begin
                        Break;
                    end;
                    
                    //
                    // G(j) = G(j-1)*( 1 + CNORM(j) )
                    //
                    GROW := GROW*(1/(1+CNORM[J]));
                    J := J+JINC;
                end;
            end;
        end;
    end
    else
    begin
        
        //
        // Compute the growth in A' * x = b.
        //
        if UPPER then
        begin
            JFIRST := 1;
            JLAST := N;
            JINC := 1;
        end
        else
        begin
            JFIRST := N;
            JLAST := 1;
            JINC := -1;
        end;
        if TSCAL<>1 then
        begin
            GROW := 0;
        end
        else
        begin
            if NOUNIT then
            begin
                
                //
                // A is non-unit triangular.
                //
                // Compute GROW = 1/G(j) and XBND = 1/M(j).
                // Initially, M(0) = max{x(i), i=1,...,n}.
                //
                GROW := 1/Max(XBND, SMLNUM);
                XBND := GROW;
                J := JFIRST;
                while (JINC>0) and (J<=JLAST) or (JINC<0) and (J>=JLAST) do
                begin
                    
                    //
                    // Exit the loop if the growth factor is too small.
                    //
                    if GROW<=SMLNUM then
                    begin
                        Break;
                    end;
                    
                    //
                    // G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) )
                    //
                    XJ := 1+CNORM[J];
                    GROW := Min(GROW, XBND/XJ);
                    
                    //
                    // M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j))
                    //
                    TJJ := ABSReal(A[J,J]);
                    if XJ>TJJ then
                    begin
                        XBND := XBND*(TJJ/XJ);
                    end;
                    if J=JLAST then
                    begin
                        GROW := Min(GROW, XBND);
                    end;
                    J := J+JINC;
                end;
            end
            else
            begin
                
                //
                // A is unit triangular.
                //
                // Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.
                //
                GROW := Min(1, 1/Max(XBND, SMLNUM));
                J := JFIRST;
                while (JINC>0) and (J<=JLAST) or (JINC<0) and (J>=JLAST) do
                begin
                    
                    //
                    // Exit the loop if the growth factor is too small.
                    //
                    if GROW<=SMLNUM then
                    begin
                        Break;
                    end;
                    
                    //
                    // G(j) = ( 1 + CNORM(j) )*G(j-1)
                    //
                    XJ := 1+CNORM[J];
                    GROW := GROW/XJ;
                    J := J+JINC;
                end;
            end;
        end;
    end;
    if GROW*TSCAL>SMLNUM then
    begin
        
        //
        // Use the Level 2 BLAS solve if the reciprocal of the bound on
        // elements of X is not too small.
        //
        if UPPER and NOTRAN or  not UPPER and  not NOTRAN then
        begin
            if NOUNIT then
            begin
                VD := A[N,N];
            end
            else
            begin
                VD := 1;
            end;
            X[N] := X[N]/VD;
            I:=N-1;
            while I>=1 do
            begin
                IP1 := I+1;
                if Upper then
                begin
                    V := APVDotProduct(@A[I][0], IP1, N, @X[0], IP1, N);
                end
                else
                begin
                    V := 0.0;
                    for i_ := IP1 to N do
                    begin
                        V := V + A[i_,I]*X[i_];
                    end;
                end;
                if NOUNIT then
                begin
                    VD := A[I,I];
                end
                else
                begin
                    VD := 1;
                end;
                X[I] := (X[I]-V)/VD;
                Dec(I);
            end;
        end
        else
        begin
            if NOUNIT then
            begin
                VD := A[1,1];
            end
            else
            begin
                VD := 1;
            end;
            X[1] := X[1]/VD;
            I:=2;
            while I<=N do
            begin
                IM1 := I-1;
                if Upper then
                begin
                    V := 0.0;
                    for i_ := 1 to IM1 do
                    begin
                        V := V + A[i_,I]*X[i_];
                    end;
                end
                else
                begin
                    V := APVDotProduct(@A[I][0], 1, IM1, @X[0], 1, IM1);
                end;
                if NOUNIT then
                begin
                    VD := A[I,I];
                end
                else
                begin
                    VD := 1;
                end;
                X[I] := (X[I]-V)/VD;
                Inc(I);
            end;
        end;
    end
    else
    begin
        
        //
        // Use a Level 1 BLAS solve, scaling intermediate results.
        //
        if XMAX>BIGNUM then
        begin
            
            //
            // Scale X so that its components are less than or equal to
            // BIGNUM in absolute value.
            //
            S := BIGNUM/XMAX;
            APVMul(@X[0], 1, N, S);
            XMAX := BIGNUM;
        end;
        if NOTRAN then
        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) / A(j,j), scaling x if necessary.
                //
                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 := 100;
                    end;
                end;
                if Flg<>100 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/b(j).
                                //
                                REC := 1/XJ;
                                APVMul(@X[0], 1, N, REC);
                                S := S*REC;
                                XMAX := XMAX*REC;
                            end;
                        end;
                        X[J] := X[J]/TJJS;
                        XJ := ABSReal(X[J]);
                    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
                                // to avoid overflow when dividing by A(j,j).
                                //
                                REC := TJJ*BIGNUM/XJ;
                                if CNORM[J]>1 then
                                begin
                                    
                                    //
                                    // Scale by 1/CNORM(j) to avoid overflow when
                                    // multiplying x(j) times column j.
                                    //
                                    REC := REC/CNORM[J];
                                end;
                                APVMul(@X[0], 1, N, REC);
                                S := S*REC;
                                XMAX := XMAX*REC;
                            end;
                            X[J] := X[J]/TJJS;
                            XJ := ABSReal(X[J]);
                        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;