Skip to content

Conversation

@jip
Copy link
Contributor

@jip jip commented Jun 14, 2021

These typos were overlooked while working on the PR #581. I hope they are the last.

@codecov
Copy link

codecov bot commented Jun 14, 2021

Codecov Report

Merging #582 (2eeed4b) into master (34c8291) will not change coverage.
The diff coverage is n/a.

Impacted file tree graph

@@           Coverage Diff           @@
##           master     #582   +/-   ##
=======================================
  Coverage   82.37%   82.37%           
=======================================
  Files        1894     1894           
  Lines      190679   190679           
=======================================
  Hits       157065   157065           
  Misses      33614    33614           

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 34c8291...2eeed4b. Read the comment docs.

*> triangular system of linear equations op(A)*X = B, when A is a
*> triangular band matrix. The test ratio is the maximum over the
*> number of right hand sides of
*> norm(op(A)*X - B) / ( norm(op(A)) * norm(X) * EPS ),
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this is still incorrect. The expression norm(op(A)*X - B) / ( norm(op(A)) * norm(X) * EPS ) is a number, so it does not have right-hand sides. I suggest something like:

The test ratio is the maximum over
norm(b - op(A)*x) / ( norm(op(A)) * norm(x) * EPS )
for every column b of B,
where op(A) = A, A**T, or A**H, x is the solution vector, and EPS is the machine epsilon.

What do you think?

Copy link
Contributor Author

@jip jip Jun 15, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, this looks better. Two adjoining "where" sections (the first one - "for every column b of B" - is implicit, and the second one is explicit) could be merged:

The test ratio is the maximum over
  norm(b - op(A)*x) / ( norm(op(A)) * norm(x) * EPS ),
where op(A) = A, A**T, or A**H, b is the column of B, x is
the solution vector, and EPS is the machine epsilon.

but this is a matter of taste.

I believe this variant is incorrect, too. The original comment states "x and b are N by NRHS matrices". So will be the residual R := b - op(A)*x. The code computes taxicab-based norm || R(j) || for each j-th column R(j) from R. Then it takes a max(j) || R(j) || which means in toto a computing of kinda 1-norm, but based on taxicab-based norm |Re(z)| + |Im(z)|, not magnitude-based norm |z|. So, the phrase "the maximum over" is excessive.

And this variant may be improved further:

  1. Since x and b are matrices not vectors here, capital letters X and B are used.
  2. I'm using norm(op(A)*X - B) not norm(B - op(A)*X) since the first is the exact expression which is computed. I know they give identical results. I don't insist on this summands swapping here.

So, the result may look like this:

The test ratio is set to
  norm(op(A)*X - B) / ( norm(op(A)) * norm(X) * EPS ),
where op(A) = A, A**T, or A**H, B is the N by NRHS matrix of
right-hand sides, X is the N by NRHS matrix of solution vectors, and
EPS is the machine epsilon.

When "kinda 1-norms" in both norm(op(A)*X - B) and norm(X) will be replaced by "true 1-norm", then the following statement could be added: "The norm used is the 1-norm.".

I'm ready to update this PR with any of variants discussed above.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, this looks better. Two adjoining "where" sections (the first one - "for every column b of B" - is implicit, and the second one is explicit) could be merged:

The test ratio is the maximum over
  norm(b - op(A)*x) / ( norm(op(A)) * norm(x) * EPS ),
where op(A) = A, A**T, or A**H, b is the column of B, x is
the solution vector, and EPS is the machine epsilon.

but this is a matter of taste.

I would go with either this one or the version I proposed. Maybe you can modify it a bit:

The test ratio is the maximum over
  norm(b - op(A)*x) / ( ||op(A)||_1 * norm(x) * EPS ),
where op(A) = A, A**T, or A**H, b is the column of B, x is
the solution vector, and EPS is the machine epsilon.

As you pointed out, the norm operator is that "kinda 1-norm" for vectors of size N, in which |Re(z)| + |Im(z)| replaces the magnitude |z| for all complex scalars z. This is what the subroutine computes if I understand the algorithm correctly.

I believe this variant is incorrect, too. The original comment states "x and b are N by NRHS matrices". So will be the residual R := b - op(A)*x. The code computes taxicab-based norm || R(j) || for each j-th column R(j) from R. Then it takes a max(j) || R(j) || which means in toto a computing of kinda 1-norm, but based on taxicab-based norm |Re(z)| + |Im(z)|, not magnitude-based norm |z|. So, the phrase "the maximum over" is excessive.

And this variant may be improved further:

1. Since x and b are matrices not vectors here, capital letters X and B are used.

2. I'm using `norm(op(A)*X - B)` not `norm(B - op(A)*X)` since the first is the exact expression which is computed. I know they give identical results. I don't insist on this summands swapping here.

So, the result may look like this:

The test ratio is set to
  norm(op(A)*X - B) / ( norm(op(A)) * norm(X) * EPS ),
where op(A) = A, A**T, or A**H, B is the N by NRHS matrix of
right-hand sides, X is the N by NRHS matrix of solution vectors, and
EPS is the machine epsilon.

I don't agree with this proposal. The algorithm does not compute norm(X), where "capital" X is the original input matrix. Maybe I am missing something, please correct me if I am wrong.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

(Sorry for delayed answer)

I think you are right. I'll use the last variant you've proposed:

I would go with either this one or the version I proposed. Maybe you can modify it a bit:

The test ratio is the maximum over
  norm(b - op(A)*x) / ( ||op(A)||_1 * norm(x) * EPS ),
where op(A) = A, A**T, or A**H, b is the column of B, x is
the solution vector, and EPS is the machine epsilon.

I'll try either to update this PR or to create a new one.

@weslleyspereira
Copy link
Collaborator

I just noticed another thing that is not related to your changes, but maybe you can help me. Is there a bug here?

lapack/TESTING/LIN/ctbt02.f

Lines 231 to 235 in 2eeed4b

IF( XNORM.LE.ZERO ) THEN
RESID = ONE / EPS
ELSE
RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) / EPS )
END IF

I think so. Suppose that the last column of X is full of zeros. Then, XNORM = ZERO in the last iteration ( J = NRHS ).
Then, the subroutine computes RESID = ONE / EPS. There are many files with the same pattern.

I can propose a PR to solve this (if this is actually a bug).

@langou langou merged commit 9e647a9 into Reference-LAPACK:master Aug 30, 2021
@jip
Copy link
Contributor Author

jip commented Sep 6, 2021

I just noticed another thing that is not related to your changes, but maybe you can help me. Is there a bug here?
[...]
I can propose a PR to solve this (if this is actually a bug).

I don't see any problems here. If any column of X is full of zeros, then the RESID takes a maximal value, and the rest columns of X doesn't matter. IMHO.

jip added a commit to jip/lapack that referenced this pull request Sep 6, 2021
@jip jip mentioned this pull request Sep 6, 2021
@julielangou julielangou added this to the LAPACK 3.10.1 milestone Nov 12, 2022
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants