Carl Love

Carl Love

28055 Reputation

25 Badges

12 years, 364 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are replies submitted by Carl Love

@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.

@JohnS 

Inplace operation is always possible without using LinearAlgebra and without using global. Though I do wonder why you'd be restricted to not using any command from LinearAlgebra. Perhaps it's an exercise.

One is allowed to assign to a parameter subscripted. Thus, the example that you mentioned can be coded like this:

MyProc:= proc(
     M::Matrix, x, y,
     i::depends(And(posint, satisfies(i-> i <= op([1,1], M)))),
     j::depends(And(posint, satisfies(j-> j <= op([1,2], M))))
)
     M[i,j]:= M[i,j]*x^i*y^j;
     [][]
end proc:

And used like this:

A:= LinearAlgebra:-RandomMatrix(2,2);

MyProc(A, x, y, 2, 2);
A;

If a procedure performs inplace operations, then there's no need for it to have a return value. So it's best to explicitly make the return value NULL. The symbol that I like to use, [][], is equivalent to NULL.

 

@JohnS Yes, it's possible to do inplace operations, and there's no need to make global. I didn't mention this possibility because perhaps it's a bit too sophisticated for this OP at this time. All matrix arithmetic commands can be handled by commands in the LinearAlgebra package, and all such commands have an inplace option. An example of such a command that can be used in this case is MatrixMatrixMultiply

Note that it's usually not appropriate for a procedure to overwrite a parameter unless the caller specifically requests inplace operation.

@Østerbro Here's the meaning of the error message: A name is an object that can be assigned to. Think of it as a variable name. Anything else is an object. So, if you try to assign to something other than a name, it's an "illegal use of an object as a name."

First 427 428 429 430 431 432 433 Last Page 429 of 709