ПОЯСНИТЕЛЬНАЯ ЗАПИСКА

к курсовому проекту по “Системному анализу”
Тема: “Решение задач линейного программирования большой размерности”

Выполнил студент гр. Э-282:
	Богдановский А. А.
Проверил преподаватель:
Тихненко Е. В.

________________________

Дата:
“___” __________ 1996 г.

СОДЕРЖАНИЕ

1. Общая постановка задачи	
2. Конкретизация задачи	
3. Математическая модель симплекс-метода с мультипликативным представлением обратной матрицы	
3.1. Модифицированный симплекс-метод	
3.2. Мультипликативная форма обратной матрицы	
3.3. Преимущества метода	
4. Алгоритм метода, реализованного в программе	
5. Текст программы SASIMPL	
6. Интерфейс пользователя	
6.1. Работа с программой	
6.2. Формат файлов, содержащих постановку задач линейного программирования	
7. Результаты работы программы	
СПИСОК ЛИТЕРАТУРЫ	
Приложение I. Текст программы SASIMPL	

1.	
Общая постановка задачи
	Реализовать на произвольной вычислительной технике с помощью любого программного средства один из методов решения задач линейного программирования большой размерности.
2.	Конкретизация задачи
	В соответствии с общей постановкой задачи, возможностями студента, доступной литературой и другими факторами, студентом была конкретизирована и сформулирована следующая задача:

•	для решения задач линейного программирования большой размерности принять “модифицированный симплекс-метод с мультипликативным представлением обратной матрицы”;
•	разработать алгоритм для этого метода;
•	описать этот алгоритм на языке программирования Borland C++ 3.1;
•	написать программу-оболочку на языке Borland C++ 3.1, реализующую примитивный пользовательский интерфейс.
3.	Математическая модель симплекс-метода с мультипликативным представлением обратной матрицы
	Решаемая задача линейного программирования представлена в канонической форме и имеет следующий вид:
min F = cx, при Ax = b, x  0 ,
где	c	-	вектор коэффициентов целевой функции (c1,c2,...,cn);
	A	-	матрица ограничений размера mn ранга m, может быть представлена также как вектора [P1, P2, ..., Pn];
	b	-	m-вектор правой части ограничений (b1,b2,...,bm).	

	Таким образом, поставленная задача имеет n - переменных и m - ограничений.
	До рассмотрения мультипликативной формы кратко опишем этапы модифицированного симплекс-метода с хранением обратной матрицы в явной форме.
3.1.	Модифицированный симплекс-метод
	В  литературе  этот  метод встречается также под названием метода обратной матрицы.
	При  решении задач линейного программирования, в которых n (количество  переменных)  существенно  больше  m (количество ограничений),  модифицированный  симплекс-метод  требует  по сравнению   с   другими   значительно   меньшего  количества вычислительных операций и объема памяти ЭВМ. 
	В   модифицированном  симплекс-методе  реализуется  та  же основная  идея, что и в обычном симплекс-методе, но здесь на каждой  итерации  пересчитывается  не  вся матрица A-1, обратная матрице  ограничений A, а лишь та часть, которая относится к текущему базису Ax.
	Рассмотрим поэтапно шаги решения задачи линейного программирования модифицированным симплекс-методом:

1.	В начале первого цикла нам известны обратная матрица  (единичная матрица), базисное решение xb = b.
2.	Образуем для каждой небазисной переменной характеристическую разность j, используя уравнение:
j = cj — = cj — Pj ,
где		-	двойственные переменные, которые можно найти следующим образом:
 = cx * ,
где	cx	-	вектор коэффициентов целевой функции при базисных переменных.
3.	Предполагая, что используется стандартное правило выбора вводимого столбца, находим:
 .
4.	Если s  0 - процедура останавливается. Текущее базисное решение является оптимальным.
5.	Если s  0, вычисляем преобразованный столбец:
 = Ps
6.	Пусть
= (, , ..., ) .
	Если все  0 - процедура останавливается: оптимум неограничен.
7.	В противном случае находим выводимую из базиса переменную:
 =  .
8.	Строим увеличенную матрицу:

	и трансформируем ее с ведущим элементом . Первые m столбцов результат дают матрицу, обратную новому базису. Преобразуем базисное решение:
xb i  xb i —  * , i  r,
xb r  
	и переходим к этапу 2.
3.2.	Мультипликативная форма обратной матрицы
	Назовем элементарной матрицей матрицу, отличающуюся от единичной только одной строкой или столбцом. При мультипликативном представлении матрица  не дается явно, а представляется в виде произведения накапливаемых элементарных матриц. Для того, чтобы показать, как это делается, примем, что  - матрица, обратная текущей базисной, и предположим, что новая обратная вычисляется с помощью ведущей операции, опирающейся на ведущий элемент . Тогда над   следует произвести следующие операции:
1.	Для i=1,...,m, ir заменить строку i на
(строка i) —  * (строка r) .
2.	Заменить строку r на  * (строка r).
	Легко доказать непосредственным перемножением матриц, что эти операции соответствуют умножению  слева на элементарную матрицу:
E =  ,	(1)
где
	(2)
т. е.
= E *  ,
где		-	новая обратная матрица.
	Если начальный базис был единичным и если осуществлены k ведущих операций, то на цикле k обратная матрица  дается формулой:
= Ek * Ek-1 * ... * E1 ,
где каждая из матриц Ei является элементарной матрицей типа (1).
	Полученное представление  в виде произведения элементарных матриц называется мультипликативной формой обратной матрицы.
	Важным свойством элементарных матриц является то, что они могут сохраняться в памяти путем записи только элементов неединичного столбца и его положения в матрице (практически необходимо запоминать только ненулевые элементы и их положение). Поскольку обратная матрица теперь в явной форме недоступна, все вычисления, требующие ее использования, сводятся к последовательности умножений на элементарные матрицы. Рассмотрим эти вычисления в том порядке, в котором они производятся в симплекс-методе:
Оценка двойственных переменных:
	Двойственные переменные даются формулой:
 = cx * ,	(3)
так что, если  дана в мультипликативной форме, то приходится оценить следующее матричное произведение:
 = (...((cx * Ek) * Ek-1) * ...) * E1 .
Операции производятся в том же порядке, как указано скобками, то есть сначала вычисляется строка cx * Ek . Затем (cx * Ek) * Ek-1 и т. д. Каждая из этих операций требует умножения элементарной матрицы слева на вектор-строку. Пусть
E = [u1, u2, ..., ur-1, , ur+1, ..., um ],	(4)
где ui есть i-ый единичный вектор, и пусть
 = (1, ..., m)	(5)
- произвольный вектор-строка. Тогда
E = (1, ..., r-1, , r+1, ..., m),
где
 =  *  =  ,
то есть вектор E есть строка. отличающаяся от  только одним элементом , который равен скалярному произведению  на неединичный столбец . Таким образом,  вычисление  при задании обратной матрицы в виде произведения k элементарных матриц требует вычисления k скалярных произведений.
Вычисление ведущего столбца :
	Вектор дается формулой
= Ek * (... * (E2 * (E1 * Ps))...) .	(6)
Здесь произведения вычисляются в прямой последовательности: сначала E1 * Ps , затем E2 * (E1 * Ps) и т. д. Каждая операция требует вычисления произведения вида E. Если E и  таковы, как в (4), (5), но  теперь вектор-столбец, то
E = (1, ..., m),
где

