The attached file, in which I play around with arrays, using them to solve
systems of linear equations, will help. It includes a really good function
by Andy Martin.
Russell [EMAIL PROTECTED]
----- Original Message -----
From: <[EMAIL PROTECTED]>
To: <[EMAIL PROTECTED]>
Sent: Thursday, February 03, 2000 6:26 AM
Subject: [REBOL] Multidimensional arrays ...
> Could anyone elaborate on following, please?
>
> ->> ar: array/initial [3 4] copy ""
> == [["" "" "" ""] ["" "" "" ""] ["" "" "" ""]]
> ->> change ar/1/1 "adfaf"
> == ""
> ->> ar
> == [["adfaf" "adfaf" "adfaf" "adfaf"] ["adfaf" "adfaf" "adfaf" "adfaf"]
> ["adfaf" "adfaf" "adfaf" "adfaf"]]
>
> How to simply write at any multidimensional array position? ar/:i/:j:
> value doesn't work ...
>
> Thanks,
>
> -pekr-
>
>
REBOL []
comment { Includes a demonstration of solving n linear equations for n
unknowns of the form
a11*x1 + a12*x2 + ... a1n*xn = b1
a21*x1 + a22*x2 + ... a2n*x2 = b2
...
an1*x1 + an2*x2 + ... ann*xn = bn
The solution approach is the congruential-gradient method
described in 20.3-2(e) of
Korn & Korn, "mathematical Handbook for Scientists and
Engineers", McGraw -Hill, NY NY 1961.
To run the demo, type "test-nxn-its" without the ""'s. In
response to subsequent questions, enter the number of unknowns
and the number of iterations. Except for the effects of
round-off, n iterations give exact
solutions. You'll see this for n's of 10 or less. For larger
n's, such as 20, more
iterations improve the results.
In the demo, 9 digit random coefficients, aij, are generated in
the range
0..1.0, and the same for the "right sides, bj's. These aij's
are converted to a symmetric
(aij = aji) positive definite set and the b's are similarly
transformed, to yeild the same
solutions xj's as for the original set of aij's.
The demo displays the transformed coefficients, the transformed
"right sides",
the "residuals" (the amounts by which the left sides fail to
equal the right sides,
using the solution xj's), and finally, the solution xj's to 9
decimal places.
Try n = 20 with 20, 30, 60, and 400 iterations. You'll see the
effect of round-off
errors and how more iterations reduce them. The 9 place results
don't change between
60 and 400 iterations, but the residuals are smaller, so the
solutions must be improving
in the non-displayed digits!
}
comment { Arrays will be taken to be generated by
z: array [r c]
which results in a block containing r "row" blocks, each conatining
c column items. Array is a REBOL
function. The first argument block element is the number of rows,
the second the number of
columns.
}
inlzarrayint: func [" Initializes the array with integers 1, 2, ..." a
/local x
][
x: 0
foreach row a [
for col 1 length? row 1 [
poke row col ( x: + x 1)
]
]
return a
]
pokearyrow: func [a r blk /local ar]
[ar: copy/deep a
insert/only remove skip ar (subtract r 1) blk
return ar
]
display: func [a][
foreach row a [print row]
exit
]
pokeij: func [a i j n][
; print ["i,j= " i j]
poke (pick a i) j n
]
comment { The following definition is a broader application, where the
path is determined
by the values in ixlist, which may be integers or path names. It
was designed by
Andrew Martin.
It is not used herein.
}
arPoke: func [ar ixlist x /local temp comd][
temp: :ar
comd: "temp"
foreach i ixlist [
either integer? i [comd: rejoin [comd "/" i]]
[comd: rejoin [comd "/" get i]]
]
do bind load rejoin [comd ": " x ] 'comd
]
Transpose: func [{ Transposes the rectangular MATRIX, a,} a
/local r trow row t s][
t: array [(length? first a) (length? a)]
r: 0
foreach row a [
r: + r 1
for c 1 length? row 1 [
pokeij t c r a/:r/:c
]
]
return t
]
comment {The following three functions are used to create a matrix equal
to the matrix product of Atranspose and A. But they are erroneous.
They are not used further herein; they are commented out by being
part of this comment.
AAhi: function [a h i ] [s][
s: 0
for k 1 length? a 1 [
s: + s (* a/:h/:k a/:k/:i )
]
]
AArow: function [a h ][ row ][
row: make block! []
for i 1 length? a 1 [
append row AAhi a h i
]
return row
]
AA: function [a][ blk][
blk: make block! []
for h 1 length? a 1 [
insert/only tail blk ( AArow a h)
]
return head blk
]
} ; end of comment
inlzarrand: func [{Initializes the rectangular array a with random
"integers" from 1 to n, of the same type as n.}a n][
random/seed now
foreach row a [
for col 1 length? row 1 [
poke row col (random n)
]
]
return a
]
inlzarrand01: func [a n][
print [ type? n]
if = (type? n) integer! [n: to-decimal n]
random/seed now
foreach row a [
for col 1 length? row 1 [
poke row col (/ random n n)
]
]
return a
]
sumdotprod: func [a b /local dpv sum][
dpv: make block! []
sum: 0
for i 1 length? a 1 [
append dpv (* (pick a i) (pick b i))
]
foreach elt dpv [
sum: + sum elt
]
return sum
]
matprod: func [a b /local ca ra cb rb p s ][
ra: length? a
rb: ca: length? first a
comment { ca = no cols of a and also equals num rows of b}
cb: length? first b
p: array reduce [ra cb]
for i 1 ra 1[
for j 1 cb 1 [
s: 0
for k 1 ca 1 [
s: + s (* a/:i/:k b/:k/:j)
]
pokeij p i j s
]
]
return p
]
det3: func [a /local d ][
d: 0
d: + d (* a/1/1 ((* a/2/2 a/3/3) - (* a/2/3 a/3/2)))
d: d - (* a/1/2 ((* a/2/1 a/3/3) - (* a/3/1 a/2/3)))
d: + d (* a/1/3 ((* a/2/1 a/3/2) - (* a/3/1 a/2/2)))
return d
]
comment { Xi = / det3 cramer a b i det3 a }
cramer: func [a b i /local c][
comment { a is nxn matrix of coeffs - b is nx1 matrix of right sides}
c: copy/deep a
c: transpose (pokearyrow (transpose c) i (first transpose b))
return c
]
matsum: func [a b /local s ][
s: copy/deep a
for i 1 length? a 1[
for j 1 length? first a [
pokeij s i j (+ a/:i/:j b/:i/:j)
]
]
return s
]
matsub: func [a b /local s ][
s: copy/deep a
for i 1 length? a 1[
for j 1 length? first a [
pokeij s i j ( a/:i/:j - b/:i/:j)
]
]
return s
]
comment { In the next function, a is nxn symmetric pos definite coeff
matrix for system of linear equations.
b is a nx1 col matrix of "right sides". On output x is nx1 col
matrix of solutions.
}
cgsoln: func [a b /local x v f n den num1 num2 nod1 nod2 saikxk saikvk]
[
n: length? a
f: array reduce [n 1]
for i 1 n 1 [
pokeij f i 1 (subtract 0 b/:i/1)
]
x: array/initial reduce [n 1] 0.0
v: copy/deep f
for j 1 (+ n 1) 1 [ ; jth iteration
den: 0
for k 1 n 1[
for h 1 n 1[
den: + den (* a/:k/:h (* v/:k/1 v/:h/1))
]
]
num1: 0
for k 1 n 1[
num1: + num1 (* f/:k/1 v/:k/1)
]
nod1: / num1 den
for i 1 n 1 [
pokeij x i 1 (subtract x/:i/1 (* nod1 v/:i/1))
]
for i 1 n 1 [
saikvk: 0
for k 1 n 1[
saikvk: + saikvk (* a/:i/:k v/:k/1)
]
pokeij f i 1 (subtract f/:i/1 (* nod1 saikvk))
]
num2: 0
for k 1 n 1[
for h 1 n 1[
num2: + num2 (* a/:k/:h (* f/:k/1 v/:h/1))
]
]
nod2: / num2 den
for i 1 n 1 [
pokeij v i 1 (subtract f/:i/1 (* nod2 v/:i/1))
]
]
display f
return x
]
cgsoln-its: func [a b m /local x v f n den num1 num2 nod1 nod2 saikxk
saikvk]
[
n: length? a
f: array reduce [n 1]
for i 1 n 1 [
pokeij f i 1 (subtract 0 b/:i/1)
]
x: array/initial reduce [n 1] 0.0
v: copy/deep f
for j 1 m 1 [ ; jth iteration
den: 0
for k 1 n 1[
for h 1 n 1[
den: + den (* a/:k/:h (* v/:k/1 v/:h/1))
]
]
num1: 0
for k 1 n 1[
num1: + num1 (* f/:k/1 v/:k/1)
]
nod1: / num1 den
for i 1 n 1 [
pokeij x i 1 (subtract x/:i/1 (* nod1 v/:i/1))
]
for i 1 n 1 [
saikvk: 0
for k 1 n 1[
saikvk: + saikvk (* a/:i/:k v/:k/1)
]
pokeij f i 1 (subtract f/:i/1 (* nod1 saikvk))
]
num2: 0
for k 1 n 1[
for h 1 n 1[
num2: + num2 (* a/:k/:h (* f/:k/1 v/:h/1))
]
]
nod2: / num2 den
for i 1 n 1 [
pokeij v i 1 (subtract f/:i/1 (* nod2 v/:i/1))
]
]
print "Residuals"
display f
Print "Solution"
display x
exit
]
comment { The next function is a demonstration of creating a set of n
linear equations in n unknowns with randomized coefficients and then
solving them by the congruential gradient method using m iterations.
}
test-nxn-its: func [ /local n m a b at ata atb xs][
n: to-integer ask "Number of unknowns? "
m: to-integer ask "Number of iterations (not less than unknowns)? "
a: array reduce [n n]
b: array reduce [n 1]
a: inlzarrand01 a 1e9
b: inlzarrand01 b 1e9
at: transpose a
ata: matprod at a
print "SymPosDef-A"
display ata
atb: matprod at b
print "Right sides"
display atb
cgsoln-its ata atb m
]
