@Thomas Dean That's completely separate from this Question. Please ask it as a new Question, not as an irrelevant Reply to another Question. A Reply is meant to be a response to the material that precedes it.

@roman_pearce How'd you get Maple to say that the system is inconsistent? And if it is, why did I get a result from the command in my Answer (using Maple 16)? Are you sure that you used the correct variable names? It was a little tricky for me to get them right.

Axel and Acer, please see my Answer in the original thread, where I use the `property/object` table to get the intersection of intervals.

@sand15 That was my Answer. I deleted it because on a subsequent reading I decided that it didn't properly address the Question. That's because it only determined whether the intervals have nonempty intersection rather than finding what the intersection was.

I'll put the Answer back on the other Question. Sorry for the confusion.

Please post your worksheet as an attached file. The parts of your post that aren't plaintext can't be copy-and-pasted.

Also, distance and distancevector are used in your code and are never defined.

@Kitonum 

I excluded 0 from the allowed input for the multiplier in addition thinking that that case was never useful (the output being the identity). It's a waste to do the matrix multiplication when the elementary matrix is the identity, so the call to addition should not be done when the multiplier is 0. That makes the code

reduced:= proc(B::Matrix)
uses LA= LinearAlgebra;
local
     M:= B, l:= 1, #l is current column.
     m:= LA:-RowDimension(M), n:= LA:-ColumnDimension(M), i, j
;
     for i to m do   #going through every row item
          #l needs to be less than column number n.
          if n < l then return M end if;
          j:= i;   #Initialize current row number.
          while M[j,l]=0 do   #Search for 1st row item <> 0.
               j:= j+1;
               if m < j then   #End of row: Go to next column.
                    j:= i;
                    l:= l+1;
                    if n < l then return M fi   #end of column and row
               end if
          end do;
          if j<>i then M:= perm(m,j,i).M end if;   #Permute rows j and i
          #Multiply row i with 1/M[i,l], if it's not 0.
          if M[i,l] <> 0 then M:= multiplikation(m,i,1/M[i,l]).M fi;
          #Subtract each row j with row i for M[j,l]-times.
          for j to m do if j<>i and M[j,l]<>0 then M:= addition(m,j,i,-M[j,l]).M fi od;
          l:= l+1   #Increase l by 1; next iteration i increase either.
     end do;
     return M
end proc:

I also restricted the coefficient domain to {float, rational}. This can be relaxed as needed.

@Kitonum It's not my procedure, nor did I claim that it was correct. I only made minor improvements to the OP's main procedure. When I said "here is the corrected reduced," I was refering only to the minor corrections that had already been discussed; and I recommended further testing.

I do claim correctness of the three elementary-matrix procedures.

@JohnS Note, however, that the way that felixp is using it doesn't cause the overwriting of the original argument. In other words, his code doesn't work inplace. His code generates a great many intermediate matrices that need to be garbage collected.

To get inplace operation, every statement of the form

M:= ...

(except for the initialization M:= B!) should be replaced with the equivalent statement of the form

M[..,..]:= ....

(Note that M[..,..] is the literal code; the dots are supposed to be there.)

All the elementary matrices would still need to be garbage collected.

@felixp Please see my code improvements in the Reply immediately above your most-recent Reply. I think that we may have been entering both Replies at the same time, so you may have missed it.

@acer Aha! I failed to test perm in the a=b case. Those zero rows are what causes the bug in reduced. And I forgot to mention that I'd already changed the typing of the last parameters of addition and multiplikation.

Below, I've taken the three elementary-matrix procedures and simplified the code, while maintaining the same logic as the originals, and added stringent type checking. I also changed them to only do square matrices, so only one dimension needs to be passed.

multiplikation:= proc(
     m::posint,
     a::depends(And(posint, satisfies(a-> a <= m))),
     b::And({float, rational}, Not(identical(0,0.)))
)
     Matrix((m,m), (i,j)-> `if`(i=j, `if`(i=a, b, 1), 0))
end proc:

addition:= proc(
     m::posint,
     a::depends(And(posint, satisfies(a-> a <= m))),
     b::depends(And(posint, satisfies(b-> b <= m))),
     c::And({float, rational}, Not(identical(0,0.)))
)
     Matrix((m,m), (i,j)-> `if`(i=a and j=b, c, `if`(i=j, 1, 0)))
end proc:

perm:= proc(
     m::posint,
     a::depends(And(posint, satisfies(a-> a <= m))),
     b::depends(And(posint, satisfies(b-> b <= m and a<>b)))
)
     Matrix((m,m), (i,j)-> `if`({i,j}={a,b} or i=j and not i in {a,b}, 1, 0))
end proc:

And here is the corrected reduced. In this, I only made minor corrections discussed here and added indenting and spacing.

reduced:= proc(B::Matrix)
uses LA= LinearAlgebra;
local
     M:= B, l:= 1, #l is current column.
     m:= LA:-RowDimension(M), n:= LA:-ColumnDimension(M), i, j
;
     for i to m do   #going through every row item
          #l needs to be less than column number n.
          if n < l then return M end if;
          j:= i;   #Initialize current row number.
          while M[j,l]=0 do   #Search for 1st row item <> 0.
               j:= j+1;
               if m < j then   #End of row: Go to next column.
                    j:= i;
                    l:= l+1;
                    if n < l then return M fi   #end of column and row
               end if
          end do;
          if j<>i then M:= perm(m,j,i).M end if;   #Permute rows j and i
          #Multiply row i with 1/M[i,l], if it's not 0.
          if M[i,l] <> 0 then M:= multiplikation(m,i,1/M[i,l]).M fi;
          #Subtract each row j with row i for M[j,l]-times.
          for j to m do if j<>i then M:= addition(m,j,i,-M[j,l]).M fi od;
          l:= l+1   #Increase l by 1; next iteration i increase either.
     end do;
     return M
end proc:

This works on the two test matrices given earlier. In particular,

M:= Matrix(2,3,[1,2,5,4,5,6]):
reduced(M);

You should test it on some extreme cases such as a zero matrix and matrices with some rows or some column all zeros.

@felixp The last line of procedure reduced should be

return M

However, there's still a bug in your procedure, because it returns the zero matrix. I haven't tracked this bug down yet. I've confirmed that your three elementary-matrix procedures are correct. To track down the bug, I recommend that you give the command

trace(reduced):

before running the procedure. This will show you a lot of what's happening inside the procedure as it's running.

Thanks for the very useful information. The timing of the for loops is quite a shame because the nested loops are hard to generalize to arbitrary loop depth.

@felixp Please post your code as it exists now.

@JohnS I should've mentioned this earlier: There's a simpler way that the OP felixp can achieve inplace operation. Rather than doing

M1:= copy(M);

he can do simply

M1:= M;

with local M1, and then continue assigning updates to M1.

This makes an address-only copy, not a true copy, so there's no doubling of memory usage.

I merely mentioning this as a possibility, since you brought up the memory issue; I'm NOT recommending it in this particular case. Inplace operation isn't something for newbies. Its inadvertent use is a major source of bugs in newbie code.

@vahid65 I thank you for your interest in my work, but I'm not currently associated with any university or other academic or research institution. I am a private tutor of mathematics and related subjects. I tutor both online and in person. My usage of Maple comes from my pure love of Maple and MaplePrimes and currently serves no other purpose. If you'd like to arrange an online session, contact me privately (use the More -> Contact Author menu at the bottom of this message).

Best regards,
Carl Love, formerly known as Carl Devore