|
Текст подпрограммы и версий asb2r_c.zip , asb2d_c.zip , asb2c_c.zip |
Тексты тестовых примеров tasb2r_c.zip , tasb2d_c.zip , tasb2c_c.zip |
Решение вещественной системы линейных алгебраических уравнений A*X = b или AT*X = b методом Гаусса с выбором ведущего элемента по столбцу и оценка числа обусловленности ленточной матрицы A, заданной в компактной форме.
Для заданной в компактной форме ленточной вещественной матрицы А порядка N выполняется треугольная факторизация А = L*U, где U - верхняя треугольная ленточная матрица, вычисляется величина, обратная числу обусловленности матрицы А:
RCOND = 1 / ( || A ||1 * || A-1 ||1 ) ,
где
|| A ||1 = maxj = 1,...,N ( a1j + a2j + ... + aNj ) ,
и затем решается система А * Х = b или АT*Х = b.
Вычисления проводятся в два этапа:
1) решается система L * Y = b (для AT решается UT*Y=b) ;
2) решается система U * X = Y (для AT решается LT*X=Y) .
Дж.Форсайт, М.Малькольм, К.Моулер. Машинные методы математических вычислений. М.: Мир, 1980.
int asb2r_c (real *a, integer *ma, integer *n, integer *ml,
integer *mu, integer *nlead, real *b, integer *ltr, integer *l,
real *rcond, real *z, integer *ierr)
Параметры
| a - | вещественный двумерный массив размера ma на n, в первых ml + mu + 1 столбцах которого задается в компактном виде исходная ленточная матрица порядка n; на выходе в первых ml столбцах массива находятся нижние кодиагонали ленточной матрицы L1*...*Ln - 1 (Li, i = 1, ..., n - 1, суть элементарные матрицы исключения метода Гаусса), в следующих ml + mu + 1 столбцах содержится в компактном виде матрица U (см. замечания по использованию); |
| ma - | первая размерность массива в вызывающей программе (тип: целый); |
| n - | порядок матрицы A (тип: целый); |
| ml - | число нижних кодиагоналей матрицы A (тип: целый); |
| mu - | число верхних кодиагоналей матрицы A (тип: целый); |
| nlead - | целый вектор длины n, содержащий на выходе информацию о выполненных в процессе исключения перестановках (см. замечания по использованию); |
| b - | вещественный вектор длины n, в котором задается правая часть системы; на выходе содержит вычисленное решение системы (см. замечания по использованию); |
| ltr - | признак решаемой системы (тип: целый), причем |
| ltr = 0 - | если решается система A*X = b |
| ltr ≠ 0 - | если решается система AT*X = b; |
| l - | признак решаемой системы (тип: целый), причем |
| l = 0 - | если система с данной матрицей решается впервые, |
| l ≠ 0 - | если система с данной матрицей решается повторно (см. замечания по использованию); |
| rcond - | вещественная переменная, содержащая на выходе вычисленное значение величины 1/(|| A ||1*|| A- 1||1) (см. замечания по использованию); |
| z - | вещественный рабочий вектор длины n; |
| ierr - | целая переменная, содержащая на выходе информацию о прохождении счета, при этом |
| ierr=65 - | если ma ≤ 0 или n ≤ 0; |
| ierr=66 - | если в процессе работы произошло переполнение (это говорит о том, что либо ||A||1, либо некоторые элементы матрицы U, либо некоторые компоненты решения системы превосходят по абсолютной величине максимально представимое на данной машине число); |
| ierr=-k - | если в результате факторизации диагональный элемент в k - й строке матрицы U равен нулю (это свидетельствует о вырожденности матрицы A). Если таких строк у матрицы U несколько, то значение k полагается равным номеру последней из них; (см. замечания по использованию); |
| ierr=67 - | если система несовместна. |
Версии
| asb2d_c - | решение системы линейных алгебраических уравнений A*X = b или AT*X = b методом Гаусса с выбором ведущего элемента по столбцу и оценка числа обусловленности для вещественных A и b, заданных с удвоенной точностью, где ленточная матрица A задана в компактной форме. |
| asb2c_c - | решение системы линейных алгебраических уравнений A*X = b или AT*X = b методом Гаусса с выбором ведущего элемента по столбцу и оценка числа обусловленности для комплексных A и b, где ленточная матрица A задана в компактной форме. |
Вызываемые подпрограммы
| afb2r_c - | подпрограмма треугольной факторизации и оценки числа обусловленности ленточной матрицы A. |
| utafsi_c - | подпрограмма выдачи диагностических сообщений. |
Замечания по использованию
| 1. |
В подпрограмме asb2d_c массивы a, b, z и переменная rcond имеют тип double, для треугольного разложения и оценки числа обусловленности ленточной матрицы A вызывается подпрограмма afb2d_c. | |
| 2. |
В подпрограмме asb2c_c массивы a, b, z имеют тип complex, для треугольного разложения и оценки числа обусловленности ленточной матрицы A вызывается подпрограмма afb2c_c. | |
| 3. |
На выходе k - ый элемент вектора nlead равен номеру строки, переставленной на k - м шаге факторизации с k - ой строкой матрицы A. Поскольку факторизация Гаусса требует n - 1 шагов, то nlead (n) = n. | |
| 4. |
Так как в результате выполненных в ходе факторизации перестановок число верхних кодиагоналей матрицы U равно mu + ml, а также в силу некоторых конструктивных особенностей подпрограммы, для правильной ее работы необходимо выполнение условия ma ≥ n > mu + 2*ml + 1. Если же mu + 2*ml + 1 ≥ n, то более целесообразно, задав матрицу A не в компактной, а в полной форме, обратиться к подпрограмме asg8r_c. | |
| 5. |
Если задано l ≠ 0, то есть система решается повторно, то перед выполнением подпрограммы нужно запомнить, если требуется, вычисленную ранее оценку числа обусловленности, а в качестве матрицы A при обращении к подпрограмме нужно взять результат предыдущего обращения к asb2r_c либо к afb2r_c, т.е. матрица A должна быть уже факторизована, т.к. при l ≠ 0 обращения к afb2r_c не происходит и полагается rcond = 0.0. | |
| 6. | Если вырабатывается значение переменной ierr, отличное от нуля, то выдается соответствующее диагностическое сообщение, и, если ierr > 0, то происходит выход из подпрограммы. Если система совместна, то матрица A вырождена, то есть для некоторых номеров k U (k, k) = 0.0, то полагается b (k) = 1.0. |
int main(void)
{
/* Local variables */
static int ierr;
extern int asb2r_c(float *, int *, int *, int *, int *, int *,
float *, int *, int *, float *, float *, int *);
static float a[25] /* was [5][5] */, b[5];
static int i__, j, k, l, n, nlead[5];
static float z__[5], rcond;
static int j0, j1, ma, ml, mu, ltr;
int i__1, i__2, i__3;
#define a_ref(a_1,a_2) a[(a_2)*5 + a_1 - 6]
ma = 5;
n = 5;
ml = 1;
mu = 1;
l = 0;
ltr = 0;
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;
/* 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);
/* l3: */
}
b[i__ - 1] = 7.f;
/* 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));
}
asb2r_c(a, &ma, &n, &ml, &mu, nlead, b, <r, &l, &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 \n", b[i__-1]);
}
return 0;
} /* main */
Результаты:
| 0 21.0 22.0 23.0 0 |
| -0.524 32.0 33.0 34.0 0 |
a = | -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 ,
b = (-7.0625, 7.0598, 0.0024, -6.4409, 6.4511) .