Текст подпрограммы и версий
asp0r_c.zip , asp0d_c.zip , asp0c_c.zip 		Тексты тестовых примеров
tasp0r_c.zip , tasp0d_c.zip , tasp0c_c.zip

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

Назначение

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

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

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

 (1)     QMQM-1...Q1AS = R ,

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

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

        RY = QMQM-1...Q1B

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

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

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

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

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

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

  2. 

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

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

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

int main(void)
{
    /* Initialized data */
    static float a[20] /* was [5][4] */ = { 2520.f,1260.f,840.f,630.f,504.f,
            1260.f,840.f,630.f,504.f,420.f,840.f,630.f,504.f,420.f,360.f,
            630.f,504.f,420.f,360.f,315.f };
    static float b[5] = { 840.f,630.f,504.f,420.f,360.f };
    static float b1[5] = { 5250.f,3234.f,2394.f,1914.f,1599.f };

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

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

    for (i__ = 1; i__ <= 5; ++i__) {
        for (j = 1; j <= 4; ++j) {
/* l92: */
            a1_ref(i__, j) = a_ref(i__, j);
        }
    }
    for (k = 1; k <= 2; ++k) {
      for (i__ = 1; i__ <= 5; ++i__) {
           printf("\n  %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));
      }
      for (i__ = 1; i__ <= 5; ++i__) {
           printf("\n  %12.4e \n", b[i__-1]);
      }
      p = 1;
      asp0r_c(a, b, x, t, s, &c__5, &c__4, &p);

      for (i__ = 1; i__ <= 5; ++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));
      }
      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__ <= 5; ++i__) {
           printf("\n  %16.7e \n", b[i__-1]);
      }
      for (i__ = 1; i__ <= 4; ++i__) {
           printf("\n  %16.7e \n", x[i__-1]);
      }
        for (i__ = 1; i__ <= 5; ++i__) {
            b[i__ - 1] = b1[i__ - 1];
            for (j = 1; j <= 4; ++j) {
/* l91: */
                a_ref(i__, j) = a1_ref(i__, j);
            }
        }
/* l9: */
    }
    return 0;
} /* main */


Результат:

      x  =   (1.00000000005, 0.500000000016, 0.3333333334, 0.25000000003, 0.20000000001)