DECLARE SUB rowk (n%, k%, r#()) DECLARE SUB dple (n%, b#(), c#(), ierr%) '*********************************************************************** '* Solve AX = B using a partial pivoting algorithm and reduced storage * '* ------------------------------------------------------------------- * '* SAMPLE RUN: * '* * '* System to solve: * '* 2.0000 -1.0000 1.0000 7.0000 -12.5400 5.0000 * '* 1.0000 5.0000 -2.0000 -8.0000 100.0000 1.0000 * '* 3.0000 -2.0000 3.0000 45.0000 27.3333 3.0000 * '* 11.0000 0.5500 -2.0000 -4.0000 1.0000 4.0000 * '* 33.0000 2.0000 -3.0000 5.0000 7.3333 -10.0000 * '* * '* Solution is: * '* 2.111496 * '* -25.829027 * '* 8.174232 * '* -2.521467 * '* 1.242104 * '* * '* ------------------------------------------------------------------- * '* Ref.: "Wassyng, A. - Solving Ax = b: A method with reduced storage * '* requirements, SIAM J. Numerical Analysis, vol.19 (1982), * '* pp. 197-204". * '* * '* Basic Release By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '*********************************************************************** 'PROGRAM Test_DPLE DEFDBL A-H, O-Z DEFINT I-N DIM a(25), b(25), soln(25) DIM row AS INTEGER ' Define the right-hand side (n=5) b(1) = 5#: b(2) = 1#: b(3) = 3#: b(4) = 4#: b(5) = -10# CLS PRINT PRINT " System to solve:" FOR i = 1 TO 5 rowk 5, i, soln() FOR j = 1 TO 5 PRINT USING " ###.####"; soln(j); NEXT j PRINT USING " ###.####"; b(i) NEXT i PRINT dple 5, b(), soln(), ierr IF ierr = 0 THEN PRINT " Solution is:" FOR i = 1 TO 5 PRINT USING "####.######"; soln(i) NEXT i ELSE PRINT " Error = "; ierr END IF END 'of main program SUB dple (n, b(), c(), ierr) '****************************************************************** ' SOLUTION OF LINEAR EQUATIONS WITH REDUCED STORAGE '****************************************************************** ' Uses the Henderson-Wassyng partial pivot algorithm. ' Wassyng, A. "Solving Ax = b: A method with reduced storage requirements", ' SIAM J. Numerical Analysis, vol.19 (1982), pp. 197-204. ' The user must provide a routine ROWK to return the requested row of the ' matrix A. zero = 0# DIM wk(n * n / 4 + n + 3) DIM iwk(25) ' SET THE NECESSARY CONSTANTS ierr = 0 maxwk = n * n / 4 + n + 3 np1 = n + 1 k = 1 iflag = -1 ' GET THE FIRST COLUMN OF THE TRANSPOSED SYSTEM rowk n, 1, c() bk = b(1) IF n <= 1 THEN IF c(1) = zero THEN GOTO 130 c(1) = bk / c(1) EXIT SUB END IF ' FIND THE PIVOT FOR COLUMN 1 m = 1 FOR i = 2 TO n IF ABS(c(m)) < ABS(c(i)) THEN m = i NEXT i iwk(1) = m c1 = c(m) c(m) = c(1) c(1) = c1 IF c(1) <> zero THEN ' FIND THE FIRST ELEMENTARY MATRIX AND STORE IT IN D FOR i = 2 TO n wk(i - 1) = -c(i) / c(1) NEXT i wk(n) = bk / c(1) ' K LOOP - EACH K FOR A NEW COLUMN OF THE TRANSPOSED SYSTEM FOR k = 2 TO n kp1 = k + 1 km1 = k - 1 ' GET COLUMN K rowk n, k, c() FOR j = 1 TO km1 m = iwk(j) cj = c(j) c(j) = c(m) c(m) = cj NEXT j bk = b(k) iflag = -iflag lcol = np1 - k lcolp1 = lcol + 1 lastm1 = 1 last = maxwk - n + k IF k <> 2 THEN lastm1 = maxwk - n + km1 IF iflag < 0 THEN last = last - n + k - 2 IF iflag > 0 THEN lastm1 = lastm1 - n + k - 3 END IF ' J LOOP - EFFECT OF COLUMNS 1 TO K-1 OF L-INVERSE FOR j = 1 TO km1 cj = c(j) ij = (j - 1) * lcolp1 IF j = km1 THEN ij = lastm1 - 1 ' I LOOP - EFFECT OF L-INVERSE ON ROWS K TO N+1 FOR i = k TO n ij = ij + 1 c(i) = c(i) + wk(ij) * cj NEXT i bk = bk - wk(ij + 1) * cj NEXT j ' K=N CASE m = k IF k >= n THEN IF c(k) = zero THEN GOTO 130 wk(last) = bk / c(k) ELSE ' FIND THE PIVOT FOR i = kp1 TO n IF ABS(c(m)) < ABS(c(i)) THEN m = i NEXT i iwk(k) = m ck = c(m) c(m) = c(k) c(k) = ck IF c(k) = zero THEN GOTO 130 ' FIND THE K-TH ELEMENTARY MATRIX ik = last FOR i = kp1 TO n wk(ik) = -c(i) / c(k) ik = ik + 1 NEXT i wk(ik) = bk / c(k) END IF ' FORM THE PRODUCT OF THE ELEMENTARY MATRICES FOR j = 1 TO km1 kjold = j * lcolp1 + k - np1 mjold = kjold + m - k ij = (j - 1) * lcol ijold = ij + j IF j = km1 THEN kjold = lastm1 mjold = lastm1 + m - k ijold = lastm1 END IF ik = last - 1 dkj = wk(mjold) wk(mjold) = wk(kjold) FOR i = kp1 TO np1 ij = ij + 1 ijold = ijold + 1 ik = ik + 1 wk(ij) = wk(ijold) + wk(ik) * dkj NEXT i NEXT j NEXT k last = maxwk IF iflag < 0 THEN last = maxwk - 2 wk(n) = wk(last) ' INSERT THE SOLUTION IN C FOR i = 1 TO n c(i) = wk(i) NEXT i nm1 = n - 1 FOR i = 1 TO nm1 k = n - i m = iwk(k) ck = c(k) c(k) = c(m) c(m) = ck NEXT i EXIT SUB END IF ' THE SYSTEM IS SINGULAR 130 ierr = k END SUB 'dple SUB rowk (n, k, r()) IF k = 1 THEN r(1) = 2#: r(2) = -1#: r(3) = 1#: r(4) = 7#: r(5) = -12.54 ELSEIF k = 2 THEN r(1) = 1#: r(2) = 5#: r(3) = -2#: r(4) = -8#: r(5) = 100# ELSEIF k = 3 THEN r(1) = 3#: r(2) = -2#: r(3) = 3#: r(4) = 45#: r(5) = 27.3333 ELSEIF k = 4 THEN r(1) = 11#: r(2) = .55: r(3) = -2#: r(4) = -4#: r(5) = 1# ELSEIF k = 5 THEN r(1) = 33#: r(2) = 2#: r(3) = -3#: r(4) = 5#: r(5) = 7.3333 END IF END SUB 'rowk