Таким образом,  оценка E требует m умножений и m сложений.
Изменение обратной матрицы:
	Здесь требуется только добавление к списку, сохраняемому в памяти, нового вектора  вида (1), (2), элементы которого вычисляются по текущему .
3.3.	Преимущества метода
	Модифицированный симплекс-метод, в особенности в мультипликативной форме, обладает значительными преимуществами по сравнению со стандартной формой. Это относится к точности, скорости и требованиям к памяти. Большая часть этих преимуществ определяется тем фактором, что, как правило, матрицы больших линейных задач (то есть с n>m>100) являются слабозаполненными, содержат малый процент ненулевых элементов. Обычной является плотность 5% или менее. Модифицированная форма симплекс-метода в большей степени способна использовать преимущества, вытекающие из этого факта. В этой форме характеристические разности и ведущий вектор вычисляются непосредственно по исходным данным. Поскольку исходная матрица слабозаполнена, а перемножение следует производить только тогда, когда оба сомножителя отличны от нуля, то время вычислений значительно сокращается. В дополнение к этому использование только исходных данных приводит к тому, что уменьшается возможность накопления ошибок округления. Наоборот, стандартные симплексные таблицы, даже если они первоначально являются слабозаполненными, в ходе итеративного процесса быстро заполняются ненулевыми элементами. Таким образом,  время вычислений увеличивается, и, поскольку каждая таблица вычисляется из предшествующей, накопление ошибок может начать играть более серьезную роль.
4.	Алгоритм метода, реализованного в программе
	Исходные данные: 
Количество переменных:	n
Количество ограничений:	m
Вектор коэффициентов целевой функции:	c 
Матрица ограничений:	A
Правая часть ограничений:	b
	Предполагается, что задача линейного программирования задана в канонической форме для задачи минимизации.
Алгоритм:
1.	Нахождение начального опорного плана:
a)	перебираются все строки ограничений и ищется переменная с коэффициентом 1;
b)	проверяется, чтобы в остальных ограничениях у этой переменной был коэффициент 0; 
c)	если условие выполняется, то эта переменная заносится в базис;
d)	так ищем переменные до тех пор, пока не найдем m базисных переменных.
В результате получаем:
•	массив NXBasis[m], в котором хранятся индексы базисных переменных;
•	массив X[m], с самими значениями базисных переменных, которые на начальном этапе должны равняться правой части ограничений b;
2.	Счетчик итераций iteration приравниваем 1: iteration := 1.
3.	Расчет двойственных переменных Pi () :
a)	создаем массив коэффициентов целевой функции при базисных переменных Cx (cx) и записываем его в массив двойственных переменных Pi: Pi := Cx;
b)	организуем последовательное умножение Pi на E[i] (Ek) в соответствии с (3), где i изменяется от iteration до 1 (если iteration = 1, то данный пункт пропускается).
4.	Поиск ведущего столбца (вводимой переменной):
a)	расчет характеристических оценок всех небазисных переменных и поиск минимальной среди них c_min (c_min  s, s_min  s ) ;
b)	если c_min  0, то алгоритм заканчивает работу и текущие базисные значения переменных являются оптимальным решением задачи.
c)	иначе, определяем, что в базис будем вводить переменную с индексом s_min.
5.	Вычисление ведущего столбца AA ():
a)	Заносим в массив AA[m] столбец s_min из матрицы A;
b)	организуем последовательное умножение E[i] (Ek) на AA в соответствии с (6), где i изменяется от 1 до iteration 1 (если iteration = 1, то данный пункт пропускается).
c)	проверяем, если все элементы AA положительны или равны нулю, то текущее решение является оптимальным, но не единственным, и алгоритм заканчивает работу.
6.	Поиск выводимой переменной из базиса:
a)	ищем минимальное значение th = X[i]/AA[i], где i от 1 до m
b)	найденное th_min (  th_min, r  r_min) будет соответствовать переменной с индексом в базисе r, выводимой из базиса;
7.	Добавляем в список элементарных матриц еще одну:
a)	вычисляем новую элементарную матрицу Eiteration в соответствии с (1), (2) и добавляем ее в список элементарных матриц;
8.	Делаем замену переменных в базисе:
a)	присваиваем NXBasis[r_min] = s_min
b)	перевычисляем значения базисных переменных:
X[i] = X[i] - th_min * AA[i], i  r_min
X[r_min] = th_min .
9.	Увеличиваем счетчик итераций на единицу: iteration := iteration + 1.
10.	Переходим к пункту 3.
5.	Текст программы SASIMPL
	Программа была написана на языке программирования Borland C++ 3.1 с использованием библиотеки работы с матрицами MATRIX и библиотеки интерфейса пользователя TSWM.
	Программа SASIMPL состоит из 6 модулей:
1.	SA.CPP	главный модуль программы (интерфейс)
2.	MSIMPLEX.CPP	реализация метода
3.	MATRIX.CPP	библиотека работы с матрицами
4.	GERROR.CPP	обработка ошибок
5.	TSWM.LIB	библиотека интерфейса пользователя
6.	GETFILE.LIB	то же

В приложении I представлены тексты модулей MSIMPLEX.CPP, SA.CPP, MATRIX.CPP и GERROR.CPP.
Интерфейс пользователя
	Программа SASIMPL имеет примитивный пользовательский интерфейс, который позволяет загружать данные о задачах из внешних файлов, решать их и просматривать результаты. Кроме того, программа предоставляет скромную справочную систему для удобства работы.
6.1.	Работа с программой
Главный экран программы, который появляется при ее запуске изображен на рисунке 1:

Рисунок 1
Здесь программа предоставляет пользователю выбрать необходимый режим работы:
1. "Решить задачу линейного программирования"
  в  этом  разделе  программы можно  задать данные по задаче   линейного  программирования  и  решить  ее; данные нужно  будет  загрузить  из  внешнего  файла; формат файлов, содержащих  постановку  задачи  линейного  программирования, можно узнать, если выбрать режим "Узнать о программе" главного   меню   данной  программы;  после  решения  задачи программа  выведет  на  экран  результаты  вычисления, кроме того,   эти  результаты  будут сохранены  в  текущей  директории в файле SASIMPL.RES.
2. "Почитать краткую теорию метода"
  в  этом  разделе  описывается  краткая  теория модифицированного      симплекс-метода с мультипликативным представлением  обратной  матрицы;
3. "Узнать о программе"
 здесь можно узнать о программе (назначение, автор, дата, ...)  и,  кроме того,  формат файлов с постановками задач линейного программирования.
4. "Выйти из программы"
 выход из данной программы.
	Если пользователь выбрал первый пункт меню, на экране появится изображение похожее на рисунок 2:

Здесь пользователю предлагается выбрать файл, содержащий данные о задаче, которую ему необходимо решить (формат таких файлов см. в 6.2 ).
	После того, как пользователь укажет программе необходимый файл, программа произведет вычисления по заданному алгоритму метода решения задач линейного программирования и на экране появятся результаты вычислений в виде типа рис. 3:

6.2.	Формат файлов, содержащих постановку задач линейного программирования
	Программа SASIMPL понимает только определенные файлы с постановкой задач. Рассмотрим один из них:
