(*~\Глагол\Отделы\ЛинАл\32\к\Набор\~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
(**)                         ОТДЕЛ кнСистЛАУ;
(*============================================================================*
 * НАЗНАЧЕНИЕ: решение Систем Линейных алгебраических Уравнений
 *  (для комплексных в виде набора из 2-х чисел с плавающей точкой в 32 разряда)
 *                                   Матричный Анализ для комплексных
 *============================================================================*)

ИСПОЛЬЗУЕТ
 дМатем   ИЗ "..\..\..\..\..\Отделы\ЧисАн\32\д\",

 дВектор  ИЗ "..\..\..\..\..\Отделы\ЧисАн\32\д\",
 дМатрица ИЗ "..\..\..\..\..\Отделы\ЧисАн\32\д\",

 дСистЛАУ ИЗ "..\..\..\..\..\Отделы\ЛинАл\32\д\",

кнМатем   ИЗ "..\..\..\..\..\Отделы\ЧисАн\32\к\Набор\",
кнВектор  ИЗ "..\..\..\..\..\Отделы\ЧисАн\32\к\Набор\",
кнМатрица ИЗ "..\..\..\..\..\Отделы\ЧисАн\32\к\Набор\";
(*
кнВектВВ  ИЗ "..\..\..\..\..\Отделы\ЧисАн\32\к\Набор\ВводВывод",

  НаЭкран ИЗ "..\..\..\..\..\Отделы\Обмен\ВводВывод\";
*)

ВИД
        (* вид Вещ - это удвоенная точность 32 разряда (6 десятичных знаков)  *)
  Вещ = дМатем.Вещ;                   (* локально *)


(*-оооооооооооооо- опережающие определения/заготовки задач -ооооооооооооооооо-*)


(*----------------------------------------------------------------------------*)
ЗАДАЧА^ РазложитьНаLU-(      A+:кнМатрица.Вид;
                       строкиД+:дСистЛАУ.ДкРядуПерестановок):ЦЕЛ;
(* Цель:  LU-разложение (A=L*U) квадратной матрицы A на два сомножителя L и U,
 *    где L - это нижняя треугольная матрица с единичной главной диагональю,
 *       а U - верхняя треугольная матрица. Эта задача используется
 *        совместно с задачей РешитьИзLU() для решения систем уравнений.
 *)
(*----------------------------------------------------------------------------*)
ЗАДАЧА^ РешитьИзLU-(A-:кнМатрица.Вид;
                    b-: кнВектор.Вид; строкиД:дСистЛАУ.ДкРядуПерестановок;
                    x+: кнВектор.Вид);
(* Цель:  решение системы уравнений после преобразования исходной матрицы A
 *        с помощью задачи РазложитьНаLU()
 *)


(*=рррррррррррррррррррррррр= реализация задач =рррррррррррррррррррррррррррррр=*)


(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА РазложитьНаLU-(      A+:кнМатрица.Вид;
                      строкиД+:дСистЛАУ.ДкРядуПерестановок):ЦЕЛ;
(* Цель:  LU-разложение (A=L*U) квадратной матрицы A на два сомножителя L и U,
 *    где L - это нижняя треугольная матрица с единичной главной диагональю,
 *       а U - верхняя треугольная матрица. Эта задача используется
 *        совместно с задачей РешитьИзLU() для решения систем уравнений.
 *----------------------------------------------------------------------------*
 * До:    <A>    - исходная комплексная матрица
 * После: <A>    - содержит элементы двух треугольных матриц - L и U
 *        <строкиД> - доступ к ряду из (действительных) номеров, т.е.
 *                   (вновь созданные записи о перестановках строк)
 * Ответ: 0 или дМатрица.ВЫРОЖДЕНА *)
ПЕР
             N,                                    (* длина вектора N - любая *)
  последний,       (* =N-1 *)
  i,maxI,j,k:ЦЕЛ;
  c,вед,t:Вещ;
  z,
  кСумма:кнМатем.Вид;
  дВекД:дВектор.Доступ;
УКАЗ
             N:=РАЗМЕР(A);
  последний:=N-1;

  СОЗДАТЬ(  дВекД,N);
  СОЗДАТЬ(строкиД,N);

  ОТ i:=0 ДО последний ВЫП
    вед:=0;
    ОТ j:=0 ДО последний ВЫП
      c:=кнМатем.мод(A[i,j]);
      ЕСЛИ c > вед ТО
        вед:=c 
      КОН 
    КОН;
    ЕСЛИ вед < дМатем.МашинEps ТО
      ВОЗВРАТ дМатрица.ВЫРОЖДЕНА
    КОН;
    дВекД[i]:=1/вед
  КОН;
  ОТ j:=0 ДО последний ВЫП
    ОТ i:=0 ДО j-1 ВЫП
      кСумма:=A[i,j];
      ОТ k:=0 ДО i-1 ВЫП
        кнМатем.умножить(A[i,k],A[k,j],z);
   (*   кнМатем.вычесть(кСумма,z,кСумма)   кСумма =кСумма-z *)
        кнМатем.отнять(z,кСумма)        (* кСумма-=z        *)
      КОН;
      A[i,j]:=кСумма
    КОН;
    maxI:=0;
    вед:=0;
    ОТ i:=j ДО последний ВЫП
      кСумма:=A[i,j];
      ОТ k:=0 ДО j-1 ВЫП
        кнМатем.умножить(A[i,k],A[k,j],z);
   (*   кнМатем.вычесть(кСумма,z,кСумма)   кСумма =кСумма-z *)
        кнМатем.отнять(z,кСумма)        (* кСумма-=z        *)
      КОН;
      A[i,j]:=кСумма;
      t:=дВекД[i]*кнМатем.мод(кСумма);
      ЕСЛИ t > вед ТО
        вед:=t;
        maxI:=i
      КОН
    КОН;
    ЕСЛИ j # maxI ТО
      ОТ k:=0 ДО последний ВЫП
        кнМатем.обмен(A[maxI,k],A[j,k]) 
      КОН;
      дВекД[maxI]:=дВекД[j]
    КОН;
    строкиД[j]:=maxI;
    ЕСЛИ кнМатем.мод(A[j,j]) = 0 ТО
      кнМатем.алгебрВид(дМатем.МАЛОЕ,дМатем.МАЛОЕ,A[j,j])
    КОН;
    ЕСЛИ j # последний ТО
      ОТ i:=j+1 ДО последний ВЫП
 (*     кнМатем.делить(A[i,j],A[j,j],A[i,j])   *)
         кнМатем.уменьшВкомпл(A[j,j],A[i,j])  
      КОН
    КОН
  КОН;
  ВОЗВРАТ 0
КОН РазложитьНаLU;

(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА РешитьИзLU-(A-:кнМатрица.Вид;
                   b-: кнВектор.Вид; строкиД:дСистЛАУ.ДкРядуПерестановок;
                   x+: кнВектор.Вид);
(* Цель:  решение системы уравнений после преобразования исходной матрицы A
 *        с помощью задачи РазложитьНаLU()
 *----------------------------------------------------------------------------*
 * До:    <A>   - исходная комплексная матрица A после задачи РазложитьНаLU()
 *        <b>   - комплексный вектор свободных членов
 *    <строкиД> - доступ к ряду из (действительных) номеров, т.е.
 *                   (записи о перестановках строк из РазложитьНаLU)
 * После: <x>   - комплексный вектор решения                                  *)
ПЕР
             N,                                    (* длина вектора N - любая *)
  последний,       (* =N-1 *)
  i,Ip,j,k:ЦЕЛ;
  z,
  кСумма:кнМатем.Вид;
УКАЗ
             N:=РАЗМЕР(A);
  последний:=N-1;

  k:=-1;
  ОТ i:=0 ДО последний ВЫП
    x[i]:=b[i]
  КОН;
  ОТ i:=0 ДО последний ВЫП
    Ip:=строкиД[i];
    кСумма:=x[Ip];
    x[Ip]:=x[i];
    ЕСЛИ k >= 0 ТО
      ОТ j:=k ДО i-1 ВЫП
        кнМатем.умножить(A[i,j],x[j],z);
   (*   кнМатем.вычесть(кСумма,z,кСумма)   кСумма=кСумма-z *)
        кнМатем.отнять(z,кСумма)        (* кСумма=кСумма-z *)
      КОН 
    ИНАЧЕ 
      ЕСЛИ кнМатем.мод(кСумма) # 0 ТО
        k:=i 
      КОН 
    КОН;
    x[i]:=кСумма
  КОН;
  ОТ i:=последний ДО 0 ПО -1 ВЫП
    кСумма:=x[i];
    ЕСЛИ i < последний ТО
      ОТ j:=i+1 ДО последний ВЫП
        кнМатем.умножить(A[i,j],x[j],z);
   (*   кнМатем.вычесть(кСумма,z,кСумма)   кСумма=кСумма-z *)
        кнМатем.отнять(z,кСумма)        (* кСумма=кСумма-z *)
      КОН 
    КОН;
    кнМатем.делить(кСумма,A[i,i],x[i])
  КОН
КОН РешитьИзLU;

КОН кнСистЛАУ.

