БЧА НИВЦ МГУ. ASS2R_C. Решение систем с матрицами специального вида (невырожденными)

Текст подпрограммы и версий
ass2r_c.zip , ass2d_c.zip

Тексты тестовых примеров
tass2r_c.zip , tass2d_c.zip

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

Назначение

Решение невырожденной разреженной линейной системы итерационным методом Гаусса - Зейделя (матрица системы представлена в формате RR (LU) U) .

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

Сокращенное название формата RR (LU) U происходит от английского словосочетания "Row - wise Representation, Lower - Upper, Unordered" (строчное представление, нижний треугольник - верхний треугольник, неупорядоченное).

Данный формат используется для прямоугольных (квадратных) матриц, у которых все или большинство диагональных элементов не равны нулю. Матрица A в этом формате представляется в виде суммы L + D + U, где L - нижняя треугольная, D - диагональная, U - верхняя треугольная матрицы. Диагональные элементы хранятся в отдельном одномерном массиве AD, а элементы матриц L и U содержатся в одномерном массиве AN и связаны между собой списками (портретами) IA и JA.

Например, пусть задана матрица

               |    4    0    1    0    2    |
               |    1    2    0    0    0    |
   A   =    |    0    0    2    1    0    |
               |    1    1    0    1    1    |
               |    0    0    0    0    16  |

Тогда в формате RR (LU) U данная матрица с точностью до упорядоченности ее ненулевых внедиагональных элементов имеет вид:

                IA  =  ( 1; 3; 4; 5; 8; 8 )
                JA  =  ( 5, 3; 1; 4; 5, 1, 2 )
               AN  =  ( 2, 1, 1, 1, 1, 1, 1 )
               AD  =  ( 4, 2, 2, 1, 16 )

Пусть квадратная матрица A порядка n линейной системы Ax = b задана в формате RR (LU) U и пусть известно, что все диагональные элементы матрицы A не равны нулю.

Метод Гаусса - Зейделя заключается в следующем. Перепишем i - е уравнение системы в виде:

         _{i -1}                                  _n
          ∑   a_{i j} x_j  +  a_{i i} x_i  +   ∑     a_{i j} x_j  =  b_i 
         ^{j =1}                               ^{j= i +1}

Поскольку по предположению все a_i i ≠ 0, то это уравнение можно записать в виде:

                             _{i -1}                  _n
          x  =  ( b_i -   ∑   a_{i j} x_j -    ∑     a_{i j} x_j ) / a_{i i}
                            ^{j =1}               ^{j= i +1}

Это уравнение решается итерационно, если задать начальное приближение к решению x_j ⁽¹⁾ ( j = 1, 2, ..., n). В случае, когда начальное приближение x_j ⁽¹⁾ заранее неизвестно, то можно положить

        x_j⁽¹⁾  =  b_j / a_{j j} ,   j = 1, 2, ..., n .

Следует отметить, что скорость сходимости итерационного процесса Гаусса - Зейделя может быть значительно увеличена, если начальное приближение к решению выбрано удачно.

Далее на m - м итерационном шаге процесса Гаусса - Зейделя имеем (m = 2, 3, 4, ...)

                                _{i -1}
     x_i^(m)  =  ( b_i  -     ∑    a_{i j} x_j^(m) -
                               ^{j =1}
                                                     _n
                                                 -  ∑     a_{i j} x_j^{(m -1)} ) / a_{i i}
                                                   ^{j= i +1}

Заметим, что здесь вновь вычисленные компоненты решения на m - м шаге сразу же используются для вычисления следующих компонент.

В некоторых случаях сходимость метода Гаусса - Зейделя можно ускорить, если использовать прием верхней релаксации:

               x_i^(m)  =  x_i^{(m -1)}  +  Q( x_i^(m)  -  x_i^{(m -1)} )

Релаксационный множитель Q обычно выбирается в пределах от 1 до 2. Если матрица положительно определенная, то Q выбирается в пределах от 0 < Q < 1. Иногда в качестве Q используется диагональная матрица, позволяющая подбирать релаксационные множители индивидуально для каждого уравнения системы.

Метод Гаусса - Зейделя хорошо сходится, если матрица A является нижней или "почти нижней" треугольной. Для его сходимости не требуется диагонального преобладания в матрице A. Необходимым и достаточным условием этого метода является следующее: все корни уравнения

               |   a_{1 1}λ    a_{1 2}      a_{1 3}    ...  a_{1 n}    |   
               |   a_{2 1}λ    a_{2 2}λ    a_{2 3}    ...  a_{2 n}    |
       det   |   .  .  .  .  .  .  .  .  .  .  .  .  .  .  .         =   0
               |   a_{n 1}λ    a_{n 2}λ    a_{n 3}λ  ... a_{n n}λ  |

