| ?- sudoku(X, Y).
X = [_#3(1..4),_#24(1..4),...etc.]
Y = [_#3(1..4),_#24(1..4),...etc.]
In the article, a more concrete query is also shown, as well as its solution. This solution is found via constraint propagation alone. This means that the system has deduced the unique solution purely by reasoning about the available constraints, which it does automatically every time a constraint is posted. This is one example of an algorithm that does not guess.
This also shows that even supposing you cannot make an algorithm that "never guesses", you can definitely make an algorithm that has to do much less guessing than searching over the whole tree. In Prolog, such algorithms are implemented and encapsulated in powerful combinatorial constraints like all_distinct/1, which are available in many Prolog systems and which you can use in your applications. Internally, the associated constraint propagators prune the search tree by applying various algorithms for you behind the scene. The ability to use such constraints and their strong pruning are major attractions of using Prolog for combinatorial tasks.
You are right that there may be cases where such propagation, albeit quite strong, may not suffice to fully solve a concrete combinatorial task. For this reason, you have to apply a concrete enumeration of remaining variables. In Constraint Logic Programming over Integers, this search is called labeling and provided by predicates like fd_labeling/2 or similar, depending on your Prolog system. The article does not use them though, and even if it did, the search could be guided by various heuristics by simply supplying a few options, which together with the pruning applied by constraints distinguish such a search from more uninformed search strategies.
Note also that the posted solution is quite hard to generalize elegantly. A shorter and equivalent Prolog program that describes 4x4 Sudoku puzzles is:
sudoku(Rows) :-
Rows = [Row1,Row2,Row3,Row4],
maplist(same_length(Rows), Rows),
blocks(Row1, Row2), blocks(Row3, Row4),
append(Rows, Vs),
fd_domain(Vs, 1, 4),
maplist(fd_all_different, Rows),
transpose(Rows, Cols),
maplist(fd_all_different, Cols).
blocks([], []).
blocks([A1,A2|As], [B1,B2|Bs]) :-
fd_all_different([A1,A2,B1,B2]),
blocks(As, Bs).