Текст подпрограммы и версий
asg8r_c.zip , asg8d_c.zip , asg8c_c.zip 		Тексты тестовых примеров
tasg8r_c.zip , tasg8d_c.zip , tasg8c_c.zip

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

Назначение

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

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

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

     RCOND = 1/(|| A ||1*|| A-1||1) ,
 где
     || A ||1 = maxj=1,...,N   ( | a1j | + | a2j | + ... + | aNj | ) .
 Затем решается система
       L*y = b
 (для системы AT*x = b решается  UT*y = b) ,
 и затем решение  x находится как решение системы 
       U*x = y
 (для AT*x = b решается LT*x= y).

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

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

    int asg8r_c (real *a, integer *m, integer *n, integer *nlead,
            real *b, integer *ltr, integer *l, real *rcond, real *x,
            integer *ierr)

   Параметры
a - вещественный двумерный массив размера m на n, в котором задается исходная матрица; на выходе из подпрограммы на соответствующих местах находятся элементы матрицы U и поддиагональные элементы матриц исключения метода Гаусса  Li , i = 1, ..., n - 1;
m - первая размерность массива a в вызывающей программе (тип: целый);
n - порядок матрицы системы (тип: целый);
nlead - целый вектор длины n, содержащий на выходе информацию о выполненных в ходе факторизации перестановках (см. замечания по использованию);
b - вещественный вектор длины n, в котором задается правая часть системы;
ltr - задает вид решаемой системы (тип: целый):
ltr = 0 - если решается система A*x = b,
ltr ≠ 0 - если решается система AT*x = b;
l - признак решаемой системы (тип: целый);
l = 0 - если ситема с данной матрицей решается впервые,
l ≠ 0 - если система с заданной матрицей решается повторно;
rcond - вещественная переменная, содержащая на выходе вычисленное значение величины, обратной числу обусловленности матрицы a (см. замечания по использованию);
x - вещественный вектор длины n, содержащий на выходе решение системы (см. замечания по использованию);
ierr - целая переменная, содержащая на выходе информацию о прохождении счета; при этом:
ierr=65 - если m ≤ 0 или n ≤ 0 ;
ierr=66 - если в процессе счета возникло переполнение (это говорит о том, что либо ||A||1, либо некоторые элементы матрицы U, либо некоторые компоненты решения системы превосходят по абсолютной величине максимальное представимое на данной машине число);
ierr=-k - если в результате факторизации в k - ой строке матрицы U диагональный элемент равен нулю (это свидетельствует о вырожденности матрицы A). Если таких строк у матрицы U несколько, то значение k полагается равным номеру последней из них;
ierr=67 - если система несовместна (см. замечания по использованию). 
   Версии
asg8d_c - решение системы линейных алгебраических уравнений A*x = b методом Гаусса с выбором ведущего элемента по столбцу и оценка числа обусловленности матрицы A для вещественных A и b, заданных с удвоенной точностью.
asg8c_c - решение системы линейных алгебраических уравнений A*x = b методом Гаусса с выбором ведущего элемента по столбцу и оценка числа обусловленности матрицы A для комплексных A и b. 
   Вызываемые подпрограммы
afg4r_c - подпрограмма треугольной факторизации и оценки числа обусловленности матрицы A.
utafsi_c - подпрограмма выдачи диагностических сообщений. 
   Замечания по использованию
  1. 

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

  2. 

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

  3. 

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

  4. 

Если задано L ≠ 0, т.е. система решается повторно, то перед выполнением подпрограммы нужно запомнить, если требуется, вычисленную ранее оценку числа обусловленности, а в качестве матрицы A при обращении к подпрограмме нужно взять результат предыдущего обращения к asg8r_c или к afg4r_c, т.е. матрица A должна быть уже факторизована, т.к. при L ≠ 0 треугольная факторизация не выполняется, а переменной rcond присваивается значение ноль.

  5.  Если вырабатывается значение переменной ierr, отличное от нуля, то выдается соответствующее диагностическое сообщение, и, если ierr ≥ 65, то происходит выход из подпрограммы. Если система совместна, но матрица A вырождена (ierr < 0), т.е. для некоторых номеров k  U (k, k) = 0, то x (k) полагается равным 1.

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

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 };
    static float b[4] = { 7.f,7.f,7.f,7.f };

    /* Local variables */
    static int ierr;
    extern int asg8r_c(float *, int *, int *, int *, float *, int *,
                       int *, float *, float *, int *);
    static int i__, l, m, n, nlead[4];
    static float x[4], rcond;
    static int ltr;

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

    m = 4;
    n = m;
    ltr = 0;
    l = 0;
    for (i__ = 1; i__ <= 4; ++i__) {
         printf("\n %15.7e %15.7e %15.7e %15.7e \n",
                a_ref(i__, 1), a_ref(i__, 2), a_ref(i__, 3), a_ref(i__, 4));
    }
    asg8r_c(a, &m, &n, nlead, b, <r, &l, &rcond, x, &ierr);

    printf("\n  %5i %5i %5i %5i \n",
           nlead[0], nlead[1], nlead[2], nlead[3]);
    for (i__ = 1; i__ <= 4; ++i__) {
         printf("\n %15.7e %15.7e %15.7e %15.7e \n",
                a_ref(i__, 1), a_ref(i__, 2), a_ref(i__, 3), a_ref(i__, 4));
    }
    printf("\n  %5i \n", ierr);
    printf("\n  %15.7e  \n", rcond);
    for (i__ = 1; i__ <= 4; ++i__) {
         printf("\n  %15.7e  \n", b[i__-1]);
    }
    for (i__ = 1; i__ <= 4; ++i__) {
         printf("\n  %15.7e  \n", x[i__-1]);
    }
    return 0;
} /* main */


Результат:

               |   8.5           -4.8            0.8               3.5        |
               | -0.92941    10.06118     4.95647     -10.45294 |
      a  =  | -0.50588    -0.65879    -9.40171       2.40526  |
               | -0.37647    -0.04046    -0.73072      12.65817 |

      nlead  =   (2, 2, 4, 4)

      rcond  =  0.41764

               |  1.023054 |
               |  0.273777 |
      x  =  | -0.462957 |
               | -0.003275 |