From: "Vitaly N. Golovach"
<vitaly.n.golovach@Domain.Removed<mailto:vitaly.n.golovach@Domain.Removed>>
Date: December 13, 2011 7:25:22 AM MST
To: "Langou, Julien"
<Julien.Langou@Domain.Removed<mailto:Julien.Langou@Domain.Removed>>
Subject: Re: [Lapack] mistake in BLAS/scabs1.f
Hi Julien,
Thank you for the detailed answer! I enjoyed reading the references that you
sent me. I agree, it makes sense to define a norm like  Re( z )  +  Im ( z )
, because it is so easy to compute and that such a norm can be useful for
certain routines where it is not crucial to have the exact Euclidean norm. I
imagine that ZGETF2 is such a routine, although I am not familiar with the
algorithm (partial pivoting with row interchanges), so it is difficult to judge
whether replacing  Re( z )  +  Im ( z )  by a properly coded (  Re( z )
^2 +  Im ( z ) ^2 )^(1/2) should improve the results. The two norms differ
only by a factor of 2^(1/2), at most.
However, "absolute value" is not a good name for this measure of length when
you refer to complex numbers. Most people will only read the description of the
routine IZAMAX and trust that it "finds the index of element having max.
absolute value". If you have the power to do it, you might want to consider
adding another sentence to the description saying what is actually meant by
"absolute value". I also liked the term "pseudo 1norm" used at the end of the
paragraph in Nick Higham's book, p.500. But making a too drastic change is also
not good, because the people in your field are, probably, used to seeing
"absolute value" in that place.
A simple example where IZAMAX would give different results depending on which
norm is used is as follows. Consider an array consisting of two complex numbers
z_1 = (6,0) and z_2 = (3,4). The Euclidean version would compare 6>5, whereas
the "pseudo 1norm" would compare 6<7.
Thanks again for your explanation and wish you good luck!
Best,
Vitaly
*******************************
Vitaly N. Golovach
LPMMC, Maison des Magist?res CNRS
25 avenue des Martyrs, BP166
38042 Grenoble Cedex, FRANCE
Tel.: +33 4 76 88 79 82
Fax.: +33 4 76 88 79 83
Email: vitaly.n.golovach@Domain.Removed<mailto:vitaly.n.golovach@Domain.Removed>
*******************************
On Mon, Dec 12, 2011 at 9:39 PM, Langou, Julien
<Julien.Langou@Domain.Removed<mailto:Julien.Langou@Domain.Removed>> wrote:
Hi Vitaly,
The definition of absolute value, abs( z ), of a complex number, z, in the
BLAS has always been
abs( z ) =  Re( z )  +  Im ( z ) .
( Always means since the seventies. )
I guess the reason is that a reliable computation of the modulus:
 z  = (  Re( z ) ^2 +  Im ( z ) ^2 )^(1/2)
is complicated (see for example dlapy2 to do this operation without
unneccessary overflows or underflows)
and since it is complicate, it consumes time.
There is a nice paragraph in Nick Higham's book "Accuracy and Stability of
Numerical Algorithms", 2nd edition p.500, Sec 27.8.
If you look at the original Level 1 BLAS paper from 1979 form Lawson, Hanson,
Kincaid, and Krogh.
( available from Kincaid webpage
http://www.cs.utexas.edu/users/kincaid/blas.pdf )
( page 311) you would read some justification on this choice.
So
* a subroutine like ICAMAX would look like the largest element in absolute
value with the definition:
abs( z ) =  Re( z )  +  Im ( z ) .
* a subroutine like ICAMAX does not compute the infinite norm of a complex
vector.
* a subroutine like ICASUM does not compute the 1norm of a complex vector.
Such routine do exist in the new BLAS standard.
http://www.netlib.org/blas/blastforum/chapter2.pdf
page 19, table 2.1, under NORM.
Anyway, you do not think "absolute value" is appropriate to reference:
abs( z ) =  Re( z )  +  Im ( z ) .
Fair enough. Suggestions? "taxi modulus"?
One thing I have always wondered is whether GETRF which uses GETF2 which uses
ICAMAX is more or less stable than one who would use  z  = (  Re( z ) ^2 +
 Im ( z ) ^2 )^(1/2) ....
Cheers,
Julien.
On Dec 12, 2011, at 8:40 AM, Vitaly N. Golovach wrote:
Dear LPACK development team,
There is a mistake in the function SCABS1 in BLAS, see file BLAS/scabs1.f.
The absolute value of a complex number is not the sum of the absolute values of
its real and imaginary parts.
Either the documentation is wrong or the code is wrong.
The function SCABS1 is used only in BLAS/icamax.f. Actually, it is also used in
BLAS/caxpy.f, but in a trivial way, where it cannot cause errors.
Best,
Vitaly
*******************************
Vitaly N. Golovach
LPMMC, Maison des Magist?res CNRS
25 avenue des Martyrs, BP166
38042 Grenoble Cedex, FRANCE
Tel.: +33 4 76 88 79 82<tel:%2B33%204%2076%2088%2079%2082>
Fax.: +33 4 76 88 79 83<tel:%2B33%204%2076%2088%2079%2083>
Email: vitaly.n.golovach@Domain.Removed<mailto:vitaly.n.golovach@Domain.Removed>
*******************************
REAL FUNCTION SCABS1(Z)
* .. Scalar Arguments ..
COMPLEX Z
* ..
*
* Purpose
* =======
*
* SCABS1 computes absolute value of a complex number
*
* =====================================================================
*
* .. Intrinsic Functions ..
INTRINSIC ABS,AIMAG,REAL
* ..
SCABS1 = ABS(REAL(Z)) + ABS(AIMAG(Z))
RETURN
END
_______________________________________________
Lapack mailing list
Lapack@Domain.Removed<mailto:Lapack@Domain.Removed>
http://lists.eecs.utk.edu/mailman/listinfo/lapack
 next part 
An HTML attachment was scrubbed...
URL:
http://lists.eecs.utk.edu/mailman/private/lapack/attachments/20111213/e28ce26e/attachment0001.html
