/*
* Apologies for the layout but this following compiles as a Chapel
* program, just in case you want to mess with it.
*
* In an SVD decomposition, there are lots of Planar Rotations on
* adjacent Columns of two 2-dimensional (2D) matrices, each of a
* size of order N*N (say). These operations do not parallelize, as
* the Chapel column-major storage means the elements of a column
* are not adjacent in memory, causing cache hit problems when the
* column is being accessed or rewritten in memory.
*
* At a cost of storing the Givens planar rotations in two (2) complex
* 1-dimensional (1D) arrays, each of size N, the parameter 'g' noted
* below, these operations can be parallelized as also seen there. The
* procedure 'qr' will be applied to both the original 2D- matrices,
* each of which will have a different planar rotation array full of
* complex numbers.
*
* Questions:
*
* a) Are there better ways to pull an array's domain apart? I
* not only need the upper and lower limits, but also the
* dimensions (as is) and then a slice thereof, i.e. the
* slice being the dimensions but shrunk by 1 at the start.
*
* b) The recursive inner loop looks like it will cause grief
* to an optimizer. Is there a smarter way to write it?
*
* c) Is there an optimal way of tupelizing the components of a
* complex number? See my grubby attempt!
*
* In the generic routine 'qr', if R is real(w), then C is complex(2*w)
*
* The planar rotation is called 'g' because it is a Given's rotation.
*/
inline proc complexAsTuple(type T, t : ?C) : (T, T) return (t.re, t.im);
proc qr(ref z : [?zD] ?R, ref g : [?gD] ?C)
{
const (allRows, _) = zD.dims(); // only want the rows
const (j, k) = (gD.first, gD.last); // pull gD apart
const Z = j..k; // put gD back together (sort of)
const G = j+1 .. k; // ignore the first element
forall r in allRows do
{
ref m = z[r, Z]; // grab the whole row
var t = z[r, j]; // grab the first element
for i in G do // highly recursive
{
(t, m[i - 1]) = complexAsTuple(R, g[i] * (m[i],
t):C);
}
m[k] = t;
writeln(t);
}
}
proc main
{
var a : [1..10, 1..8] real(64);
var r : [2..7] complex(128);
a = 1.0;
r = (0.6, 0.8):complex;
qr(a, r);
}
Regards - Damian
Pacific Engineering Systems International, 277-279 Broadway, Glebe NSW 2037
Ph:+61-2-8571-0847 .. Fx:+61-2-9692-9623 | unsolicited email not wanted here
Views & opinions here are mine and not those of any past or present employer
_______________________________________________
Chapel-developers mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/chapel-developers
