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

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

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

ВИД
  ДкРядуПерестановок- = ДОСТУП К цВектор.Вид;
(*
  ДкРядуПерестановок- = ДОСТУП К РЯД ИЗ ЦЕЛ;
*)
        (* вид Вещ - это удвоенная точность 32 разряда (6 десятичных знаков)  *)
  Вещ = дМатем.Вещ;                   (* локально *)
 
(*=
(*----------------------------------------------------------------------------*
 * Код ошибки после выполнения функций
 *----------------------------------------------------------------------------*)
ПОСТ
  ВЫРОЖДЕНА-           =-1;  (* вырожденная матрица *)
  НЕТ_СХОДИМОСТИ-      =-2;  (* нет сходимости *)
  НЕ_ПоложитОПРЕДЕЛЕНа-=-3;  (* матрица не является положительно определенной *)
ПОСТ
  МАЛОЕ = 1.D-20;               (* вместо 0. с точностью 32 разряда *)
=*)

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


(*----------------------------------------------------------------------------*)
ЗАДАЧА^ МетодГаусса-(A- :дМатрица.Вид; b-:дВектор.Вид;
                     Ao+:дМатрица.Вид; x+:дВектор.Вид):ЦЕЛ;
(*
 * Цель:  решение системы линейных алгебраических уравнений методом
 *        исключения Гаусса *)
(*----------------------------------------------------------------------------*)
ЗАДАЧА^ Обратить-(A-,Ao+:дМатрица.Вид):ЦЕЛ;
(*
   Цель:  обращение квадратной матрицы методом Гаусса
 *)
(*----------------------------------------------------------------------------*)
ЗАДАЧА^ Определитель-(A+:дМатрица.Вид):Вещ;
(*
   Цель:  вычисление определителя квадратной матрицы
 *)
