@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

@JohnS There's no inconsistency, and this behavior is the same as most other programming languages. Let me try to explain. M is the name (address of) a Matrix, which is a container; each M[i,j] is a name of (an address of) some object in the container. Changing the contents of the container does not change the container. Here's an analogy: M is a house, and each M[i,j] is a person living in the house. Even if every person is replaced by another person, it's still the same house, and it still has the same address.

But assigning to M is changing the house. Let's say that M stands for Mary's house (meaning the house owned by Mary---whether she actually lives there is irrelevant). That house has an address. Let's say that Mary sells this house and buys another house at a different address. She names the new house M. That's equivalent to assigning to M.

If the system weren't set up this way, then inplace operations would be impossible; and, as you noted earlier, that would be a major waste of memory.

@adrian5213 The command iquo(N,n,r) returns the integer quotient of N divided by n and returns the remainder in the variable r. If you need more information than that, you'll have to look at the source code, which is in the GNU Multi Precision Arithmetic Library and is open source.

Yes, of course you can use variables with &B as in

&B QR(z, rn);

There's no way to both represent something as an integer and for it to have leading zeros. If you want leading zeros, you should use a symbol or string representation.

P.S. A list representation, as suggested by John, also works. It's the option with the most flexibility, but also the least-compact representation, both internally and as displayed. It's because of this lack of compactness that I chose not to mention it earlier.

@Kitonum The order of equations may be irrelevant mathematically, but computationally it can have a major effect on the speed with which a solution is found.

@J Rod As a workaround, you can use character style Maple Input (aka 1D input) rather than the default 2D Input. The character style is the leftmost pull-down on the context toolbar (the bottommost toolbar). You can make this change permanently or for the session by using the Tools -> Options menu.

A great many long-time users, me included, use exclusively Maple Input.