|
Текст подпрограммы и версий aig3r_c.zip aig3d_c.zip aig3c_c.zip |
Тексты тестовых примеров taig3r_c.zip taig3d_c.zip taig3c_c.zip |
Обращение вещественной матрицы общего вида методом Гаусса с выбором ведущего элемента по столбцу и оценка числа обусловленности матрицы.
Для заданной вещественной квадратной матрицы А порядка N выполняется треугольная факторизация L- 1А = U, где L - нижняя треугольная матрица, U - верхняя треугольная матрица, вычисляется величина, обратная числу обусловленности матрицы А
rcond_c = 1/(|| A ||1*|| A-1 ||1) ,
где || A ||1 = maxj = 1, ..., N{ | a1 j | + | a2 j | + ...+ | aN j | }.
Далее вычисляется матрица U-1, и затем вычисляется А-1 = U-1L-1.
Дж. Форсайт, М. Малькольм, К. Моулер. Машинные методы математических вычислений. М.: Мир, 1980.
int aig3r_c (real *a, integer *m, integer *n, integer *nlead,
real *rcond, real *z__, integer *ierr)
Параметры
| a - | вещественный двумерный массив размера m на n, в котором задается исходная матрица; на выходе на месте исходной матрицы находится вычисленная обращенная матрица; |
| m - | первая размерность массива a в вызывающей программе (тип: целый); |
| n - | порядок обращаемой матрицы (тип: целый); |
| nlead - | целый вектор длины n, содержащий на выходе информацию о выполненных в ходе факторизации перестановках (см. замечания по использованию); |
| rcond - | вещественная переменная, содержащая на выходе вычисленное значение величины, обратной числу обусловленности матрицы; |
| z - | вещественный рабочий вектор длины n; |
| ierr - | целая переменная, содержащая на выходе информацию о прохождении счета; при этом: |
| ierr=65 - | если m ≤ 0 или n ≤ 0; |
| ierr=66 - | если в процессе счета возникло переполнение (это говорит о том, что либо ||A||1, либо некоторые элементы матрицы U, либо некоторые элементы обращенной матрицы превосходят по абсолютной величине максимальное представимое на данной машине число); |
| ierr=-k - | если в результате факторизации в к-й строке матрицы U диагональный элемент равен нулю (это свидетельствует о вырожденности матрицы A). Если таких строк у матрицы U несколько, то значение k полагается равным номеру последней из них. |
Версии
| aig3d_c - | обращение методом Гаусса с выбором ведущего элемента по столбцу и оценка числа обусловленности вещественной матрицы A, заданной с удвоенной точностью. |
| aig3c_c - | обращение методом Гаусса с выбором ведущего элемента по столбцу и оценка числа обусловленности комплексной матрицы A. |
Вызываемые подпрограммы
| afg4r_c - | подпрограмма треугольной факторизации и оценки числа обусловленности матрицы A; |
| utafsi_c - | подпрограмма выдачи диагностических сообщений. |
Замечания по использованию
| 1. |
В подпрограмме aig3d_c массивы a и z и переменная rcond имеют тип double, для треугольного разложения и оценки числа обусловленности матрицы A вызывается подпрограмма afg4d_c. | |
| 2. |
В подпрограмме aig3c_c массивы a и z имеют тип complex, для треугольного разложения матрицы A вызывается подпрограмма afg4c_c. | |
| 3. |
На выходе k - й элемент вектора nlead равен номеру строки, переставленной на k - м шаге факторизации с k - й строкой матрицы A. Так как факторизация Гаусса требует n - 1 шагов, то nlead (n) = n. | |
| 4. |
Если вырабатывается значение переменной ierr, отличное от нуля, то выдается соответствующее диагностическое сообщение и происходит выход из подпрограммы. |
int main(void)
{
/* Initialized data */
static float a[16] /* was [4][4] */ = { 7.9f,8.5f,4.3f,3.2f,5.6f,-4.8f,
4.2f,-1.4f,5.7f,.8f,-3.2f,-8.9f,-7.2f,3.5f,9.3f,3.3f };
/* System generated locals */
int i__1, i__2, i__3;
/* Local variables */
static int ierr;
extern int aig3r_c(float *, int *, int *, int *, float *, float *, int *);
static float b[16] /* was [4][4] */,
c__[16] /* was [4][4] */;
static int i__, j, k, m, n, nlead[4];
static float z__[4], rcond;
#define a_ref(a_1,a_2) a[(a_2)*4 + a_1 - 5]
#define b_ref(a_1,a_2) b[(a_2)*4 + a_1 - 5]
#define c___ref(a_1,a_2) c__[(a_2)*4 + a_1 - 5]
m = 4;
n = m;
i__1 = m;
for (i__ = 1; i__ <= i__1; ++i__) {
i__2 = n;
for (j = 1; j <= i__2; ++j) {
b_ref(i__, j) = a_ref(i__, j);
/* l1: */
}
/* l2: */
}
for (i__ = 1; i__ <= 4; ++i__) {
printf("\n %15.5e %15.5e %15.5e %15.5e \n",
a_ref(i__,1), a_ref(i__,2), a_ref(i__,3), a_ref(i__,4));
}
aig3r_c(a, &m, &n, nlead, &rcond, z__, &ierr);
printf("\n %5i %5i %5i %5i \n",
nlead[0], nlead[1], nlead[2], nlead[3]);
for (i__ = 1; i__ <= 4; ++i__) {
printf("\n %15.5e %15.5e %15.5e %15.5e \n",
a_ref(i__,1), a_ref(i__,2), a_ref(i__,3), a_ref(i__,4));
}
printf("\n %5i %20.12e \n", ierr, rcond);
i__1 = m;
for (i__ = 1; i__ <= i__1; ++i__) {
i__2 = n;
for (j = 1; j <= i__2; ++j) {
c___ref(i__, j) = 0.f;
i__3 = n;
for (k = 1; k <= i__3; ++k) {
c___ref(i__, j) = c___ref(i__, j) + b_ref(i__, k) * a_ref(k,
j);
/* l3: */
}
/* l4: */
}
/* l5: */
}
for (i__ = 1; i__ <= 4; ++i__) {
printf("\n %15.5e %15.5e %15.5e %15.5e \n",
c___ref(i__,1), c___ref(i__,2), c___ref(i__,3), c___ref(i__,4));
}
i__2 = m;
for (i__ = 1; i__ <= i__2; ++i__) {
i__1 = n;
for (j = 1; j <= i__1; ++j) {
c___ref(i__, j) = 0.f;
i__3 = n;
for (k = 1; k <= i__3; ++k) {
c___ref(i__, j) = c___ref(i__, j) + a_ref(i__, k) * b_ref(k,
j);
/* l6: */
}
/* l7: */
}
/* l8: */
}
for (i__ = 1; i__ <= 4; ++i__) {
printf("\n %15.5e %15.5e %15.5e %15.5e \n",
c___ref(i__,1), c___ref(i__,2), c___ref(i__,3), c___ref(i__,4));
}
return 0;
} /* main */
Результаты:
| 0.05056 0.05429 0.00629 0.03500 |
a = | 0.05189 -0.08460 0.07212 -0.00030 |
| -0.00841 0.04319 0.02021 -0.12113 |
| -0.04971 0.02797 0.07900 -0.05773 |
nlead = (2, 2, 4, 4)
rcond = 0.41764