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

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

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

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

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

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


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


(*----------------------------------------------------------------------------*)
ЗАДАЧА^ Якоби-(                    A+:дМатрица.Вид;
                 проходов:ЦЕЛ;
                             точность:Вещ;
               СобстВектСтроки_дМатр+:дМатрица.Вид;
                   собЗначений_дВект+: дВектор.Вид):ЦЕЛ;
(* Цель:  вычисление с. значений и с. векторов методом Якоби
 *----------------------------------------------------------------------------*
 * До:    <A>        - исходная матрица
 *        <проходов> - максимальное число проходов
 *        <точность> - требуемая точность
 * После: <СобстВектСтроки_дМатр> - матрица нормированных с. векторов (по строкам), 
 *               у которых первый элемент > 0
 *        <собЗначений_дВект> - с. значения
 *        <A>  - разрушенная исходная матрица
 * Ответ: 0 или дМатрица.НЕТ_СХОДИМОСТИ *)

(*----------------------------------------------------------------------------*)
ЗАДАЧА^ Значений_кВектор-(               A+:дМатрица.Вид;
                          дСобстЗнач_кВект+,
                          мСобстЗнач_кВект+: дВектор.Вид):ЦЕЛ;
(* Цель:  вычисление собственных значений квадратной матрицы
 *----------------------------------------------------------------------------*
 * До:    <A>     - исходная матрица
 * После: <дСобстЗнач_кВект> - действительная часть собственных значений
 *        <мСобстЗнач_кВект> -   мнимая       часть собственных значений
 *        <A>     - разрушенная исходная матрица
 * Ответ: 0 или дМатрица.НЕТ_СХОДИМОСТИ *)

(*----------------------------------------------------------------------------*)
ЗАДАЧА^ Вектор_1собЗначения-(          A- :дМатрица.Вид;
                              дСобЗначение,
                                  точность:Вещ;
                             дСобстВектор+: дВектор.Вид):ЦЕЛ;
(* Цель:  вычисление одного собственного вектора для соответствующего
 *        вещественного (кратного) собственного значения
 *----------------------------------------------------------------------------*
 * До:    <A>        - исходная матрица
 *        <сз>       - собственное значениеЭ для которого вычисляется <дСобстВектор+>
 *        <точность> - требуемая точность
 * После: <дСобстВектор> - нормированный собственный вектор, 
 *                         у которого первый элемент > 0
 * Ответ: 0 или дМатрица.НЕТ_СХОДИМОСТИ *)

(*----------------------------------------------------------------------------*)
ЗАДАЧА^ КорниПолинома-(         Коэфф-:дВектор.Вид; степень:ЦЕЛ;
                       дКорни+,мКорни+:дВектор.Вид):ЦЕЛ;
(* Цель:  вычисление вектора действительных и мнимых корней
 *        действительного полинома методом companion матрицы
 *----------------------------------------------------------------------------*
 * До:    <Коэфф>   - вектор коэффициентов действительного полинома
 *        <степень> - степень этого полинома
 * После: <дКорни> - действительные части корней (в порядке возрастания)
 *        <мКорни> - мнимые         части корней
 * Ответ: 0 или дМатрица.НЕТ_СХОДИМОСТИ *)


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


