This should go into lapack/{gonum,testlapack}, it's not a BLAS routine.

The test should follow the lapack/testlapack style, for example dpbtf2.go is kind of related and relatively recent. dtrtrs.go is not so related but also quite recent and shows the style on which we've recently settled in lapack.

Got confused since the page says "BLAS Level 2 Algorithm". will fix.

Thanks for the suggestions!

Got confused since the page says "BLAS Level 2 Algorithm". will fix.

No worries. Some Lapack routines come in 2 variants, a simpler one that uses BLAS Level 2 and operates on columns/rows, and a more complex one that uses BLAS Level 3, operates on blocks and should be more performant. The comment here means that this routine uses BLAS Level 2.

Let's see what Github does with these review comments done when the code was still in the blas directory.

Ok, it shows them as outdated but they are of course all still valid.

I'm at a loss as to how to construct the permutation matrices from pivot indices. Cholesky factorization function does not have a similar implementation.

Why do you need to construct a permutation matrix?

@kortschak I'm way out of my element here, but I'd thought that a good way to test Dgetc2 was to reconstruct the matrix A from the results of the LU factorization A = P * L * U * Q where P and Q are permutation matrices, L is lower triangular with unit diagonal elements and U is upper triangular.

What my test is doing thus far is checking if the permutation indices yield unique combinations, which I'm not even sure is an expected result. It is failing for now.

I'd thought that a good way to test Dgetc2 was to reconstruct the matrix A from the results of the LU factorization A = P * L * U * Q where P and Q are permutation matrices, L is lower triangular with unit diagonal elements and U is upper triangular.

Yes, reconstructing A is the way to test Dgetc2 but you don't need to form P and Q explicitly and multiply, you can use ipiv and jpiv to swap rows/columns to achieve the same effect.

What my test is doing thus far is checking if the permutation indices yield unique combinations, which I'm not even sure is an expected result. It is failing for now.

Testing the content of ipiv and jpiv is easy and will not hurt. I would fill them with negative numbers before calling Dgetc2, then make their copies, sort the copies and then check that ipivCopy[i] == i for all i. Then I would use ipiv, jpiv to reconstruct A.

Testing the content of ipiv and jpiv is easy and will not hurt. I would fill them with negative numbers before calling Dgetc2, then make their copies, sort the copies and then check that ipivCopy[i] == i for all i. Then I would use ipiv, jpiv to reconstruct A.

I got around to implementing this and letting it fail. I've read the code and compared it with reference but nothing seems out of place. Any ideas on what I should be looking out for?

For the time being ipiv and jpiv are being filled (no negative numbers) with repeats i.e. ipiv:[1 1], jpiv:[0 1] or :

ipiv:[2 9 9 9 10 10 12 13 14 14 14 15 16 16 16 17 18 18 19 19]
jpiv:[1 5 9 12 13 13 14 14 14 15 16 16 16 16 17 18 18 19 19 19]

That's my fault, I didn't pay enough attention to how the permutations are represented in ipiv and jpiv. Now I see there is not much that can be tested except that all elements of ipiv, jpiv are set, in your case that no element is -1. After that you should focus on checking that A can be reconstructed from the factorization.

you can use ipiv and jpiv to swap rows/columns to achieve the same effect.

There's definitely more to reconstructing A than swapping columns and rows. I did just that in my test but noticed the lines

for j := i + 1; j < n; j++ {
	a[j*lda+i] /= a[i*lda+i]
}
bi.Dger(n-i-1, n-i-1, -1.0, a[(i+1)*lda+i:], lda, a[i*lda+i+1:], 1, a[(i+1)*lda+i+1:], lda)

modify a in every case. There's gotta be a better way of testing this than programming inverseDgetc2.

There's definitely more to reconstructing A than swapping columns and rows.

Yes, Dgetc2 stores the L and U triangular factors in-place into a (and P and Q into ipiv and jpiv).

Yes, Dgetc2 stores the L and U triangular factors in-place into a (and P and Q into ipiv and jpiv).

First: thanks for the hint, I was not aware of permutation matrices. I'm feeling good about where the test is at, though I'm not sure what is causing it to fail. The first column remains zero for a 2x2 matrix. If it were not much to ask, I could use a tidbit of help on this one! It's been less than a week since I've learned what a LAPACK is! Figured it out: was missing an = on the greater-or-equal. Comparing with expected output it seems Dgetc2 works.

@kortschak This branch is ready to merge. Dgetc2 passes all tests.

Codecov Report

Merging #1655 (13cc21c) into master (2b838f4) will decrease coverage by 0.00%.
The diff coverage is 63.46%.

@@            Coverage Diff             @@
##           master    #1655      +/-   ##
==========================================
- Coverage   73.22%   73.21%   -0.01%     
==========================================
  Files         504      505       +1     
  Lines       59811    59867      +56     
==========================================
+ Hits        43794    43831      +37     
- Misses      13537    13553      +16     
- Partials     2480     2483       +3     

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 2b838f4...13cc21c. Read the comment docs.

Changes to tests done as per your suggestion. Removed those cases generated by the reference since they may fail depending on implementation details.

I found it CONTRIBUTING.md:

you will need to be added to the CONTRIBUTORS and AUTHORS file; after your contribution has been accepted you will be asked to open a pull request adding yourself to them.

@soypat , can you submit a new PR where you add yourself to those two files?

Make the commit message "A+C: add Patricio Whittingslow" in the other PR.