(*----------------------------------------------------------------------------*)
ЗАДАЧА^ РазложитьНаLL-(A-,L+:дМатрица.Вид):ЦЕЛ;
(*
 * Цель:  Разложение Холецкого (A=L*L') симметрической положительно
 *        определенной матрицы A на две треугольные матрицы L и L', где L - это
 *        нижняя треугольная матрица, а L' - матрица, транспонированная к L.
 *)
(*----------------------------------------------------------------------------*)
ЗАДАЧА^ РазложитьНаLU-(      A+:дМатрица.Вид;
                       строкиД+:дСистЛАУ.ДкРядуПерестановок):ЦЕЛ;
(* Цель:  LU-разложение (A=L*U) квадратной матрицы A на два сомножителя L и U, 
 *        где L - это нижняя треугольная матрица с единичной главной диагональю,
 *        а U - верхняя треугольная матрица. Эта задача используется
 *        совместно с РешитьИзLU для решения систем уравнений.                *)
(*----------------------------------------------------------------------------*)
ЗАДАЧА^ РешитьИзLU-(A-:дМатрица.Вид;
                    b-: дВектор.Вид; строкиД:дСистЛАУ.ДкРядуПерестановок;
                    x+: дВектор.Вид);
(* Цель:  решение системы уравнений после преобразования исходной матрицы
 *        с помощью РазложитьНаLU *)
(*----------------------------------------------------------------------------*)
ЗАДАЧА^ РазложитьНаSV-(A+:дМатрица.Вид;
                       s+: дВектор.Вид;
                       V+:дМатрица.Вид):ЦЕЛ;
(* Цель:  SV разложение прямоугольной матрицы A (N x M, с N >= M) на три
 *        сомножителя (A=U*S*V'), где U (N x M) состоит из столбцов - собственных
 *        векторов, S (M x M) - диагональная матрица с элементами >= 0 (собственные
 *        значения) и V (M x M) - ортогональная матрица. Эта задача используется
 *        совместно с РешитьИзSV для решения систем уравнений. *)
(*----------------------------------------------------------------------------*)
ЗАДАЧА^ ОбнулитьS-(s+:дВектор.Вид; откл:Вещ);
(*
   Цель:  обнуляет наименьшие собственные значения по заданной границе
 *)
(*----------------------------------------------------------------------------*)
ЗАДАЧА^ РешитьИзSV-(U-:дМатрица.Вид; s-   :дВектор.Вид;
                    V-:дМатрица.Вид; b-,x+:дВектор.Вид);
(* Цель:  решение системы уравнений после преобразования исходной матрицы
 *        с помощью РазложитьНаSV и обнуления наименьших собственных значений
 *        с помощью ОбнулитьS *)
(*----------------------------------------------------------------------------*)
ЗАДАЧА^ ВосстановитьИзSV-(U-:дМатрица.Вид;
                          s-: дВектор.Вид;
                       V-,A+:дМатрица.Вид);
(* Цель:  восстанавливает матрицу A перемножая U, S и V' после
 *        обнуления наименьших собственных значений задачей ОбнулитьS *)
(*----------------------------------------------------------------------------*)
ЗАДАЧА^ РазложитьНаQR-(A+,R+:дМатрица.Вид):ЦЕЛ;
(*
 * Цель:  QR разложение прямоугольной матрицы <A> (N x M, с N >= M) на два
 *        сомножителя (A=Q*R), где Q (N x M) - ортогональная матрица, а
 *        R (M x M) - верхняя треугольная матрица. Эта задача используется
 *        совместно с РешитьИзQR для решения систем уравнений.
 *)
(*----------------------------------------------------------------------------*)
ЗАДАЧА^ РешитьИзQR-(q-,R-:дМатрица.Вид;
                    b-,x+: дВектор.Вид);
(*
 * Цель:  решение системы уравнений после преобразования исходной матрицы
 *        с помощью РазложитьНаQR
 *)
(*----------------------------------------------------------------------------*)
ЗАДАЧА^ ГильбертМатрица-(М+:дМатрица.Вид); 
(*
 * Построение матрицы Гильберта <порядок>-го порядка
 *)
(*----------------------------------------------------------------------------*)
ЗАДАЧА^ ГильбертСистема-(A+:дМатрица.Вид;
                         b+: дВектор.Вид); 
(*
   Построение системы Гильберта n-го порядка
 *)


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


(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА МетодГаусса-(A- :дМатрица.Вид; b-:дВектор.Вид;
                    Ao+:дМатрица.Вид; x+:дВектор.Вид):ЦЕЛ;
(* Цель:  решение системы линейных алгебраических уравнений методом
 *        исключения Гаусса
 *----------------------------------------------------------------------------*
 * До:    <A>  - матрица системы
 *        <b>  - вектор свободных членов
 * После: <Ao> - обратная матрица
 *        <x>  - вектор решения
 * Ответ: 0 или дМатрица.ВЫРОЖДЕНА *)
ПЕР
  i,j,k,посл:ЦЕЛ;
  вед,t:Вещ;
  ведСтр,ведСтлб:ДОСТУП К РЯД ИЗ ЦЕЛ; (* строки и столбцы ведущего элемента *)
УКАЗ
  посл:=РАЗМЕР(A)-1;
  СОЗДАТЬ(ведСтр,посл+1);
  СОЗДАТЬ(ведСтлб,посл+1);
   дМатрица.Списать(A,Ao); (* Ao=A *)
(* дМатрица.СписатьВектор(b,x);  (* x=b     СписатьВектор(b,x) *) *)
    дВектор.Списать(b,x);        (* x=b     СписатьВектор(b,x) *)

  ОТ k:=0 ДО посл ВЫП
    (* поиск наибольшего ведущего элемента в подматрице Ao[k..посл,k..посл] *)
    вед:=Ao[k,k];
    ведСтр[k]:=k;
    ведСтлб[k]:=k;
    ОТ i:=k ДО посл ВЫП
      ОТ j:=k ДО посл ВЫП
        ЕСЛИ МОДУЛЬ(Ao[i,j]) > МОДУЛЬ(вед) ТО
          вед:=Ao[i,j];
          ведСтр[k]:=i;
          ведСтлб[k]:=j
        КОН 
      КОН 
    КОН;

    (* ведущий элемент слабоват => псевдо-вырожденная матрица *)
    ЕСЛИ МОДУЛЬ(вед) < дМатем.МашинEps ТО
      ВОЗВРАТ дМатрица.ВЫРОЖДЕНА
    КОН;

    (* обмен текущей строки (k) со строкой ведущего элемента *)
    ЕСЛИ ведСтр[k] # k ТО
      дМатрица.ОбменСтрок(Ao,ведСтр[k],k);
      дМатем.обмен(x[ведСтр[k]],x[k])
    КОН;

    (* обмен текущего столбца (k) со столбцом ведущего элемента *)
    ЕСЛИ ведСтлб[k] # k ТО
      дМатрица.ОбменСтолбцов(Ao,ведСтлб[k],k) 
    КОН;

    (* преобразование строки ведущего элемента *)
    Ao[k,k]:=1;
    ОТ j:=0 ДО посл ВЫП
      Ao[k,j]:=Ao[k,j]/вед 
    КОН;
    x[k]:=x[k]/вед;

    (* преобразование оставшихся строк *)
    ОТ i:=0 ДО посл ВЫП
      ЕСЛИ i # k ТО
        t:=Ao[i,k];
        Ao[i,k]:=0;
        ОТ j:=0 ДО посл ВЫП
          Ao[i,j]:=Ao[i,j]-t*Ao[k,j]
        КОН;
        x[i]:=x[i]-t*x[k]
      КОН 
    КОН
  КОН;

  (* перестановка обратной матрицы *)
  ОТ i:=посл ДО 0 ПО -1 ВЫП
    ЕСЛИ ведСтлб[i] # i ТО
      дМатрица.ОбменСтрок(Ao,ведСтлб[i],i);
      дМатем.обмен(x[ведСтлб[i]],x[i])
    КОН 
  КОН;
  ОТ j:=посл ДО 0 ПО -1 ВЫП
    ЕСЛИ ведСтр[j] # j ТО
      дМатрица.ОбменСтолбцов(Ao,ведСтр[j],j) 
    КОН 
  КОН;

  ВОЗВРАТ 0
КОН МетодГаусса;

(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА Обратить-(A- :дМатрица.Вид;
                 Ao+:дМатрица.Вид):ЦЕЛ;
(*
 * Цель:  обращение квадратной матрицы методом Гаусса
 *----------------------------------------------------------------------------*
 * До:    <A>  - исходная матрица
 * После: <Ao> - обратная матрица
 * Ответ: 0 или дМатрица.ВЫРОЖДЕНА
 *)
ПЕР
  i,j,k,посл:ЦЕЛ;
  вед,t:Вещ;
  ведСтр,ведСтлб:ДОСТУП К РЯД ИЗ ЦЕЛ;   (* строки и столбцы ведущего элемента *)
УКАЗ
  посл:=РАЗМЕР(A)-1;
  СОЗДАТЬ(ведСтр,посл+1);
  СОЗДАТЬ(ведСтлб,посл+1);
   дМатрица.Списать(A,Ao);

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

  ВОЗВРАТ 0
КОН Обратить;

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

    ЕСЛИ НЕ опред0 ТО
      ОТ i:=k+1 ДО посл ВЫП
        ЕСЛИ МОДУЛЬ(A[i,k]) > дМатем.МашинEps ТО
          (* обнуляем k элемент этого ряда *)
          t:=-A[i,k]/A[k,k];
          ОТ j:=1 ДО посл ВЫП
            A[i,j]:=A[i,j]+t*A[k,j]
          КОН
        КОН 
      КОН 
    КОН;

    d:=d*A[k,k]; 
    УВЕЛИЧИТЬ(k)
  КОН;

  ЕСЛИ опред0 ТО
    ВОЗВРАТ 0
  КОН;
  ВОЗВРАТ d*A[посл,посл]
КОН Определитель;

(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА РазложитьНаLL-(A-,L+:дМатрица.Вид):ЦЕЛ;
(*
 * Цель:  Разложение Холецкого (A=L*L') симметрической положительно
 *        определенной матрицы A на две треугольные матрицы L и L', где L - это
 *        нижняя треугольная матрица, а L' - матрица, транспонированная к L.
 *----------------------------------------------------------------------------*
 * Прим:  Эту задачу можно использовать для испытания на положительно
 *        определенность.
 * До:    <A>    - исходная матрица A
 * После: <L>    - матрица L
 * Ответ: 0 или дМатрица.НЕ_ПоложитОПРЕДЕЛЕНа
 *)
ПЕР
  i,j,k,посл:ЦЕЛ;
  сумма:Вещ;
УКАЗ
  посл:=РАЗМЕР(A)-1;
  ОТ k:=0 ДО посл ВЫП
    сумма:=A[k,k];
    ОТ j:=0 ДО k-1 ВЫП
      сумма:=сумма-дМатем.кв(L[k,j]) 
    КОН;

    ЕСЛИ сумма <= 0 ТО
      ВОЗВРАТ дМатрица.НЕ_ПоложитОПРЕДЕЛЕНа
    КОН;

    L[k,k]:=дМатем.квкор(сумма);
    ОТ i:=k+1 ДО посл ВЫП
      сумма:=A[i,k];
      ОТ j:=0 ДО k-1 ВЫП
        сумма:=сумма-L[i,j]*L[k,j]
      КОН;
      L[i,k]:=сумма/L[k,k]
    КОН
  КОН;
  ВОЗВРАТ 0
КОН РазложитьНаLL;

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

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

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

  ОТ i:=0 ДО посл ВЫП
    Ip:=строкиД[i];
    сумма:=x[Ip];
    x[Ip]:=x[i];
    ЕСЛИ k >= 0 ТО
      ОТ j:=k ДО i-1 ВЫП
        сумма:=сумма-A[i,j]*x[j]
      КОН 
    ИНАЧЕ 
      ЕСЛИ сумма # 0 ТО
        k:=i 
      КОН 
    КОН;
    x[i]:=сумма
  КОН;
  ОТ i:=посл ДО 0 ПО -1 ВЫП
    сумма:=x[i];
    ЕСЛИ i < посл ТО
      ОТ j:=i+1 ДО посл ВЫП
        сумма:=сумма-A[i,j]*x[j]
      КОН 
    КОН;
    x[i]:=сумма/A[i,i]
  КОН
КОН РешитьИзLU;

(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА РазложитьНаSV-(A+:дМатрица.Вид;
                      s+: дВектор.Вид;
                      V+:дМатрица.Вид):ЦЕЛ;
(*
 * Цель:  SV разложение прямоугольной матрицы A (N x M, с N >= M) на три
 *        сомножителя (A=U*S*V'), где U (N x M) состоит из столбцов - собственных
 *        векторов, S (M x M) - диагональная матрица с элементами >= 0 (собственные
 *        значения) и V (M x M) ортогональная матрица. Эта задача используется
 *        совместно с РешитьИзSV для решения систем уравнений;
 *----------------------------------------------------------------------------*
 * Прим:  текст этой задачи переписан из подпрограммы SVD пакета EISPACK
 * До:    <A>      - исходная матрица
 * После: <A>      - содержит элементы U
 *        <s>      - вектор собственных значений
 *        <V>      - ортогональная матрица       
 * Ответ: 0 или дМатрица.НЕТ_СХОДИМОСТИ
 *)
ПЕР
  i,Its,j,i1,k,l,n,посл1,посл2:ЦЕЛ;
  Anorm,c,f,g,h,Scale,t,x,y,z:Вещ;
  r:ДОСТУП К дВектор.Вид;
  label1:КЛЮЧ;
УКАЗ
  посл1:=РАЗМЕР(A,0)-1;
  посл2:=РАЗМЕР(A,1)-1;
  g:=0;
  Scale:=0;
  Anorm:=0;
  СОЗДАТЬ(r,посл2+1);
                     (* приведение к бидиагональному виду методом Хаусхолдера *)
  ОТ i:=0 ДО посл2 ВЫП
    l:=i+1;
    r[i]:=Scale*g;
    g:=0;
    t:=0;
    Scale:=0;
    ЕСЛИ i <= посл1 ТО
      ОТ k:=i ДО посл1 ВЫП
        Scale:=Scale+МОДУЛЬ(A[k,i]) 
      КОН;
      ЕСЛИ Scale # 0 ТО
        ОТ k:=i ДО посл1 ВЫП
          A[k,i]:=A[k,i]/Scale;
          t:=t+дМатем.кв(A[k,i])
        КОН;
        f:=A[i,i];
        g:=-дМатем.знакЧ(f)*дМатем.квкор(t);
        h:=f*g-t;
        A[i,i]:=f-g;
        ЕСЛИ i # посл2 ТО
          ОТ j:=l ДО посл2 ВЫП
            t:=0;
            ОТ k:=i ДО посл1 ВЫП
              t:=t+A[k,i]*A[k,j]
            КОН;
            f:=t/h;
            ОТ k:=i ДО посл1 ВЫП
              A[k,j]:=A[k,j]+f*A[k,i]
            КОН
          КОН
        КОН;
        ОТ k:=i ДО посл1 ВЫП
          A[k,i]:=Scale*A[k,i]
        КОН
      КОН
    КОН;
    s[i]:=Scale*g;
    g:=0;
    t:=0;
    Scale:=0;
    ЕСЛИ (i <= посл1) И (i # посл2) ТО
      ОТ k:=l ДО посл2 ВЫП
        Scale:=Scale+МОДУЛЬ(A[i,k]) 
      КОН;
      ЕСЛИ Scale # 0 ТО
        ОТ k:=l ДО посл2 ВЫП
          A[i,k]:=A[i,k]/Scale;
          t:=t+дМатем.кв(A[i,k])
        КОН;
        f:=A[i,l];
        g:=-дМатем.знакЧ(f)*дМатем.квкор(t);
        h:=f*g-t;
        A[i,l]:=f-g;
        ОТ k:=l ДО посл2 ВЫП
          r[k]:=A[i,k]/h 
        КОН;
        ЕСЛИ i # посл1 ТО
          ОТ j:=l ДО посл1 ВЫП
            t:=0;
            ОТ k:=l ДО посл2 ВЫП
              t:=t+A[j,k]*A[i,k]
            КОН;
            ОТ k:=l ДО посл2 ВЫП
              A[j,k]:=A[j,k]+t*r[k]
            КОН;
          КОН 
        КОН;
        ОТ k:=l ДО посл2 ВЫП
          A[i,k]:=Scale*A[i,k]
        КОН
      КОН
    КОН;
    Anorm:=дМатем.макс(Anorm,МОДУЛЬ(s[i])+МОДУЛЬ(r[i]));
  КОН;
                                          (* накопление правых преобразований *)
  ОТ i:=посл2 ДО 0 ПО -1 ВЫП
    ЕСЛИ i < посл2 ТО
      ЕСЛИ g # 0 ТО
        ОТ j:=l ДО посл2 ВЫП
                 (* двойное деление предохраняет от возможной потери точности *)
          V[j,i]:=(A[i,j]/A[i,l])/g 
        КОН;
        ОТ j:=l ДО посл2 ВЫП
          t:=0;
          ОТ k:=l ДО посл2 ВЫП
            t:=t+A[i,k]*V[k,j]
          КОН;
          ОТ k:=l ДО посл2 ВЫП
            V[k,j]:=V[k,j]+t*V[k,i]
          КОН
        КОН
      КОН;
      ОТ j:=l ДО посл2 ВЫП
        V[i,j]:=0;
        V[j,i]:=0
      КОН
    КОН;
    V[i,i]:=1;
    g:=r[i];
    l:=i
  КОН;
                                           (* накопление левых преобразований *)
  ОТ i:=посл2 ДО 0 ПО -1 ВЫП
    l:=i+1;
    g:=s[i];
    ЕСЛИ i < посл2 ТО
      ОТ j:=l ДО посл2 ВЫП
        A[i,j]:=0 
      КОН 
    КОН;
    ЕСЛИ g # 0 ТО
      ЕСЛИ i # посл2 ТО
        ОТ j:=l ДО посл2 ВЫП
          t:=0;
          ОТ k:=l ДО посл1 ВЫП
            t:=t+A[k,i]*A[k,j]
          КОН;
                 (* двойное деление предохраняет от возможной потери точности *)
          f:=(t/A[i,i])/g;
          ОТ k:=i ДО посл1 ВЫП
            A[k,j]:=A[k,j]+f*A[k,i]
          КОН;
        КОН 
      КОН;
      ОТ j:=i ДО посл1 ВЫП
        A[j,i]:=A[j,i]/g 
      КОН;
    ИНАЧЕ
      ОТ j:=i ДО посл1 ВЫП
        A[j,i]:=0 
      КОН 
    КОН;
    A[i,i]:=A[i,i]+1
  КОН;
                                       (* диагонализация бидиагональной формы *)
  ОТ k:=посл2 ДО 0 ПО -1 ВЫП
    Its:=0;
    КОЛЬЦО
      l:=k;
                                                   (* проверка на расщепление *)
      КОЛЬЦО
        ЕСЛИ l < 0 ТО
          label1:=ВКЛ;
          ВЫХОД
        КОН;
        n:=l-1;
        ЕСЛИ (МОДУЛЬ(r[l])+Anorm) = Anorm ТО 
          label1:=ОТКЛ; 
          ВЫХОД
        КОН;
                (* r[0] всегда =0, поэтому s[-1] никогда не будет вычисляться *)
        ЕСЛИ (МОДУЛЬ(s[n])+Anorm) = Anorm ТО 
          label1:=ВКЛ;
          ВЫХОД
        КОН;
        УМЕНЬШИТЬ(l)
      КОН; (* кольцо *)

      ЕСЛИ label1 ТО 
        c:=0;
        t:=1;
        i:=l; 
        КОЛЬЦО
          ЕСЛИ i > k ТО
            ВЫХОД
          КОН;
          f:=t*r[i];
          r[i]:=c*r[i];
          ЕСЛИ (МОДУЛЬ(f)+Anorm) = Anorm ТО
            ВЫХОД
          КОН;
          g:=s[i];
          h:=дМатем.Пифагор(f,g);
          s[i]:=h;
          c:=g/h;
          t:=-f/h;
          ОТ j:=0 ДО посл1 ВЫП
            y:=A[j,n];
            z:=A[j,i];
            A[j,n]:=y*c+z*t;
            A[j,i]:=-y*t+z*c
          КОН;
          УВЕЛИЧИТЬ(i)
        КОН (* кольцо *)
      КОН;
                                                       (* проверка сходимости *)
      z:=s[k];
      ЕСЛИ l = k ТО
        ЕСЛИ z < 0 ТО
                                              (* сделаем s[k] неотрицательным *)
          s[k]:=-z;
          ОТ j:=0 ДО посл2 ВЫП
            V[j,k]:=-V[j,k]
          КОН
        КОН;
                               (* сходимость *)
        ВЫХОД
      КОН;
      ЕСЛИ Its = 30 ТО
        ВОЗВРАТ дМатрица.НЕТ_СХОДИМОСТИ
      КОН;
      УВЕЛИЧИТЬ(Its);
                               (* сдвиг снизу на минор 2 х 2 *)
      x:=s[l];
      n:=k-1;
      y:=s[n];
      g:=r[n];
      h:=r[k];
      f:=0.5*(((g+z)/h)*((g-z)/y)+y/h-h/y);
      g:=дМатем.Пифагор(f,1);
      f:=x-(z/x)*z+(h/x)*(y/(f+дМатем.знакЧ(f)*МОДУЛЬ(g))-h);
                                               (* следующее QR преобразование *)
      c:=1;
      t:=1;
      ОТ i1:=l ДО n ВЫП
        i:=i1+1;
        g:=r[i];
        y:=s[i];
        h:=t*g;
        g:=c*g;
        z:=дМатем.Пифагор(f,h);
        r[i1]:=z;
        c:=f/z;
        t:=h/z;
        f:=x*c+g*t;
        g:=-x*t+g*c;
        h:=y*t;
        y:=y*c;
        ОТ j:=0 ДО посл2 ВЫП
          x:=V[j,i1];
          z:=V[j,i];
          V[j,i1]:=x*c+z*t;
          V[j,i]:=-x*t+z*c
        КОН;
        z:=дМатем.Пифагор(f,h);
        s[i1]:=z;
                                  (* вращение может быть произвольным при z=0 *)
        ЕСЛИ z # 0 ТО
          c:=f/z;
          t:=h/z
        КОН;
        f:=c*g+t*y;
        x:=-t*g+c*y;
        ОТ j:=0 ДО посл1 ВЫП
          y:=A[j,i1];
          z:=A[j,i];
          A[j,i1]:=y*c+z*t;
          A[j,i]:=-y*t+z*c
        КОН 
      КОН;
      r[l]:=0;
      r[k]:=f;
      s[k]:=x
    КОН (* кольцо *)
  КОН;
  ВОЗВРАТ 0
КОН РазложитьНаSV;

(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА ОбнулитьS-(s+:дВектор.Вид; откл:Вещ);
(*
 * Цель:  обнуляет наименьшие собственные значения по заданной границе
 *----------------------------------------------------------------------------*
 * До:    <s>    - вектор собственных значений S
 *        <откл> - относительное отклонение (граница будет равна откл*макс(S) )
 * После: <s>    - измененный вектор собственных значений
 *)
ПЕР
  граница:Вещ;
  i:ЦЕЛ;
УКАЗ
  граница:=откл*дВектор.макс(s);
  ОТ i:=0 ДО РАЗМЕР(s)-1 ВЫП
    ЕСЛИ s[i] < граница ТО 
      s[i]:=0 
    КОН 
  КОН
КОН ОбнулитьS;

(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА РешитьИзSV-(U-:дМатрица.Вид; s-   :дВектор.Вид;
                   V-:дМатрица.Вид; b-,x+:дВектор.Вид);
(*
 * Цель:  решение системы уравнений после преобразования исходной матрицы
 *        с помощью РазложитьНаSV и обнуления наименьших собственных значений
 *        с помощью ОбнулитьS
 *----------------------------------------------------------------------------*
 * До:    <U> - матрица U после РазложитьНаSV
 *        <s> - вектор  S после РазложитьНаSV
 *        <V> - матрица V после РазложитьНаSV
 *        <b> - вектор свободных членов
 * После: <x> - вектор решения (V*Diag(1/S(i))*U'*b, для S(i) # 0)
 *)
ПЕР
  i,j,посл1,посл2:ЦЕЛ;
  сумма:Вещ;
  врем:ДОСТУП К дВектор.Вид;
УКАЗ
  посл1:=РАЗМЕР(U,0)-1;
  посл2:=РАЗМЕР(U,1)-1;
  СОЗДАТЬ(врем,посл2+1);
  ОТ j:=0 ДО посл2 ВЫП
    сумма:=0;
    ЕСЛИ s[j] > 0 ТО
      ОТ i:=0 ДО посл1 ВЫП
        сумма:=сумма+U[i,j]*b[i]
      КОН;
      сумма:=сумма/s[j]
    КОН;
    врем[j]:=сумма
  КОН;
  ОТ j:=0 ДО посл2 ВЫП
    сумма:=0;
    ОТ i:=0 ДО посл2 ВЫП
      сумма:=сумма+V[j,i]*врем[i]
    КОН;
    x[j]:=сумма
  КОН;
КОН РешитьИзSV;

(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА ВосстановитьИзSV-(U-:дМатрица.Вид;
                         s-: дВектор.Вид;
                      V-,A+:дМатрица.Вид);
(*
 * Цель:  восстанавливает матрицу A перемножая U, S и V' после
 *        обнуления наименьших собственных значений задачей ОбнулитьS() 
 *----------------------------------------------------------------------------*
 * До:    <U> - матрица U после РазложитьНаSV
 *        <s> - вектор  S после РазложитьНаSV
 *        <V> - матрица V после РазложитьНаSV
 * После: <A> - восстановленная матрица
 *)
ПЕР
  i,j,k,посл1,посл2:ЦЕЛ;
УКАЗ
  посл1:=РАЗМЕР(U,0)-1;
  посл2:=РАЗМЕР(U,1)-1;
  ОТ i:=0 ДО посл1 ВЫП
    ОТ j:=0 ДО посл2 ВЫП
      A[i,j]:=0;
      ОТ k:=0 ДО посл2 ВЫП
        ЕСЛИ s[k] > 0 ТО
          A[i,j]:=A[i,j]+U[i,k]*V[j,k]
        КОН 
      КОН
    КОН 
  КОН
КОН ВосстановитьИзSV;

(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА РазложитьНаQR-(A+,R+:дМатрица.Вид):ЦЕЛ;
(*
 * Цель:  QR разложение прямоугольной матрицы <A> (N x M, с N >= M) на два
 *        сомножителя (A=Q*R), где Q (N x M) ортогональная матрица, а
 *        R (M x M) - верхняя треугольная матрица. Эта задача используется
 *        совместно с РешитьИзQR для решения систем уравнений.
 *----------------------------------------------------------------------------*
 * До:    <A> - исходная матрица
 * После: <A> - содержит элементы Q
 *        <R> - верхняя треугольная матрица
 * Ответ: 0 или ВЫРОЖДЕНА
 *)
ПЕР
  i,j,k,посл1,посл2:ЦЕЛ;
  сумма:Вещ;
УКАЗ
  посл1:=РАЗМЕР(A,0)-1;
  посл2:=РАЗМЕР(A,1)-1;
  ОТ k:=0 ДО посл2 ВЫП
                                    (* вычисляем k-й диагональный элемент <R> *)
    сумма:=0;
    ОТ i:=0 ДО посл1 ВЫП
      сумма:=сумма+дМатем.кв(A[i,k]) 
    КОН;

    ЕСЛИ сумма = 0 ТО
      ВОЗВРАТ дМатрица.ВЫРОЖДЕНА
    КОН;

    R[k,k]:=дМатем.квкор(сумма);

                    (* делим k-й столбец <A> на k-й диагональный элемент <R>
                                         * (это начало перезаписи <A> на <q>) *)
    ОТ i:=0 ДО посл1 ВЫП
      A[i,k]:=A[i,k]/R[k,k]
    КОН;

    ОТ j:=(k+1) ДО посл2 ВЫП
                               (* дополняем оставшиеся элементы в столбце <R> *)
      сумма:=0;
      ОТ i:=0 ДО посл1 ВЫП
        сумма:=сумма+A[i,k]*A[i,j]
      КОН;
      R[k,j]:=сумма;
                              (* обновляем элементы столбца матрицы <q> (<A>) *)
      ОТ i:=0 ДО посл1 ВЫП
        A[i,j]:=A[i,j]-A[i,k]*R[k,j]
      КОН
    КОН
  КОН;

  ВОЗВРАТ 0;
КОН РазложитьНаQR;

(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА РешитьИзQR-(q-,R-:дМатрица.Вид;
                   b-,x+: дВектор.Вид);
(*
 * Цель:  решение системы уравнений после преобразования исходной матрицы
 *        с помощью РазложитьНаQR
 *----------------------------------------------------------------------------*
 * До:    <q> - матрица Q после РазложитьНаQR
 *        <R> - матрица R после РазложитьНаQR
 *        <b> - вектор свободных членов
 * После: <x> - вектор решения
 *)
ПЕР
  i,j,посл1,посл2:ЦЕЛ;
  сумма:Вещ;
УКАЗ
  посл1:=РАЗМЕР(q,0)-1;
  посл2:=РАЗМЕР(q,1)-1;
                                   (* строим q'b и записываем результат в <x> *)
  ОТ j:=0 ДО посл2 ВЫП
    x[j]:=0;
    ОТ i:=0 ДО посл1 ВЫП
      x[j]:=x[j]+q[i,j]*b[i]
    КОН
  КОН;
                                            (* обновляем <x> вектором решения *)
  x[посл2]:=x[посл2]/R[посл2,посл2];
  ОТ i:=(посл2-1) ДО 0 ПО -1 ВЫП
    сумма:=0;
    ОТ j:=(i+1) ДО посл2 ВЫП
      сумма:=сумма+R[i,j]*x[j]
    КОН;
    x[i]:=(x[i]-сумма)/R[i,i]
  КОН
КОН РешитьИзQR;

(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА ГильбертМатрица-(М+:дМатрица.Вид); 
(*
   Построение матрицы Гильберта <порядок>-го порядка 
     ( 1      1/2     1/3     1/4     ... 1/порядок      ) 
     ( 1/2    1/3     1/4     1/5     ... 1/(порядок+1)  ) 
 М = ( 1/3    1/4     1/5     1/6     ... 1/(порядок+2)  ) 
     ( ........................................... ) 
     ( 1/порядок    1/(порядок+1) 1/(порядок+2) 1/(порядок+3) ... 1/(2N-1) ) 
 *) 
ПЕР 
  i,j,посл:ЦЕЛ; 
УКАЗ 
  посл:=РАЗМЕР(М)-1;
                             (* первая строка *) 
  М[0,0]:=1; 
  ОТ j:=1 ДО посл ВЫП 
    М[0,j]:=1/(j+1) 
  КОН; 
                             (* оставшиеся строки *) 
  ОТ i:=1 ДО посл ВЫП 
                             (* последний столбец *) 
    М[i,посл]:=1/(посл+i+1); 
                             (* заполнить матрицу *) 
    ОТ j:=0 ДО посл-1 ВЫП 
      М[i,j]:=М[i-1,j+1]
    КОН 
  КОН 
КОН ГильбертМатрица; 
 
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА ГильбертСистема-(A+:дМатрица.Вид; b+:дВектор.Вид); 
(*
// Построение системы Гильберта n-го порядка
// 
//  A[0..n-1,0..n-1] = матрица системы: 
// 
//        ( 1      1/2     1/3     1/4     ... 1/n      ) 
//        ( 1/2    1/3     1/4     1/5     ... 1/(n+1)  ) 
//    A = ( 1/3    1/4     1/5     1/6     ... 1/(n+2)  ) 
//        ( ........................................... ) 
//        ( 1/n    1/(n+1) 1/(n+2) 1/(n+3) ... 1/(2N-1) ) 
// 
//  b[0..n-1] = вектор свободных членов: 
// 
//            n 
//    b[i] = Sum A[i,j] 
//           j=1 
// 
//  Тогда вектор решения будет x = [1 1 1 ... 1]
*) 
ПЕР 
  i,j,посл:ЦЕЛ; 
УКАЗ 
      посл:=РАЗМЕР(A)-1;

  дСистЛАУ.ГильбертМатрица(A); 
                                (* вектор постоянных членов *) 
  ОТ i:=0 ДО посл ВЫП 
    b[i]:=0; 
    ОТ j:=0 ДО посл ВЫП 
      b[i]:=b[i]+A[i,j]
    КОН 
  КОН 
КОН ГильбертСистема; 
 
КОН дСистЛАУ.
