Документ взят из кэша поисковой машины. Адрес оригинального документа : http://num-anal.srcc.msu.su/lib_na/cat/am_htm_c/ammhr_c.htm
Дата изменения: Thu Dec 3 13:41:15 2015
Дата индексирования: Sun Apr 10 01:29:35 2016
Кодировка: Windows-1251
БЧА НИВЦ МГУ. AMMHR_C. Умножение матриц
Текст подпрограммы и версий
ammhr_c.zip , ammhd_c.zip , ammhc_c.zip
Тексты тестовых примеров
tammhr_c.zip , tammhd_c.zip , tammhc_c.zip

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

Назначение

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

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

Для заданной в компактной форме ленточной вещественной

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

      RCOND  =  1 / ( || A ||1 - || A-1 ||1 )  ,   где
                                           N
        || A ||1 =      max     {   ∑   | ai j |   }
                         j=1,...N       i=1 

и затем для заданной вещественной прямоугольной матрицы B размера N на NN вычисляется прямоугольная матрица С = А- 1 В размера N на NN путем решения NN систем линейных алгебраических уравнений А*С (J) = В(J), где С(J) и В(J) суть J - е столбцы матриц C и B,  J = 1,...,NN.

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

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

    int ammhr_c (real *a, integer *ma, integer *n, integer *ml,
        integer *mu, integer *nlead, real *b, integer *mb, integer *nn,
        real *rcond, real *z, integer *ierr)

Параметры

a - вещественный двумерный массив размера ma*n, в первых ml+mu+1 столбцах которого задается в компактном виде ленточная матрица A порядка n; на выходе в первых ml столбцах массива находятся нижние кодиагонали ленточной матрицы L1 *...* Ln - 1 (Li, i = 1,...,n - 1, суть элементарные матрицы исключения метода Гаусса), в следующих ml+mu+1 столбцах содержится в компактном виде матрица U;
ma - первая размерность массива a в вызывающей программе (тип: целый);
n - порядок матрицы A и число строк матриц B и C (тип: целый);
ml - число нижних кодиагоналей матрицы A (тип: целый);
mu - число верхних кодиагоналей матрицы A (тип: целый);
nlead - целый вектор длины n, содержащий на выходе информацию о выполненных в ходе факторизации перестановках (см. замечания по использованию);
b - вещественный двумерный массив размера mb*nn, в котором задается матрица B; на выходе на соответствующих местах находятся элементы матрицы C = A- 1 B;
mb - первая размерность массива b в вызывающей программе (тип: целый);
nn - число столбцов матрицы B и C (тип: целый);
rcond - вещественная переменная, содержащая на выходе вычисленное значение величины 1 / (|| A ||1 - || A- 1 ||1) (см. замечания по использованию);
z - вещественный рабочий вектор длины n;
ierr - целая переменная, содержащая на выходе информацию о прохождении счета, при этом:
ierr=65 - если хотя бы одна из переменных ma, n, mb, nn имеет значение, меньшее единицы;
ierr=66 - если в процессе счета произошло переполнение (это говорит о том, что либо || A ||1, либо некоторые элементы матрицы U или матрицы C превосходят по абсолютной величине максимальное представимое на данной машине число);
ierr= -k - если в результате факторизации в k - ой строке матрицы U диагональный элемент равен нулю (это свидетельствует о вырождености матрицы A). Если таких строк у матрицы U несколько, то значение k полагается равным номеру последней из них;
ierr=67 - если для некоторого j, 1 ≤ j ≤ nn, система A*C (j) = B(j) несовместна.

Версии

ammhd_c - умножение матрицы, обратной к заданной с удвоенной точностью в компактной форме вещественной ленточной матрице A, на вещественную прямоугольную матрицу, заданную с удвоенной точностью, и оценка числа обусловленности матрицы A;
ammhc_c - умножение матрицы, обратной к заданной в компактной форме комплексной ленточной матрице A, на комплексную прямоугольную матрицу и оценка числа обусловленности матрицы A.

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

afb2r_c - подпрограмма треугольной факторизации и оценки числа обусловленности ленточной матрицы A.
utafsi_c - подпрограмма выдачи диагностических сообщений.

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

  1. 

В подпрограмме ammhd_c массивы a, b, z и переменная rcond имеют тип double, для треугольного разложения и оценки числа обусловленности ленточной матрицы A вызывается подпрограмма afb2d_c.

  2. 

В подпрограмме ammhc_c массивы a, b и z имеют тип complex, для треугольного разложения и оценки числа обусловленности ленточной матрицы A вызывается подпрограмма afb2c_c.

  3. 

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

  4. 

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

  5.  Если вырабатывается значение переменной ierr, отличное от нуля, то выдается соответствующее диагностическое сообщение и если ierr > 0, то происходит выход из подпрограммы.

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

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

    /* Local variables */
    static int ierr;
    static float a[25] /* was [5][5] */,
                 b[50] /* was [5][10] */;
    static int i__, j, k, n, nlead[5];
    static float z__[5];
    extern int ammhr_c(float *, int *, int *, int *, int *, int *,
                       float *, int *, int *, float *, float *, int *);
    static float rcond;
    static int j0, j1, ma, mb, ml, nn, mu;

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

    ma = 5;
    n = 5;
    ml = 1;
    mu = 1;
    mb = 5;
    nn = 10;
    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;
            b_ref(i__, j) = 0.f;
            b_ref(i__, j + 5) = 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);
            b_ref(i__, j) = a_ref(i__, k);
            b_ref(i__, j + 5) = a_ref(i__, k);
/* L3: */
        }
/* 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));
    }
    for (i__ = 1; i__ <= 5; ++i__) {
         printf("\n  %12.4e %12.4e %12.4e %12.4e %12.4e \n",
                b_ref(i__, 1), b_ref(i__, 2), b_ref(i__, 3),
                b_ref(i__, 4), b_ref(i__, 5));
         printf("\n  %12.4e %12.4e %12.4e %12.4e %12.4e \n",
                b_ref(i__, 6), b_ref(i__, 7), b_ref(i__, 8),
                b_ref(i__, 9), b_ref(i__, 10));
    }
    ammhr_c(a, &ma, &n, &ml, &mu, nlead, b, &mb, &nn, &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 %12.4e %12.4e %12.4e %12.4e \n",
                b_ref(i__, 1), b_ref(i__, 2), b_ref(i__, 3),
                b_ref(i__, 4), b_ref(i__, 5));
         printf("\n  %12.4e %12.4e %12.4e %12.4e %12.4e \n",
                b_ref(i__, 6), b_ref(i__, 7), b_ref(i__, 8),
                b_ref(i__, 9), b_ref(i__, 10));
    }
    return 0;
} /* main */


Результаты:

                 |   0         21.0     22.0   23.0    0 |
                 | -0.524   32.0     33.0   34.0    0 |
   a_ref  =  | -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 ,

                 | 1.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0 |
                 | 0.0  1.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0 |
   b_ref  =  | 0.0  0.0  1.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0 |
                 | 0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  1.0  0.0 |
                 | 1.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  1.0 |