; *** Начало файла с данными о задаче ***
; Пример  из  Ю.  П. Зайченко, С. А. Шумилова
; "Исследование операций" N1.46
n = 8                        ; количество переменных
m = 4                        ; количество ограничений
F = -3 -1 0 1000 0 0 1000    ; целевая функция минимизации
LIMITS:                      ; задание ограничений
1 2 -1 1 = 5
2 4 0 0 1 = 16
3 1 0 0 0 -1 1 = 6
1 3 0 0 0 0 0 1 = 9
;*** Конец файла ***

	Символ ';' является символом-коментария.
	Создавая  файлы  такого  формата, обратите внимание на то, что  целевая  функция  F  должна  задаваться для минимизации (преобразование  максимизирующей  к  минимизирующей  функции производится  путем умножения ее коэффициентов на -1). Кроме того,  в ограничениях должны стоять только (!) равенства, то есть  должны  уже  быть  введены  свободные  и искусственные переменные.
	Если   в   целевой   функции  или  в  ограничениях  заданы коэффициенты  не  для  всех  переменных, то эти коэффициенты принимаются  за нулевые, но в целях проверки, в этом случае, программа выведет предупреждающее сообщение об этом.
     При  вводе искусственных переменных используйте большие коэффициенты  в  целевой  функции,  чтобы  они  были на два-три порядка больше, чем для обычных переменных. 
7.	Результаты работы программы
	В качестве контрольных примеров бралось множество примеров из задачника /2/.
	Рассмотрим следующий пример:
Fmax = 3x1 + x2
x1   + 2x2  5,
2x1 + 4x2  16,
3x1 + x2  6,
x1 + 3x2  9,
x1, x2  0.

Так как программе SASIMPL требуется постановка задачи в канонической форме, то приведем ее к канонической форме:
Fmax = 3x1 + x2 + 0x3 + Mx4 + 0x5 + 0x6 + Mx7 + 0x8
x1   + 2x2 - x3 + x4 = 5,
2x1 + 4x2 + x5 = 16,
3x1 + x2  - x6 + x7 = 6,
x1 + 3x2 + x8 = 9,
x1, x2  0.
Таким образом,  мы ввели две свободные переменные x3, x6, и две искусственные - x4 и x7 .
	Составим входной файл для программы, содержащий данную задачу:
n = 8                        ; количество переменных
m = 4                        ; количество ограничений
F = -3 -1 0 1000 0 0 1000    ; целевая функция минимизации
LIMITS:                      ; задание ограничений
1 2 -1 1 = 5
2 4 0 0 1 = 16
3 1 0 0 0 -1 1 = 6
1 3 0 0 0 0 0 1 = 9

	Как можно увидеть, функция Fmax преобразована к Fmin путем умножения на -1. В качестве больших коэффициентов M взято число 1000, которое во много раз превосходит соседние коэффициенты в целевой функции.
	В результате вычислений программы SASIMPL при подаче на вход указанного файла, получили следующие результаты:
Входные данные:
Количество переменных  n = 9
Количество ограничений m = 5
Вектор коэффициентов целевой функции C:
-1.000 -1.000 -2.000 -1.000  0.000  0.000  0.000  0.000  0.000

Матрица коэффициентов в ограничениях A:
 1.000  2.000  1.000 -1.000  1.000  0.000  0.000  0.000  0.000
 2.000  1.000  0.000  0.000  0.000  1.000  0.000  0.000  0.000
 1.000  2.000  0.000  0.000  0.000  0.000  1.000  0.000  0.000
 0.000  0.000  1.000  1.000  0.000  0.000  0.000  1.000  0.000
 0.000  0.000 -2.000  3.000  0.000  0.000  0.000  0.000  1.000

Вектор правой части ограничений B:
10.000  4.000  6.000  5.000  6.000

********************************************************
Итерация N0
Значения базисных переменных:
X5 = 10
X6 = 4
X7 = 6
X8 = 5
X9 = 6
Значение целевой функции:  0.000

Итерация N1
Значения базисных переменных:
X5 = 5
X6 = 4
X7 = 6
X3 = 5
X9 = 16
Значение целевой функции: -10.000

Итерация N2
Значения базисных переменных:
X5 = 3
X1 = 2
X7 = 4
X3 = 5
X9 = 16
Значение целевой функции: -12.000

Итерация N3
Значения базисных переменных:
X2 = 2
X1 = 1
X7 = 1
X3 = 5
X9 = 16
Значение целевой функции: -13.000

Данное решение является оптимальным!

1.	Оптимизация больших систем. Л. С. Лэсдон
2.	Исследование операций. Сборник задач. Ю. П. Зайченко, С. А. Шумилова. - Киев: Вища школа. Изд. при Киев. ун-те, 1984. - 224 с.
3.	Исследование операций: В 2-х томах. Пер. с англ./Под ред. Дж. Моудера, С. Элмаграби.-М.: Мир, 1981. Т. 1. 712 с., ил.
4.	Вычислительные методы линейной алгебры. Фадеев Д. К., Фадеева В. Н.
5.	Вычислительные основы линейной алгебры. В. В. Воеводин
6.	Исследование операций в задачах, алгоритмах и программах. Е. М. Кудрявцев
7.	Исследование операций. Дегтярев Ю. И. -М.: Высшая школа, 1986-320 с., ил.

Файл SA.CPP:
/**********************************************************************

                 Главный модуль к программе SASIMPL

   Реализован модифицированный симплекс-метод с мультипликативным
   представлением обратной матрицы

В данном модуле описывается интерфейсная часть программы (считывание
файла с данными)

В модуле MSIMPLEX.CPP реализован непосредственно сам метод

Эта программа является курсовым проектом по системному анализу

Язык написания: C++
Компилятор: Borland C++ 3.1
Автор: студент ДВГТУ группы Э-282 Богдановский А.
Дата: осенний семестр 1995-96 гг.
Преподаватель: Тихненко Е. В.

**********************************************************************/

// Нашенской программы описания
#include "sa.h"

#include <ctype.h>
#include <alloc.h>
#include <dir.h>
#include <process.h>
#include <stdio.h>
#include <fstream.h>

RUSSIAN;
KEYBOARD;

#define HELPFILE "SASIMPL.HLP"

char *Copyright = "Богдановский Александр, ДВГТУ, Э-282, 1996 г.";

char *err_torun = "(Прочитайте информацию о формате входного файла\n"
                  "в разделе главного меню \"Узнать о программе\")\n";

char *ViewFile = "README.COM";

// Входной файл с данными
char DATFile[101] = "";

// Эта функция находится в библиотеке GETFILE.LIB и нужна для
// выбора VEC-файла ( написана она на очень старой библиотеке
// WMT.LIB, т.к. в этом дурацком TSWM'е нет такой функции ).
// Первый параметр - переменная куда запишется файл (с путем)
// Второй параметр - маска (char [15]).  Возвращает  1 - если
// файл выбрали и 0 - otherwise. ПредПоследний параметр - это
// функция просмотра MAP- или VEC-файлов по клавише F3 (функ-
// ция определена в этом модуле - char View( char * );  ).  И
// последний параметр - это функция подсказки, которая встра-
// ивается в меню файлов по клавише F1
extern int SelFile( unsigned char*, unsigned char*,
                    unsigned char (*)(unsigned char*),
                    unsigned char (*)(unsigned char*)=NULL );

