Пусть верхний левый угол матрицы имеет индекс (1,1).

Будем для каждой строки i формировать вектор p[1..M] такой, что p[j] есть число последовательных единичных элементов в столбце j, начиная с элемента (i,j) и выше его. Таким образом, если p[j]=0, то A[i,j]=0, а если p[j]=u>0, то

A[i,j]=A[i,j+1]= ... =A[i,j-u+1]=1,

а элемент A[i,j-u] нулевой (если, конечно, такой элемент есть в матрице, т.е. если j-u>0).

Тогда площадь (сумма элементов) единичного прямоугольника S_i(L,R) с нижними правым и левым углами в элементах (i,R) и (i,L) соответственно есть площадь основания (R-L+1) умноженная на высоту этого прямоугольника

h(L,R)=минимальное p[i] по всем j, изменяющимся от L до R.

Для каждой строки i надо найти максимум величины S_i(L,R) при 1<=L<=R<=M, а максимум по всем строкам и даст искомую величину.

Очевидно, что для первой строки

p[j]=A[1,j].

Если мы знаем вектор l для строки i, то мы можем вычислить его для строки (i+1) по следующей формуле

if A[i+1,j]=0 then p[j]:=0
    else p[j]:=p[j]+1;

Более коротко это же можно записать и так: p[j]:=(p[j]+1)*A[i+1,j];

Будем рассматривать вектор p, соответствующий строке i, и для каждого индекса j, 1<=j<=M, определим максимальный размер единичного прямоугольника с высотой p(j), располагающегося в столбцах с номерами от L(j) до R(j), L(j)<=j<=R(j), нижняя граница которого строка i. Очевидно, что L(j) есть увеличенный на единицу индекс первого меньшего чем p[j] элемента вектора p при просмотре p от j-го элемента влево, или L(j)=1, если такого меньшего элемента нет. Аналогично, R(j) есть уменьшенный на 1 индекс первого меньшего чем p[j] элемента вектора p при просмотре p от j-го элемента вправо, или R(j)=M, если такого элемента нет.

Как быстро вычислить L(j) и R(j)? Используя алгоритм 8 Главы "Структуры данных" мы можем найти все L и R за два прохода по вектору p:

Будем заполнять массив R. Вектор p просматриваем слева направо. Организуем стек для позиций элементов. Для каждого текущего p[j] будем выталкивать из стека все позиции, на которых стоят элементы большие текущего, и заносить в соответствующие этим позициям места массива R число (j-1). Замет позицию текущего элемента j поместить в стек. После просмотра всех элементов в стеке будут стоять индексы позиций массива R, в которые необходимо занести число M.

var Stack: array [0..M] of byte;
...
S[0]:=0; {стек пуст, S[0] есть указатель на последнюю
помещенную в стек позицию}
for j:=1 to M do
begin
while p[j]<p[S[S[0]]] do
{S[0] - это индекс последней занятой позиции в стеке, 
на которой находится число S[S[0]] - индекс элемента массива p, 
а сам этот элемент - это p[S[S[0]]]}
begin
R[S[S[0]]]:=j-1; {для элемента массива p с индексом S[S[0]] 
нашли координату правой стенки}
S[0]:=S[0]-1;{убрать элемент из стека}
end;
S[0]:=S[0]+1;{индекс - в стек}
S[S[0]]:=j;
end;
for j:=1 to S[0] do R[S[S[0]]]:=M;

Для заполнения массива L необходимо проделать те же самые операции, но вектор p будет просматриваться справа налево:

...

S[0]:=0; {стек пуст, S[0] есть указатель на последнюю
помещенную в стек позицию}
for j:=M downto 1 do
begin
while p[j]<p[S[S[0]]] do
begin
L[S[S[0]]]:=j+1; {для элемента массива p с индексом S[S[0]]
 нашли координату левой стенки}
S[0]:=S[0]-1;{убрать элемент из стека}
end;
S[0]:=S[0]+1;{индекс - в стек}
S[S[0]]:=j;
end;
for j:=1 to S[0] do L[S[S[0]]]:=1;

Последнее, что остается сделать - это за один проход вычислить максимум по всем j выражения

p[j]*(R[j]-L[j]+1).

Этот максимум есть размер максимального единичного прямоугольника с нижней гранью в строке i.

Максимум по всем i и даст решение.

Сложность алгоритма O(N*M), т.е. количество операции линейно зависит от числа элементов матрA!