Текст подпрограммы и версий
ass1r_c.zip , ass1d_c.zip 		Тексты тестовых примеров
tass1r_c.zip , tass1d_c.zip

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

Назначение

Решение систем линейных алгебраических уравнений с разреженными матрицами методом сопряженных градиентов.

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

ass1r_c вычисляет решение системы линейных алгебраических уравнений Ax = b (где  A - квадратная матрица общего вида, в том числе разреженная, порядка  N) с заданной точностью EPS методом сопряженных градиентов, а также сумму квадратов компонент вектора невязки Ax - b .

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

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

    int ass1r_c (real *b, integer *n, S_fp suba, S_fp subat,
                 real *x, real *rsq, integer *itmax, real *eps, real *g, real *h,
                 real *xi, real *xj, integer *ierr)

   Параметры
b - вещественный вектор длины  n, содержащий компоненты вектора правой части системы;
n - порядок матрицы системы (тип: целый);
suba - подпрограмма вычисления вектора  ax, первый оператор которой имеет вид:
int suba (float *x, float *ax),
где  x - вещественный вектор длины  n, содержащий текущий вектор решения  x;  ax - вещественный вектор длины  n, содержащий результирующий вектор  ax;
subat - подпрограмма вычисления вектора atx, первый оператор которой имеет вид:
int subat (float *x, float *atx),
где  x - вещественный вектор длины  n, содержащий текущий вектор решения  x;  atx - вещественный вектор длины  n, содержащий результирующий вектор atx;
x - вещественный вектор длины  n, содержащий на входе в подпрограмму начальное приближение к решению системы, а на выходе - вычисленное решение системы;
rsq - вещественная переменная, содержащая вычисленную сумму квадратов компонент вектора невязки ax - b;
itmax - заданное максимальное допустимое количество итераций метода сопряженных градиентов (тип: целый);
eps - заданная точность, с которой необходимо вычислить решение системы (тип: вещественный);
        g, h, -
      xi, xj   		вещественные векторы длины  n, используемые в подпрограмме в качестве рабочих;
ierr - целая переменная, служащая для сообщения об ошибках, обнаруженных в ходе работы подпрограммы; при этом
ierr=65 - когда матрица системы вырождена;
ierr=66 - когда заданная точность не может быть достигнута за максимальное количество итераций. 
   Версии
ass1d_c - решение систем линейных алгебраических уравнений с разреженными матрицами методом сопряженных градиентов в режиме удвоенной точности. При этом параметры b, x, rsq, eps, g, h, xi, xj должны иметь тип double. 
   Вызываемые подпрограммы
utas10_c - подпрограмма выдачи диагностических сообщений при работе подпрограммы ass1r_c.
utas11_c - подпрограмма выдачи диагностических сообщений при работе подпрограммы ass1d_c. 
   Замечания по использованию
 

В качестве начального приближения к решению системы может быть взят нулевой вектор. Все особенности структуры матрицы  a системы должны учитываться в подпрограммах suba и subat.

Если значение rsq велико, то это означает, что матрица системы близка к вырожденной и полученный вектор  x представляет собой наилучшее приближение к решению в смысле наименьших квадратов.

В подпрограммах ass1r_c и ass1d_c имеется внешняя структура с именем ass1rr_ , содержащая элемент целого типа iter. По окончании работы этих подпрограмм значение переменной iter полагается равной количеству выполненных итераций метода сопряженных градиентов.

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

struct {
    int iter;
} ass1rr_;

#define ass1rr_1 ass1rr_

struct {
    float a[100] /* was [10][10] */;
} _BLNK__;

#define _BLNK__1 _BLNK__

int main(void)
{
    /* System generated locals */
    int i__1, i__2;

    /* Local variables */
    extern int suba_c();
    static int ierr;
    extern int ass1r_c(float *, int *, U_fp, U_fp, float *, float *, int *,
                       float *, float *, float *, float *, float *, int *);
    static float b[10], g[10], h__[10];
    static int i__, j, n;
    static float r__, x[10];
    extern int subat_c();
    static int itmax;
    static float xi[10], xj[10], eps, rsq;

#define a_ref(a_1,a_2) _BLNK__1.a[(a_2)*10 + a_1 - 11]

    n = 10;
    i__1 = n;
    for (i__ = 1; i__ <= i__1; ++i__) {
        x[i__ - 1] = 0.f;
/* l1: */
        b[i__ - 1] = 1.f;
    }
    r__ = .1f;
    i__1 = n;
    for (i__ = 1; i__ <= i__1; ++i__) {
        i__2 = n;
        for (j = 1; j <= i__2; ++j) {
/* l2: */
            a_ref(i__, j) = 0.f;
        }
    }
    i__2 = n;
    for (i__ = 1; i__ <= i__2; ++i__) {
        a_ref(i__, i__) = r__;
/* l3: */
        r__ += .1f;
    }
    eps = 1e-5f;
    itmax = n * 10;
    ass1r_c(b, &n, (U_fp)suba_c, (U_fp)subat_c, x, &rsq, &itmax, &eps, g,
            h__, xi, xj, &ierr);

    for (i__ = 0; i__ <= 5; i__+=5) {
         printf("\n  %15.5e %15.5e %15.5e %15.5e %15.5e \n",
                x[i__], x[i__+1], x[i__+2], x[i__+3], x[i__+4]);
    }
    printf("\n  %15.5e \n", rsq);
    printf("\n  %5i \n", ass1rr_1.iter);
    return 0;
} /* main */


int suba_c(float *x, float *v)
{
    /* System generated locals */
    int i__1, i__2;

    /* Local variables */
    static int i__, j, n;

#define a_ref(a_1,a_2) _BLNK__1.a[(a_2)*10 + a_1 - 11]

    /* Parameter adjustments */
    --v;
    --x;

    /* Function Body */
    n = 10;
    i__1 = n;
    for (i__ = 1; i__ <= i__1; ++i__) {
        v[i__] = 0.f;
        i__2 = n;
        for (j = 1; j <= i__2; ++j) {
/* l1: */
            v[i__] += a_ref(i__, j) * x[j];
        }
    }
    return 0;
} /* suba_c */


int subat_c(float *x, float *v)
{
    /* System generated locals */
    int i__1, i__2;

    /* Local variables */
    static int i__, j, n;

#define a_ref(a_1,a_2) _BLNK__1.a[(a_2)*10 + a_1 - 11]

    /* Parameter adjustments */
    --v;
    --x;

    /* Function Body */
    n = 10;
    i__1 = n;
    for (i__ = 1; i__ <= i__1; ++i__) {
        v[i__] = 0.f;
        i__2 = n;
        for (j = 1; j <= i__2; ++j) {
/* l1: */
            v[i__] += a_ref(j, i__) * x[j];
        }
    }
    return 0;
} /* subat_c */


Результаты: 

      x  =  ( 10 , 5 , 10/3 , 2.5 , 2 , 5/3 , 10/7 , 1.25 , 10/9 , 1 )

      rsq  = 0.19785e - 8 
      iter = 11