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

Назначение

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

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

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

           rcond = 1/( || A ||1 * || A-1 ||1 ) ,

 где   || A ||1 = maxj = 1, ..., N { | a1 j | + | a2 j | + ...+ | aN j | } ,

а затем вычисляется определитель матрицы А как произведение диагональных элементов матрицы U, умноженное на (-1)^I, где I - число выполненных в процессе факторизации перестановок.

Определитель записывается в виде:

     det A = D1*10D2 ,   где 1.0 ≤ D1 < 10.0 

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

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

    int adb1r_c(real *a, integer *ma, integer *n, integer *ml,
        integer *mu, integer *nlead, real *det1, real *det2, real *rcond, 
        real *z__, integer *ierr)

   Параметры 

   Версии 

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

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

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

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

#define a_ref(a_1,a_2) a[(a_2)*5 + a_1 - 6]

    ma = 5;
    n = 5;
    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;
/* l2: */
        }
/* l1: */
    }
    ml = 1;
    mu = 1;
    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 (j = 1; j <= 5; ++j) {
        printf("\n  %12.4e %12.4e %12.4e %12.4e %12.4e \n",
               a_ref(j, 1), a_ref(j, 2), a_ref(j, 3),
               a_ref(j, 4), a_ref(j, 5));
    }
    adb1r_c(a, &ma, &n, &ml, &mu, nlead, &det1, &det2, &rcond, z__, &ierr);

        printf("\n  %5i %5i %5i %5i %5i \n",
                   nlead[0], nlead[1], nlead[2], nlead[3], nlead[4]);
    for (j = 1; j <= 5; ++j) {
        printf("\n  %12.4e %12.4e %12.4e %12.4e %12.4e \n",
               a_ref(j, 1), a_ref(j, 2), a_ref(j, 3),
               a_ref(j, 4), a_ref(j, 5));
    }
        printf("\n  %5i %16.7e \n", ierr, rcond);
        printf("\n  %16.7e %16.7e \n", det1, det2);
    return 0;
} /* main */


Результаты:

           |    0         21.0     22.0     23.0     0 |
           | -0.524    32.0     33.0     34.0     0 |
    a = | -0.015    43.0     44.0     45.0     0 |
           |  0.292    54.0     55.0        0       0 |
           | -0.228    0.569      0         0       0 |

       nlead = (2, 3, 4, 5, 5)
       rcond = 1.47362E-03
       det1 = 8.88360
       det2 = 5.0