// Обработка входного файла
// считывание из него входных данных и заполнение нужных переменных
int LoadInputData( char *filename )
{
  ifstream ifile( filename );
  if( !ifile )
  {
    GError << "\nОшибка открытия входного файла: " << filename << endl
           << "(запустите программу без параметров для подсказки)" << endl;
    GlobalError();
    return 0;
  }

  static char buf[201];
  static char *s;
  static int *NorM;
  static char fl1, fl2;
  char fl_jnen = 0;

  fl1 = fl2 = 0;

  // считывание входного файла
  do
  {
    ifile.getline( buf, 200 );   // ... по строчке
    strupr( buf );  // поднимаем все символы строки
    s = strchr( buf, ';' ); if( s != NULL ) *s = 0;   // убираем комментарии

// считываем директиву LIMITS ***

    s = strstr( buf, "LIM" );
    if( s != NULL )
    {
      if( n == 0 || m == 0 )
      {
        GError << "\nЗначения n и m должны определяться до определения\n"
                  "ограничений !\n" << err_torun;
        GlobalError();
        return 0;
      }

      A.SetSize( m, n ); A = 0.;
      if( B != NULL ) delete B; B = new double [m];
      if( B == NULL )
      {
        GError << "\nНе хватает памяти под вектор B в правой части\n"
                  "ограничений!\n";
        GlobalError();
        return 0;
      }
      for( int i = 0; i < m && !ifile.eof(); i++ )
      {
        ifile.getline( buf, 200 );
        s = strchr( buf, ';' ); if( s != NULL ) *s = 0;

        s = buf;
        for( int j = 0; j < n+1 && *s; j++ )
        {
          while( !isdigit( *s ) && *s != '-' && *s != '.'
                 && *s != '=' && *s != 'E' && *s ) s++;
          if( *s )
          {
            if( *s == '=' )
            {
              if( j != n )
              {
                if( fl_jnen == 0 )
                {
                  GError << "\nВнимание: в ограничении заданы не все коэффициенты!\n"
                            "          они принимаются нулевыми...\n";
                  GlobalError();
                  fl_jnen = 1;
                }
              }
              while( !isdigit( *s ) && *s != '-' && *s != '.'&& *s ) s++;
              if( *s )
                B[i] = atof( s );
              else
              {
                GError << "\nОшибка задания правой части ограничений\n"
                          "в строке " << (i+1) << endl << err_torun;
                GlobalError(); return 0;
              }
              break;
            }
            if( j == n )
            {
              GError << "\nКоличество переменных в ограничении на строке"
                     << (i+1) << "\nбольше, чем n = " << n << endl
                     << err_torun;
              GlobalError(); return 0;
            }
            A[i][j] = atof( s );
          }
          else
          {
            if( j == 0 )
            {
              GError << "\nОшибка в ограничениях LIMITS: строка " << (i+1)
                     << endl << err_torun;
              GlobalError(); return 0;
            }
          }
          while( isdigit( *s ) || *s == '-' || *s == '.' || *s == 'E' ) s++;
        }
      }
      if( i != m )
      {
        GError << "\nЗадано мало ограничений или большое число m !\n"
               << err_torun;
        GlobalError(); return 0;
      }

      fl2 = 1;
      continue;
    }

// *** считываем директивы n и m ***

    s = strchr( buf, 'N' );
    if( s != NULL )
    {
      NorM = &n; // n - глобальная переменная
    }
    else
    {
      s = strchr( buf, 'M' );
      if( s != NULL )
      {
        NorM = &m; // m - глобальная переменная
      }
    }

    if( s != NULL )
    {
      while( !isdigit(*s) && *s ) s++;
      if( *s )
      {
        *NorM = atoi( s );
        if( *NorM <= 0 )
        {
          GError << "\nОшибка во входном файле: " << endl
                 << " неправильное значение n или m ...\n" << err_torun;
          GlobalError(); return 0;
        }
        continue;
      }
      else
      {
        GError << "\nДиректива n или m без продолжения !\n" << err_torun;
        GlobalError(); return 0;
      }
    }

//*** считываем директиву F ***

    s = strchr( buf, 'F' );
    if( s != NULL )
    {
      if( n == 0 )
      {
        GError << "\nКоличество переменных n должно определяться\n"
               << "до определения целевой функции во входном файле!\n"
               << err_torun;
        GlobalError(); return 0;
      }
      static char *err_f = "\nОшибка в определении вектора C"
                           "коэффициентов при x\nв целевой функции F...\n";
      s = strchr( buf, '=' ); if( s == NULL ) s = strchr( buf, ':' );
      if( s != NULL )
      {
        if( C != NULL ) delete C;
        C = new double [ n ]; memset( (void*)C, 0, sizeof(double) * n );
        if( C == NULL )
        {
          GError << "\nНе хватает памяти под вектор C коэффициентов при x" << endl
                 << "в целевой функции F...\n";
          GlobalError(); return 0;
        }
        for( int i = 0; i < n && *s; i++ )
        {
          while( !isdigit(*s) && *s != '-' && *s != '.' && *s != 'E'
                 && *s ) s++;
          if( *s )
            C[i] = atof( s );
          else
          {
            if( i == 0 )
            {
              GError << err_f << err_torun;
              GlobalError(); return 0;
            }
            break;
          }
          while( isdigit( *s ) || *s == '-' || *s == '.' || *s == 'E' ) s++;
        }
        if( i != n )
        {
          GError << "\nВнимание: в целевой функции были заданы не все\n"
                    "          коэффициенты для x !\n";
          GlobalError();
        }
        fl1 = 1;
        continue;
      }
      else
      {
        GError << err_f << err_torun;
        GlobalError(); return 0;
      }
    }

  } while( !ifile.eof() );

  if( !fl1 )
  {
    GError << "\nВо входном файле не задана целевая функция!\n" << err_torun;
    GlobalError(); return 0;
  }
  if( !fl2 )
  {
    GError << "\nВо входном файле не заданы ограничения!\n" << err_torun;
    GlobalError(); return 0;
  }
  return 1;
}

// Уничтожение переменных, под которые была запрошена память
void UnLoadData( void )
{
  if( C != NULL ){ delete C; C = NULL; }
  if( B != NULL ){ delete B; B = NULL; }
}

// Проверить, что нам дали в параметры
int TestArg( int argc, char **argv )
{
  if( argc == 2 )
  {
    ifstream ifile( argv[1] );
    if( !ifile ); else{ strcpy( DATFile, argv[1] ); return 1; }
  }
  return 0;
}

// Просмотрщик MAP- и VEC- файлов,  встраиваемый в меню файлов
// SelFile, по клавише F3
unsigned char View( unsigned char *file )
{

  strupr( file );

  // запомнить экран
  StoreScreen();

  int i = spawnlp( P_WAIT, ViewFile, ViewFile, file, NULL );

  // если с первого раза не запустилась,  то пробуем запустить
  // без путя
  if( i == -1 )
  {
    char *s = strrchr( ViewFile, '\\' );
    if( s != NULL )
      i = spawnlp( P_WAIT, s+1, s+1, file, NULL );
  }

  // восстановить старый текстовый режим
  textmode( T.currmode );
  SetBackGround( 16 ); // немигающие фоновые цвета

  // восстановить экран
  RestoreScreen();

  // зачем включить курсор, а потом выключить ???
  cursorON;
  cursorOFF;
  // я сам не знаю,  но если убрать  cursorON,  то  курсор  не
  // OFFнется

  if( i == -1 )
  {
    StartMonitor();
    Warning( NULL, "Извините, но программа (%s), с помощью кот"
                   "орой\nможно посмотреть файл %s,\nпочему-то"
                   " не запускается.\nЖелательно, чтобы эта пр"
                   "ограмма находилась\nв той же дирректории, "
                   "из которой Вы запустили MAKEMAP.EXE",
                   ViewFile, file );
    FinishMonitor();
  }

  return 15;

}

char *all = "Все, кому нужны исходники этой программы, обращайтесь"
            "к автору этой программы.  Он запросто подарит вам их!";

