Текст подпрограммы и версий
ast3r_c.zip , ast3d_c.zip , ast3c_c.zip 		Тексты тестовых примеров
tast3r_c.zip , tast3d_c.zip

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

Назначение

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

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

Заданная в компактной формe cимметричная положительно определенная трехдиагональная матрица А порядка N с помощью модифицированного алгоритма Гаусса без выбора ведущего элемента приводится к специальному виду L- 1А = С, где трехдиагональная матрица c имеет нулевые элементы в верхней половине нижней кодиагонали и в нижней половине верхней кодиагонали, а в строке с номером [N/2] + 1 отличен от нуля лишь диагональный элемент.

Затем решается система Сx = L- 1b.

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

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

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

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

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

  2. 

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

  3. 

Если вырабатывается значение ierr ≠ 0, то выдается соответствующее диагностическое сообщение и, если ierr > 0, то происходит выход из подпрограммы.

  4. 

Так как факторизация происходит без выбора ведущего элемента, то подпрограмму рекомендуется использовать для решения систем с положительно определенными матрицами, так как в этом случае процесс факторизации будет устойчивым относительно погрешностей данных и ошибок округления.

Для матриц, не являющихся положительно определенными, процесс будет, вообще говоря, неустойчивым; кроме того, в процессе факторизации может появиться нулевой ведущий элемент (в этом случае также происходит выход из подпрограммы).

Поэтому для решения систем с такими матрицами лучше обратиться к подпрограмме ast2r_c, задав матрицу A двумерным массивом размера m на 3.

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

int main(void)
{
    /* Initialized data */
    static float a[10] /* was [5][2] */ = { 0.f,-1.f,-1.f,-1.f,-1.f,1.f,2.f,
                                            2.f,2.f,2.f };
    static float b[5] = { 1.f,0.f,0.f,0.f,0.f };

    /* System generated locals */
    int i__1;

    /* Local variables */
    static int ierr;
    extern int ast3r_c(float *, int *, int *, float *, int *);
    static int i__, j, m, n;

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

    m = 5;
    n = 5;
    for (i__ = 1; i__ <= 5; ++i__) {
         printf("\n  %12.4e %12.4e \n", a_ref(i__, 1), a_ref(i__, 2));
    }
    printf("\n  %12.4e %12.4e %12.4e %12.4e %12.4e \n",
               b[0], b[1], b[2], b[3], b[4]);
    ast3r_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  =   ( 5.0,  4.0,  3.0,  2.0,  1.0 )