My source

home

Допустимо свободное использование материалов данного раздела, в том числе и для создания законченных коммерческих программ с указанием источника. http://drobotenko.com/ . Для написания платных книг, библиотек и плагинов требуется разрешение автора.

Метод формальных манипуляций для обращения матриц.

Линейные системы основа всех численных методов. Решение линейных систем проходят еще в школе. Предлагаемый алгоритм обращения матриц отличается просто более тщательной группировкой членов и совмещением операций. Как и LU разложение он обращает матрицу “на месте” , без дополнительной памяти, но более быстрый. Все операции выполняются в одном цикле, нет декомпозиций, прямого-обратного прохода. Число операций ровно N³ - столько же сколько и умножение матриц. Коды много лет используются в коммерческих проектах.

Примерный алгоритм

Заводятся 2 целых вектора длинной N. Векторы расписываются номерами строк, столбцов.
Среди всех непомеченных строк, столбцов ищется максимальный элемент M_pq
Вычисления
- Для всех элементов
  M_jk-= M_pk*M_jq/M_pq
- Для строк и столбцов ведущего элемета
  M_jq/= M_pq
  M_pj/= M_pq
- сам ведущий
  M_pq= -1/ M_pq
строка и столбец ведущего помечаются как обработанные (значение -1) и производятся перестановки
Если остались необработанные - возврат к новому поиску ведущего

В результате, на месте матрицы, получается обратная матрица с обратным знаком. Помимо отсутствия обратного хода этот алгоритм отличается еще и тем, что не делает различия между строками и столбцами. При перестановках строка меняется с номером столбца, а столбец с номером строки, при этом используются исходные номера. Номера приписываются строкам/столбцам изначально в двух массивах. Строка идет при перестановке на истинный номер столбца и остается там ,получая номер -1. Вытесняемая строка/столбец идут на ее место и пока располагаются по несвоему номеру. Чтоб пояснить перестановки возьмем предельный случай, в каждой строке и столбце по одному ненулевому элементу, все они ведущие. Для прямого преобразования i->a_ji*j . Обратное j->1/a_ji*i . Значит 1/a_ji в обратной матрице стоит на месте ij. Отсюда ясно где, в конечном счете, должна оказаться строка и столбец ведущего элемента. Вычисления в текстах эквивалентны, но несколько соптимизированы.

В этих кодах выполняется полный выбор максимального элемента. Исключительно из за этого он несколько медленнее LU с частичным выбором ведущего, иначе был бы быстрее. Есть примеры, где частичный выбор ведущего срабатывает плохо. И я не вижу универсального критерия для частичного выбора, например при симметричной матрице ведущий нужно искать на диагонали, и поэтому частичный выбор ведущего - это уже специализация зависящая от задачи. Примерно половина времени уходит на поиск ведущего при полном переборе, и при этом скорость процентов на 30 медленнее чем у алгоритмов с LU с частичным выбором.

Обоснование методов – элементарные преобразования + приписывание параллельной единичной матрицы. При элементарных преобразованиях зануляется элемент в исходной матрице и образуется новый в единичной. Это позволяет использовать только одну матрицу. До сих пор известный вариант совмещения этого LU разложение. Но имеются 2 стратегии – либо занулять только ниже диагонали как в LU, либо занулять сразу весь столбец. Как оказывается, мой вариант более производителен в случае обращения матриц. Для решения линейных уравнений рациональней двухпроходный метод. Но LU-разложение проигрывает обычному методу Гаусса с обратным ходом из-за абстрактности метода. Производятся операции над коэффициентами, вместо конкретных чисел. Соотношение 1/3N³ к 2/3N³ на LU-разложение (только LU, без обращения матрицы).

Здесь классическое LU-обращение матриц ( netlib.org ) - помимо использования нерационального LU, эти коды весьма громоздки и обращаются к огромному количеству низкоуровневых процедур написанных на ассемблере. Но современные компиляторы обычно не нуждаются в этом. Мат библиотеки от Intel и AMD отличаются большей скоростью. Они неплохо работают с большими матрицами, которые не помещаются в кэш. Основная операция это вычитание из строки другой строки, помноженной на число, которое в регистре с обратной запись результата. Выпадение строк из кэша, а это несколько мегабайт в зависимости от процессора приводит к замедлению на порядок. Это на одно ядро, тоесть для четырехядерных можно сказать, что производительность падает в сорок раз. Библиотеки от Intel и AMD неплохо справляются с этой проблемой. Чтоб уменьшить замедление нужно за один проход по матрице делать как можно больше операций над одной строкой. Можно искать следующий разрешающий элемент параллельно с преобразованиями и вычитать сразу несколько строк. Правда, при вычитание нескольких строк за один проход невозможен поиск ведущего по всей матрице. Возможно в каких-то случаях использовать сравнение кодов вместо сравнения плавающих чисел при поиске ведущего. Если удалить знаковый бит, то поскольку порядок в старшей части результат плавающего сравнения такой же как и целого, но вряд ли это актуально. Вот оптимизированное обращение матрицы корреляций, которая симметрична а ведущие элементы положительны и на главной диагонали. Используется верхняя часть, потом переписывается в нижнюю CovarianceMtrix