void TheoryReading( void )
{
  TSWIN w( ATTR( YELLOW, GREEN ) );
  w.Locate( T.screenwidth/2, T.screenheight/2, 62, 18, CENTER );
  w.Shadow( RIGHT_BOTTOM );
  w.AttrBor( ATTR( WHITE, GREEN ) );
  w.Open();
  w.PutTitle( " Краткое описание метода " );
  w.PutTitle(CENTER_BOTTOM,0," Esc, , PgUp, PgDn, Home, End ");

  TextShow( HELPFILE, "Описание_Метода" );
}

void AboutReading( void )
{
  TSWIN w( ATTR( YELLOW, GREEN ) );
  w.Locate( T.screenwidth/2, T.screenheight/2, 62, 18, CENTER );
  w.Shadow( RIGHT_BOTTOM );
  w.AttrBor( ATTR( WHITE, GREEN ) );
  w.Open();
  w.PutTitle( " О Программе " );
  w.PutTitle(CENTER_BOTTOM,0," Esc, , PgUp, PgDn, Home, End ");

  TextShow( HELPFILE, "About" );
}

// Подсказка встроенная в меню файлов SelFile, которое  не  из
// этой библиотеки
unsigned char SelFileHelp( unsigned char * )
{
  StartMonitor();
  HelpHotKey();
  FinishMonitor();
  return 15;
}

int GetDATFile( void )
{
  TSWIN w( ATTR( BLACK, YELLOW ) );
  w.Locate( 1, 6, 18, 6, LEFT_TOP );
  w.Border( "        " );
  w.Shadow( RIGHT_BOTTOM );
  w.Open();

  Ccputs( 1, "Выберите файл" );
  Ccputs( 2, "с данными задачи" );
  Ccputs( 3, "линейного" );
  Ccputs( 4, "программирования" );

  BL bl( "\x1 F1\x2 - подсказка\x1 F3\x2 - просмотреть файл\x1 Enter\x2 - выбрать\x1"
         "_Esc\x2 - выход " );

  // то берем имя из меню файлов
  unsigned char path[100] ;
  unsigned char mask[15] = "*.*";
  getcwd( path, 99 );  // взяли текущий путь

  // запускаем меню файлов
  FinishMonitor();     // т.к. SelFile из старой библиотеки
  SetHelp( "SelFile - взять файл с данными" );
  int i = SelFile( path, mask, View, SelFileHelp );
  SetHelp( "Nothing" );
  StartMonitor();

  if( i == 0 ) return 0;

  strcpy( DATFile, path );
  strupr( DATFile );
  return 1;
}

int GetDATA( void )
{
  if( DATFile[0] == 0 )
  {
    if( !GetDATFile() ) return 0;
  }
  if( !LoadInputData( DATFile ) ) return 0;
  DATFile[0] = 0;
  return 1;
}

void OutResults( void )
{
  TSWIN w( ATTR( WHITE, MAGENTA ) );
  w.Locate( T.screenwidth/2, T.screenheight/2, 64, 18, CENTER );
  w.Shadow( RIGHT_BOTTOM );
  w.AttrBor( ATTR( YELLOW, MAGENTA ) );
  w.Open();
  w.PutTitle( " Результаты вычислений " );
  w.PutTitle(CENTER_BOTTOM,0," Esc, , PgUp, PgDn, Home, End ");

  TextShow( OutFile, NULL );
}

void GeneralWork( void )
{
  // загрузка данных в память
  if( !GetDATA() ) return;

  TSWIN w( ATTR( WHITE, BROWN ) );
  w.Locate( T.screenwidth/2, T.screenheight/2, 40, 7, CENTER );
  w.Shadow( RIGHT_BOTTOM );
  w.Open();
  Ccputs( 1, "Тихо!" );
  Ccputs( 2, "Идут вычисления..." );
  Ccputs( 3, "Интерация номер:" );

  // решение ЛП нашим диким методом
  if( !MSimplex() ) return;
  w.Close();

  // вывод результатов
  OutResults();
  // очистка памяти от всякой ерунды
  UnLoadData();
}

// Начало программы
void main( int argc, char **argv )
{

  StartMonitor();

  SetBackGround( 16 ); // немигающие фоновые цвета
  cursorOFF;

  // Установка подсказки
  char HelpColor[3] = { ATTR(BLACK,CYAN), ATTR(YELLOW,CYAN), 0 };
  SetHelp( HELPFILE, "Главное меню", HelpColor );
  HotKey hk( 0, F0+1, HelpHotKey );

  {
    TSWIN wmain( ATTR( GREEN, BLACK ) );
    wmain.Locate( 1, 1, T.screenwidth, T.screenheight - 1 );
    wmain.FillChar( 0xB1 );
    wmain.Border( "\xB1\xB1\xB1\xB1\xB1\xB1\xB1\xB1" );
    wmain.Open();
    PutString( "Bog                       ПРОГРАММА SASIMPL                            Версия 1 ", 1, 1, 80, ATTR( BLACK, LIGHTGREEN ) );

    BL_Color[0] = ATTR( RED, WHITE );
    BL_Color[1] = ATTR( BLACK, WHITE );

    BL bl( "\x1_ _", NULL );

    // если аргумент не VEC-файл, то начало-тягомотина
    if( TestArg( argc, argv ) == 0 )
    {
      BL bl( "\x1_ ,\x2 - \"ходить\" по меню \x1 Enter\x2 - "
             "выбрать пункт меню \x1_ Esc\x2 - выход_", NULL );

      TSWIN w( ATTR( WHITE, BLUE ) );
      w.Locate( T.screenwidth/2, T.screenheight/2 + 1, 70, 19, CENTER );
      w.Open( CENTER, 25 );

      textcolor( YELLOW );

      Ccputs(  1, "РЕШЕНИЕ ЗАДАЧ, ЛИНЕЙНОГО ПРОГРАММИРОВАНИЯ" );
      Ccputs(  2, "БОЛЬШОЙ РАЗМЕРНОСТИ МОДИФИЦИРОВАННЫМ СИМПЛЕКС-МЕТОДОМ" );
      Ccputs(  3, "С МУЛЬТИПЛИКАТИВНЫМ ПРЕДСТАВЛЕНИЕМ ОБРАТНОЙ МАТРИЦЫ" );

      Ccputs( 5, "  Программа  работает  с  линейными  задачами минимизации,");
      Ccputs( 6, "поставленными в канонической форме, то есть в ограничениях");
      Ccputs( 7, "стоят   только   равенства  (введены  все  искусственные и");
      Ccputs( 8, "свободные перемнные).                                     ");

      textcolor( LIGHTMAGENTA );
      Ccputs( 10, "Выберите, что вы желаете:" );
      textcolor( WHITE );

      TSMENU mmain( 4 );
      mmain.Color( ATTR( BLACK, CYAN ) );
      mmain.CurrentShow( Attr );
      mmain.Item( 12, 12, "Решать задачу линейного программирования" );
      mmain.Item( 12, 14, "     Почитать краткую теорию метода     " );
      mmain.Item( 12, 15, "          Узнать о программе            " );
      mmain.Item( 12, 16, "          Выйти из программы            " );

      while( 1 )
      {
        int i = mmain.Exec();
        if( i )
          switch( mmain.Current() )
          {
            case 0: GeneralWork(); continue;
            case 1: TheoryReading(); continue;
            case 2: AboutReading(); continue;
            default: break;
          }
        break;
      }
    }
    else GeneralWork();
  }
  cursorON;
  FinishMonitor();
}

