|
Текст подпрограммы и версий ammhr_c.zip , ammhd_c.zip , ammhc_c.zip |
Тексты тестовых примеров tammhr_c.zip , tammhd_c.zip , tammhc_c.zip |
Умножение матрицы, обратной к заданной в компактной форме вещественной ленточной матрице A, на прямоугольную матрицу и оценка числа обусловленности матрицы A.
Для заданной в компактной форме ленточной вещественной
Для заданной в компактной форме ленточной вещественной матрицы А порядка N выполняется треугольная факторизация L- 1 А = U, где U - верхняя треугольная ленточная матрица, вычисляется величина, обратная числу обусловленности матрицы A:
RCOND = 1 / ( || A ||1 - || A-1 ||1 ) , где
N
|| A ||1 = max { ∑ | ai j | }
j=1,...N i=1
и затем для заданной вещественной прямоугольной матрицы B размера N на NN вычисляется прямоугольная матрица С = А- 1 В размера N на NN путем решения NN систем линейных алгебраических уравнений А*С (J) = В(J), где С(J) и В(J) суть J - е столбцы матриц C и B, J = 1,...,NN.
Дж.Форсайт, М.Малькольм, К.Моулер. Машинные методы математических вычислений. М.: Мир, 1980.
int ammhr_c (real *a, integer *ma, integer *n, integer *ml,
integer *mu, integer *nlead, real *b, integer *mb, integer *nn,
real *rcond, real *z, integer *ierr)
Параметры
| a - | вещественный двумерный массив размера ma*n, в первых ml+mu+1 столбцах которого задается в компактном виде ленточная матрица A порядка n; на выходе в первых ml столбцах массива находятся нижние кодиагонали ленточной матрицы L1 *...* Ln - 1 (Li, i = 1,...,n - 1, суть элементарные матрицы исключения метода Гаусса), в следующих ml+mu+1 столбцах содержится в компактном виде матрица U; |
| ma - | первая размерность массива a в вызывающей программе (тип: целый); |
| n - | порядок матрицы A и число строк матриц B и C (тип: целый); |
| ml - | число нижних кодиагоналей матрицы A (тип: целый); |
| mu - | число верхних кодиагоналей матрицы A (тип: целый); |
| nlead - | целый вектор длины n, содержащий на выходе информацию о выполненных в ходе факторизации перестановках (см. замечания по использованию); |
| b - | вещественный двумерный массив размера mb*nn, в котором задается матрица B; на выходе на соответствующих местах находятся элементы матрицы C = A- 1 B; |
| mb - | первая размерность массива b в вызывающей программе (тип: целый); |
| nn - | число столбцов матрицы B и C (тип: целый); |
| rcond - | вещественная переменная, содержащая на выходе вычисленное значение величины 1 / (|| A ||1 - || A- 1 ||1) (см. замечания по использованию); |
| z - | вещественный рабочий вектор длины n; |
| ierr - | целая переменная, содержащая на выходе информацию о прохождении счета, при этом: |
| ierr=65 - | если хотя бы одна из переменных ma, n, mb, nn имеет значение, меньшее единицы; |
| ierr=66 - | если в процессе счета произошло переполнение (это говорит о том, что либо || A ||1, либо некоторые элементы матрицы U или матрицы C превосходят по абсолютной величине максимальное представимое на данной машине число); |
| ierr= -k - | если в результате факторизации в k - ой строке матрицы U диагональный элемент равен нулю (это свидетельствует о вырождености матрицы A). Если таких строк у матрицы U несколько, то значение k полагается равным номеру последней из них; |
| ierr=67 - | если для некоторого j, 1 ≤ j ≤ nn, система A*C (j) = B(j) несовместна. |
Версии
| ammhd_c - | умножение матрицы, обратной к заданной с удвоенной точностью в компактной форме вещественной ленточной матрице A, на вещественную прямоугольную матрицу, заданную с удвоенной точностью, и оценка числа обусловленности матрицы A; |
| ammhc_c - | умножение матрицы, обратной к заданной в компактной форме комплексной ленточной матрице A, на комплексную прямоугольную матрицу и оценка числа обусловленности матрицы A. |
Вызываемые подпрограммы
| afb2r_c - | подпрограмма треугольной факторизации и оценки числа обусловленности ленточной матрицы A. |
| utafsi_c - | подпрограмма выдачи диагностических сообщений. |
Замечания по использованию
| 1. |
В подпрограмме ammhd_c массивы a, b, z и переменная rcond имеют тип double, для треугольного разложения и оценки числа обусловленности ленточной матрицы A вызывается подпрограмма afb2d_c. | |
| 2. |
В подпрограмме ammhc_c массивы a, b и z имеют тип complex, для треугольного разложения и оценки числа обусловленности ленточной матрицы A вызывается подпрограмма afb2c_c. | |
| 3. |
На выходе k - й элемент вектора nlead равен номеру строки, переставленной на k - м шаге факторизации с k - ой строкой матрицы A. Поскольку факторизация Гаусса требует n - 1 шагов, то nlead(n) = n. | |
| 4. |
При обращении к подпрограмме необходимо выполнение условия mb ≥ n. Кроме того, так как в результате выполненных в ходе факторизации перестановок число верхних кодиагоналей матрицы U равно mu+ml, а также в силу некоторых конструктивных особенностей подпрограммы, для правильной ее работы необходимо выполнение условия ma ≥ n > mu + 2 ml + 1. Если же mu + 2 ml + 1 ≥ n, то более целесообразно, задав матрицу A не в компактной, а в полной форме, обратиться к подпрограмме ammgr_c. | |
| 5. | Если вырабатывается значение переменной ierr, отличное от нуля, то выдается соответствующее диагностическое сообщение и если ierr > 0, то происходит выход из подпрограммы. |
int main(void)
{
/* System generated locals */
int i__1, i__2, i__3;
/* Local variables */
static int ierr;
static float a[25] /* was [5][5] */,
b[50] /* was [5][10] */;
static int i__, j, k, n, nlead[5];
static float z__[5];
extern int ammhr_c(float *, int *, int *, int *, int *, int *,
float *, int *, int *, float *, float *, int *);
static float rcond;
static int j0, j1, ma, mb, ml, nn, mu;
#define a_ref(a_1,a_2) a[(a_2)*5 + a_1 - 6]
#define b_ref(a_1,a_2) b[(a_2)*5 + a_1 - 6]
ma = 5;
n = 5;
ml = 1;
mu = 1;
mb = 5;
nn = 10;
i__1 = ma;
for (i__ = 1; i__ <= i__1; ++i__) {
i__2 = n;
for (j = 1; j <= i__2; ++j) {
a_ref(i__, j) = 0.f;
b_ref(i__, j) = 0.f;
b_ref(i__, j + 5) = 0.f;
/* L1: */
}
/* L2: */
}
i__1 = ma;
for (i__ = 1; i__ <= i__1; ++i__) {
/* Computing max */
i__2 = 1, i__3 = i__ - ml;
j0 = max(i__2,i__3);
/* Computing min */
i__2 = n, i__3 = i__ + mu;
j1 = min(i__2,i__3);
i__2 = j1;
for (j = j0; j <= i__2; ++j) {
k = j - i__ + ml + 1;
a_ref(i__, k) = (float) (i__ * 10 + j);
b_ref(i__, j) = a_ref(i__, k);
b_ref(i__, j + 5) = a_ref(i__, k);
/* L3: */
}
/* L4: */
}
for (i__ = 1; i__ <= 5; ++i__) {
printf("\n %12.4e %12.4e %12.4e %12.4e %12.4e \n",
a_ref(i__, 1), a_ref(i__, 2), a_ref(i__, 3),
a_ref(i__, 4), a_ref(i__, 5));
}
for (i__ = 1; i__ <= 5; ++i__) {
printf("\n %12.4e %12.4e %12.4e %12.4e %12.4e \n",
b_ref(i__, 1), b_ref(i__, 2), b_ref(i__, 3),
b_ref(i__, 4), b_ref(i__, 5));
printf("\n %12.4e %12.4e %12.4e %12.4e %12.4e \n",
b_ref(i__, 6), b_ref(i__, 7), b_ref(i__, 8),
b_ref(i__, 9), b_ref(i__, 10));
}
ammhr_c(a, &ma, &n, &ml, &mu, nlead, b, &mb, &nn, &rcond, z__, &ierr);
printf("\n %5i %5i %5i %5i %5i \n",
nlead[0], nlead[1], nlead[2], nlead[3], nlead[4]);
for (i__ = 1; i__ <= 5; ++i__) {
printf("\n %12.4e %12.4e %12.4e %12.4e %12.4e \n",
a_ref(i__, 1), a_ref(i__, 2), a_ref(i__, 3),
a_ref(i__, 4), a_ref(i__, 5));
}
printf("\n %5i \n", ierr);
printf("\n %16.7e \n", rcond);
for (i__ = 1; i__ <= 5; ++i__) {
printf("\n %12.4e %12.4e %12.4e %12.4e %12.4e \n",
b_ref(i__, 1), b_ref(i__, 2), b_ref(i__, 3),
b_ref(i__, 4), b_ref(i__, 5));
printf("\n %12.4e %12.4e %12.4e %12.4e %12.4e \n",
b_ref(i__, 6), b_ref(i__, 7), b_ref(i__, 8),
b_ref(i__, 9), b_ref(i__, 10));
}
return 0;
} /* main */
Результаты:
| 0 21.0 22.0 23.0 0 |
| -0.524 32.0 33.0 34.0 0 |
a_ref = | -0.015 43.0 44.0 45.0 0 |
| 0.292 54.0 55.0 0 0 |
| -0.228 0.569 0 0 0 |
nlead = (2, 3, 4, 5, 5) , rcond = 1.47362E-03 ,
| 1.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 |
| 0.0 1.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 |
b_ref = | 0.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 |
| 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0 0.0 |
| 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0 |