Текст подпрограммы и версий
ast2r_c.zip , ast2d_c.zip , ast2c_c.zip 		Тексты тестовых примеров
tast2r_c.zip , tast2d_c.zip , tast2c_c.zip

Подпрограмма:  ast2r_c

Назначение

Решение вещественной системы линейных алгебраических уравнений Ax = b с трехдиагональной матрицей A.

Математическое описание

Для заданной в компактной фopмe вещественной трехдиагональной матрицы А порядка N выполняется треугольная факторизация Гаусса L- 1А = U, где U - верхняя треугольная ленточная матрица, и затем решается система Ux = L- 1b.

Дж.Форсайт, М.Малькольм, К.Моулер. Машинные методы математических вычислений. М.: Мир, 1980.

Использование 

    int ast2r_c (real *a, integer *m, integer *n, real *b,
            integer *ierr)

   Параметры
a - вещественный двумерный массив размера m на 3, в котором задается в компактном виде исходная трехдиагональная матрица порядка n; на выходе содержит в компактном виде верхнюю треугольную ленточную матрицу U;
m - первая размерность массива a в вызывающей программе (тип: целый);
n - порядок матрицы A (тип: целый);
b - вещественный вектор длины n, в котором задается правая часть системы; на выходе содержит вычисленное решение системы (см. замечания по использованию);
ierr - целая переменная, содержащая на выходе информацию о прохождении счета, при этом:
ierr=65 - если m ≤ 0 или n ≤ 0;
ierr=66 - если в процессе работы произошло переполнение (это говорит о том, что некоторые элементы матрицы U, либо некоторые компоненты решения системы превосходят по абсолютной величине максимальное представимое на данной машине число);
ierr=-k - если в результате факторизации в k - ой строке матрицы U диагональный элемент равен нулю (это свидетельствует о вырожденности матрицы A); если таких строк у матрицы U несколько, то значение k полагается равным номеру последней из них (см. замечания по использованию);
ierr=67 - если система несовместна. 
   Версии
ast2d_c - решение системы линейных алгебраических уравнений Ax = b с трехдиагональной матрицей A, заданной в компактной форме, для вещественных a и b, заданных с удвоенной точностью.
ast2c_c - решение системы линейных алгебраических уравнений Ax = b с трехдиагональной матрицей A, заданной в компактной форме, для комплексных a и b. 
   Вызываемые подпрограммы
utafsi_c - подпрограмма выдачи диагностических сообщений. 
   Замечания по использованию
  1. 

В подпрограмме ast2d_c массивы a и b имеют тип double.

  2. 

В подпрограмме ast2c_c массивы a и b имеют тип complex.

  3.  Если вырабатывается значение переменной ierr, отличное от нуля, то выдается соответствующее диагностическое сообщение, и если ierr > 0, то происходит выход из подпрограммы. Если система совместна, но матрица A вырождена, т.е. для некоторых номеров k существуют U (k, k) = 0., то полагается b (k) = 1.

Пример использования 

int main(void)
{

    /* System generated locals */
    int i__1, i__2, i__3;

    /* Local variables */
    static int ierr;
    extern int ast2r_c(float *, int *, int *, float *, int *);
    static float a[15] /* was [5][3] */, b[5];
    static int i__, j, k, m, n;
    static float z__[5];
    static int j0, j1;

#define a_ref(a_1,a_2) a[(a_2)*5 + a_1 - 6]

    m = 5;
    n = 5;
    i__1 = m;
    for (i__ = 1; i__ <= i__1; ++i__) {
        z__[i__ - 1] = (float) i__;
/* Computing max */
        i__2 = 1, i__3 = i__ - 1;
        j0 = max(i__2,i__3);
/* Computing min */
        i__2 = n, i__3 = i__ + 1;
        j1 = min(i__2,i__3);
        i__2 = j1;
        for (j = j0; j <= i__2; ++j) {
            k = j - i__ + 2;
            a_ref(i__, k) = (float) (i__ * 10 + j);
/* l3: */
        }
/* l4: */
    }
    a_ref(1, 1) = 0.f;
    a_ref(5, 3) = 0.f;
    i__1 = m;
    for (i__ = 1; i__ <= i__1; ++i__) {
        b[i__ - 1] = 0.f;
        for (k = 1; k <= 3; ++k) {
            j = k + i__ - 2;
            b[i__ - 1] += a_ref(i__, k) * z__[j - 1];
/* l5: */
        }
/* l6: */
    }
    for (i__ = 1; i__ <= 5; ++i__) {
         printf("\n  %12.4e %12.4e %12.4e \n",
                a_ref(i__, 1), a_ref(i__, 2), a_ref(i__, 3));
    }
    printf("\n  %12.4e %12.4e %12.4e %12.4e %12.4e \n",
               z__[0], z__[1], z__[2], z__[3], z__[4]);
    printf("\n  %12.4e %12.4e %12.4e %12.4e %12.4e \n",
               b[0], b[1], b[2], b[3], b[4]);
    ast2r_c(a, &m, &n, b, &ierr);

    printf("\n  %12.4e %12.4e %12.4e %12.4e %12.4e \n",
               b[0], b[1], b[2], b[3], b[4]);
    printf("\n  %5i \n", ierr);
    return 0;
} /* main */


Результат:

      b  =   ( 1., 2., 3., 4., 5. )