Текст подпрограммы и версий
afb2r_c.zip , afb2d_c.zip , afb2c_c.zip 		Тексты тестовых примеров
tafb2r_c.zip , tafb2d_c.zip , tafb2c_c.zip

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

Назначение

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

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

Для заданной в компактной форме ленточной вещественной матрицы А порядка N выполняется треугольная факторизация А = LU, где U - верхняя треугольная ленточная матрица, причем матрица L- 1 имеет вид: 

     L-1 = LN-1 * PN-1 *...* L1 * P1 ,

где Рi,  i = 1, ..., N - 1, суть матрицы перестановок, обеспечивающие стратегию выбора ведущего элемента по столбцам; Li,  i = 1, ..., N - 1, суть элементарные матрицы исключения в методе Гаусса. Все матрицы Li являются нижними треугольными ленточными матрицами с единичными диагональными элементами. После выполнения факторизации вычисляется величина RСОND, обратная числу обусловленности матрицы А: 

 rcond = 1 / (|| A ||1 * || A-1 ||1) ,   где  || A ||1 = max j = 1,...,N  Sj ,
здесь     N
     Sj =  ∑ | ai j |
             i=1

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

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

    int afb2r_c (real *a, integer *ma, integer *n, integer *ml,
            integer *mu, integer *nlead, real *rcond, real *z, integer *ierr)

   Параметры
a - вещественный двумерный массив размера ma*n, в первых ml + mu + 1 столбцах которого задается в компактном виде исходная ленточная матрица порядка N; на выходе в первых ml столбцах массива находятся нижние кодиагонали ленточной матрицы L1 *...* LN - 1, в следующих mu + ml + 1 столбцах содержится в компактном виде матрица U;
ma - первая размерность массива a в вызывающей программе (тип: целый);
n - порядок матрицы A (тип: целый);
ml - число нижних кодиагоналей матрицы A (тип: целый);
mu - число верхних кодиагоналей матрицы A (тип: целый);
nlead - целый вектор длины n, содержащий на выходе информацию о выполненных в процессе исключения перестановках (см. замечания по использованию);
rcond - вещественная переменная, содержащая на выходе вычисленное значение величины, обратной числу обусловленности матрицы A (см. замечания по использованию);
z - вещественный рабочий вектор длины n;
ierr - целая переменная, содержащая на выходе информацию о прохождении счета, при этом:
ierr=65 - если ma ≤ 0 или n ≤ 0;
ierr=66 - если в процессе работы произошло переполнение (это говорит о том, что либо || A ||1, либо некоторые элементы матрицы U превосходят по абсолютной величине максимальное представимое на данной машине число);
ierr=-k - если в результате факторизации диагональный элемент в K - й строке матрицы U равен нулю (это свидетельствует о вырожденности матрицы A). Если таких строк у матрицы u несколько, то значение K полагается равным номеру последней из них. 
   Версии
afb2d_c - треугольное разложение методом Гаусса с выбором ведущего элемента по столбцу и оценка числа обусловленности ленточной вещественной матрицы, заданной с удвоенной точностью в компактной форме.
afb2c_c - треугольное разложение методом Гаусса с выбором ведущего элемента по столбцу и оценка числа обусловленности ленточной комплексной матрицы, заданной в компактной форме. 
   Вызываемые подпрограммы
utafsi_c - подпрограмма выдачи диагностических сообщений. 
   Замечания по использованию
  1. 

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

  2. 

В подпрограмме afb2d_c массивы a, z и переменная rcond имеют тип double.

  3. 

На выходе K - й элемент вектора nlead равен номеру строки, переставленной на K - м шаге факторизации с K - й строкой матрицы A. Поскольку факторизация Гаусса требует n - 1 шагов, то nlead (n) = n.

  4. 

Так как в результате выполненных в ходе факторизации перестановок число верхних кодиагоналей матрицы U равно mu + ml, а также в силу некоторых конструктивных особенностей подпрограммы, для правильной ее работы необходимо выполнение условия ma ≥ n >mu + 2*ml + 1. Если mu + 2*ml + 1 ≥ n, то более целесообразно, задав матрицу A не в компактной, а в полной форме, обратиться к подпрограмме afg4r_c.

  5.  Если переменной ierr присвоено значение, отличное от нуля, то выдается соответствующее диагностическое сообщение, полагается RCOND = 0.0 и происходит выход из подпрограммы (если ierr < 0, то выход происходит по окончании факторизации).

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

