The problem of calculating the inverse of the matrix type float

It is necessary to calculate the condition number of the matrix by definition and give a lower estimate for this number (that is, an experimental estimate).

The condition number by definition: condA = || A ^ (- 1) || * || A ||. The fact is that before this task, the program implemented the Gauss method for solving the SLAE and the qr-decomposition method. To reverse the matrix (and a 4x4 matrix is ​​given), then I simply solve four equations, where the right side is first (1,0,0,0), (0,1,0,0), (0,0,1,0) , (0,0,0,1). I get the inverse matrix, multiply by the original and get the unit. Everything would be fine, but we picked up a matrix in which it is with the type of "float" that the multiplication is not single. With the type of "double" everything is good, but it is with "float" - nonsense.

The matrix itself (also in the readable file):

0.2910 1.8100 9.3110 9.1100 4.2280 1.4500 8.5790 44.1950 42.9950 20.4290 -0.2900 -1.7980 -9.2500 -9.0500 -4.2080 0.0000 0.0820 0.4100 0.4500 0.1220 

By multiplying the original matrix (see above) and the inverse, it turns out that the matrix is ​​not exactly the same (in some places the numbers in the spirit of 0.06, etc.) emerge. Once again, this happens only with the type "float". For the code, I ask you not to scold too much, I know that it is heavily cluttered.

 #include <iostream> #include <iomanip> #include <cmath> #include <fstream> #include <float.h> using namespace std; typedef float tdouble; const tdouble eps = 1e-16; //при 1e-8 DATA9 не работает const unsigned n = 4; void product(tdouble**, tdouble*, tdouble*); void product(tdouble**, tdouble**, tdouble**); void output(tdouble**, tdouble*); void output(tdouble*); void output(tdouble **); void changing(tdouble**, tdouble*, int, int); bool gausss(tdouble**, tdouble*); void qr(tdouble**, tdouble*, tdouble**); void qr(tdouble**, tdouble*, tdouble*, tdouble*, tdouble*, tdouble**); void reverse(tdouble**, tdouble*, tdouble*); void spherical(tdouble*, tdouble*)ж tdouble octahedron(tdouble*); tdouble octahedron(tdouble*, tdouble*); tdouble cubic(tdouble*, tdouble*); tdouble cubic(tdouble *); tdouble cubic_matrix(tdouble **); tdouble octah_matrix(tdouble **); int main() { setlocale(LC_ALL, "Russian"); tdouble **mtx, **mtxun, *vect, *vectun; mtx = new tdouble *[n]; //массив указателей mtxun = new tdouble *[n]; //неизменяемая матрица vect = new tdouble[n]; vectun = new tdouble[n]; ifstream filemat; filemat.open("U:\\FS13\\FS2-x3\\Melikhov\\Metody_vych\\LLaba1 \\DATA9.txt"); for (int i = 0; i < n; i++) { mtxun[i] = new tdouble[n]; mtx[i] = new tdouble[n]; for (int j = 0; j < n + 1; j++) { filemat >> (j == 4 ? vect[i] : mtx[i][j]); //если j == 4 ? vectun[i] = vect[i] : mtxun[i][j] = mtx[i][j]; } } filemat.close(); cout << "\t\tВывод исходной матрицы и вектора правой части" << endl << endl; output(mtx, vect); //вывод на экран матрицы и правого вектор //Проверка на вырожденность bool h = gausss(mtx, vect); tdouble product_diagonal_elements = mtx[0][0]; //произведение диагональных элементов (для проверки вырожденности) for (int i = 1; i < n; i++) product_diagonal_elements *= mtx[i][i]; if (!h || abs(product_diagonal_elements) < eps) //если модуль произведения диагональных элементов меньше эпсилон, то вывод сообщения { cout << "Матрица вырожденная. До свидания!"; } else //иначе { tdouble *sol = new tdouble[n], *check = new tdouble[n]; tdouble *sin1 = new tdouble[n], *sin2 = new tdouble[n], *sin3 = new tdouble[n], *sin4 = new tdouble[n]; //единичная правая часть tdouble *solr1 = new tdouble[n], *solr2 = new tdouble[n], *solr3 = new tdouble[n], *solr4 = new tdouble[n]; //решение с единичной правой частью tdouble *solv1 = new tdouble[n], *solv2 = new tdouble[n], *solv3 = new tdouble[n], *solv4 = new tdouble[n]; //решение с единичной правой частью tdouble **mtxrev = new tdouble*[n]; //обратная матрица tdouble nun_cub, nrev_cub, condex_cub; //норма исходной, норма обратной, число обусловленности tdouble nun_oct, nrev_oct, condex_oct; for (int i = 0; i < n; i++) { (i == 0 ? sin1[i] = 1 : sin1[i] = 0); (i == 1 ? sin2[i] = 1 : sin2[i] = 0); (i == 2 ? sin3[i] = 1 : sin3[i] = 0); (i == 3 ? sin4[i] = 1 : sin4[i] = 0); } //приводим к первоначальному виду mtx и vect tdouble **T = new tdouble*[n], **mtxcheck = new tdouble*[n]; //создаем новый массив T - единичную матрицу. Понадобится для получения матрицы Q for (int i = 0; i < n; i++) { T[i] = new tdouble[n]; mtxrev[i] = new tdouble[n]; mtxcheck[i] = new tdouble[n]; for (int j = 0; j < n; j++) { T[i][j] = (i == j ? 1 : 0); //если i=j, тогда присваем 1, иначе - 0 mtx[i][j] = mtxun[i][j]; //в mtx записываем неизмененную матрицу mtxun для восстановления первой } vect[i] = vectun[i]; //аналогично с вектором } //QR-метод qr(mtx, vect, sin1, sin2, sin3, sin4, T); ////Транспонируем T в Q, создаем массив A tdouble **Q = new tdouble*[n], **A = new tdouble*[n]; for (int i = 0; i < n; i++) { Q[i] = new tdouble[n]; A[i] = new tdouble[n]; for (int j = 0; j < n; j++) Q[i][j] = T[j][i]; } //Находим обратную матрицу reverse(mtx, sin1, solr1); reverse(mtx, sin2, solr2); reverse(mtx, sin3, solr3); reverse(mtx, sin4, solr4); for (int j = 0; j < n; j++) { for (int i = 0; i < n; i++) { if (j == 0) mtxrev[i][j] = solr1[i]; else if (j == 1) mtxrev[i][j] = solr2[i]; else if (j == 2) mtxrev[i][j] = solr3[i]; else if (j == 3) mtxrev[i][j] = solr4[i]; } } cout << "Обратная матрица" << endl; output(mtxrev); product(mtxun, mtxrev, mtxcheck); cout << "Проверка" << endl; output(mtxcheck); cout << endl << endl << "Кубическая норма" << endl; nun_cub = cubic_matrix(mtxun); nrev_cub = cubic_matrix(mtxrev); condex_cub = nun_cub * nrev_cub; cout << endl << "Норма исходной матрицы = " << nun_cub << endl << "Норма обратной матрицы = " << nrev_cub << endl << "Произведение норм = " << condex_cub << endl; cout << endl << endl << "Октаэдрическая норма" << endl; nun_oct = octah_matrix(mtxun); nrev_oct = octah_matrix(mtxrev); condex_oct = nun_oct * nrev_oct; cout << endl << "Норма исходной матрицы = " << nun_oct << endl << "Норма обратной матрицы = " << nrev_oct << endl << "Произведение норм = " << condex_oct << endl; //Обратный ход reverse(mtx, vect, sol); tdouble n_dx, n_db, //нормы разностей dx, db nx_d, nb_d, //нормы деления х, b n_x, n_b; //нормы x и b tdouble max_div, norm_div; //максимум деления, норма деления tdouble o_dx, o_db, ox_d, ob_d, o_x, o_b, maxo_div, normo_div; //октаэдрическая n_x = cubic(sol); o_x = octahedron(sol); //cout << endl << "n_x " << n_x << endl; n_b = cubic(vectun); o_b = octahedron(vectun); //cout << endl << "n_b " << n_b << endl; max_div = 0; maxo_div = 0; //Приводим к исходному виду vect и mtx for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) mtx[i][j] = mtxun[i][j]; vect[i] = vectun[i]; } for (int i = 0; i < n; i++) { vect[i] += 0.01; //возмущаем i-ую компоненту вектора vect //cout << "VECT" << endl; //output(vect); n_db = cubic(vect, vectun); o_db = octahedron(vect, vectun); //cout << endl << "n_db " << n_db << endl; qr(mtx, vect, T); reverse(mtx, vect, solv1); //cout << endl; output(solv1); cout << endl; for (int k = 0; k < n; k++) { for (int j = 0; j < n; j++) mtx[k][j] = mtxun[k][j]; vect[k] = vectun[k]; } n_dx = cubic(solv1, sol); o_dx = octahedron(solv1, sol); //cout << endl << "n_dx " << n_dx << endl; nx_d = n_dx / n_x; ox_d = o_dx / o_x; //cout << endl << "nx_d " << nx_d << endl; nb_d = n_db / n_b; ob_d = o_db / o_b; //cout << endl << "nb_d " << nb_d << endl; norm_div = nx_d / nb_d; normo_div = ox_d / ob_d; //cout << endl << "norm_div " << norm_div << endl; if (max_div < norm_div) max_div = norm_div; if (maxo_div < norm_div) maxo_div = normo_div; } cout << endl << "Оценка снизу (кубическая норма): " << max_div << endl; cout << endl << "Оценка снизу (октаэдрическая норма): " << maxo_div << endl; delete[] sin1; delete[] sin2; delete[] sin3; delete[] sin4; delete[] check; delete[] sol; for (int i = 0; i < n; i++) { delete[] A[i]; delete[] Q[i]; delete[] T[i]; } delete[] A; delete[] Q; delete[] T; } cin.get(); delete[] vect; delete[] vectun; for (int i = 0; i < n; i++) { delete[] mtx[i]; delete[] mtxun[i]; } delete[] mtx; delete[] mtxun; } void reverse(tdouble **mtx, tdouble *vect, tdouble *sol) { //Обратный ход/ for (int i = n - 1; i >= 0; i--) { tdouble temp = 0; for (int j = i + 1; j < n; j++) //цикл с первого раза не пройдет, поэтому пустой массив sol[j] себя не обнаружит сначала, а обнаружит только после первой итерации внешнего цикла temp += mtx[i][j] * sol[j]; sol[i] = (vect[i] - temp) / mtx[i][i]; } } bool gausss(tdouble **mtx, tdouble *vect) { tdouble c; //коэффициент масштабирования for (int j = 0; j < n; j++) //двигаемся по столбцам { tdouble temp = abs(mtx[j][j]); //ищем максимум в столбце, запоминая j,j-й элемент (ищем всегда на диагонали и ниже) int mem = j; //переменная для запоминания номера строки, где ищется максимум for (int i = j + 1; i < n; i++) //от j+1, потому что под диагональю ищем максимальный элемент if (temp < abs(mtx[i][j])) { temp = abs(mtx[i][j]); //запоминаем новый максимум mem = i; //запоминаем соответствующую строку } if (temp < eps) return false; //если максимум сравним с нулем, то возвращаем false, то есть дальше работа программы бессмысленна changing(mtx, vect, mem, j); //меняем местами строки for (int i = j + 1; i < n; i++) //цикл для вычисления коэффициентов масштабирования и преобразования { c = mtx[i][j] / mtx[j][j]; for (int k = j; k < n; k++) mtx[i][k] -= c * mtx[j][k]; //пересчитываем i-ю строку в матрице (линейная комбинация)... vect[i] -= c * vect[j]; //... и в векторе } } return true; //возвращаем истину } void qr(tdouble **mtx, tdouble *vect, tdouble **T) //принцип: меняем местами строки так, чтобы не было деления на ноль. Затем вычисляем коэффициенты и заменяем строки на линейные комбинации { //QR-метод for (int k = 0; k < n; k++) //идем по строке { //int z = k; //запоминаем номер строки //while (abs(mtx[z][k]) < eps && z < n) //пока нулевой элемент, увеличиваем z, запоминающую строку. Также проверяется выход за пределы массива. // z++; //if (z != k) // changing(mtx, vect, z, k); tdouble c, s, temp_sqr; //c,s - коэффициенты, temp_sqr - временный корень for (int i = k + 1; i < n; i++) //сдвигаемся на строчку вниз и считаем: { //подсчет коэффициентов if (abs(mtx[i][k]) < eps) continue; temp_sqr = sqrt(mtx[k][k] * mtx[k][k] + mtx[i][k] * mtx[i][k]); c = mtx[k][k] / temp_sqr; s = mtx[i][k] / temp_sqr; tdouble temp, tempT; for (int j = 0; j < n; j++)//k { //замена на линейные комбинации tempT = T[k][j]; T[k][j] = c * T[k][j] + s * T[i][j]; T[i][j] = -s * tempT + c * T[i][j]; if (j >= k) // чтобы лишний раз не пересчитывать нули { temp = mtx[k][j]; //сохраняем первую строку во временную переменную mtx[k][j] = c * mtx[k][j] + s * mtx[i][j]; //заменяем k-ю строку на лк с коэффициентами c и s mtx[i][j] = -s * temp + c * mtx[i][j]; //заменяем j-ю строку на лк с использованием первоначальной k-й строки } } //аналогично с вектором правой части temp = vect[k]; vect[k] = c * vect[k] + s * vect[i]; vect[i] = -s * temp + c * vect[i]; } } } void qr(tdouble **mtx, tdouble *vect, tdouble *sin1, tdouble *sin2, tdouble *sin3, tdouble *sin4, tdouble **T) //принцип: меняем местами строки так, чтобы не было деления на ноль. Затем вычисляем коэффициенты и заменяем строки на линейные комбинации { //QR-метод for (int k = 0; k < n; k++) //идем по строке { //int z = k; //запоминаем номер строки //while (abs(mtx[z][k]) < eps && z < n) //пока нулевой элемент, увеличиваем z, запоминающую строку. Также проверяется выход за пределы массива. // z++; //if (z != k) // changing(mtx, vect, z, k); tdouble c, s, temp_sqr; //c,s - коэффициенты, temp_sqr - временный корень for (int i = k + 1; i < n; i++) //сдвигаемся на строчку вниз и считаем: { //подсчет коэффициентов if (abs(mtx[i][k]) < eps) continue; temp_sqr = sqrt(mtx[k][k] * mtx[k][k] + mtx[i][k] * mtx[i][k]); c = mtx[k][k] / temp_sqr; s = mtx[i][k] / temp_sqr; tdouble temp, tempT; for (int j = 0; j < n; j++)//k { //замена на линейные комбинации tempT = T[k][j]; T[k][j] = c * T[k][j] + s * T[i][j]; T[i][j] = -s * tempT + c * T[i][j]; if (j >= k) // чтобы лишний раз не пересчитывать нули { temp = mtx[k][j]; //сохраняем первую строку во временную переменную mtx[k][j] = c * mtx[k][j] + s * mtx[i][j]; //заменяем k-ю строку на лк с коэффициентами c и s mtx[i][j] = -s * temp + c * mtx[i][j]; //заменяем j-ю строку на лк с использованием первоначальной k-й строки } } //аналогично с вектором правой части temp = vect[k]; vect[k] = c * vect[k] + s * vect[i]; vect[i] = -s * temp + c * vect[i]; temp = sin1[k]; sin1[k] = c * sin1[k] + s * sin1[i]; sin1[i] = -s * temp + c * sin1[i]; temp = sin2[k]; sin2[k] = c * sin2[k] + s * sin2[i]; sin2[i] = -s * temp + c * sin2[i]; temp = sin3[k]; sin3[k] = c * sin3[k] + s * sin3[i]; sin3[i] = -s * temp + c * sin3[i]; temp = sin4[k]; sin4[k] = c * sin4[k] + s * sin4[i]; sin4[i] = -s * temp + c * sin4[i]; } } } void product(tdouble **mtx, tdouble *sol, tdouble *check) { for (int i = 0; i < n; i++) //фиксировали строку { check[i] = 0; for (int j = 0; j < n; j++) //двигаемся по столбцам и.... check[i] += mtx[i][j] * sol[j]; //перемножаем } } void product(tdouble **mtxun, tdouble **mtxrev, tdouble **mtxcheck) { for (int i = 0; i < n; i++) for (int j = 0; j < n; j++) mtxcheck[i][j] = 0; for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { for (int m = 0; m < n; m++) { mtxcheck[i][j] += mtxun[i][m] * mtxrev[m][j]; } } } } void output(tdouble *sol) { for (int i = 0; i < n; i++) cout << fixed << setprecision(4) << setw(42) << sol[i] << endl; } void output(tdouble **mtx, tdouble *vect) { for (int i = 0; i < n; i++) //фиксировали строку { for (int j = 0; j < n + 1; j++) //поехали по столбцам cout << fixed << setprecision(4) << setw(14) << (j == n ? vect[i] : mtx[i][j]); //если j==n, то вывести вектор (последний столбец), иначе - матрицу cout << endl; } } void output(tdouble **mtx) { for (int i = 0; i < n; i++) //фиксировали строку { for (int j = 0; j < n; j++) //поехали по столбцам cout << fixed << setprecision(4) << setw(14) << mtx[i][j]; //если j==n, то вывести вектор (последний столбец), иначе - матрицу cout << endl; } } void changing(tdouble **mtx, tdouble *vect, int line1, int line2) { tdouble *t = mtx[line2]; //запоминаем адрес mtx[line2] = mtx[line1]; mtx[line1] = t; tdouble r = vect[line2]; vect[line2] = vect[line1]; vect[line1] = r; } void spherical(tdouble *vectun, tdouble *check) { tdouble diff, diff2 = diff = 0.0; for (int i = 0; i < n; i++) { diff += (vectun[i] - check[i]) * (vectun[i] - check[i]); diff2 += vectun[i] * vectun[i] - 2 * vectun[i] * check[i] + check[i] * check[i]; } cout << scientific << endl << "\t\t\t\tНевязка сферическая" << endl << setw(14) << diff << ' ' << setw(14) << diff2 << endl; //os << scientific << endl << "\t\t\t\tНевязка сферическая" << endl << setw(14) << diff << ' ' << setw(14) << diff2 << endl; } tdouble octahedron(tdouble *vectun) { tdouble diffo_2 = 0.0; for (int i = 0; i < n; i++) diffo_2 += abs(vectun[i]); //cout << scientific << endl << "\t\t\t\tНевязка октаэдрическая" << endl << setw(14) << diff << endl; return diffo_2; // os << scientific << endl << "\t\t\t\tНевязка октаэдрическая" << endl << setw(14) << diff << endl; } tdouble octahedron(tdouble *vectun, tdouble *check) { tdouble diffo_1 = 0.0; for (int i = 0; i < n; i++) diffo_1 += abs(vectun[i] - check[i]); //cout << scientific << endl << "\t\t\t\tНевязка октаэдрическая" << endl << setw(14) << diff << endl; return diffo_1; // os << scientific << endl << "\t\t\t\tНевязка октаэдрическая" << endl << setw(14) << diff << endl; } tdouble cubic(tdouble *vectun, tdouble *check) { tdouble diffc_1 = abs(vectun[0] - check[0]); //максимум модуля ищется => запоминаем самую первую разность, а затем с ней сравниваем for (int i = 1; i < n; i++) { if (diffc_1 < abs(vectun[i] - check[i])) diffc_1 = abs(vectun[i] - check[i]); } // cout << scientific << endl << "\t\t\t\tНевязка кубическая" << endl << setw(14) << diff << endl; return diffc_1; //os << scientific << endl << "\t\t\t\tНевязка кубическая" << endl << setw(14) << diff << endl; } tdouble cubic(tdouble *vectun) { tdouble diffc_2 = abs(vectun[0]); //максимум модуля ищется => запоминаем самую первую разность, а затем с ней сравниваем for (int i = 1; i < n; i++) { if (diffc_2 < abs(vectun[i])) diffc_2 = abs(vectun[i]); } // cout << scientific << endl << "\t\t\t\tНевязка кубическая" << endl << setw(14) << diff << endl; return diffc_2; //os << scientific << endl << "\t\t\t\tНевязка кубическая" << endl << setw(14) << diff << endl; } tdouble cubic_matrix(tdouble **mtxun) { tdouble norm = 0; tdouble maxn = norm; for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) norm += abs(mtxun[i][j]); if (norm > maxn) maxn = norm; norm = 0; } norm= maxn; return norm; } tdouble octah_matrix(tdouble **mtxun) { tdouble norm = 0; tdouble maxn = norm; for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) norm += abs(mtxun[j][i]); if (norm > maxn) maxn = norm; norm = 0; } norm = maxn; return norm; }