Текст подпрограммы и версий
ammgr_c.zip , ammgd_c.zip , ammgc_c.zip
Тексты тестовых примеров
tammgr_c.zip , tammgd_c.zip , tammgc_c.zip

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

Назначение

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

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

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

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

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

    int ammgr_c (real *a, integer *m, integer *n, integer *nlead,
        real *b, integer *mm, integer *nn, integer *ierr)

Параметры

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

Версии

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

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

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

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

  1. 

В подпрограмме ammgd_c массивы a и b имеют тип double, для треугольного разложения матрицы A вызывается подпрограмма afg5d_c.

  2. 

В подпрограмме ammgc_c массивы a и b имеют тип complex, для треугольного разложения матрицы A вызывается подпрограмма afg5c_c.

  3. 

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

  4. 

При обращении к подпрограмме необходимо соблюдение условия mm ≥ n.

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

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

int main(void)
{
    /* Initialized data */
    static float a[16] /* was [4][4] */ = { 7.9f,8.5f,4.3f,3.2f,5.6f,-4.8f,
                       4.2f,-1.4f,5.7f,.8f,-3.2f,-8.9f,-7.2f,3.5f,9.3f,3.3f };
    /* Local variables */
    static int ierr;
    static float b[32] /* was [4][8] */;
    static int i__, j, m, n, nlead[4];
    extern int ammgr_c(float *, int *, int *, int *, float *, int *,
                       int *, int *);
    static int mm, nn;

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

    m = 4;
    n = m;
    mm = m;
    nn = 8;
    for (i__ = 1; i__ <= 4; ++i__) {
        for (j = 1; j <= 4; ++j) {
            b_ref(i__, j) = a_ref(i__, j);
            b_ref(i__, j + 4) = a_ref(i__, j);
/* L1: */
        }
/* L2: */
    }
    for (i__ = 1; i__ <= 4; ++i__) {
         printf("\n  %16.7e %16.7e %16.7e %16.7e \n",
                a_ref(i__, 1), a_ref(i__, 2), a_ref(i__, 3), a_ref(i__, 4));
    }
    ammgr_c(a, &m, &n, nlead, b, &mm, &nn, &ierr);

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


Результат:
   
   nlead  =  (2, 2, 4, 4)
   
                 |  8.5            -4.8            0.8             3.5        |
   a_ref  =  | -0.92941   10.06118    4.95647  -10.45294 |
                 | -0.50588    -0.65879   -9.40171    2.40526  |
                 | -0.37647    -0.04046   -0.73072   12.65817 |
 
                 | 1  0  0  0  1  0  0  0 |
   b_ref  =  | 0  1  0  0  0  1  0  0 |
                 | 0  0  1  0  0  0  1  0 |
                 | 0  0  0  1  0  0  0  1 |