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

Назначение

Решение вещественной системы линейных алгебраических уравнений A*X = b или A^T*X = b методом Гаусса с выбором ведущего элемента по столбцу и оценка числа обусловленности ленточной матрицы A, заданной в компактной форме.

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

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

     RCOND = 1 / ( || A ||1 * || A-1 ||1 ) ,
  где
     || A ||1 = maxj = 1,...,N  ( a1j + a2j + ... + aNj ) , 

и затем решается система А * Х = b или А^T*Х = b.

Вычисления проводятся в два этапа:

     1) решается система  L * Y = b  (для AT решается UT*Y=b) ;
     2) решается система U * X = Y  (для AT решается LT*X=Y) . 

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

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

    int asb2r_c (real *a, integer *ma, integer *n, integer *ml,
            integer *mu, integer *nlead, real *b, integer *ltr, integer *l,
            real *rcond, real *z, integer *ierr)

   Параметры 

   Версии 

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

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

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

int main(void)
{
    /* Local variables */
    static int ierr;
    extern int asb2r_c(float *, int *, int *, int *, int *, int *,
                       float *, int *, int *, float *, float *, int *);
    static float a[25] /* was [5][5] */, b[5];
    static int i__, j, k, l, n, nlead[5];
    static float z__[5], rcond;
    static int j0, j1, ma, ml, mu, ltr;
    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;
    ml = 1;
    mu = 1;
    l = 0;
    ltr = 0;
    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: */
    }
    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: */
        }
        b[i__ - 1] = 7.f;
/* l4: */
    }
    for (i__ = 1; i__ <= 5; ++i__) {
         printf("\n  %12.4e %12.4e %12.4e %12.4e %12.4e \n",
                a_ref(i__, 1), a_ref(i__, 2), a_ref(i__, 3),
                a_ref(i__, 4), a_ref(i__, 5));
    }
    asb2r_c(a, &ma, &n, &ml, &mu, nlead, b, <r, &l, &rcond, z__, &ierr);

    printf("\n  %5i %5i %5i %5i %5i \n",
           nlead[0], nlead[1], nlead[2], nlead[3], nlead[4]);
    for (i__ = 1; i__ <= 5; ++i__) {
         printf("\n  %12.4e %12.4e %12.4e %12.4e %12.4e \n",
                a_ref(i__, 1), a_ref(i__, 2), a_ref(i__, 3),
                a_ref(i__, 4), a_ref(i__, 5));
    }
    printf("\n  %5i \n", ierr);
    printf("\n  %16.7e \n", rcond);
    for (i__ = 1; i__ <= 5; ++i__) {
         printf("\n  %12.4e \n", b[i__-1]);
    }
    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 , 

      b  =  (-7.0625,  7.0598,  0.0024,  -6.4409,  6.4511) .