Разрешите не описывать программу полностью, это позволит сосредоточится на важных моментах. Организовать ввод данных может каждый программист с маломальским опытом. Все примеры и куски программных кодов написаны на языке MASM-64 для 64-разрядной Windows. Сразу разберем пример и методы его решения. Имеем следующее:
Nmatr – размерность матрицы, точнее количество строк. Количество столбцов на 1 больше.
lpMatr– указатель на область памяти, где хранятся элементы матрицы. Объем памяти = Nmatr * (Nmatr+1) * 10 байт.
indRow– массив смещений строк.
indCol– массив смещений столбцов.
iStep – номер выполняемого шага, от 0 до Nmatr-1
cntElem – начальное состояние счетчиков
Для наглядности примем Nmatr = 4, а саму матрицу заполним следующими числами:
241 0.0178 32 5879 23853.0356
7 254 0.125 845 3895.3750
94.2 154 748 14 2702.2000
158 158 45 45 789.0000
Чтобы не было сомнений, я заранее задал корни решения и посчитал свободные члены в Excel, так чтобы в конце работы можно было проверить результат. Массив индексов для столбцов заполняется, как смещения элементов от начала строки, а массив индексов для строк заполняется, как смещения начала строк от начала матрицы:
Столбцы = 0, 10, 20, 30, 40;
Строки = 0, 40, 80, 120;
Как мы уже говорили, будем искать число с максимальным порядком вместо реально максимального числа, так как на точность работы сопроцессора мантисса почти не влияет. В шестнадцатиричном виде первые четыре коэффициенты выглядят так
4006 F100000000000000
3FF9 91D14E3BCD35A858
4004 8000000000000000
400B B7B8000000000000
Первые два байта – порядок, далее мантисса. Правда для процессора они последние, ведь в памяти компьютера числовые переменные хранятся с обратным следованием байт. Сравнивая десятичную форму и шестнадцатиричную, легко видеть, что достаточно анализировать только порядок, чтобы выбрать одно из наибольших чисел. Выполнение поиска максимального числа сводится к следующей подпрограмме:
mov rbx, lpMatr ; базовый указатель на матрицу
mov r15, lpCol ; базовый указатель на смещения столбцов
mov r14, lpRow ; базовый указатель на смещения строк
mov iStep, 0 ; счетчик шагов.
; Будет отрезать обработанные строки и столбцы
mov rax, Nmatr
mov cntElem, rax ; размерность обрабатываемой области
; Будет укорачиваться по мере обработки строк и столбцов
ciclStepDn: ; цикл прямого хода метода Гаусса
mov r12, iStep ; индекс верхней строки и левого столбца,
; после обрезания обработанных строк и столбцов.
xor rdx, rdx ; будущий максимальный порядок
xor r10, r10 ; индекс строки с максимумом
xor rbx, rbx ; индекс столбца с максимумом
mov r8, cntElem ; счетчик строк
Cicl1:
mov rcx, cntElem ; счетчик столбцов
mov r11, iStep ; индекс текущего столбца
Cicl2:
mov rax, qword ptr [r14 + r12*8] ; смещение строки
add rax, qword ptr [r15 + r11*8] ; плюс смещение столбца
movzx rax, word ptr [rbx+rax+8] ; взять порядок
and rax, 07FFFh ; отбросить знаковый бит
cmp rax, rdx ; сравнить с текущим максимумом
jle neMax ; если меньше или равно - обойти
mov rdx, rax ; обновить текущий максимум
mov r10, r11 ; обновить индекс столбца максимума
mov rbx, r12 ; обновить индекс строки максимума
neMax:
inc r11 ; переходим к следующему элементу
loop Cicl2 ; переходим на начало цикла Cicl2
; и уменьшаем счетчик столбцов
inc r12 ; увеличиваем индекс текущей строки
dec r8 ; уменьшаем счетчик строк
jnz Cicl1 ; переходим на начало цикла Cicl1
test rdx, rdx
jz Avaria ; порядок максимального числа равен нулю,
; дальнейший расчет нецелесообразен
mov r12, iStep ; переход к угловому элементу
mov r11, iStep
; mov r8, qword ptr [r14 + r12*8] ; смещение строки угла
; add r8, qword ptr [r15 + r11*8] ; плюс смещение столбца угла
; mov r9, qword ptr [r14 + rbx*8] ; смещение строки максимума
; add r9, qword ptr [r15 + r10*8] ; плюс смещ столбца максимума
; xchg
inc iStep ; переход к следующему шагу
dec cntElem ; уменьшаем число строк и столбцов на следующем шаге
jne ciclStepDn ; продолжение расчета или выход
Как видим, самый глубокий цикл Cicl2 занимает всего 11 машинных команд, значит выполняться он будет менее чем за 70 наносекунд. Если бы проводили сравнение с помощью команд сопроцессора время было бы в десятки раз больше.
Теперь нужно переставить смещения наибольшего числа в верхний угол рассматриваемого участка матрицы и можно проводить вычисления.
Мне кажется, проще предоставить на ваше рассмотрение отлаженный фрагмент программы, чтобы Вы могли сами оценить его возможности. Если появятся непонятные моменты, всегда рад помочь.
Ниже приведен текст полностью отлаженной программы, решающей заданную систему уравнений.
; ¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤
.NOLIST
.NOCREF
include \masm32\include64\masm64rt.inc
.LIST
.data
Nmatr dq 4
lpMatr dq Matrica ; в реальной программе это будет указатель
lpCol dq indCol
lpRow dq indRow
indCol dq 0,10,20,30,40
indRow dq 0,50,100,150
Matrica dt 241.0, 0.0178, 32.0, 5879.0, 23853.0356
dt 7.0, 254.0, 0.125, 845.0, 3895.3750
dt 94.2, 154.0, 748.0, 14.0, 2702.2000
dt 158.0, 158.0, 45.0, 45.0, 789.0000
.data?
iStep dq ?
cntElem dq ?
.code
; ¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤
entry_point proc
invoke CalcGauss
conout "Howdy, your new console template here.",lf,lf
waitkey
invoke ExitProcess,0
ret
entry_point endp
; ¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤
CalcGauss proc
LOCAL srsi:QWORD, srdi:QWORD, srbx:QWORD
LOCAL sr12:QWORD, sr13:QWORD, sr14:QWORD, sr15:QWORD
mov srsi, rsi
mov srdi, rdi
mov srbx, rbx
mov sr12, r12
mov sr13, r13
mov sr14, r14
mov sr15, r15
mov rbx, lpMatr ; базовый указатель на матрицу
mov r15, lpCol ; базовый указатель на смещения столбцов
mov r14, lpRow ; базовый указатель на смещения строк
mov iStep, 0 ; счетчик шагов.
; Будет отрезать обработанные строки и столбцы
mov rax, Nmatr
mov cntElem, rax ; размерность обрабатываемой области
; Будет укорачиваться по мере обработки строк и столбцов
ciclStepDn: ; цикл прямого хода метода Гаусса
mov r12, iStep ; индекс верхней строки и левого столбца,
; после обрезания обработанных строк и столбцов.
xor rdx, rdx ; будущий максимальный порядок
xor r10, r10 ; индекс строки с максимумом
xor r9, r9 ; индекс столбца с максимумом
mov r8, cntElem ; счетчик строк
Cicl1:
mov rcx, cntElem ; счетчик столбцов
mov r11, iStep ; индекс текущего столбца
Cicl2:
mov rax, qword ptr [r14 + r12*8] ; смещение строки
add rax, qword ptr [r15 + r11*8] ; плюс смещение столбца
movzx rax, word ptr [rbx+rax+8] ; взять порядок
and rax, 07FFFh ; отбросить знаковый бит
cmp rax, rdx ; сравнить с текущим максимумом
jle neMax ; если меньше или равно - обойти
mov rdx, rax ; обновить текущий максимум
mov r10, r11 ; обновить индекс столбца максимума
mov r9, r12 ; обновить индекс строки максимума
neMax:
inc r11 ; переходим к следующему элементу
loop Cicl2 ; переходим на начало цикла Cicl2
; и уменьшаем счетчик столбцов
inc r12 ; увеличиваем индекс текущей строки
dec r8 ; уменьшаем счетчик строк
jnz Cicl1 ; переходим на начало цикла Cicl1
test rdx, rdx
jz Avaria ; порядок максимального числа равен нулю,
; дальнейший расчет нецелесообразен
mov r12, iStep ; переход к угловому элементу
mov rax, qword ptr [r14 + r9*8] ; смещение строки
xchg qword ptr [r14+r12*8], rax
mov qword ptr [r14 + r9*8], rax ; поменяли смещения строк
mov rax, qword ptr [r15 + r10*8] ; смещение столбцов
xchg qword ptr [r15+r12*8], rax
mov qword ptr [r15 + r10*8], rax ; поменяли смещения столбцов
; надо разделить элементы первой строки на угловой элемент
mov r12, iStep ; индекс текущей строки
mov r11, iStep ; индекс текущего столбца
mov rcx, cntElem ; счетчик столбцов
mov rax, qword ptr [r14 + r12*8] ; смещение строки
add rax, qword ptr [r15 + r11*8] ; плюс смещение столбца
fld tbyte ptr [rbx + rax] ; угловой элемент
CiclDivUgol:
inc r11
mov rax, qword ptr [r14 + r12*8] ; смещение строки
add rax, qword ptr [r15 + r11*8] ; плюс смещение столбца
fld tbyte ptr [rbx + rax] ; элемент строки
fdiv st, st(1) ; разделить
fstp tbyte ptr [rbx + rax] ; сохранить
loop CiclDivUgol
finit
mov r12, iStep ; индекс текущей строки
lea r13, [r12 + 1] ; индекс следующей строки
mov r11, iStep ; индекс текущего столбца
mov r8, cntElem
dec r8 ; счетчик строк, из которых будем вычитать
jz NextStepDwn
CiclSubRows:
lea r10, [r11 + 1] ; индекс следующего столбца
mov rcx, cntElem ; счетчик столбцов
mov rax, qword ptr [r14 + r13*8] ; смещение строки
add rax, qword ptr [r15 + r11*8] ; плюс смещение столбца
fld tbyte ptr [rbx + rax] ; множитель
CiclSubCols:
mov rdx, qword ptr [r14 + r13*8] ; смещение строки
add rdx, qword ptr [r15 + r10*8] ; плюс смещение столбца,
; адрес запомним, чтоб назад положить результат
fld tbyte ptr [rbx + rdx] ; элемент вычитания
mov rax, qword ptr [r14 + r12*8] ; смещение строки
add rax, qword ptr [r15 + r10*8] ; плюс смещение столбца множимого
fld tbyte ptr [rbx + rax] ; элемент над текущим
fmul st, st(2) ; умножаем
fsubp st(1), st ; вычитаем
fstp tbyte ptr [rbx + rdx] ; кладем на место
inc r10 ; следующий столбец
loop CiclSubCols ; возврат цикла по столбцам
finit
inc r13 ; следующая строка
dec r8
jnz CiclSubRows ; возврат цикла по строкам
NextStepDwn:
; шаг завершен
inc iStep ; переход к следующему шагу
dec cntElem ; уменьшаем число строк и столбцов на следующем шаге
jne ciclStepDn ; продолжение расчета или на обратный ход метода
mov r12, Nmatr ; индекс свободного члена
lea r8, [r12-2] ; индекс необработанных строк
ciclObrXod:
mov rdx, qword ptr [r14 + r8*8] ; смещение необработанной строки
add rdx, qword ptr [r15 + r12*8] ; смещение свободного члена
fld tbyte ptr [rbx + rdx] ; загружен свободный член
lea r10, [r12-1] ; индекс последнего коэффициента
lea rcx, [r12-1]
sub rcx, r8 ; счетчик обратных вычислений
ciclObrCalc:
mov rax, qword ptr [r14 + r10*8] ; смещение обработанной строки
add rax, qword ptr [r15 + r12*8] ; смещение готового корня
fld tbyte ptr [rbx + rax] ; загружен готовый корень
mov rax, qword ptr [r14 + r8*8] ; смещение необработанной строки
add rax, qword ptr [r15 + r10*8] ; смещение множителя
fld tbyte ptr [rbx + rax] ; загружен множитель
fmulp
fsubp st(1), st ; вычитаем произведение из свободного члена
dec r10
loop ciclObrCalc
fstp tbyte ptr [rbx + rdx] ; сохранить вычисленный корень
dec r8
jge ciclObrXod ; повторить для следующего корня. r8==0 тоже обрабатывать
; далее идет загрузка корней в обратном порядке, чтобы в стеке выглядело правильно
mov rax, qword ptr [r14]
add rax, qword ptr [r15 + r12*8]
fld tbyte ptr [rbx + rax]
mov rax, qword ptr [r14+8]
add rax, qword ptr [r15 + r12*8]
fld tbyte ptr [rbx + rax]
mov rax, qword ptr [r14+16]
add rax, qword ptr [r15 + r12*8]
fld tbyte ptr [rbx + rax]
mov rax, qword ptr [r14+24]
add rax, qword ptr [r15 + r12*8]
fld tbyte ptr [rbx + rax]
nop
nop ; можно взглянуть на результат проделанной работы
nop
finit
Avaria:
mov rsi, srsi
mov rdi, srdi
mov rbx, srbx
mov r12, sr12
mov r13, sr13
mov r14, sr14
mov r15, sr15
ret
CalcGauss endp
;======================================================
End
Открою наконец секрет, какие же корни были заложены в эти уравнения. Это 1, 2, 3, 4. А содержимое стека сопроцессора после проведения вычислений приведено в заставке к статье.
Спасибо, что дочитали программу до конца.
