LAPACK Archives

[Lapack] xLARFP and scaling.

Dear LAPACKers:

Also while attempting to upgrade to LAPACK 3.2, more eigs tests began to fail, 
but this time due to differences from DGEQR2 used by the ARPACK routine DSEUPD. 
 It may be that there is something peculiar going on in ARPACK (which I'm 
discussing with Lehoucq and Sorensen) but this encouraged me to look closely at 
xLARFP.

We're somewhat troubled by the fact that xLARFP can produce reflectors I - tau 
* v * v' with tau being arbitrarily small, and the elements in v being 
arbitrarily large.  The reflectors produced by xLARFG have 1 <= tau <= 2 and 
each element in v bounded above by 1 in absolute value, both of which seem to 
be desirable properties.

Here is a very specific example:

Consider x = [1e100; 5e-62].  Calling DLARFP to construct a reflector gives tau 
\approx 1.48e-323!  Applying that reflector to y = [1; 1] gives Hy \approx [1; 
-1.3175] although the reflector should reflect over something quite nearly 
parallel to the x-axis.  By contrast, DLARFG constructs a reflector with tau 
\approx 2 and then maps y to Hy = [-1; 1].  Also, making x(2) smaller (but 
nonzero) will cause DLARFG to continue to produce a reflector H such that Hy 
(as applied by DLARF) will be [-1;1], whereas DLARFP will stop producing such a 
reflector once tau underflows.  Note that the attached DLARFP2 (described 
below) will produce reflector that maps y to [1; -1] when constructed with x as 
above, and will continue to do so as x(2) -> 0.

We don't particularly want to patch LAPACK for our usage, but again, we are a 
bit concerned about the use of xLARFP underneath xGEQR2.  Our customers very 
frequently have poorly scaled data and will probably run into examples like 
that shown above.  At this point, we are inclined to replace usages of xLARFP 
with xLARFG.

As for general usage within LAPACK, there seem to be a few ways to deal with 
this behavior:
1.  Don't use xLARFP underneath xGEQR2 and introduce a new suite of routines 
which form QR using xLARFP.

2.  Change xLARFP to construct tau and v with certain bounds.  One possible 
approach that occurred to me is the following:
When alpha > 0, DLARFP computes tau = sin^2(t)/(1 + cos(t)) where t is the 
angle between the vector to be reflected and the target.  Here, t \in [0, pi/2) 
(as alpha > 0) and it is apparent that tau may become arbitrarily small.  One 
could factor sin(t) out of tau and instead set tau = 1/(1+cos(t)) which is now 
bounded (0.5 <= tau <= 1) and one could scale v such that v(1) = sin(t) instead 
of 1.  This scaling then makes v such that |v_k| <= 2.  Applying the reflectors 
then necessitates setting v(1) to sin(t) instead of 1 prior to calling DLARF.  
This requires either storing sin(t) and carrying it around, or attempting to 
recover it from tau as v(1) = sin(t) = sqrt(2*tau - 1)/tau.  While the latter 
does not require extra memory, it looks sort of dodgy thanks to a subtraction 
and a square root.  I'm attaching a version of DLARFP which has been modified 
to take the above approach and which stores v(1) in workspace.  I confess that 
I haven't thought deeply about how this generalizes to complex.

3.  Do nothing except please correct the comments in xLARFP which state that 1 
<= tau <= 2.

Any advice you can provide is welcome.

Thanks.
Pat.


-------------- next part --------------
A non-text attachment was scrubbed...
Name: dlarfp2.f
Type: application/octet-stream
Size: 4688 bytes
Desc: dlarfp2.f
Url : 
http://lists.eecs.utk.edu/mailman/private/lapack/attachments/20090729/cce8f52f/dlarfp2.obj
 

<Prev in Thread] Current Thread [Next in Thread>


For additional information you may use the LAPACK/ScaLAPACK Forum.
Or one of the mailing lists, or