Beefy Boxes and Bandwidth Generously Provided by pair Networks
Perl: the Markov chain saw
 
PerlMonks  

comment on

( [id://3333]=superdoc: print w/replies, xml ) Need Help??
The problem is that the 'x' operator is threading over the internal (real-valued) matrices rather than computing the correct complex products.

That, in turn, is because you are using the pdl() constructor around your matrices. Thus they are not blessed into the PDL::Complex class, but only into the PDL class. The pdl() constructor is pretty permissive and turns just about anything into a PDL, but can be confused on occasion. In this case, pdl() recognizes a collection of nested list refs, some of which contain PDLs (or at least PDL subclasses). It is trying to do the Right Thing by putting your inner PDL dims in the deepest (fastest-running) indices of the final constructed PDL. That looks almost right, as PDL::Complex is implemented with an extra 0th dim that runs over {real|imaginary}.

Because your matrices are PDLs and not PDL::Complexes, they are being multiplied differently than you expect: MatrixM is a 2x2x2, which is treated as a collection of two 2x2 matrices in (column, row, thread) order rather than a single complex matrix in ({r|i},column,row) order. MatrixX is a 2x2, which is being treated as a single matrix in (column, row) order rather than a complex row vector in ({r|i},column) order. You have a couple of ways around the complex issue:

  • (1) Use the correct constructor, "new PDL::Complex", rather than the incorrect one "pdl";
  • (2) Manually bless your PDLs back into PDL::Complex with 'bless($matrixM,"PDL::Complex");' (this is what PDL::Complex::new does anyway); or
  • (3) Assemble your matrices from independently constructed PDLs, as in '$matrix=pdl([5,4])+pdl([1,2])*i;'

I think that you are about to hit another problem, which is that $matrixX is a row vector rather than (as you probably think, given that you put it on the right hand side of a matrix) a column vector.

The row vector/column vector problem arises from the need to choose between (row,column) indexing order, which renders incorrectly on screen, or (column, row) indexing order, which is backwards but renders correctly on screen. PDL uses (column,row) order by default.

To make a row vector, you say '$row = pdl(1,2,3)' or (if you prefer) '$row = pdl [1,2,3]'.

To make a column vector you must insert a trivial 0th dim (of size 1), which you can do thusly:

  • (1) by nesting: $col = pdl([1],[2],[3])
  • (2) by transposition: $col = pdl(1,2,3)->transpose
  • (3) by dummy index: $col = pdl(1,2,3)->(*1)

Incidentally, there's one more index pitfall to watch for (though it isn't yet in your example): sometimes you may want to thread across row vectors (for example, if you want to multiply a collection of N row vectors by a collection of N matrices that are stored in a CxRxN PDL). Then you should remember to insert the dummy dimension explicitly. For example, if your vectors are (1,2), (4,5), and (7,8) you should say:

$matrices = sequence(2,2,3); $threadable_rows = pdl([1,2],[4,5],[7,8])->(:,*1,:); print $threadable_rows x $matrices;
which will yield a different answer than
$matrices = sequence(2,2,3); $wrong_rows = pdl([1,2],[4,5],[7,8]); print $wrong_rows x $matrices;

In reply to Re: PDL::Complex question by drzowie
in thread PDL::Complex question by gmacfadden

Title:
Use:  <p> text here (a paragraph) </p>
and:  <code> code here </code>
to format your post; it's "PerlMonks-approved HTML":



  • Are you posting in the right place? Check out Where do I post X? to know for sure.
  • Posts may use any of the Perl Monks Approved HTML tags. Currently these include the following:
    <code> <a> <b> <big> <blockquote> <br /> <dd> <dl> <dt> <em> <font> <h1> <h2> <h3> <h4> <h5> <h6> <hr /> <i> <li> <nbsp> <ol> <p> <small> <strike> <strong> <sub> <sup> <table> <td> <th> <tr> <tt> <u> <ul>
  • Snippets of code should be wrapped in <code> tags not <pre> tags. In fact, <pre> tags should generally be avoided. If they must be used, extreme care should be taken to ensure that their contents do not have long lines (<70 chars), in order to prevent horizontal scrolling (and possible janitor intervention).
  • Want more info? How to link or How to display code and escape characters are good places to start.
Log In?
Username:
Password:

What's my password?
Create A New User
Domain Nodelet?
Chatterbox?
and the web crawler heard nothing...

How do I use this?Last hourOther CB clients
Other Users?
Others sharing their wisdom with the Monastery: (4)
As of 2024-03-29 12:23 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?

    No recent polls found