Подпрограмма: ast1r_c
Назначение
Решение вещественной системы линейных алгебраических
уравнений Ax = b или ATx = b
с треугольной матрицей A.
Математическое описание
Для заданной вещественной тpeyгoльной матрицы А порядка N
решается система Аx = b (АTx = b),
причем для нахождения компонент решения xi,
i = 1, ..., N система
рассматриваeтся как одно векторнoe уравнeние
x1a1 + ... + xNaN = b ,
где векторы ai, i = 1, ..., N
суть столбцы (строки) матрицы А.
Дж.Форсайт, М.Малькольм, К.Моулер. Машинные методы
математических вычислений. М., Мир, 1980.
Использование
int ast1r_c (real *a, integer *m, integer *n, real *b,
integer *ltr, integer *low, integer *ierr)
Параметры
|
a -
|
двумерный массив размера m на n, в котором
задается матрица системы (тип: вещественный);
|
|
m -
|
первая размерность массива a в вызывающей
подпрограмме (тип: целый);
|
|
n -
|
порядок матрицы системы (тип: целый);
|
|
b -
|
вещественный вектор длины n в котором задается
правая часть системы; на выходе содержит
вычисленное решение системы (см. замечания по
использованию);
|
|
ltr -
|
признак решаемой системы (тип: целый), причем
|
|
ltr = 0 -
|
если решается система Ax = b,
|
|
ltr ≤ 0 -
|
если решается система ATx = b;
|
|
low -
|
признак решаемой системы (тип: целый), причем
|
|
low = 0 -
|
если матрица A верхняя треугольная,
|
|
low ≠ 0 -
|
если матрица A нижняя треугольная;
|
|
ierr -
|
целая переменная, содержащая на выходе
информацию о прохождении счета; при этом
|
|
ierr=65 -
|
если m ≤ 0 или
n ≤ 0;
|
|
ierr=66 -
|
если в процессе работы произошло
переполнение (это говорит о том, что
некоторые компоненты решения системы
превосходят по абсолютной величине
максимальное представимое на данной машине число);
|
|
ierr=-k -
|
если в k - й строке матрицы A
диагональный элемент равен нулю (это
свидетельствует о вырожденности матрицы). Если
таких строк несколько, то значение k
полагается равным номеру последней из
них (см. замечания по использованию);
|
|
ierr=67 -
|
если система несовместна.
|
Версии
|
ast1d_c -
|
решение системы линейных алгебраических
уравнений Ax = b или ATx = b с
треугольной матрицей A для вещественных a и b, заданных
с удвоенной точностью.
|
|
ast1c_c -
|
решение системы линейных алгебраических
уравнений Ax =b или ATx = b с
треугольной матрицей A для комплексных a и b.
|
Вызываемые подпрограммы
|
utafsi_c -
|
подпрограмма выдачи диагностических сообщений.
|
Замечания по использованию
|
| 1.
|
В подпрограмме ast1d_c массивы a и b имеют тип
double.
|
|
| 2.
|
В подпрограмме ast1c_c массивы a и b имеют тип
complex.
|
|
| 3.
|
Если вырабатывается значение переменной ierr,
отличное от нуля, то выдается соответствующее
диагностическое сообщение, и если ierr > 0, то происходит
выход из подпрограммы. Если система совместна, но
матрица A вырождена, т.е. для некоторых номеров k
a (k, k) = 0., то полагается
x (k) = 1.
|
Пример использования
int main(void)
{
/* System generated locals */
int i__1, i__2;
float r__1;
/* Local variables */
static int ierr;
extern int ast1r_c(float *, int *, int *, float *, int *, int *, int *);
static float a[25] /* was [5][5] */, b[5];
static int i__, j, m, n;
static float s, w, z__[5];
static int ltr, low;
#define a_ref(a_1,a_2) a[(a_2)*5 + a_1 - 6]
m = 5;
n = 5;
ltr = 0;
low = 0;
i__1 = n;
for (j = 1; j <= i__1; ++j) {
i__2 = n;
for (i__ = j; i__ <= i__2; ++i__) {
a_ref(i__, j) = 0.f;
/* l2: */
}
i__2 = j;
for (i__ = 1; i__ <= i__2; ++i__) {
a_ref(i__, j) = (float) (i__ * 10 + j);
/* l3: */
}
z__[j - 1] = (float) j;
/* l4: */
}
i__1 = m;
for (i__ = 1; i__ <= i__1; ++i__) {
b[i__ - 1] = 0.f;
i__2 = n;
for (j = i__; j <= i__2; ++j) {
b[i__ - 1] += a_ref(i__, j) * z__[j - 1];
/* l5: */
}
/* l6: */
}
for (i__ = 1; i__ <= 5; ++i__) {
printf("\n %10.4e %10.4e %10.4e %10.4e %10.4e \n",
a_ref(i__, 1), a_ref(i__, 2), a_ref(i__, 3),
a_ref(i__, 4), a_ref(i__, 5));
}
printf("\n %10.4e %10.4e %10.4e %10.4e %10.4e \n",
z__[0], z__[1], z__[2], z__[3], z__[4]);
ast1r_c(a, &m, &n, b, <r, &low, &ierr);
w = 0.f;
i__2 = n;
for (j = 1; j <= i__2; ++j) {
s = (r__1 = b[j - 1] - z__[j - 1], abs(r__1));
if (s > w) {
w = s;
}
/* l7: */
}
printf("\n %10.4e %10.4e %10.4e %10.4e %10.4e \n",
b[0], b[1], b[2], b[3], b[4]);
printf("\n %16.7e \n", w);
printf("\n %5i \n", ierr);
return 0;
} /* main */
Результат:
b = ( 1.000, 2.000, 3.000, 4.000, 5.000 ) .