Текст подпрограммы и версий
asp1r_c.zip , asp1d_c.zip , asp1c_c.zip 		Тексты тестовых примеров
tasp1r_c.zip , tasp1d_c.zip , tasp1c_c.zip

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

Назначение

Нахождение нормального решения неопределенной системы линейных алгебраических уравнений полного ранга методом отражений.

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

Находится нормальное решениe cиcтемы АХ = В, где А - матрица полного ранга размера N на М (N ≤ М), В - заданный вектор длины N. Для решения используется нормализованное приведение матрицы системы к нижней треугольной форме L с помощью последовательности преобразований отражения 

 (1)     SAQ1Q2...QN = L ,

где Qi - соответствующие матрицы отражения, 1 ≤ i ≤ N , S - результирующая матрица перестановок, L - нижняя треугольная матрица размера N на М.

Из полученной треугольной системы 

        LY = SB

находится нормальное решение Y, по которому затем определяется искомое решение Х = Q1Q2 ... QNY .

Д.К.Фадеев, В.Н.Кублановская, В.Н.Фадеева, О решении линейных алгебраических систем с прямоугольными матрицами. Тp. Мат. ин-та АН СССР, 1968, 96.

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

    int asp1r_c (real *a, real *b, real *x, real *t, integer *s,
            integer *n, integer *m, integer *l)

   Параметры
a - вещественный двумерный массив размера n на m (n ≤ m), в котором задается исходная матрица A; в результате работы подпрограммы в массиве a на соответствуюших местах запоминаются поддиагональные элементы вычисленной нижней треугольной матрицы L; в остальной части массива в последовательных строках запоминаются векторы, порождающие соответствующие матрицы отражения;
b - вещественный вектор длины n, в котором задается правая часть системы;
x - вещественный вектор длины m , в котором запоминается найденное решение системы;
t - вещественный вектор длины m, используемый подпрограммой как рабочий; в результате работы подпрограммы в первых n компонентах массива t запоминаются диагональные элементы нижней треугольной матрицы L;
is - целый вектор длины n, в котором запоминаются перестановки строк при проведении нормализованного процесса (1); в is (k) запоминается номер строки, переставленной с k - ой строкой на k - ом шаге преобразования;
n, m - число строк и столбцов исходной матрицы A, n ≤ m (тип: целый);
l - задает режим работы подпрограммы (тип: целый); при этом:
l = 1 - если система с данной матрицей A решается в первый раз;
l ≠ 1 - если система решается повторно с той же матрицей A, но с другой правой частью; в этом случае разложение (1) матрицы A не выполняется. 
   Версии
asp1d_c - нахождение нормального решения недоопределенной системы линейных алгебраических уравнений полного ранга, заданной с удвоенной точностью, методом отражений.
asp1c_c - нахождение нормального решения комплексной недоопределенной системы линейных алгебраических уравнений полного ранга методом отражений. 
   Вызываемые подпрограммы:  нет 
   Замечания по использованию
  1. 

В подпрограмме asp1c_c параметры a, b, x, t имеют тип complex.

  2. 

В подпрограмме asp1d_c параметры a, b, x, t имеют тип double.

  3.  При повторном решении системы с той же матрицей информация, полученная ранее в массивах a, t, is не должна портиться.

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

int main(void)
{
    /* Initialized data */
    static float a[20] /* was [4][5] */ = { 1.f,-2.f,3.f,-4.f,-5.f,6.f,-7.f,
            8.f,9.f,10.f,-11.f,12.f,0.f,1.f,-1.f,4.f,-5.f,-15.f,16.f,-20.f };
    static float b[12] = { 1.f,-24.f,27.f,-36.f,6.f,13.f,-13.f,16.f,6.f,-8.f,
                          10.f,-12.f };

    /* Local variables */
    extern int asp1r_c(float *, float *, float *, float *, int *, int *,
                       int *, int *);
    static int i__, j, p, s[4];
    static float t[5], x[5], a1[20] /* was [4][5] */, b1[4];
    static int i1, j1;

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

    for (i__ = 1; i__ <= 4; ++i__) {
        for (j = 1; j <= 5; ++j) {
/* l100: */
            a1_ref(i__, j) = a_ref(i__, j);
        }
    }
    j1 = 0;
    for (i1 = 1; i1 <= 3; ++i1) {
        for (i__ = 1; i__ <= 4; ++i__) {
/* l110: */
            b1[i__ - 1] = b[j1 + i__ - 1];
        }
        for (i__ = 1; i__ <= 4; ++i__) {
             printf("\n  %16.7e %16.7e %16.7e \n",
                    a_ref(i__, 1), a_ref(i__, 2), a_ref(i__, 3));
             printf("\n  %16.7e %16.7e \n",
                    a_ref(i__, 4), a_ref(i__, 5));
        }
        for (i__ = 1; i__ <=4; ++i__) {
             printf("\n  %16.7e \n", b1[i__-1]);
        }
        p = 1;
        asp1r_c(a, b1, x, t, s, &c__4, &c__5, &p);

        for (i__ = 1; i__ <= 4; ++i__) {
             printf("\n  %16.7e %16.7e %16.7e \n",
                    a_ref(i__, 1), a_ref(i__, 2), a_ref(i__, 3));
             printf("\n  %16.7e %16.7e \n",
                    a_ref(i__, 4), a_ref(i__, 5));
        }
        for (i__ = 1; i__ <= 5; ++i__) {
             printf("\n  %16.7e \n", t[i__-1]);
        }
        printf("\n  %5i %5i %5i %5i \n", s[0], s[1], s[2], s[3]);
        for (i__ = 1; i__ <= 4; ++i__) {
             printf("\n  %16.7e \n", b1[i__-1]);
        }
        for (i__ = 1; i__ <= 5; ++i__) {
             printf("\n  %16.7e \n", x[i__-1]);
        }
        for (i__ = 1; i__ <= 4; ++i__) {
            for (j = 1; j <= 5; ++j) {
/* l101: */
                a_ref(i__, j) = a1_ref(i__, j);
            }
        }
        j1 += 4;
    }
/* l102: */
    return 0;
} /* main */



Результат:

      x  =   ( 0.00000000044,  -0.069999999928,  0.799999999989,  -2.7000000002,
                  3.4999999997,  -1.5400000002 )