должны быть по модулю меньше 1. Если матрица A является вещественной симметричной положительно определенной матрицей, то метод Гаусса - Зейделя всегда сходится. Скорость сходимости будет быстрее, если матрица A близка к диагональной.

Описанный метод Гаусса - Зейделя реализован в виде подпрограммы ass2r_c. В этой подпрограмме не требуется задавать на входе начальное приближение к решению: оно вычисляется в ass2r_c по формуле x_i ⁽¹⁾ = b_i /a_i i. Подпрограмма предполагает задание максимально допустимого количества итераций. Контроль точности ведется по абсолютной погрешности EPS: текущая итерация считается приемлемым приближением к решению, если | x_i ^(m) - x_i ^(m - 1) | < EPS для всех i.

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

Н.С.Бахвалов. Численные методы. Изд - во "Наука", 1973.

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

    int ass2r_c (integer *ia, integer *ja, real *an, real *ad,
            integer *n, real *b, real *x, real *q, real *eps, integer *itmax, 
            integer *iflag)

   Параметры

ia, ja, - an	заданные портрет и ненулевые внедиагональные элементы матрицы A в формате RR (LU) U;
ad -	вещественный одномерный массив длины n, содержащий диагональные элементы матрицы A; требуется, чтобы ни один диагональный элемент не был равен нулю;
n -	заданный порядок матрицы A (тип: целый);
b -	вещественный одномерный массив длины n, содержащий вектор правой части системы;
x -	вещественный одномерный массив длины n, на выходе из подпрограммы содержащий вычисленные компоненты вектора решения системы;
q -	заданный релаксационный множитель (тип: вещественный);
eps -	заданная абсолютная погрешность, с которой требуется вычислить решение;
itmax -	заданное максимальное количество итераций по методу Гаусса - Зейделя;
iflag -	целая переменная, служащая для сообщения о режиме окончания работы подпрограммы:

iflag=0 -	решение системы вычислено с заданной точностью;
iflag=1 -	решение не получено с заданной точностью за заданное максимальное количество итераций.

   Версии

ass2d_c -

решение невырожденной разреженной линейной системы итерационным методом Гаусса - Зейделя в режиме удвоенной точности; при этом параметры an, ad, b, x, q и eps должны иметь тип double.

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

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

Подпрограммы ass2r_c и ass2d_c имеют внешнюю структуру с именем ass2rr_ , содержащую элемент целого типа iter. Значение переменной iter полагается равным количеству в действительности выполненных итераций, если iflag = 0. Если iflag = 1, то iter = itmax.

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

struct {
    int iter;
} ass2rr_;

#define ass2rr_1 ass2rr_

int main(void)
{
    /* Initialized data */
    static int ia[6] = { 1,2,3,5,6,8 };
    static int ja[7] = { 5,1,2,1,2,3,1 };
    static float an[7] = { 1.f,1.f,1.f,1.f,1.f,1.f,2.f };
    static float ad[5] = { 4.f,2.f,2.f,8.f,16.f };
    static float b[5] = { 1.f,1.f,1.f,1.f,1.f };

    /* System generated locals */
    int i__1;

    /* Local variables */
    extern int ass2r_c(int *, int *, float *, float *, int *, float *,
                       float *, float *, float *, int *, int *);
    static int i__, n;
    static float q;
    static int iflag;
    static float x[5];
    static int itmax;
    static float eps;

    n = 5;
    q = 1.5f;
    i__1 = n;
    for (i__ = 1; i__ <= i__1; ++i__) {
/* l1: */
        x[i__ - 1] = b[i__ - 1] / ad[i__ - 1];
    }
    eps = .001f;
    itmax = 500;
    ass2r_c(ia, ja, an, ad, &n, b, x, &q, &eps, &itmax, &iflag);

    printf("\n  %5i \n", iflag);
    printf("\n  %5i \n", ass2rr_1.iter);
    printf("\n  %15.6e %15.6e %15.6e %15.6e %15.6e \n",
           x[0], x[1], x[2], x[3], x[4]);
    return 0;
} /* main */


Результаты: 

       iflag = 0 ,   iter = 7 ,

       x = (0.245396, 0.377041, 0.188364, 0.0778308, 0.0203379)