Файл MSIMPLEX.CPP:
/**********************************************************************

           Один из вспомогательных модулей программы SASIMPL

                   (Главный модуль программы SA.CPP)

 Модуль, содержащий непосредственно модифицированный сиплекс-метод
 с мультипликативным представлением обратной матрицы

 Предполагается, что входные данные уже заданы ...
 Вывод результатов производится в выходной файл ...

 Линейная задача задается в каноническом виде (это значит, что
 в ограничениях никаких неравенств и введены уже все свободные
 и исскуственные переменные) для задачи минимизации

 Соглашения: m x n означает матрицу высотой в m строк и шириной
             n столбцов ;-)

**********************************************************************/

// Нашенской программы описания
#include "sa.h"

#include <setjmp.h>

#define Empty ;

#define NoMem { GError << "\nНет памяти!!!"; GlobalError(); longjmp(Gjmpbuf,1); }

// Выходной файл с результатами работы
char *OutFile = "SASIMPL.RES";

// Тот же файл, но в виде указателя для системы
ofstream OFile;

//************************ ВХОДНЫЕ ДАННЫЕ *****************************
// ( должны задаваться до вызова функции MultiSimplex() )

// Самые главные значения: количество переменных
// и количество ограничений
int n, m;

// Вектор коэффициентов при x в ЦФ (n элементов)
double *C = NULL;

// Матрица коэффициентов ограничений (m x n)
Matrix A;

// Вектор правой части ограничений (m элементов)
double *B = NULL;

//************************ ПЕРЕМЕННЫЕ *********************************
int *NXBasis = NULL;     // индексы базисных переменных
double *X = NULL;        // значения базисных переменных
jmp_buf Gjmpbuf;         // всякая ерунда

//******************* КЛАСС - СПИСОК ЭЛЕМЕНТАРНЫХ МАТРИЦ **************
class EMatrix
{
public:
  static int N;  // размер элементарных матриц NxN
  int fl;        // как представлены данные  (0 - в массиве, n - в
                 // структуре, где n - ко-во элементов в ней)
  // представление столбца элементарной матрицы, который содержит данные
  union E_Data
  {
    double *E;     // либо в виде массива (если много не 0-вых элементов)
    struct E_Ef    // либо в виде структуры, хранящей только не 0-вые эл.
    { int n; double En; } *Eef;
  } Edata;

  int Column;      // номер столбца элементарной матрицы, хранящего данные
  EMatrix *Next;   // следующая элементарная матрица списка

  void Init( void );    // функция инициализации списка
  void destroy( void ); // функция уничтожения списка

public:
  EMatrix(){ Init(); }                // конструктор по умолчания (списка)
  EMatrix( int h ){ Init(); N = h; }  // конструктор одной элементарной мат.
  ~EMatrix(){ destroy(); }            // деструктор
  void Add( double *e, int column );  // добавить в список элемент. матрицу
  EMatrix &operator[]( int n );       // взять из списка элемент. матрицу
  double Get( int n );                // взять из матрицы конкрет. элемент
};

// *** ОПИСАНИЕ ФУНКЦИЙ КЛАССА СПИСКА ЭЛЕМЕНТАРНЫХ МАТРИЦ ***
void EMatrix::Init( void )
{
  Next = NULL;
  fl = 0;
  Edata.E = NULL;
}

void EMatrix::Add( double *e, int column )
{
  EMatrix *s = this;
  if( this->Edata.E != NULL || this->Edata.Eef != NULL )
  {
    while( s->Next != NULL ) s = s->Next;
    s->Next = new EMatrix; if(s->Next == NULL ) NoMem;
    s = s->Next;
  }
  // считаем ко-во ненулевых элементов
  int sk = 0;
  for( int i = 0; i < N; i++ ) if( e[i] ) sk++;

  if( N * sizeof(double) < sk * sizeof(E_Data::E_Ef) )
  {
    s->fl = sk;
    s->Edata.Eef = new E_Data::E_Ef [sk]; if( s->Edata.Eef == NULL ) NoMem;
    for( i = 0, sk = 0; i < N; i++ )
      if( e[i] ){ s->Edata.Eef[sk].n = i; s->Edata.Eef[sk].En = e[i]; sk++; }
  }
  else
  {
    s->fl = 0;
    s->Edata.E = new double [N]; if( s->Edata.E == NULL ) NoMem;
    for( i = 0; i < N; i++ ) s->Edata.E[i] = e[i];
  }
  s->Column = column;
}

void EMatrix::destroy( void )
{
  if( Next != NULL ) delete Next;
  if( fl )
    if( Edata.Eef != NULL ) delete Edata.Eef;
  else
    if( Edata.E != NULL ) delete Edata.E;
  Next = NULL; Edata.E = NULL; fl = 0;
}

EMatrix &EMatrix::operator[]( int n )
{
  EMatrix *s = this;
  for( int i = 0; i != n; i++, s = s -> Next ) Empty;
  return *s;
}

double EMatrix::Get( int n )
{
  if( fl )
  {
    for( int i = 0; i < fl; i++ )
    {
      if( Edata.Eef[i].n == n ) return Edata.Eef[i].En;
      if( Edata.Eef[i].n > n ) break;
    }
    return 0.;
  }
  return Edata.E[n];
}

int EMatrix::N = 0;

//**************** ФУНКЦИИ МОДУЛЯ ***********************************

// Поиск базиса в начальных данных
void FindBasis( void )
{
  for( int j = 0; j < n; j++ )
  {
    for( int i = 0; i < m; i++ )
    {
      // ищем переменную с коэф. 1 в ограничениях (индекс переменной j)
      if( A[i][j] == 1 )
      {
        // проверяем нет ли ее уже в базисе
        for( int k = 0; NXBasis[k] != -1 &&
                        NXBasis[k] != j && k < m; k++ ) Empty;
        if( NXBasis[k] != j )
        {
          // ...если нет, то проверяем коэффициенты в столбце
          // ограничений этой переменной
          for( int k = 0; k < m; k++ ) if( k != i ) if( A[k][j] ) break;
          if( k == m )
          {
            // если все Ob, то заносим эту переменную в базис
            for( int l = 0; NXBasis[l] != -1 && l < m; l++ ) Empty;
            NXBasis[l] = j;
            if( l == m - 1 ) return;  // все переменные нашли
            break; // переходим на другую строку ограничений
          }
        }
      }
    }
    if( j == n )
    {
      GError << "\nЗадача  линейного  программирования должна быть\n"
             << "задана в  каноническом  виде.  Это значит,  что\n"
             << "в ограничениях уже должен присутствовать базис!\n"
             << "А программа его что-то не находит...\n"
             << err_torun;
      GlobalError(); longjmp( Gjmpbuf, 1);
    }
  }
}

// Подготовительная работа
void BeforeTimeWork( void )
{
  if( NXBasis != NULL ) delete NXBasis;
  NXBasis = new int [ m ]; if( NXBasis == NULL ) NoMem;
  for( int i = 0; i < m; i++ ) NXBasis[i] = -1;
  FindBasis();    // ищем начальный базис
  // заполняем значения базисных переменных из правой части B
  if( X != NULL ) delete X; X = new double [ m ]; if( X == NULL ) NoMem;
  for( i = 0; i < m; i++ ) X[i] = B[i];
}

// Завершающая работа
void AfterTimeWork( void )
{
  if( NXBasis != NULL ){ delete NXBasis; NXBasis = NULL; }
  if( X != NULL ){ delete X; X = NULL; }
}