(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА Якоби-(                   A+ :дМатрица.Вид;
                проходов:ЦЕЛ;
                            точность:Вещ;
              СобстВектСтроки_дМатр+:дМатрица.Вид;
                  собЗначений_дВект+: дВектор.Вид):ЦЕЛ;
(* Цель:  вычисление с. значений и с. векторов методом Якоби
 *----------------------------------------------------------------------------*
 * До:    <A>        - исходная матрица
 *        <проходов> - наибольшее число проходов
 *        <точность> - требуемая точность
 * После: <СобстВектСтроки_дМатр> - матрица нормированных с. векторов (по строкам), 
 *               у которых первый элемент > 0
 *        <собЗначений_дВект> - с. значения
 *        <A>  - разрушенная исходная матрица
 * Ответ: 0 или дМатрица.НЕТ_СХОДИМОСТИ *)
ПЕР
  синус,косинус,тангенс,tg2тета:Вещ;
  квCos,квSin,СинКос,сумКвДиаг:Вещ;
  Aii,Ajj,Aij,Aik,Ajk,
  СобстВектСтроки_дМатрik,
  СобстВектСтроки_дМатрjk,d:Вещ;
  i,j,k,проход,посл:ЦЕЛ;
  готово:КЛЮЧ;
УКАЗ
  посл:=РАЗМЕР(A)-1;
  проход:=0;
  ОТ i:=0 ДО посл ВЫП
    ОТ j:=0 ДО посл ВЫП
      ЕСЛИ i = j ТО
        СобстВектСтроки_дМатр[i,j]:=1
      ИНАЧЕ
        СобстВектСтроки_дМатр[i,j]:=0 
      КОН 
    КОН 
  КОН;

  ПОВТОРЯТЬ
    готово:=ВКЛ;
    УВЕЛИЧИТЬ(проход);
    сумКвДиаг:=0;
    ОТ i:=0 ДО посл ВЫП
      сумКвДиаг:=сумКвДиаг+дМатем.кв(A[i,i])
    КОН;

    ОТ i:=0 ДО посл-1 ВЫП
      ОТ j:=i+1 ДО посл ВЫП
        ЕСЛИ МОДУЛЬ(A[i,j]) > точность*сумКвДиаг ТО
          готово:=ОТКЛ;
          (* вычисляем поворот *)
          d:=A[i,i]-A[j,j];
          ЕСЛИ МОДУЛЬ(d) > дМатем.МашинEps ТО
            tg2тета:=d/(2*A[i,j]);
            тангенс:=дМатем.знакЧ(tg2тета) * дМатем.квкор(1+дМатем.кв(tg2тета)) - tg2тета;
            косинус:=1/дМатем.квкор(1+дМатем.кв(тангенс));
              синус:=тангенс*косинус
          ИНАЧЕ
            косинус:=дМатем.ПолКвКОРиз2;                   (* квкор(2)/2 *)
              синус:=дМатем.знакЧ(A[i,j])*дМатем.ПолКвКОРиз2
          КОН;
          (* поворачиваем матрицу *)
          квCos:=дМатем.кв(косинус);
          квSin:=дМатем.кв(  синус);
          СинКос:=синус*косинус;
          Aii:=A[i,i]*квCos+2*A[i,j]*СинКос+A[j,j]*квSin;
          Ajj:=A[i,i]*квSin-2*A[i,j]*СинКос+A[j,j]*квCos;
          Aij:=(A[j,j]-A[i,i])*СинКос + A[i,j]*(квCos-квSin);
          ОТ k:=0 ДО посл ВЫП
            ЕСЛИ (k # i) И (k # j) ТО
              Aik:= A[i,k]*косинус + A[j,k]*синус;
              Ajk:=-A[i,k]*синус   + A[j,k]*косинус;
              A[i,k]:=Aik;
              A[k,i]:=Aik;
              A[j,k]:=Ajk;
              A[k,j]:=Ajk
            КОН 
          КОН;
          A[i,i]:=Aii;
          A[j,j]:=Ajj;
          A[i,j]:=Aij;
          A[j,i]:=Aij;

          (* поворачиваем с. вектора *)
          ОТ k:=0 ДО посл ВЫП
            СобстВектСтроки_дМатрik:= косинус*СобстВектСтроки_дМатр[i,k] +   синус*СобстВектСтроки_дМатр[j,k];
            СобстВектСтроки_дМатрjk:=  -синус*СобстВектСтроки_дМатр[i,k] + косинус*СобстВектСтроки_дМатр[j,k];
            СобстВектСтроки_дМатр[i,k]:=СобстВектСтроки_дМатрik;
            СобстВектСтроки_дМатр[j,k]:=СобстВектСтроки_дМатрjk
          КОН
        КОН
      КОН
    КОН
  ДО готово ИЛИ (проход > проходов);

  (* на диагонали преобразованной матрицы находятся с. значения *)
  ОТ i:=0 ДО посл ВЫП
    собЗначений_дВект[i]:=A[i,i]
  КОН;

  ЕСЛИ проход > проходов ТО
    ВОЗВРАТ дМатрица.НЕТ_СХОДИМОСТИ
  КОН;

  (* упорядочивание с. значений и с. векторов *)
  ОТ i:=0 ДО посл-1 ВЫП
    k:=i;
    d:=собЗначений_дВект[i];
    ОТ j:=i+1 ДО посл ВЫП
      ЕСЛИ собЗначений_дВект[j] > d ТО
        k:=j;
        d:=собЗначений_дВект[j]
      КОН 
    КОН;
    дМатем.обмен(собЗначений_дВект[i],собЗначений_дВект[k]);
    дМатрица.ОбменСтрок(СобстВектСтроки_дМатр,i,k)
  КОН;

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

  ВОЗВРАТ 0
КОН Якоби;

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

(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА ElmHes(A+:дМатрица.Вид);
(* исключение матрицы до верхней формы Гезенберга *)
ПЕР
  i,j,m,посл:ЦЕЛ;
  x,y:Вещ;
УКАЗ
  посл:=РАЗМЕР(A)-1;
  ОТ m:=1 ДО посл-1 ВЫП
    x:=0;
    i:=m;
    ОТ j:=m ДО посл ВЫП
      ЕСЛИ МОДУЛЬ(A[j,m-1]) > МОДУЛЬ(x) ТО
        x:=A[j,m-1];
        i:=j
      КОН 
    КОН;
    ЕСЛИ i # m ТО
      ОТ j:=m-1 ДО посл ВЫП
        дМатем.обмен(A[i,j],A[m,j]) 
      КОН;
      ОТ j:=0 ДО посл ВЫП
        дМатем.обмен(A[j,i],A[j,m]) 
      КОН
    КОН;
    ЕСЛИ x # 0 ТО
      ОТ i:=m+1 ДО посл ВЫП
        y:=A[i,m-1];
        ЕСЛИ y # 0 ТО
          y:=y/x;
          A[i,m-1]:=y;
          ОТ j:=m ДО посл ВЫП
            A[i,j]:=A[i,j]-y*A[m,j]
          КОН;
          ОТ j:=0 ДО посл ВЫП
            A[j,m]:=A[j,m]+y*A[j,i]
          КОН
        КОН
      КОН 
    КОН
  КОН;
  ОТ i:=2 ДО посл ВЫП
    ОТ j:=0 ДО i-2 ВЫП
      A[i,j]:=0 
    КОН 
  КОН
КОН ElmHes;

(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА Hqr(               A+:дМатрица.Вид;
           дСобстЗнач_кВект+,
           мСобстЗнач_кВект+: дВектор.Вид):ЦЕЛ;
(* находит собственные значения у верхней матрицы Гезенберга *)
ПЕР
  i,Its,j,k,l,m,n,посл:ЦЕЛ;
  Anorm,p,q,r,s,t,u,v,w,x,y,z:Вещ;

  ЗАДАЧА Sign(a,b:Вещ):Вещ;
  УКАЗ
    ЕСЛИ b < 0 ТО
      ВОЗВРАТ -МОДУЛЬ(a)
    ИНАЧЕ
      ВОЗВРАТ  МОДУЛЬ(a)
    КОН
  КОН Sign;

УКАЗ
  посл:=РАЗМЕР(A)-1;
  (* вычисление нормы *)
  k:=0;
  Anorm:=0;
  ОТ i:=0 ДО посл ВЫП
    ОТ j:=k ДО посл ВЫП
      Anorm:=Anorm+МОДУЛЬ(A[i,j]) 
    КОН;
    k:=i 
  КОН;
  n:=посл;
  t:=0;
  Its:=0;
  КОЛЬЦО
(* 60: *)
(* 70: *)
    (* поиск малого одиночного поддиагонального элемента *)
    l:=n;
    КОЛЬЦО
      ЕСЛИ l < 1 ТО
        ВЫХОД
      КОН;
      s:=МОДУЛЬ(A[l-1,l-1])+МОДУЛЬ(A[l,l]);
      ЕСЛИ s = 0 ТО s:=Anorm КОН;
      ЕСЛИ s+МОДУЛЬ(A[l,l-1]) = s ТО 
        ВЫХОД
      КОН;
      УМЕНЬШИТЬ(l)
    КОН;
(* 100: *)
    ЕСЛИ l = n ТО
(* 270: *)
      (* нашли однократный корень *)
      дСобстЗнач_кВект[n]:=A[n,n]+t;
      мСобстЗнач_кВект[n]:=0;
      n:=n-1;
      ЕСЛИ n < 0 ТО
        ВОЗВРАТ 0
      КОН;
      Its:=0
      (* goto 60 *)
    ИНАЧЕ
      (* формируем сдвиг *)
      x:=A[n,n];
      y:=A[n-1,n-1];
      w:=A[n,n-1]*A[n-1,n];
      ЕСЛИ l = n-1 ТО
(* 280: *)
        (* нашли кратный корень *)
        p:=0.5*(y-x);
        q:=дМатем.кв(p)+w;
        z:=дМатем.квкор(МОДУЛЬ(q));
        x:=x+t;
        ЕСЛИ q >= 0 ТО
          (* действительноя пара *)
          z:=p+Sign(z,p);
          дСобстЗнач_кВект[n]:=x+z;
          дСобстЗнач_кВект[n-1]:=дСобстЗнач_кВект[n];
          ЕСЛИ z # 0 ТО дСобстЗнач_кВект[n]:=x-w/z КОН;
          мСобстЗнач_кВект[n]:=0;
          мСобстЗнач_кВект[n-1]:=0
        ИНАЧЕ
(* 320: *)
          (* комплексная пара *)
          дСобстЗнач_кВект[n]  :=x+p;
          дСобстЗнач_кВект[n-1]:=дСобстЗнач_кВект[n];
          мСобстЗнач_кВект[n]  :=-z;
          мСобстЗнач_кВект[n-1]:=z
        КОН;
(* 330: *)
        n:=n-2;
        ЕСЛИ n < 0 ТО
          ВОЗВРАТ 0
        КОН;
        Its:=0
        (* goto 60 *)
      ИНАЧЕ
        УВЕЛИЧИТЬ(Its);
        ЕСЛИ (Its = 10) ИЛИ (Its = 20) ТО
          (* формируем исключительный сдвиг *)
          t:=t+x;
          ОТ i:=0 ДО n ВЫП
            A[i,i]:=A[i,i]-x 
          КОН;
          s:=МОДУЛЬ(A[n,n-1])+МОДУЛЬ(A[n-1,n-2]);
          x:=0.75*s;
          y:=x;
          w:=-0.4375*дМатем.кв(s)
        АЕСЛИ Its = 30 ТО
          ВОЗВРАТ дМатрица.НЕТ_СХОДИМОСТИ
        КОН;
(* 130: *)
        (* поиск двух малых смежных поддиагональных элементов *)
        m:=n-2;
        КОЛЬЦО
          ЕСЛИ m < l ТО
            ВЫХОД
          КОН;
          z:=A[m,m];
          r:=x-z;
          s:=y-z;
          p:=(r*s-w)/A[m+1,m]+A[m,m+1];
          q:=A[m+1,m+1]-z-r-s;
          r:=A[m+2,m+1];
          s:=МОДУЛЬ(p)+МОДУЛЬ(q)+МОДУЛЬ(r);
          p:=p/s;
          q:=q/s;
          r:=r/s;
          ЕСЛИ m = l ТО 
            ВЫХОД
          КОН;
          v:=МОДУЛЬ(p)*(МОДУЛЬ(A[m-1,m-1])+МОДУЛЬ(z)+МОДУЛЬ(A[m+1,m+1]));
          u:=МОДУЛЬ(A[m,m-1])*(МОДУЛЬ(q)+МОДУЛЬ(r));
          ЕСЛИ u+v=v ТО 
            ВЫХОД
          КОН;
          УМЕНЬШИТЬ(m)
        КОН;
(* 150: *)
        ОТ i:=m+2 ДО n ВЫП
          A[i,i-2]:=0;
          ЕСЛИ i # (m+2) ТО 
            A[i,i-3]:=0
          КОН 
        КОН;
        (* двойной QR шаг разрешая строки l на n и столбцы m на n *) 
        ОТ k:=m ДО n-1 ВЫП
          ЕСЛИ k # m ТО
            p:=A[k,k-1];
            q:=A[k+1,k-1];
            r:=0;
            ЕСЛИ k # (n-1) ТО
              r:=A[k+2,k-1]
            КОН;
            x:=МОДУЛЬ(p)+МОДУЛЬ(q)+МОДУЛЬ(r);
            ЕСЛИ x # 0 ТО
              p:=p/x;
              q:=q/x;
              r:=r/x
            КОН 
          КОН;
(* 170: *)
          s:=Sign(дМатем.квкор(дМатем.кв(p)+дМатем.кв(q)+дМатем.кв(r)),p);
          ЕСЛИ s # 0 ТО
            ЕСЛИ k = m ТО
              ЕСЛИ l # m ТО
                A[k,k-1]:=-A[k,k-1]
              КОН
            ИНАЧЕ
              A[k,k-1]:=-s*x
            КОН;
            p:=p+s;
            x:=p/s;
            y:=q/s;
            z:=r/s;
            q:=q/p;
            r:=r/p;
            (* изменение строк *)
            ОТ j:=k ДО n ВЫП
              p:=A[k,j]+q*A[k+1,j];
              ЕСЛИ k # (n-1) ТО
                p:=p+r*A[k+2,j];
                A[k+2,j]:=A[k+2,j]-p*z
              КОН;
              A[k+1,j]:=A[k+1,j]-p*y;
              A[k,j]:=A[k,j]-p*x
            КОН;
            (* изменение столбцов *)
            ОТ i:=l ДО цМатем.мин(n,k+3) ВЫП
              p:=x*A[i,k]+y*A[i,k+1];
              ЕСЛИ k # (n-1) ТО
                p:=p+z*A[i,k+2];
                A[i,k+2]:=A[i,k+2]-p*r
              КОН;
              A[i,k+1]:=A[i,k+1]-p*q;
              A[i,k]:=A[i,k]-p
            КОН 
          КОН 
        КОН
      КОН 
    КОН 
  КОН
КОН Hqr;

(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА Значений_кВектор-(               A+:дМатрица.Вид;
                         дСобстЗнач_кВект+,
                         мСобстЗнач_кВект+: дВектор.Вид):ЦЕЛ;
(* Цель:  вычисление собственных значений квадратной матрицы
 *----------------------------------------------------------------------------*
 * До:    <A>       - исходная матрица
 * После: <дСобстЗнач_кВект> - действительная часть собственных значений
 *        <мСобстЗнач_кВект> -   мнимая       часть собственных значений
 *        <A>       - разрушенная исходная матрица
 * Ответ: 0 или дМатрица.НЕТ_СХОДИМОСТИ *)
УКАЗ
  Balance(A);
  ElmHes(A);
  ВОЗВРАТ Hqr(A,дСобстЗнач_кВект,мСобстЗнач_кВект)
КОН Значений_кВектор;

(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА Вектор_1собЗначения-(          A- :дМатрица.Вид;
                             дСобЗначение,
                                 точность:Вещ;
                            дСобстВектор+: дВектор.Вид):ЦЕЛ;
(* Цель:  вычисление одного собственного вектора для соответствующего
 *        вещественного (кратного) собственных значения
 *----------------------------------------------------------------------------*
 * До:    <A>        - исходная матрица
 *        <дСобЗначение> - собственное значение квадратной матрицы
 *        <точность> - требуемая точность
 * После: <дСобстВектор> - нормированный собственный вектор, 
 *                         у которого первый элемент > 0
 * Ответ: 0 или дМатрица.НЕТ_СХОДИМОСТИ *)
ПЕР
  код,i,посл:ЦЕЛ;
  A1:дМатрица.Доступ;

(*===*)
  ЗАДАЧА SetMatrix(A-,A1+:дМатрица.Вид; дСобЗначение:Вещ);
  (* строит A1 = A-дСобЗначение*i *)
  ПЕР
    i:ЦЕЛ;
  УКАЗ
    дМатрица.Списать(A,A1);
    ОТ i:=0 ДО РАЗМЕР(A)-1 ВЫП
      A1[i,i]:=A[i,i]-дСобЗначение 
    КОН
  КОН SetMatrix;

(*===*)
  ЗАДАЧА Solve(          A- :дМатрица.Вид;
               n:ЦЕЛ;
                    точность:Вещ;
               дСобстВектор+: дВектор.Вид):ЦЕЛ;
  (* решает систему A*x = 0 после закрепления за n-м неизвестным 1 *)
  ПЕР
    A1, W:дМатрица.Доступ;
    b,s,x: дВектор.Доступ;
    код,i,i1,j,j1,посл:ЦЕЛ;
  УКАЗ
    посл:=РАЗМЕР(A)-1;
    СОЗДАТЬ(A1,посл,посл);
    СОЗДАТЬ(W,посл,посл);
    СОЗДАТЬ(b,посл);
    СОЗДАТЬ(s,посл);
    СОЗДАТЬ(x,посл);

    i1:=0;
    ОТ i:=0 ДО посл ВЫП
      ЕСЛИ i # n ТО
        j1:=0;
        ОТ j:=0 ДО посл ВЫП
          ЕСЛИ j # n ТО
            A1[i1,j1]:=A[i,j];
            УВЕЛИЧИТЬ(j1)
          КОН
        КОН;
        b[i1]:=-A[i,n];
        УВЕЛИЧИТЬ(i1)
      КОН 
    КОН;

    код:=дСистЛАУ.РазложитьНаSV(A1^,s^,W^);

    ЕСЛИ код = 0 ТО
      дСистЛАУ.ОбнулитьS(s^,точность);
      дСистЛАУ.РешитьИзSV(A1^,s^,W^,b^,x^);

      (* обновить собственный вектор *)
      i1:=0;
      ОТ i:=0 ДО посл ВЫП
        ЕСЛИ i = n ТО
          дСобстВектор[i]:=1
        ИНАЧЕ
          дСобстВектор[i]:=x[i1];
          УВЕЛИЧИТЬ(i1)
        КОН 
      КОН
    КОН;
    ВОЗВРАТ код;
  КОН Solve;

(*===*)
  ЗАДАЧА ZeroVector(b-:дВектор.Вид; точность:Вещ):КЛЮЧ;
  (* ВКЛ, если вектор b 0-й *)
  ПЕР
    i:ЦЕЛ;
  УКАЗ
    ОТ i:=0 ДО РАЗМЕР(b)-1 ВЫП
      ЕСЛИ МОДУЛЬ(b[i]) >= точность ТО
        ВОЗВРАТ ОТКЛ
      КОН
    КОН;
    ВОЗВРАТ ВКЛ
  КОН ZeroVector;

(*===*)
  ЗАДАЧА CheckEigenVector(          A1-:дМатрица.Вид;
                          дСобстВектор-: дВектор.Вид; точность:Вещ):КЛЮЧ;
  (* ВКЛ, если выполняется равенство A1*дСобстВектор = 0 *)
  ПЕР
    i,k,посл:ЦЕЛ;
    bД:дВектор.Доступ;
  УКАЗ
    посл:=РАЗМЕР(A1)-1;
    СОЗДАТЬ(bД,посл+1);
    (* строим bД = A1*дСобстВектор *)
    ОТ i:=0 ДО посл ВЫП
      ОТ k:=0 ДО посл ВЫП
        bД[i]:=bД[i]+A1[i,k]*дСобстВектор[k]
      КОН 
    КОН;
    ВОЗВРАТ ZeroVector(bД^,точность)
  КОН CheckEigenVector;

(*===*)
  ЗАДАЧА Normalize(дСобстВектор+:дВектор.Вид);
  (* нормируем собственный вектор и устанавливаем первый элемент >= 0 *)
  ПЕР
    сумма,Norm:Вещ;
    i,посл:ЦЕЛ;
  УКАЗ
    посл:=РАЗМЕР(дСобстВектор)-1;
    сумма:=0;
    ОТ i:=0 ДО посл ВЫП
      сумма:=сумма+дМатем.кв(дСобстВектор[i]) 
    КОН;
    Norm:=дМатем.квкор(сумма);
    ОТ i:=0 ДО посл ВЫП
      дСобстВектор[i]:=дСобстВектор[i]/Norm 
    КОН;
    ЕСЛИ дСобстВектор[0] < 0 ТО
      ОТ i:=0 ДО посл ВЫП
        дСобстВектор[i]:=-дСобстВектор[i]
      КОН 
    КОН
  КОН Normalize;

(*===*)
УКАЗ
  посл:=РАЗМЕР(A)-1;
  СОЗДАТЬ(A1,посл+1,посл+1);

 (* строим A1 = A-дСобЗначение*i *)
  SetMatrix(A,A1^,дСобЗначение);

  (* попытаться решить систему A1*дСобстВектор=0 исключая 1 уравнение *)
  ОТ i:=0 ДО посл ВЫП
    ЕСЛИ (Solve(A1^,i,точность,дСобстВектор) = 0)
        И
           CheckEigenVector(A1^,дСобстВектор,точность)
    ТО
      Normalize(дСобстВектор);
      ВОЗВРАТ 0
    КОН
  КОН;
  ВОЗВРАТ дМатрица.НЕТ_СХОДИМОСТИ 
КОН Вектор_1собЗначения;

(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА КорниПолинома-(         Коэфф-:дВектор.Вид; степень:ЦЕЛ;
                      дКорни+,мКорни+:дВектор.Вид):ЦЕЛ;
(* Цель:  вычисление вектора действительных и мнимых корней
 *        действительного полинома методом companion матрицы
 *----------------------------------------------------------------------------*
 * До:    <Коэфф>   - вектор коэффициентов действительного полинома
 *        <степень> - степень этого полинома
 * После: <дКорни> - действительные части корней (в порядке возрастания)
 *        <мКорни> - мнимые         части корней
 * Ответ: 0 или дМатрица.НЕТ_СХОДИМОСТИ *)
ПЕР
  n:ЦЕЛ;                (* размер матрицы    *)
  i,j,k:ЦЕЛ;  
  код:ЦЕЛ;
  врем:Вещ;
  companionМатД:дМатрица.Доступ;   (* companion матрица nxn *)
УКАЗ
  n:=степень-1;
  СОЗДАТЬ(companionМатД,n+1,n+1);

  (* заполняем companion матрицу *)
  ОТ j:=0 ДО n ВЫП
    companionМатД[0,j]:=-Коэфф[степень-j-2]/Коэфф[степень-1]
  КОН;
  ОТ j:=0 ДО n-1 ВЫП
    companionМатД[j+1,j]:=1 
  КОН;

  (* корни полинома являются и с. значениями companion матрицы *)
  Balance(companionМатД^);
  код := Hqr(companionМатД^,дКорни,мКорни);

  ЕСЛИ код = 0 ТО
    (* упорядочивание корней в порядке возрастания действительных частей *)
    ОТ i:=0 ДО n-1 ВЫП
      k:=i;
      врем:=дКорни[i];
      ОТ j:=i+1 ДО n ВЫП
        ЕСЛИ дКорни[j] < врем ТО
          k:=j;
          врем:=дКорни[j]
        КОН 
      КОН;
      дМатем.обмен(дКорни[i],дКорни[k]);
      дМатем.обмен(мКорни[i],мКорни[k])
    КОН
  КОН;
  ВОЗВРАТ код
КОН КорниПолинома;

КОН дСобств.
