Текст подпрограммы и версий ( Фортран )
mnp1r.zip 		Тексты тестовых примеров ( Фортран )
tmnp1r.zip
Текст подпрограммы и версий ( Си )
mnp1r_c.zip 		Тексты тестовых примеров ( Си )
tmnp1r_c.zip
Текст подпрограммы и версий ( Паскаль )
mnp1r_p.zip 		Тексты тестовых примеров ( Паскаль )
tmnp1r_p.zip

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

Назначение

Решение задачи квадратичного программирования при наличии линейных ограничений на переменные методом покоординатного спуска.

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

Решается задача квадратичного программирования 

     min  { (p, x) + (Cx, x)  |  Ax ≤ b } ,

где  p и  x - векторы длины  n,  C - симметричная положительно определенная матрица размерности  n * n,  A - матрица размерности  m * n,  b - вектоp длины  m.

Матрица  C задается в компактной форме, т.е. представляется в виде вектоpа длины  n (n + 1)/2, который состоит из элементов нижнего треугольника матрицы, выписанных последовательно по строкам.

Для получения решения исходной задачи строится двойственная ей задача: найти 

     min  { φ(u) = (h, u) + (Gu, u)  |  u ≥ 0 } ,

где  h = (1/2)AC - 1 p + b,  G = (1/4)AC - 1 A*.

Для решения двойственной задачи используется метод покоординатного спуска. Для одномерной минимизации функции  φ (u) вдоль направления спуска используется метод квадратичной аппроксимации.

Некоторая вычисленная точка  uk ≥ 0 считается точкой минимума функции  φ (u), если выполнено хотя бы одно из следующих условий:

1.  | uik - uik - 1 | ≤ XE для всех  i =1, ..., m,  где  uk = (u1k, ..., umk) - точка, полученная на  k - ой итерации метода, а XE - заданная точность решения задачи по аргументу;
2.  | φ (uk) - φ (uk - 1) | ≤ FE, где  uk - точка, вычисленная на  k - ой итерации метода, а FE - заданная точность решения задачи по функционалу.

Решение исходной задачи вычисляется по формуле 

     x(u) = - (1/2) C -1 (A*u + p) .

Кюнци Г.П., Крелле B., Нелинейное программирование, Изд - во "Советское радио", M., 1965.

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

    SUBROUTINE  MNP1R (M, N, A, B, U, P, FSTEP, IPAR, MAXK, XE,
                                             FE, KOUNT, CO, G, I0, IERR) 
   Параметры
M - число стpок матрицы  A (тип: целый);
N - число столбцов матрицы  A (тип: целый);
A - вещественный двумерный массив размерности  M * N;
B - вещественный вектоp длины  M;
U - вещественный вектоp длины  max {6M + N + 11; N (N + 1)/2 }; на входе первые N (N + 1)/2 компонент вектоpа содержат компактную запись матрицы  C;
P - вещественный вектоp длины  N; на входе содержит компоненты вектоpа  p, а на выходе содержит компоненты решения исходной задачи;
FSTEP - начальная длина шага (тип: вещественный);
IPAR - параметр, управляющий вариантом покоординатного спуска: если IPAR = 1, то используется вариант с ускоренным дроблением шага (тип: целый);
MAXK - целая переменная, на входе равная максимально допустимому числу вычислений функции, а на выходе - числу фактически выполненных вычислений функций;
XE - заданная точность решения задачи по аргументу (тип: вещественный);
FE - заданная точность решения задачи по функционалу (тип: вещественный);
KOUNT - целая переменная, содержащая число фактически выполненных итераций метода;
CO - вещественный вектоp длины N (N + 1)/2, используемый как рабочий;
G - вещественный вектоp длины M (M + 1)/2, используемый как рабочий;
I0 - целочисленный вектоp длины  M, используемый как рабочий;
IERR - целая переменная, указывающая причину окончания процесса:
IERR= 0 - если найден минимум функции с заданной точностью по аргументу или по функционалу;
IERR= 4 - если выполнено заданное максимальное число вычислений функции, но точность не была достигнута;
IERR=65 - если матрица  C не является положительно определенной. 
   Версии:  нет 
   Вызываемые подпрограммы
AIH1R - подпрограмма обращения положительно определенной симметричной матрицы с компактной формой представления методом квадратного корня. 
   Замечания по использованию
 

Матрица  C должна быть положительно определенной.

Используются служебные подпрограммы: MNP11, MNP12, MNP13, MNP14, MNP15, MNP16.

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

   Найти
         min  [ Q(x) = 0.5x12 + 0.5x22 - x1 - 2x2 ]

   при ограничениях

         2x1 + 3x2 ≤ 6 , 
           x1 + 4x2 ≤ 5 , 
           x1          ≥ 0 , 
                    x2 ≥ 0 .

   Точка минимума    x*  =  (13/17, 18/17) .

       INTEGER  KOUNT, M, MAXK, N, I0, IERR, IPAR
       REAL  A, B, CO, FE, FSTEP, G, P, U, XE
       DIMENSION  A(4, 2), B(4), I0(4), P(2), U(37), CO(3), G(10)
       DATA  XE, FE, M, N /1.E - 8, 1.E - 9, 4, 2/
       DATA  FSTEP, IPAR, MAXK / 1., 1, 250/
       DATA  A /2., 1., - 1., 0., 3., 4., 0., - 1./
       DATA  B /6., 5., 0., 0./
       DATA  P / -1., - 2./
       DATA U /0.5, 0., 0.5, 34*0./
       CALL  MNP1R (M, N, A, B, U, P, FSTEP, IPAR, MAXK, XE, FE,
      *                          KOUNT, CO, G, I0, IERR)

Результаты:

      IERR  =  0

      KOUNT  =  95
      MAXK  =  117

      P(1)  =  7.646729 - 01
      P(2)  =  1.058691 + 00