J: Гаусс-Жордана

Question

Задача закодировать метод Гаусса-Жордана решения системы линейных алгебраических уравнений является упражнение, которое я выбрал, чтобы продвинуться в изучении J. системы является Ах = Ь, где является пп матрицы с размерностью матрицей, б и неизвестные х являются п -векторов. Во-первых, я начал с простейшей формы с контрольными структурами:

gj0 =: dyad :0 NB. usage is same to %. 
    y=.y,.b 
    for_d.i.#y do. 
    for_r.i.#y do. 
     if.r=d do.continue.end. NB. do not eliminate d'th row 
     t=.%/ (<"1(r,d),:(d,d)) { y 
     for_c.d}.>:i.#y do. 
     y=.(((<r,c){y)-(t*(<d,c){y)) (<r,c)} y 
     end. 
     y=.0 (<r,d)} y NB. ensure zero 
    end. 
    end. NB. now A is diagonal but not identity matrix, so: 
    x=.{:"1 y NB. x = b 
    for_r.i.#y do. 
    x=.((r{x)%(<r,r){y) r} x NB. divide by coefficients on diagonal 
    end. 
    ) 
    Ab =: (".;._2) 0 :0 
    0.25 _0.16 _0.38 0.17 
    0.19 _0.22 _0.02 0.41 
    0.13 0.08 _0.08 _0.13 
    0.13 _0.1 _0.32 0.65 
    ) 
    b =: 0.37 0.01 0.01 1.51 
    (,.".&.>)('A';'b';'gj0 A,.b';'b %. A') 
┌────────┬──────────────────────┐ 
│A  │0.25 _0.16 _0.38 0.17│ 
│  │0.19 _0.22 _0.02 0.41│ 
│  │0.13 0.08 _0.08 _0.13│ 
│  │0.13 _0.1 _0.32 0.65│ 
├────────┼──────────────────────┤ 
│b  │0.37 0.01 0.01 1.51 │ 
├────────┼──────────────────────┤ 
│b gj0 A │_1 3 _2 2    │ 
├────────┼──────────────────────┤ 
│b %. A │_1 3 _2 2    │ 
└────────┴──────────────────────┘ 

Правильно! Далее я решил избавиться от как можно большего числа структур управления, как это возможно:

gj1 =:dyad :0 
    y=.y,.b 
    for_d.i.#y do. 
    for_r.d ({.,]}.~[:>:[) i.#y do. NB. for indices without d 
     t=.%/ (<"1(r,d),:(d,d)) { y 
     y=.((r{y)-(t*d{y)) r}y NB. no need to iterate for each column 
     y=.0 (<r,d)} y 
    end. 
    end. 
    ({:"1 y)%(+/}:"1 y) NB. b divide by sum of each item of A (drop zeroes) 
    ) 
    b gj1 A 
_1 3 _2 2 

ОК, теперь я могу попробовать перевести for_r. -loop в молчаливой форме ... но мне кажется, что это будет выглядеть более громоздким и Я думаю, что я ошибаюсь, но что учится без ошибок? Я действительно хочу, чтобы закодировать метод Гаусса-Жордана молчаливо:

• упражнения в J кодирования
• увидеть, если это лучше в исполнении
• попытаться понять код через несколько недель :)

Помогите мне, пожалуйста, написать его до конца или указать лучший подход.

Answer 1

Благодаря Eelvex, который посоветовал мне посмотреть в addons/math/misc/linear.ijs, я завершил задачу с этим прекрасным кодом:

gj=: monad :0 
    I=. i.#y 
    for_i. I do. y=. y - (col - i=I) */ (i{y) % i{col=. i{"1 y end. 
    ) 
    gj Ab 
1 0 0 0 _1 
0 1 0 0 3 
0 0 1 0 _2 
0 0 0 1 2 

Это заняло некоторое время, чтобы понять глагол pivot в linear.ijs - но pencil- метод бумаги помогает.