// Непосредственная реализация симплекс-метода с мультипликатывным
// представлением обратной базисной матрицы
void MultiSimplex( void )
{
  BeforeTimeWork(); // предварительная подготовка (поиск нач. базиса,...)

  EMatrix E( m + 1 ); // инициализация списка элементарных матриц

  Matrix Cx( 1, m ); // коэффициенты ЦФ при базисных переменных

  for( int iteration = 0; ; iteration++ )
  {
    int i;
    // заполняем матрицу Cx - коэф.-ты при ЦФ базисных переменных
    for( i = 0; i < m; i++ ) Cx[0][i] = C[ NXBasis[i] ];

    { // вывод значений базисных переменных
    char str[10]; Ccputs( 4, itoa(iteration,str,10) );
    OFile << "Итерация N" << iteration << endl
          << "Значения базисных переменных:\n";
    for( i = 0; i < m; i++ )
      OFile << "X" << (NXBasis[i]+1) << " = " << X[i] << endl;
    }
    { // вывод значения целевой функции
    Matrix bb(m,1,X);
    Matrix *aaa = Cx * bb;
    OFile << "Значение целевой функции: " << (*aaa) << endl;
    delete aaa;
    }

    // расчет двойственных перемнных
    Matrix Pi( 1, m );
    if( iteration )
    {
      for( i = iteration - 1; i >= 0; i-- )
      {
        double Sum = 0;

        for( int j = 0; j < m; j++ )
        {
          double a = E[i].Get(j);
          if( a )
          {
            double b = Cx[0][j];
            if( b )
              Sum += a * b;
          }
        }
        Cx[0][E[i].Column] = Sum;
      }
    }
    Pi = Cx;

    // нахождения вводимой переменной
    double c_min = 9e100, s_min = -1;
    for( i = 0; i < n; i++ )
    {
      Matrix *P_itoe = A.GetColumn( i );
      Matrix *P = Pi * (*P_itoe);

      double c = C[i] - (*P)[0][0];
      if( c < c_min ){ c_min = c; s_min = i; }

      delete P_itoe; delete P;
    }
    // вводим переменную с индексом s_min

    if( c_min >= -0.1e-12 )        // дальше оптимизировать некуда...
    {                              // -0.1e-12 - это почти ноль
      OFile << "\nДанное решение"  //              - дань погрешности :-(
               " является оптимальным!\n";
      break;
    }

    // вычисление ведущего столбца
    Matrix *P_soe = A.GetColumn( s_min );
    Matrix AA( m, 1 );
    AA = *P_soe;
    for( i = 0; i < iteration; i++ )
    {
      int clmn = E[i].Column;
      for( int j = 0; j < m; j++ )
      {
        if( j != clmn )
          AA[j][0] += E[i].Get(j) * AA[clmn][0];
      }
      AA[clmn][0] *= E[i].Get(clmn);
    }
    // ...он у нас теперь в AA

    // проверям, если в нем все элементы положительны, то
    // множество решений неограничено...
    for( i = 0; i < m; i++ ) if( AA[i][0] > 0  ) break;
    if( i == m )
    {
      OFile << "\nОптимум неограничен...!\n";
      break;
    }

    // ищем выводимую переменную из базиса
    double th_min = 9e100;
    int r_min = -1;
    for( i = 0; i < m; i++ )
    {
      if( AA[i][0] > 0 )
      {
        double th = X[i]/AA[i][0];
        if( th < th_min ){ th_min = th; r_min = i; }
      }
    }
    // выводим из базиса r_min-ую переменную

    // добавляем в список элементарных матриц еще одну
    // кроме, того изменяем попутно значения базисных переменных
    double *EE = new double [ m ]; if( EE == NULL ) NoMem;
    for( i = 0; i < m; i++ )
    {
      if( i != r_min )
      {
        EE[i] = - AA[i][0] / AA[r_min][0];
        X[i] = X[i] - th_min * AA[i][0];
      }
    }
    EE[r_min] = 1. / AA[r_min][0];
    X[r_min] = th_min;
    NXBasis[r_min] = s_min;
    E.Add( EE, r_min );
    delete EE;

    delete P_soe;
  }
  AfterTimeWork();
}

// Самая главная функция нашего метода
//  на основе входных данных создает выходной файл с выходными рез-ми
int MSimplex( void )
{
   OFile.open( OutFile );
   if( !OFile )
   {
     GError << "\nНевозможно создать выходной файл: " << OutFile;
     GlobalError();
     return 0;
   }

   { // вывод входных данных
   OFile << "Входные данные:\n";
   OFile << "Количество переменных  n = " << n << endl
         << "Количество ограничений m = " << m << endl;

   Matrix R( 1, n, C );
   OFile << "Вектор коэффициентов целевой функции C:\n" << R << endl;
   OFile << "Матрица коэффициентов в ограничениях A:\n" << A << endl;
   Matrix RR( 1, m, B );
   OFile << "Вектор правой части ограничений B:\n" << RR << endl;
   OFile << "********************************************************\n";
   }

   int value = setjmp( Gjmpbuf );
   if( value == 1 ){ return 0; } // For Whom The Bell Tolls...?
                                 // Time Matches On!
   MultiSimplex();

   OFile.close();

   return 1;
}

Файл MATRIX.CPP:
/*
  Файл содержит определение функций класса  Matrix  и  дру-
жественных функций этого класса. Это главный файл необходи-
мый для компоновки библиотеки MATRIX.LIB (см. MATRIX.H).
*/

#include <iostream.h>
#include <stdio.h>
#include <stdlib.h>
#include <mem.h>
#include "matrix.h"

// освобождение памяти под матрицу
void Matrix::destroy_matrix( void )
{
  if( A != NULL )
  {
    for( int i = 0; i < height; i++ ) if( A[i] != NULL ) delete A[i];
    delete A;
    A = NULL; width = height = 0;
  }
}

// (пере)установка размеров матрицы
void Matrix::SetSize( int h, int w )
{
  destroy_matrix();

  A = (double **) new char [ h * sizeof(double*) ];
  if( A != NULL )
  {
    for( int i = 0; i < h; i++ )
    {
      A[i] = new double [w];
      if( A[i] == NULL ) break;
    }
    if( i == h )
    {
      height = h; width = w;
      return;
    }
  }
  cout << "Не хватило памяти для создания матрицы ! "
            "SetSize( int, int ) !\n";
  exit( 1 );
}

void Matrix::Init( void )
{
  A = NULL; width = height = 0;
}

// конструктирование матрицы с помощью массива
Matrix::Matrix( int h, int w, double *a )
{
  Init();
  SetSize( h, w );
  if( a != NULL ) operator=( a );
    else operator=( 0. );
}

// конструктирование матрицы с помощью целого (int) массива
Matrix::Matrix( int h, int w, int *a )
{
  Init();
  SetSize( h, w );
  if( a != NULL ) operator=( a );
    else operator=( 0. );
}

// конструирование единичной матрицы
Matrix::Matrix( int hw )
{
  Init();
  SetSize( hw, hw );
  for( int i=0; i<height; i++ )
    for( int j=0; j<width; j++ )
      if( i == j ) A[i][j] = 1.;
      else A[i][j] = 0.;
}

// конструктирование матрицы с помощью другой матрицы
Matrix::Matrix( Matrix &M )
{
  Init();
  operator=( M );
}

// уничтожение объекта матрицы
void Matrix::destroy( void )
{
  destroy_matrix();
}

