This should probably be found in library(lists)

Program above, with comments, as posted on STack Overflow

What happens is:

• sudoku/1 imposes constraints on a 9x9 matrix of elements
  • with the help of blocks/2
• problem/2 lists a particular problem

And we can then put both together:

?- problem(1,Matrix),
   % Matrix is now a 9x9 matrix of partially bound elements.
   % Impose constraints; these suffice to find a unique solution
   sudoku(Matrix).

Matrix = [[9,8,7,6,5,4,3,2,1],
          [2,4,6,1,7,3,9,8,5],
          [3,5,1,9,2,8,7,4,6],
          [1,2,8,5,3,7,6,9,4],
          [6,3,4,8,9,2,1,5,7],
          [7,9,5,4,6,1,8,3,2],
          [5,1,9,2,8,6,4,7,3],
          [4,7,2,3,1,9,5,6,8],
          [8,6,3,7,4,5,2,1,9]].

The constraint solver starts to work immediately (apparently it detects that there is enough information to proceed; is that always the case?). There is no need to call label/1.

Let's see how the constraints between the matrix elements are set up:

sudoku(Rows) :-

   % If "Rows" is unbound, bind it to a list of 9 unbound variables.
   % If "Rows" is bound (as in our call above), this just verifies that there
   % are indeed 9 elements in Rows. Same as:
   %
   % Rows = [R1,
   %         R2,
   %         R3,
   %         R4,
   %         R5,
   %         R6,
   %         R7,
   %         R8,
   %         R9].

   length(Rows, 9),

   % For each element in "Rows", bind the element to a list of the same
   % length as "Rows" (i.e. 9); we now have a 9 x 9 matrix of possibly
   % unbound variables. A bit too clever maybe.
   %
   % Same as:
   %
   % Rows = [[E11,E12,E13,E14,E15,E16,E17,E18,E19],
   %         [E21,E22,E23,E24,E25,E26,E27,E28,E29],
   %         [E31,E32,E33,E34,E35,E36,E37,E38,E39],
   %         [E41,E42,E43,E44,E45,E46,E47,E48,E49],
   %         [E51,E52,E53,E54,E55,E56,E57,E58,E59],
   %         [E61,E62,E63,E64,E65,E66,E67,E68,E69],
   %         [E71,E72,E73,E74,E75,E76,E77,E78,E79],
   %         [E81,E82,E83,E84,E85,E86,E87,E88,E89],
   %         [E91,E92,E93,E94,E95,E96,E97,E98,E99]].

   maplist(same_length(Rows), Rows),

   % Concatenate all the variables in all the rows into a new list Vs (that
   % predicate is badly named, it concatenates lists form a list into a list)

   append(Rows, Vs),

   % A first constraint: every element Exx in Vs must be in domain [1..9]
   % If we passed in a matrix Rows with concrete values for Exx that are not
   % in domain [1..9], this will fail.

   Vs ins 1..9,

   % Another constraint: For every row (list of 9 elements) in Rows:
   % all the elements of the row must be distinct

   maplist(all_distinct, Rows),

   % We now want to impose the constraint that for every column,
   % all the elements of the column must be distinct

   % Transpose the 9x9 matrix Rows into the 9x9 matrix Columns.
   % (the new rows of Columns are the columns of Rows).
   %
   % Columns = [[E11,E21,E31,E41,E51,E61,E71,E81,E91],
   %            [E12,E22,E32,E42,E52,E62,E72,E82,E92]
   %            [E13,E23,E33,E43,E53,E63,E73,E83,E93]
   %            [E14,E24,E34,E44,E54,E64,E74,E84,E94]
   %            [E15,E25,E35,E45,E55,E65,E75,E85,E95]
   %            [E16,E26,E36,E46,E56,E66,E76,E86,E96]
   %            [E17,E27,E37,E47,E57,E67,E77,E87,E97]
   %            [E18,E28,E38,E48,E58,E68,E78,E88,E98]
   %            [E19,E29,E39,E49,E59,E69,E79,E89,E99]].

   transpose(Rows, Columns),

   % Another constraint: For every "column of Rows"
   % all the elements of the column must be distinct

   maplist(all_distinct, Columns),

   % Now we need to impose the "all-distinct" constraint
   % on the 3x3 sub-matrices

   % Give the rows distinct names

   Rows = [As,Bs,Cs,Ds,Es,Fs,Gs,Hs,Is],

   % Impose constraint on "3-element-wide" groups of the
   % "3-row-high" group of As, Bs, Cs: namely, "all distinct"

   blocks(As, Bs, Cs),
                   
   % Same of the next "3-high" group
   
   blocks(Ds, Es, Fs),

   % Same of the next "3-high" group

   blocks(Gs, Hs, Is).