J: исключение Гаусса-Иордана

Question

Задача для кодирования метода Гаусса-Джордана для решения линейной системы алгебраических уравнений - это упражнение, которое я выбрал для продвижения в обучении J. Система - это Ax = b, где A - это n- by-n матрица, b и неизвестные x равны n-векторы. Во-первых, я начал с самой простой формы со структурами управления:

   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.- переходить в молчаливую форму... но кажется, что это будет выглядеть более громоздким, и я думаю, что я ошибаюсь - но что учиться без ошибок? Я действительно хочу кодировать метод Гаусса-Джордана, чтобы негласно:

• упражнение в 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 - но метод карандашной бумаги помогает.