Текст подпрограммы и версий
asg1r_c.zip , asg1d_c.zip , asg1c_c.zip 		Тексты тестовых примеров
tasg1r_c.zip , tasg1d_c.zip , tasg1c_c.zip

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

Назначение

Решение вещественной системы линейных алгебраических уравнений A*x = b с матрицей A и вектором правой части b методом Гаусса с выбором ведущего элемента по столбцу. При повторном решении системы с той же матрицей и другой правой частью эффективно используется выполненная ранее факторизация матрицы.

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

Матрица системы А N - го порядка приводится к верхнетреугольному виду U последовательностью элементарных преобразований Гаусса L1, L2, ..., LN - 1 и матрицей перестановок Q так, что 

     LN-1*LN-2...L1*Q*A = U ,

и затем решается система с треугольной матрицей U и правой частью
LN - 1*LN - 2...L1*Q*b . Матрицы Li ,  i = 1, 2, ..., N - 1 имеют единичные диагональные элементы, матрица Q осуществляет перестановку строк матрицы А и обеспечивает стратегию выбора ведущего элемента по столбцам.

В.В.Воеводин, Р.В.Петрина, Комплекс алгоритмов, основанных на преобразованиях типа Гаусса, в пакете линейной алгебры, Сб. "Численный анализ на ФОРТРАНе", вып,3, Изд-во МГУ, 1973.

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

    int asg1r_c (real *a, real *b, real *x, integer *s, integer *n,
             integer *l)

   Параметры
a - двумерный n на n массив, в котором задается исходная матрица (тип: real); по окончании работы подпрограммы в массиве a на сответствующих местах запоминаются элементы матрицы U и векторы-столбцы, порождающие матрицы Li ;
b - вектор длины n, в котором задается вектор правой части системы (тип: вещественный);
x - вектор длины n, в котором запоминается вектор вычисленного решения (тип: вещественный);
s - вектор длины n, в котором запоминается вектор, порождающий матрицу перестановок; при этом в s (k) запоминается номер строки переставленной на k - ом шаге с k - ой строкой (тип: целый);
n - заданный порядок матрицы системы (тип: целый);
l - признак решаемой задачи (тип: целый), причем
l = 1 - если система с данной матрицей решается впервые;
l ≠ 1 - если система с данной матрицей решается повторно. 
   Версии
asg1d_c - решение системы линейных алгебраических уравнений A*x = b методом Гаусса с выбором ведущего элемента по столбцу для вещественных A и b, заданных с удвоенной точностью.
asg1c_c - решение системы линейных алгебраических уравнений A*x = b методом Гаусса с выбором ведущего элемента по столбцу для комплексных A и b. 
   Вызываемые подпрограммы:  нет 
   Замечания по использованию
  1. 

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

  2. 

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

  3. 

Поскольку факторизация Гаусса требует n - 1 шагов, то s (n) = 0.

  4.  Подпрограммы asg1r_c, asg1d_c, asg1c_c позволяют помещать вычисленное решение x на место правой части b.

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

int main(void)
{
    /* Initialized data */
    static float pa[16] /* was [4][4] */ = { 7.9f,8.5f,4.3f,3.2f,5.6f,-4.8f,
                  4.2f,-1.4f,5.7f,.8f,-3.2f,-8.9f,-7.2f,3.5f,9.3f,3.3f };

    /* Local variables */
    extern int asg0r_c(float *, float *, float *, int *, int *),
               asg1r_c(float *, float *, float *, int *, int *, int *);
    static float b[4];
    static int i__, j;
    static float b1[4], b2[4];
    static int s2[4];
    static float x1[4], x2[4], pb[16] /* was [4][4] */;

#define pa_ref(a_1,a_2) pa[(a_2)*4 + a_1 - 5]
#define pb_ref(a_1,a_2) pb[(a_2)*4 + a_1 - 5]

    for (i__ = 1; i__ <= 4; ++i__) {
        for (j = 1; j <= 4; ++j) {
            pb_ref(i__, j) = pa_ref(i__, j);
/* l14: */
        }
    }
    for (i__ = 1; i__ <= 4; ++i__) {
/* l3: */
        b[i__ - 1] = 0.f;
    }
    b[0] = 1.f;
    for (i__ = 1; i__ <= 4; ++i__) {
        b1[i__ - 1] = b[i__ - 1];
        b2[i__ - 1] = b[i__ - 1];
/* l11: */
    }
    asg0r_c(pa, b1, x1, &c__4, &c__1);
    asg1r_c(pb, b2, x2, s2, &c__4, &c__1);

    for (i__ = 0; i__ <= 14; i__+= 2) {
         printf("\n  %16.7e %16.7e \n", pa[i__], pa[i__+1]);
    }
    for (i__ = 0; i__ <= 2; i__+= 2) {
         printf("\n  %16.7e %16.7e \n", b1[i__], b1[i__+1]);
    }
    for (i__ = 0; i__ <= 2; i__+= 2) {
         printf("\n  %16.7e %16.7e \n", x1[i__], x1[i__+1]);
    }
    for (i__ = 0; i__ <= 14; i__+= 2) {
         printf("\n  %16.7e %16.7e \n", pb[i__], pb[i__+1]);
    }
    for (i__ = 0; i__ <= 2; i__+= 2) {
         printf("\n  %16.7e %16.7e \n", b2[i__], b2[i__+1]);
    }
    for (i__ = 0; i__ <= 2; i__+= 2) {
         printf("\n  %16.7e %16.7e \n", x2[i__], x2[i__+1]);
    }
    return 0;
} /* main */


Результат:

      x  =  (0.051,  0.052,  -0.008,  0.050)