(*~\Глагол\Отделы\Числа~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
(**)                        ОТДЕЛ КомплМатр;
(*============================================================================*
 * НАЗНАЧЕНИЕ: определение матриц из комплексных чисел
 *             и основные вычисления над ними
 *
 * Источники, ссылки, библиография
Debord J., Библиотека математических подпрограмм
   tpmath.zip

Goffe B., Программа SimAnn.FOR (Simulated Annealing in Fortran)
   http://www.netlib.org/opt/simann.f

EISPACK: Библиотека Фортран функций для вычисления собственных значений и векторов
   http://www.netlib.org/eispack

Marsaglia G., Тесты для генераторов случайных чисел
   http://stat.fsu.edu/~geo/diehard.html

Moshier S., Библиотека математических подпрограмм
   http://www.moshier.net

Press W.H., Teukolsky S.A., Vetterling W.T., Flannery B.P.
   Численные рецепты на Паскале
   http://garbo.uwasa.fi/pc/turbopas.html

Пакет численных методов для Турбо Паскаля
   Borland International, 1986
 *============================================================================*)
ИСПОЛЬЗУЕТ
  Матем,
  Компл,
  Вект,
  КомплВект,
  Матр;

ВИД
  Вид-=РЯД ИЗ КомплВект.Вид;
  Доступ-=ДОСТУП К Вид;

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

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

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

КОН КомплМатр.