В этих кодах производится полный разбор всех случаев. И обращение, и решение СЛАУ возвращают дефект линейного преобразования, детерминант и частичное решение в том случае если детерминант равен 0. Частичное решение - это решение на подпространстве, где детерминант еще не равен 0. Это частичное решение, например, подходит для метода наименьших квадратов, даже в том случае если ищется приближение функциями, которые линейно зависимы.

Здесь две версии кодов. Шаблонная версия кодов отличается тем, что может работать с родными массивами c++ и любыми контейнерами. При этом допустимо их смешение в одном вызове. Для подключения любого контейнера достаточно дописать пару строчек. Реализован оператор определения типа в шаблонах – typeof. Это решение соответствует последним стандартам и проверено в VC6 – VC8. Версия double содержит удобные контейнеры, также совместима с простыми массивами.



#include <math.h>
int gaus_ostatok;
// дефект линейного преобразования
// 0 == ОК

long double gaus_deter;
//  в случае успешного обращения детерминант

double gaus_minved;
//  минимальный ведущий элемент
//  можно использовать для оценки точности


#define _NaN() (double&)*"Implementation dependent";
//  заполнитель для неразрешенных измерений
//  при невозможности обращения

     template<int N>
     inline static  void  invers(double m[N][N])
     {
       //  (c) Drobotenko    http://drobotenko.com

       int rn[N],cn[N];
       int j,k;

       for(j=N;j--;)
         rn[j]=cn[j]=j;
       gaus_minved=1e99;
       gaus_deter=1;
       for(gaus_ostatok=N;gaus_ostatok;gaus_ostatok--)
       {  int jved,kved;
          double vved=-1,t;

          // поиск ведущего
          for(j=N;j--;)
          {  if(~rn[j])
             for(k=N;k--;)
               if(~cn[k])
               if(vved<fabs(m[j][k]))
                  vved=fabs(m[j][k]),jved=j,kved=k;
          }
          
          if(gaus_minved>vved)
             gaus_minved=vved;
          gaus_deter*=m[jved][kved];
          if(vved<1e-99)
          {  for(j=N;j--;)
             {  if(~rn[j]) for(k=N;k--;)
                  m[j][k]=_NaN();
                if(~cn[j]) for(k=N;k--;)
                  m[k][j]=_NaN();
             }
             return;
          }

          int jt=rn[jved],kt=cn[kved];


          
          if(jt!=kved)
            gaus_deter=-gaus_deter;
          // меняем знак детерминанта если перестановка
          



          // перестановки
          for(j=N;j--;)
              t=m[kt][j],m[kt][j]=m[jved][j],m[jved][j]=t;
          for(j=N;j--;)
              t=m[j][jt],m[j][jt]=m[j][kved],m[j][kved]=t;

          rn[jved]=rn[kt];
          cn[kved]=cn[jt];
          rn[kt]=cn[jt]=-1;

          vved=m[kt][jt];   m[kt][jt]=1;
          for(j=N;j--;)
          {  if(j==kt)
               continue;
             double mul=m[j][jt]/vved;
             m[j][jt]=0;
             for(k=N;k--;)
               m[j][k]-=m[kt][k]*mul;
               // самый внутренний цикл ровно N^3 операций
          }
          for(k=N;k--;)
              m[kt][k]/=vved;
       }
     }



const int r=5;
double m[r][r],e[r][r];
double mval(int i,int j)
 {  // значения матрицы
    return int(1000./(.1+fabs(i-.95*j)))/1000.;
 }
int main(int, char **)
{  for(int i=r;i--;)
   for(int j=r;j--;)
     m[i][j]=mval(i,j);

   invers(m);

   for(int i=r;i--;)
   for(int j=r;j--;)
   { // проверка в Е получается единичная матрица
      e[i][j]=0;
      for(int k=r;k--;)
        e[i][j]+=m[i][k]*mval(k,j);
   }
}

Если компилятор не может справится с этим текстом, попробуйте заменить параметр на ссылку void invers(double (&m)[N][N]) , или явно укажите параметр шаблона при вызове invers<r>(m); можно наконец в параметрах шаблона задать и тип матрицы

template<int N, class C>
inline static void invers(C &m)
{
}

....................
invers<r,double [r][r]>(m);

возможно это с любым компилятором работать будет, хоть каждый вызов указывать параметры шаблона неудобно - поддержка стандарта не у всех компиляторов на высоте :(
т.к. резмерность массива можно вывести через sizeof удобно написать макрос
#define INVERS(M) invers<sizeof((M)[0])/sizeof((M)[0][0])>(M)

версия c контейнерами VC6-VC8

КОДЫ - какаято старая версия