The LAPACK code is correct as far as I can tell. One point of confusion might be the algorithm in itself. Do you understand what the pivoting strategy is?The routine is using the Bunch-Kaufman pivoting strategy from their 1977 paper.This strategy is also called 'partial pivoting' (in Higham for example).I find all this best explained in Higham's book (Accuracy and Stability of Numerical Algorithms).In the second edition, this is Section 11.1.2, Algorithm 11.2.I believe this is also explained in Golub and van Loan.

I am not going to rewrite the books. Here are some quick hints.What you are looking for with these loops at step k are:ABSAKK :: the absolute value of the diagonal element in position (K,K)IMAX :: the row-index of the largest off-diagonal element in column K, and COLMAX is its absolute valueJMAX is the column-index of the largest off-diagonal element in row IMAX, and ROWMAX is its absolute valueSo note that getting ROWMAX is tricky since you store only the lower part of A, so you need to "bounce" on the diagonal, morever here you are in packed storage.

Once you have ABSAKK, COLMAX and ROWMAX, you have some criterion to check which you can read indeed in the routine itself.Depending on this criteria, you will choose * either not to permute (i.e. keep A(K,K) as the pivot)* or pick a 2-b-2 pivot block, namely A( [ K IMAX ] , A[ K IMAX] ), (so you will need to permute it to bring in position K:K+1,K:K+1)* or pick the 1-by-1 pivot A([ IMAX, IMAX ]) (so you need to permute as well)

In these previous line, consider JMAX as a temporay variable. It is used several times for different things.The code is really hard to get sense of it ....It scans row IMAX of the symmetric matrix AP stored in packed storage from K to N and record the maximum absolute value in ROWMAX.There are two steps.

I perfectly understand how the ROW is being scanned. I have no problems with that.

For sake of clarify, let us visualize this situation. Assume A is the Upper Triangle of the matrix along with the diagonal. The lower-triangle is not required since we can get that info from the upper triangle.

Thus to read the complete "IMAX"th row of the symmetric matrix, one needs to scan the "IMAX" column and "IMAX" Row of the Triangular matrix in hand.Since the Matrix is column-packed, reading a ROW is difficult than reading a column.

The code I had originally quoted scans "IMAX"th row first and finds out the Position and Value of the largest off-diagonal element.Once this is done, it then moves on scans the "IMAX"th column and finds out the global maximum for the IMAX row.

However, I feel that the code below, where the IMAXth column (IMAX row in lower triangle)is scanned , has a bug:

In the code above, JMAX is being updated without knowing which ROWMAX got selected in the MAX function below it. So, in the case that the "MAX" function below retained the original ROWMAX (got by scanning the upper triangle), JMAX will be pointing to something that is NOT actually the MAX.

Thats my point. From a software angle, It looks like a bug -- it is non-deterministic. The way the algorithm behaves will purely dependent on the content of the Matrix.

In the code above, JMAX is being updated without knowing which ROWMAX got selected in the MAX function below it.So, in the case that the "MAX" function below retained the original ROWMAX (got by scanning the upper triangle), JMAX will be pointing to something that is NOT actually the MAX.

OK, I think I got you and you are right. It is indeed written in the comment that

JMAX is the column-index of the largest off-diagonal element in row IMAX,

Well, OK, that is not true :)check the code carefully. JMAX is not used anymore. There is no need to update JMAX to the correct value. I told you:

In these previous line, consider JMAX as a temporay variable.

However you are write: the comment is not correct and needs to be updated. The code is correct though. What we care about is ROWMAX in the pivoting strategy. JMAX is irrelevant once ROWMAX is determined.