|
Текст подпрограммы и версий adg2r_c.zip , adg2d_c.zip , adg2c_c.zip |
Тексты тестовых примеров tadg2r_c.zip , tadg2d_c.zip , tadg2c_c.zip |
Вычисление определителя методом Гаусса с выбором ведущего элемента по столбцу и оценка числа обусловленности вещественной матрицы общего вида.
Для заданной вещественной квадратной матрицы А порядка N выполняется треугольная факторизация L - 1*А = U, где U - верхняя треугольная матрица, и вычисляется величина, обратная числу обусловленности матрицы А:
rcond = 1/( || A ||1 * || A-1 ||1) ,
где || A ||1 = maxj = 1, ..., N{ | a1 j | + | a2 j | + ...+ | aN j | }.
Определитель матрицы А вычисляется как произведение диагональных элементов матрицы U, умноженное на (-1)I, где I - число выполненных в процессе факторизации перестановок, и записывается в виде:
det A = D1*10D2 , где 1. ≤ D1 < 10.
Дж. Форсайт, М. Малькольм, К. Моулер. Машинные методы математических вычислений. Изд-во "Мир", М.: 1980.
int adg2r_c (real *a, integer *m, integer *n, integer *nlead,
real *det1, real *det2, real *rcond, real *z__, integer *ierr)
Параметры
| a - | вещественный двумерный массив размера m на n, в котором задается исходная матрица; на выходе на соответствующих местах находятся элементы матрицы u и поддиагональные элементы матриц исключения метода Гаусса li , i = 1, ..., N-1; |
| m - | первая размерность массива a в вызывающей программе (тип: целый); |
| n - | порядок матрицы (тип: целый); |
| nlead - | целый вектор длины n, содержащий на выходе информацию о выполненных в ходе факторизации перестановках (см. замечания по использованию); |
| det1 - | вещественная переменная, содержащая на выходе мантиссу определителя; |
| det2 - | вещественная переменная, содержащая на выходе десятичный порядок определителя; |
| rcond - | вещественная переменная, содержащая на выходе вычисленное значение величины, обратной числу обусловленности матрицы a; |
| z - | вещественный рабочий вектор длины n; |
| ierr - | целая переменная, содержащая на выходе информацию о прохождении счета, при этом: |
| ierr=65 - | если m ≤ 0 или n ≤ 0; |
| ierr=66 - | если в процессе счета произошло переполнение (это говорит о том, что либо ||a||1, либо некоторые элементы матрицы u превосходят по абсолютной величине максимальное представимое на данной машине число); |
| ierr=-k - | если в результате факторизации в К-ой строке матрицы u диагональный элемент равен нулю (это свидетельствует о вырожденности матрицы a). Если таких строк у матрицы u несколько, то значение K полагается равным номеру последней из них. |
Версии
| adg2d_c - | вычисление определителя методом Гаусса с выбором ведущего элемента по столбцу и оценка числа обусловленности вещественной матрицы, заданной с удвоенной точностью. |
| adg2c_c - | вычисление определителя методом Гаусса с выбором ведущего элемента по столбцу и оценка числа обусловленности комплексной матрицы. |
Вызываемые подпрограммы
| afg4r_c - | подпрограмма треугольной факторизации и оценка числа обусловленности матрицы a. |
| utafsi_c - | подпрограмма выдачи диагностических сообщений. |
Замечания по использованию
| 1. |
В подпрограмме adg2d_c массивы a и z и переменные rcond, det1 и det2 имеют тип double, для треугольного разложения и оценки числа обусловленности матрицы a вызывается подпрограмма afg4d_c. | |
| 2. |
В подпрограмме adg2c_c массивы a и z и переменная det1 имеют тип complex, для треугольного разложения и оценки числа обусловленности матрицы a вызывается подпрограмма afg4c_c. | |
| 3. |
На выходе k - й элемент вектора nlead равен номеру строки, переставленной на k - ом шаге факторизации с k - ой строкой матрицы a. Так как факторизация Гаусса требует n - 1 шагов, то nlead (n) = n. | |
| 4. | Если вырабатывается значение переменной ierr, отличное от нуля, то выдается соответствующее диагностическое сообщение, полагается rcond = 0, det1 = 0 , det2 = 0 и происходит выход из подпрограммы. |
int main(void)
{
/* Initialized data */
static float a[16] /* was [4][4] */ = { 1.f,.42f,.54f,.66f,.42f,1.f,.32f,
.44f,.54f,.32f,1.f,.22f,.66f,.44f,.22f,1.f };
/* Local variables */
static int ierr;
extern int adg2r_c(float *, int *, int *, int *, float *, float *,
float *, float *, int *);
static int j, m, n, nlead[4];
static float z__[4], rcond, det1, det2;
#define a_ref(a_1,a_2) a[(a_2)*4 + a_1 - 5]
m = 4;
n = m;
for (j = 1; j <= 4; ++j) {
printf("\n %16.7e %16.7e %16.7e %16.7e \n",
a_ref(j, 1), a_ref(j, 2), a_ref(j, 3), a_ref(j, 4));
}
adg2r_c(a, &m, &n, nlead, &det1, &det2, &rcond, z__, &ierr);
printf("\n %5i %5i %5i %5i \n",nlead[0], nlead[1], nlead[2], nlead[3]);
for (j = 1; j <= 4; ++j) {
printf("\n %16.7e %16.7e %16.7e %16.7e \n",
a_ref(j, 1), a_ref(j, 2), a_ref(j, 3), a_ref(j, 4));
}
printf("\n %5i %22.14e \n",ierr, rcond);
printf("\n %22.14e %22.14e \n",det1, det2);
return 0;
} /* main */
Результаты:
nlead = (1, 2, 3, 4)
| -1.0 0.42 0.54 0.66 |
a = | -0.42 0.83260 0.09320 0.16280 |
| -0.54 -0.11316 0.69785 -0.15482 |
| -0.66 -0.19767 0.22186 0.49787 |
rcond = 0.09880
det1 = 2.86152
det2 = -1.00000