diff --git a/library/examples/matrix.yap b/library/examples/matrix.yap
index 60cfaeb39..ec33454d3 100644
--- a/library/examples/matrix.yap
+++ b/library/examples/matrix.yap
@@ -1,5 +1,6 @@
 
 :- use_module(library(matrix)).
+:- use_module(library(maplist)).
 
 t1 :-
 	X <== matrix([1,2,3,4,5,6],[dim=[3,2]]),
@@ -27,3 +28,32 @@ numbers(I0, I1, Vals) :-
 	( I0 =< I1 -> Vals = [I0|MVals], I01 is I0+1, numbers(I01, I1, MVals) ;
 	    Vals = [] ).
 
+t5 :-
+	numbers(1, 100, L),
+	X <== matrix(L, [dim=[10,10]]),
+	writeln('diagonal:'),
+	foreach([I in 0..9, J in I..I], Y^(Y <== X[I,J], writeln(Y) ) ).
+t6 :-
+	Len = 10,
+	LenSq is Len*Len,
+	Len1 is Len-1,
+	numbers(1, LenSq, L),
+	X <== matrix(L, [dim=[Len,Len]]),
+	Y <== matrix(L, [dim=[Len,Len]]),
+	Z <== matrix(L, [dim=[Len,Len]]),
+	writeln('product:'),
+	foreach([I in 0..Len1, J in 0..Len1], step(X,Y,Z,I,J) ),
+	O <== list(Z),
+	writeln(O).
+
+step(X,Y,Z,I,J) :-
+	Xs <== X[I,_],
+	Ys <== Y[_,J],
+	foldl(addprod, Xs, Ys, 0, P), % scalar product
+	Z[I,J] <== P.
+
+addprod(X, Y, S0, S) :-
+	S is S0+X*Y.
+
+
+
diff --git a/library/matrix.yap b/library/matrix.yap
index e04815204..57ff0b338 100644
--- a/library/matrix.yap
+++ b/library/matrix.yap
@@ -91,11 +91,14 @@ typedef enum {
 	    matrix_column/3,
 	    matrix_get/2,
 	    matrix_set/2,
+	    foreach/2,
 	    op(100, fy, '[]')
 	    ]).
 
 :- load_foreign_files([matrix], [], init_matrix).
 
+:- meta_predicate foreach(+,:).
+
 :- use_module(library(maplist)).
 :- use_module(library(lists)).
 
@@ -498,20 +501,21 @@ el_list([N], El, Els, NEls, I0, I1) :-
 	append(El, NEls, Els),
 	I1 is I0+1.
 
-foreach( Domain, Locals, Goal) :-
-	global_variables( V, Locals, Goal, GlobalVars ),
+foreach( Domain, M:(Locals^Goal)) :- !,
+	global_variables( Domain, Locals, Goal, GlobalVars ),
+	iterate( Domain, [], GlobalVars, M:Goal, [], [] ).
+foreach( Domain, Goal ) :-
+	global_variables( Domain, [], Goal, GlobalVars ),
 	iterate( Domain, [], GlobalVars, Goal, [], [] ).
 
-global_variables( Vs, Locals, Goal, GlobalVars ) :-
-	term_variables( [V|Locals], Pars ),
-	term_variables( Goal, GVs, Pars),
-	foldl( del, Pars, GVs, GlobalVars ).
-
-del(Ps, Vs, RVs ) :-
-	foldl(delv, Ps, Vs, RVs).
+global_variables( Domain, Locals, Goal, GlobalVars ) :-
+	term_variables( Domain+Locals, Pars ),
+	term_variables( Goal, DGVs, Pars),
+	sort( DGVs, GVs ),
+	foldl( delv, Pars, GVs, GlobalVars ).
 
 delv( V, [V1|Vs], Vs) :- V == V1, !.
-delv( V, [_|Vs], NVs) :-
+delv( V, [V1|Vs], [V1|NVs]) :-
 	delv( V, Vs, NVs).
 
 iterate( [], [], GlobalVars, Goal, Vs, Bs ) :-
@@ -525,19 +529,26 @@ iterate( [H|L], Cont, GlobalVars, Goal, Vs, Bs ) :- !,
 iterate( [] ins _A .. _B, Cont, GlobalVars, Goal ) :- !,
 	iterate(Cont, [], GlobalVars, Goal, Vs, Bs ).
 iterate( [V|Ps] ins A..B, Cont, GlobalVars, Goal, Vs, Bs  ) :-
-	eval(A, NA),
-	eval(B, NB),
-	( NA >= NB ->  iterate( Cont, [], GlobalVars, Goal, [V|Vs], [B|Bs] ) ;
+	eval(A, Vs, Bs, NA),
+	eval(B, Vs, Bs, NB),
+	( NA > NB ->  true ;
 	  A1 is NA+1,
-	  iterate( Ps ins A1..NB, GlobalVars, Goal, [V|Vs], [A|Bs] )
+	  iterate( Cont, [], GlobalVars, Goal, [V|Vs], [NA|Bs] ),
+	  iterate( Ps ins A1..NB, GlobalVars, Goal, [V|Vs], [NA|Bs] )
 	).
 iterate( V in A..B, Cont, GlobalVars, Goal, Vs, Bs) :-
 	var(V),
-	eval(A, NA),
-	eval(B, NB),
-	( NA >= NB -> iterate( Cont, [], GlobalVars, Goal, [V|Vs], [B|Bs] ) ;
+	eval(A, Vs, Bs, NA),
+	eval(B, Vs, Bs, NB),
+	( NA > NB -> true ;
 	  A1 is NA+1,
-	  iterate( V in A1..NB, Cont, GlobalVars, Goal, [V|Vs], [B|Bs] ) 
+	  iterate( Cont, [], GlobalVars, Goal, [V|Vs], [NA|Bs] ),
+	  iterate( V in A1..NB, Cont, GlobalVars, Goal, Vs, Bs ) 
 	).
 
 	 
+eval(I, _Vs, _Bs, I) :- integer(I), !.
+eval(I, Vs, Bs, NI) :-
+	copy_term(I+Vs, IA+Bs),
+	NI <== IA.
+