C++: считаем определитель матрицы

Вообще-то программировать расчёт определителей не нужно. Их умеет считать, скажем, встроенная функция МОПРЕД из Excel:

• набираем элементы матрицы в смежных ячейках, например, матрица размерностью 4*4 показана на картинке;
• в нужной ячейке вводим формулу (в нашем случае =МОПРЕД(A1:D4) и нажимаем Enter :)

Не труднее вычислить и в MathCAD - просто нажать кнопку на панели матриц...

Но иногда нужен алгоритм, а не ответ... вот немного кода на консольном C++, на совершенство он не претендует, но нули в матрице или "не тот" порядок элементов смущать функцию determinant не должны. Пример из main - 1001-й на работу с динамической матрицей средствами C++ :) Остальное закомментировано в исходнике.

#define bool int
#define true 1
#define false 0

int search (double **a, int m, int n, double what, 
 bool match, unsigned int &uI, unsigned int &uJ, unsigned int starti, unsigned int startj) {
 // Поиск в матрице a[m][n] элемента с указанным значением what
 // Возвращаеются его номер строки и столбца uI, uJ, если элемент найден.
 // match - искать равный элемент или отличный от указанного.
 // Вернёт 0 - не найдено, не 0 - найдено
 if ((!m) || (!n)) return 0;
 if ((starti >= n) || (startj >= m)) return 0;
 for (unsigned int i = starti; i < n; i++)
 for (unsigned int j = startj; j < m; j++) {
  if (match == true) {
   if (a[i][i] == what) {
    uI = i; uJ = j; return 1;
   }
  }
  else if (a[i][j] != what) {
   uI = i; uJ = j; return 1;
  }
 }
 return 0;
}

void swaprows (double **a, int n, int m, unsigned int x1, unsigned int x2) {
 //Меняет в матрице a[n][m] строки с номерами x1 и x2 местами
 if ((!n) || (!m)) return; 
 if ((x1 >= n) || (x2 >= n) || (x1 == x2)) return;
  double tmp;
  for (unsigned int x = 0; x < m; x++) {
    tmp = a[x1][x];
    a[x1][x] = a[x2][x];
    a[x2][x] = tmp;
  }
  return;
};

void swapcolumns (double **a, int n, int m, unsigned int x1, unsigned int x2) {
 //Меняет в матрице a[n][m] столбцы с номерами x1 и x2 местами
 if ((!n) || (!m)) return; 
 if ((x1 >= m) || (x2 >= m) || (x1 == x2)) return;
  double tmp;
  for (unsigned int x = 0; x < n; x++) {
    tmp = a[x][x1];
    a[x][x1] = a[x][x2];
    a[x][x2] = tmp;
  }
  return;
};

double determinant (double **a, unsigned int n) {
 //Вычисление определителя квадратной матрицы a[n][n]
 unsigned int m = n;
 if (m == 0) return 0;
 if (m == 1) return a[0][0];
 if (m == 2) return (a[0][0] * a[1][1] - a[1][0] * a[0][1]);
 bool sign = false; // смена знака определителя. по умолчанию - нет
 double det = 1; // определитель
 double tmp;
 unsigned int x, y;
 for (unsigned int i = 0; i < n; i++) { // цикл по всей главной диагонали
  if (a[i][i] == 0) { // если элемент на диагонали равен 0, то ищем ненулевой элемент в матрице
   if (!search(a,m,n,0, false, y, x, i, i)) return 0; // если все элементы нулевые, то опр. = 0
   if (i != y) { // меняем i-ую строку с y-ой
    swaprows(a,m,n,i, y);
    sign = !sign;
   }
   if (i != x) { // меняем i-ый столбец с x-ым
    swapcolumns(a,m,n,i, x);
    sign = !sign;
   }
   // таким образом, в a[i][i], теперь ненулевой элемент.
  }
  // выносим элемент a[i][i] за определитель
  det *= a[i][i];
  tmp = a[i][i];
  for (x = i; x < m; x++) {
   a[i][x] = a[i][x] / tmp;
  }
  // таким образом a[i][i] теперь равен 1
  // зануляем все элементы стоящие под (i, i)-ым,
  // при помощи вычитания с опр. коеффициентом
  for (y = i + 1; y < n; y++) {
   tmp = a[y][i];
   for (x = i; x < m; x++)
    a[y][x] -= (a[i][x]*tmp);
   }
  }
  if (sign) return det*(-1);
  return det;
};

#include <stdio.h>

int main () {
 const int n=4;
 int data[n*n] = {
  5,4,3,2,
  11,-1,2,7,
  0,1,0,4,
  -13,79,1,2
 };
 int i,j,k=0;
 double **a = new double * [n];
 for (i=0; i<n; i++) a[i] = new double [n];
 for (i=0; i<n; i++)
 for (j=0; j<n; j++) a[i][j] = data[k++];

 printf ("\nDet(A)=%.1f",determinant(a,n));

 return 0;
}

14.06.2013, 19:40 [39039 просмотров]