Текст подпрограммы и версий
am02d_c.zip , am02c_c.zip
Тексты тестовых примеров
tam02d_c.zip , tam02c_c.zip

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

Назначение

Вычисление нормированного вектора невязки системы линейных алгебраических уравнений.

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

Вычисляется вектор R = r / S, где  r = b - Ax,  b - заданный вектор длины  N,  x - заданный вектор длины  М,  A - заданная матрица размера N * М,  S = || r || ∞.

Нормировка вектора невязки  r  необходима при решении системы линейных алгебраических уравнений Аx = b  с уточнением.

В.В.Воеводин, Вычислительные основы линейной алгебры, "Наука", М., 1977.

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

    int am02d_c (doublereal *a, doublereal *b, doublereal *x,
        doublereal *r, doublereal *s, integer *n, integer *m)

Параметры

a - двумерный массив размера n * m, в котором задается матрица системы (тип: вещественный двойной точности);
b - вектор длины  n, в котором запоминается правая часть системы (тип: вещественный двойной точности);
x - вектор длины  m, в котором задается решение системы (тип: вещественный двойной точности);
r - одномерный массив длины  n, в котором запоминается вычисленный нормированный вектор невязки (тип: вещественный двойной точности);
s - переменная, в которой запоминается бесконечная норма вычисленного вектора невязки (тип: вещественный двойной точности);
n, m - число строк и столбцов матрицы системы (тип: целый).

Версии

am02c_c - вычисление нормированного вектора невязки комплексной системы линейных алгебраических уравнений.

Вызываемые подпрограммы: нет.

Замечания по использованию

  1. 

В подпрограмме am02c_c параметры a, b, x, r имеют тип complex.

  2.  Подпрограммы am02d_c, am02c_c целесообразно использовать при решении системы линейных алгебраических уравнений с уточнением (см. пример использования).

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

int main(void)
{
    /* Initialized data */
    static double b[7] = { 0.,0.,0.,0.,0.,0.,0. };

    /* System generated locals */
    int i__1, i__2;
    double d__1;

    /* Local variables */
    extern int am02d_c(double *, double *, double *, double *, double *,
                       int *, int *),
               ash0d_c(double *, double *, double *, double *, int *, int *);
    static double a[49] /* was [7][7] */, d__[7], f;
    static int i__, j, k, n;
    static double r__[7], s, x[7], z__, a1[49] /* was [7][7] */, s0,
                  s1[7], x1[7];

#define a_ref(a_1,a_2) a[(a_2)*7 + a_1 - 8]
#define a1_ref(a_1,a_2) a1[(a_2)*7 + a_1 - 8]

    n = 7;
    i__1 = n;
    for (i__ = 1; i__ <= i__1; ++i__) {
        i__2 = n;
        for (j = 1; j <= i__2; ++j) {
            a_ref(i__, j) = 360360. / (double) (i__ + j - 1);
/* l4: */
            a1_ref(i__, j) = a_ref(i__, j);
        }
    }
    b[4] = 360360.;
    ash0d_c(a, b, x, s1, &n, &c__1);
    s0 = 0.;
    i__2 = n;
    for (i__ = 1; i__ <= i__2; ++i__) {
        x1[i__ - 1] = x[i__ - 1];
        if ((d__1 = x[i__ - 1], abs(d__1)) <= s0) {
            goto l3;
        }
        s0 = (d__1 = x[i__ - 1], abs(d__1));
l3:
        ;
    }
    for (k = 1; k <= 6; ++k) {
         for (i__ = 1; i__ <= 7; ++i__) {
              printf("\n  %20.12e \n", x[i__-1]);
         }
    }
    am02d_c(a1, b, x, r__, &s, &n, &n);

    for (i__ = 1; i__ <= 7; ++i__) {
         printf("\n  %20.12e \n", r__[i__-1]);
    }
    printf("\n  %20.12e \n", s);

    ash0d_c(a, r__, d__, s1, &n, &c__2);

    for (i__ = 1; i__ <= 7; ++i__) {
         printf("\n  %20.12e \n", d__[i__-1]);
    }
        i__2 = n;
        for (i__ = 1; i__ <= i__2; ++i__) {
/* l2: */
            x[i__ - 1] += d__[i__ - 1] * s;
        }
        z__ = 0.;
        i__2 = n;
        for (i__ = 1; i__ <= i__2; ++i__) {
            f = (d__1 = x[i__ - 1] - x1[i__ - 1], abs(d__1));
            x1[i__ - 1] = x[i__ - 1];
            if (f <= z__) {
                goto l5;
            }
            z__ = f;
l5:
            ;
        }
        if (s0 + z__ == s0) {
            goto l6;
        }
/* l9: */
/*    } */
l6:
    for (i__ = 1; i__ <= 7; ++i__) {
         printf("\n  %20.12e \n", x[i__-1]);
    }
/* l21: */
    return 0;
} /* main */


Подпрограмма ash0d_c (a, b, x, s, n, p) находит решение системы ax = b, причем при повторном решении системы с той же матрицей и другой правой частью параметр  p полагают отличным от единицы. Проводится два шага уточнения решения системы с матрицей Гильберта 7 порядка и правой частью  b = (0, 0, 0, 0, 1, 0, 0).


Результат:

                x1                                x2                                x3
    4.85115043411 + 04     4.85099999532 + 04    4.85100000000 + 04
  - 1.94046009606 + 06  - 1.94039999813 + 06  - 1.94040000000 + 06
    1.87115797099 + 07     1.87109999819 + 07    1.87110000000 + 07
  - 7.27672565200 + 07  - 7.27649999300 + 07  - 7.27650000001 + 07
    1.33406641488 + 08     1.33402499871 + 08    1.33402500000 + 08
  - 1.15263342213 + 08  - 1.15259759888 + 08  - 1.15259760000 + 08
    3.78389772248 + 07     3.78377999634 + 07    3.78378000000 + 07

где   x1 - решение системы без уточнения;
        x2 - решение системы после 1-ого шага уточнения;
        x3 - решение системы после 2-ого шага уточнения.