int main(void)
{
    /* Local variables */
    static int ierr;
    extern int afb2r_c(float *, int *, int *, int *, int *, int *,
                       float *, float *, int *);
    static float a[81] /* was [9][9] */;
    static int i__, j, k, n, nlead[9];
    static float z__[9], rcond;
    static int j0, j1, ma, ml, mu;
    int i__1, i__2, i__3;

#define a_ref(a_1,a_2) a[(a_2)*9 + a_1 - 10]

    ma = 9;
    n = 9;
    i__1 = ma;
    for (i__ = 1; i__ <= i__1; ++i__) {
        i__2 = n;
        for (j = 1; j <= i__2; ++j) {
            a_ref(i__, j) = 0.f;
/* l1: */
        }
/* l2: */
    }
    ml = 2;
    mu = 3;
    i__1 = ma;
    for (i__ = 1; i__ <= i__1; ++i__) {
/* Computing max */
        i__2 = 1, i__3 = i__ - ml;
        j0 = max(i__2,i__3);
/* Computing min */
        i__2 = n, i__3 = i__ + mu;
        j1 = min(i__2,i__3);
        i__2 = j1;
        for (j = j0; j <= i__2; ++j) {
            k = j - i__ + ml + 1;
            a_ref(i__, k) = (float) (i__ * 10 + j);
/* l3: */
        }
/* l4: */
    }
    for (i__ = 1; i__ <= 9; ++i__) {
         printf("\n %13.4e %13.4e %13.4e \n %13.4e %13.4e %13.4e \n",
              a_ref(i__, 1), a_ref(i__, 2), a_ref(i__, 3),
              a_ref(i__, 4), a_ref(i__, 5), a_ref(i__, 6));
         printf("\n %13.4e %13.4e %13.4e \n",              
              a_ref(i__, 7), a_ref(i__, 8), a_ref(i__, 9));
    }
    afb2r_c(a, &ma, &n, &ml, &mu, nlead, &rcond, z__, &ierr);

    printf("\n  %5i %5i %5i %5i %5i %5i %5i %5i %5i \n",
            nlead[0], nlead[1], nlead[2], nlead[3], nlead[4], nlead[5],
            nlead[6], nlead[7], nlead[8]);
    for (i__ = 1; i__ <= 9; ++i__) {
         printf("\n %13.4e %13.4e %13.4e \n %13.4e %13.4e %13.4e \n",
              a_ref(i__, 1), a_ref(i__, 2), a_ref(i__, 3),
              a_ref(i__, 4), a_ref(i__, 5), a_ref(i__, 6));
         printf("\n %13.4e %13.4e %13.4e \n",              
              a_ref(i__, 7), a_ref(i__, 8), a_ref(i__, 9));
    }
    printf("\n  %5i %15.7e \n", ierr, rcond);
    return 0;
} /* main */


Результат:

                    |   0          0         31.0       32.0    33.0   34.0   35.0  36.0  0 |
                    |   0        -0.68     42.0       43.0    44.0   45.0   46.0  47.0  0 |
                    | -0.35    -0.02     53.0       54.0    55.0   56.0   57.0  58.0  0 |
                    | -0.008  -0.006   64.0       65.0    66.0   67.0   68.0  69.0  0 |
      a_ref  =  | -0.01    -0.01     75.0       76.0    77.0   78.0   79.0     0    0 |
                    | -0.005  -0.004   86.0       87.0    88.0   89.0     0        0    0 |
                    | -0.19      0.002   97.0       98.0    99.0     0       0        0    0 |
                    |  0.3       -0.25     0.76       15.02    0        0       0        0    0 |
                    | -0.13     -0.51     0.0003     0         0        0       0        0    0 |

      nlead  =   (3, 4, 5, 6, 7, 8, 9, 9, 9).

      rcond  =  1.2017e-07