// the same...
Matrix::~Matrix()
{
  destroy();
}

// операция присваивания матрице массива
Matrix &Matrix::operator=( double *a )
{
  for( int i=0; i<height; i++ )
     for( int j=0; j<width; j++ ) A[i][j] = a[i*width+j];
  return *this;
}

// операция присваивания матрице целого (int) массива
Matrix &Matrix::operator=( int *a )
{
  for( int i=0; i<height; i++ )
     for( int j=0; j<width; j++ ) A[i][j] = a[i*width+j];
  return *this;
}

// операция заполнения матрицы числом
Matrix &Matrix::operator=( double a )
{
  for( int i=0; i<height; i++ )
     for( int j=0; j<width; j++ ) A[i][j] = a;
  return *this;
}

// операция присваивания матрице матрицы
Matrix &Matrix::operator=( Matrix &M )
{
  SetSize( M.height, M.width );
  for( int i=0; i<height; i++ )
     for( int j=0; j<width; j++ ) A[i][j] = M.A[i][j];
  return *this;
}

// операция сложения матриц
Matrix *operator+( Matrix &M1, Matrix &M2 )
{
  if( M1.width != M2.width || M1.height != M2.height )
  {
    cout << "Ошибка сложения! Несоответствие размеров матриц.\n";
    exit( 1 );
  }
  Matrix *Temp = new Matrix( M1.height, M1.width );
  if( Temp == NULL ){ cout << "Нет памяти для создания матрицы!\n"
                              "operator+(Matrix&, Matrix&)\n";
                      exit( 1 ); }
  for( int i=0; i<M1.height; i++ )
     for( int j=0; j<M1.width; j++ )
        Temp->A[i][j] = M1.A[i][j] + M2.A[i][j];
  return Temp;
}

// операция вычитания матриц
Matrix *operator-( Matrix &M1, Matrix &M2 )
{
  if( M1.width != M2.width || M1.height != M2.height )
  {
    cout << "Ошибка вычитания! Несоответствие размеров матриц.\n";
    exit( 1 );
  }
  Matrix *Temp = new Matrix( M1.height, M1.width );
  if( Temp == NULL ){ cout << "Нет памяти для создания матрицы!\n"
                              "operator-(Matrix&,Matrix&)\n";
                      exit( 1 ); }
  for( int i=0; i<M1.height; i++ )
     for( int j=0; j<M1.width; j++ )
        Temp->A[i][j] = M1.A[i][j] - M2.A[i][j];
  return Temp;
}

// операция умножения матриц
Matrix *operator*( Matrix &M1, Matrix &M2 )
{
  if( M1.width != M2.height )
  {
    cout << "Ошибка умножения!!!\n";
    exit( 1 );
  }
  Matrix *Temp = new Matrix( M1.height, M2.width );
  if( Temp == NULL ){ cout << "Нет памяти для создания матрицы!\n"
                              "operator*(Matrix&,Matrix&)\n";
                      exit( 1 ); }
  for( int i=0; i<Temp->height; i++ )
     for( int j=0; j<Temp->width; j++ )
     {
       Temp->A[i][j] = 0;
       for( int k=0; k<M1.width; k++ )
          Temp->A[i][j] += M1.A[i][k]*M2.A[k][j];
     }
  return Temp;
}

// операция умножения матрицы на число
Matrix *operator*( Matrix &M, double a )
{
  Matrix *Temp = new Matrix( M.height, M.width );
  if( Temp == NULL ){ cout << "Нет памяти для создания матрицы!\n"
                              "operator*(Matrix&,double)\n";
                      exit( 1 ); }
  for( int i=0; i<Temp->height; i++ )
     for( int j=0; j<Temp->width; j++ )
        Temp->A[i][j] = M.A[i][j] * a;
  return Temp;
}

// то же, что и предыдущее только ...
Matrix *operator*( double a, Matrix &M )
{
  return M*a;
}

// операция вывода матрицы на поток
ostream &operator<<( ostream &OS, Matrix &M )
{
  for( int i=0; i<M.height; i++ )
  {
    for( int j=0; j<M.width; j++ )
    {
      char *s = new char [200]; if( s == NULL ) return OS;
      sprintf( s, "%6.3f ", M.A[i][j] );
      OS << s;
      delete s;
    }
    OS << "\n";
  }
  return OS;
}

// опреация ввода матрицы с потока
istream &operator>>( istream &IS, Matrix &M )
{
  for( int i=0; i<M.height; i++ )
     for( int j=0; j<M.width; j++ ) IS >> M.A[i][j];
  return IS;
}

// операция сравнения матриц (возвращает 1, если они равны)
int operator==( Matrix &M1, Matrix &M2 )
{
  if( M1.width != M2.width || M1.height != M2.height ) return 0;
  for( int i=0; i<M1.height; i++ )
     for( int j=0; j<M1.width; j++ )
        if( M1.A[i][j] != M2.A[i][j] ) return 0;
  return 1;
}

// операция сравнения матриц (возвращает 1, если они не равны)
int operator!=( Matrix &M1, Matrix &M2 )
{
  if( M1.width != M2.width || M1.height != M2.height ) return 1;
  for( int i=0; i<M1.height; i++ )
     for( int j=0; j<M1.width; j++ )
        if( M1.A[i][j] != M2.A[i][j] ) return 1;
  return 0;
}

// операция транспонирования матрицы
Matrix *operator~( Matrix &M )
{
  Matrix *Temp = new Matrix( M.width, M.height );
  if( Temp == NULL ){ cout << "Нет памяти для создания матрицы!\n"
                              "operator~(Matrix&)\n";
                      exit( 1 ); }

  for( int i=0; i< M.height; i++ )
     for( int j=0; j< M.width; j++ ) Temp->A[j][i] = M.A[i][j];
  return Temp;
}

Matrix *operator+( Matrix &M, double a )
{
  Matrix *Temp = new Matrix( M.width, M.height );
  if( Temp == NULL ){ cout << "Нет памяти для создания матрицы!\n"
                              "operator+(Matrix&,double)\n";
                      exit( 1 ); }
  for( int i = 0; i < M.height; i++ )
    for( int j = 0; j < M.width; j++ )
      Temp->A[i][j] = M.A[i][j] + a;
  return Temp;
}

Matrix *operator+( double a, Matrix &M )
{
  return operator+( M, a );
}

Matrix *operator-( Matrix &M, double a )
{
  return operator+( M, -a );
}

Matrix *Matrix::GetColumn( int n )
{
  Matrix *Temp = new Matrix( height, 1 );
  if( Temp == NULL ){ cout << "Нет памяти для создания матрицы!\n"
                              "operator~(Matrix&)\n";
                      exit( 1 ); }
  for( int i = 0; i < height; i++ ) Temp->A[i][0] = A[i][n];
  return Temp;
}
Файл GERROR.CPP:
#include <tswm.h>
#include <iostream.h>

#include <gerror.h>

static char buf[1024];

ostrstream GError( buf, 1024 );  // строка с ошибкой
                    // заполняется во время возникновения
                    // ошибки, а затем вызывается функция
                    // GloabalError()

// Стандартная обработка глобальной ошибки
// Эту функцию можно написать заново для другой программной
// среды, так как эта - только для текстового режима MS-DOS
void GlobalError( void )
{
  GError << endl << ends;   // если кто-то забыл поставить конец сторки
  Warning( NULL, GError.str() );
  GError.seekp( 0 ); // для следующего сообщения
}