Solve the chromatic puzzle Over at our friends at Puzzling.SE, the following puzzle was posted: Is this chromatic puzzle always solvable? by Edgar G. You can play it here.

Puzzle explanation

Given a m x n grid with tiles of three different colours, you may select any two adjacent tiles, if their colours are different. These two tiles are then converted to the third colour, i.e., the one colour not represented by these two tiles. The puzzle is solved if all tiles have the same colour. Apparently, one can prove that this puzzle is always solvable, if  neither m nor n are divisible by 3.



Of course, this begs for a solving algorithm. You will write a function or program that solves this puzzle. Note that functions with 'side effects' (i.e., the output is on stdout rather than in some awkward data type return value) are explicitly allowed.

Input &amp; Output

The input will be an m x n matrix consisting of the integers 1, 2 and 3 (or 0, 1, 2 if convenient). You may take this input in any sane format. Both m and n are &gt;1 and not divisible by 3. You may assume the puzzle is not solved 

You will then solve the puzzle. This will involve a repeated selection of two adjacent tiles to be 'converted' (see above). You will output the two coordinates of these tiles for each step your solving algorithm took. This may also be in any sane output format. You are free to choose between 0-based and 1-based indexing of your coordinates, and whether rows or columns are indexed first. Please mention this in your answer, however.

Your algorithm should run within reasonable time on the original 8x8 case. Brute-forcing it completely is explicitly disallowed, i.e. your algorithm should run under O(k^[m*(n-1)+(m-1)*n]) with k the number of steps needed for the solution. The solution however is not required to be optimal. The proof given in the linked question may give you an idea as to how to do this (e.g., first do all columns using only vertically adjacent tiles, and then do all rows)

Test cases

In these test cases, the coordinates are 1-based and rows are indexed first (like MATLAB/Octave and probably many others).

Input: 
[1 2]
Output: (result: all 3's)
[1 1],[1,2]


Input:
[ 1 2
  3 1 ]
Output: (result: all 1's)
[1 1],[2 1]        (turn left column into 2's)
[2 1],[2 2]        (turn right column into 3's)
[1 1],[1 2]        (turn top row into 1's)
[2 1],[2 2]        (turn bottom row into 1's)

Input:
[1 2 3 2
 3 2 1 1]

Output: (result: all 3's)
[1 1],[1 2] 
[1 3],[1 4] 
[1 2],[1 3] 
[1 1],[1 2] 
[1 2],[1 3] 
[1 1],[1 2]
[1 3],[1 4]
[2 1],[2 2]
[1 1],[2 1]
[1 2],[2 2]
[1 3],[2 3]
[1 4],[2 4]


If desired, I may post a pastebin of larger test cases, but I think this should be sufficient.
Octave, 334 313 bytes

Since the challenge may seem a bit daunting, I present my own solution. I did not formally prove that this method works (I guess that will come down to proving that the algorithm will never get stuck in a loop), but so far it works perfectly, doing 100x100 testcases within 15 seconds. Note that I chose to use a function with side effects rather than one that returns all the coordinates since that saved me a few bytes. Coordinates are row-major, 1-based, and formatted as row1 col1 row2 col2. Input colours are 0,1,2 since this works better with mod, at the cost of having to use numel rather than nnz. Golfed version: Edit: saved another few bytes by using a technique from Kevin Lau's answer.

function P(p)
k=0;[m,n]=size(p);t=@(v)mod(sum(v)*numel(v),3);for c=1:n
p(:,c)=V(p(:,c));end
k=1;for r=1:m
p(r,:)=V(p(r,:));end
function a=V(a)
while any(a-(g=t(a)))
isempty(q=find(diff(a)&amp;a(1:end-1)-g,1))&amp;&amp;(q=find(diff(a),1));
a([q q+1])=t(a([q q+1]));if k
disp([r q r q+1])
else
disp([q c q+1 c])
end;end;end;end


Example GIF of the solving algorithm:



Ungolfed version:

function solveChromaticPuzzle(p)
[m,n]=size(p);                           % Get size
t=@(v)mod(sum(v)*numel(v),3);            % Target colour function
for c=1:n                                % Loop over columns
    p(:,c)=solveVec(p(:,c));             % Solve vector
end
for r=1:m                                % Loop over rows
    p(r,:)=solveVec(p(r,:));
end
    function a=solveVec(a)               % Nested function to get globals
        g=t(a);                          % Determine target colour
        while(any(a~=g))                 % While any is diff from target...
            % Do the finding magic. Working left-to-right, we find the
            % first pair that can be flipped (nonzero diff) that does not
            % have the left colour different from our goal colour
            q=min(intersect(find(diff(a)),find(a~=g)));
            if(isempty(q))               % In case we get stuck...
                q=find(diff(a),1);       % ... just flip the first possible
            end;
            a([q q+1])=t(a([q q+1]));    % Do the actual flipping.
            if(exist('r'))               % If we're working per row
                disp([r q r q+1])        % Print the pair, using global row
            else
                disp([q c q+1 c])        % Print the pari, using global col
            end
        end
    end
end