SSPTRF implementation

[sarnath]

I was going through the SSPTR implementation at http://www.netlib.org/lapack/explore-html/ssptrf.f.html

I have a doubt in this piece of code.

Code: Select all

201: *
202: *              JMAX is the column-index of the largest off-diagonal
203: *              element in row IMAX, and ROWMAX is its absolute value
204: *
205:                ROWMAX = ZERO
206:                JMAX = IMAX
207:                KX = IMAX*( IMAX+1 ) / 2 + IMAX
208:                DO 20 J = IMAX + 1, K
209:                   IF( ABS( AP( KX ) ).GT.ROWMAX ) THEN
210:                      ROWMAX = ABS( AP( KX ) )
211:                      JMAX = J
212:                   END IF
213:                   KX = KX + J
214:    20          CONTINUE
215:                KPC = ( IMAX-1 )*IMAX / 2 + 1
216:                IF( IMAX.GT.1 ) THEN
217:                   JMAX = ISAMAX( IMAX-1, AP( KPC ), 1 )
218:                   ROWMAX = MAX( ROWMAX, ABS( AP( KPC+JMAX-1 ) ) )
219:                END IF


Everything is clear until line 214. It is as per the comment at line 202.
The entire IMAX row is scanned (except the diagonal element) to find what is the maximum value.

Starting from line 216, the "IMAX"th COLUMN is scanned fora larger value.

At the end of 219, "JMAX" stands for a ROW# (if IMAX was > 1) and ROWMAX could EITHER be on the IMAXth row or IMAXth column.

My point of confusion is that "JMAX" is updated readily inside the IF statement without actually looking at the result of the "MAX" statement below.

Can some1 please clear my doubts?

Thank you,
Best Regards,
Sarnath

[sarnath]

Can some1 answer to my post above please?

Thank you

[Julien Langou]

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 value
JMAX is the column-index of the largest off-diagonal element in row IMAX, and ROWMAX is its absolute value
So 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)

As you can see the value of JMAX becomes useless once we get ROWMAX.

--j

[sarnath]

Thanks for the detailed and quick answer.

I will go back to textbook and come back again..

As you can see, I am new to this whole LAPACK thing. Thanks!

Best REgards,
Sarnath

[sarnath]

Oh Yes! I got the point.. (but still a doubt remains..)

The "bouncing" on diagonal set me thinking on right lines...

It is a symmetric matrix...! And hence the IMAX "ROW" and IMAX "column" is also scanned. Perfect! Thank you so much!

That all sounds great except for the setting of "JMAX" -- which could potentially be wrong.

Can I ask an ignorant question at this point?

JMAX has no role in finding ROWMAX. It just records where "ROWMAX" was found.
So, I hope it is being used by the sections of code that follows it.

So, In such a case, setting of JMAX (inside the IF IMAX.GT.1) without checking the result of MAX statement that follows it ----- could potentially be a bug... No?

[Julien Langou]

Code: Select all

202: *              JMAX is the column-index of the largest off-diagonal
203: *              element in row IMAX, and ROWMAX is its absolute value
204: *
205:                ROWMAX = ZERO
206:                JMAX = IMAX
207:                KX = IMAX*( IMAX+1 ) / 2 + IMAX
208:                DO 20 J = IMAX + 1, K
209:                   IF( ABS( AP( KX ) ).GT.ROWMAX ) THEN
210:                      ROWMAX = ABS( AP( KX ) )
211:                      JMAX = J
212:                   END IF
213:                   KX = KX + J
214:    20          CONTINUE
215:                KPC = ( IMAX-1 )*IMAX / 2 + 1
216:                IF( IMAX.GT.1 ) THEN
217:                   JMAX = ISAMAX( IMAX-1, AP( KPC ), 1 )
218:                   ROWMAX = MAX( ROWMAX, ABS( AP( KPC+JMAX-1 ) ) )
219:                END IF


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.

Try to look at SSYTF2 (full format), that makes thing somewhat more readable ....

Code: Select all

219: *
220: *              JMAX is the column-index of the largest off-diagonal
221: *              element in row IMAX, and ROWMAX is its absolute value
222: *
223:                JMAX = IMAX + ISAMAX( K-IMAX, A( IMAX, IMAX+1 ), LDA )
224:                ROWMAX = ABS( A( IMAX, JMAX ) )
225:                IF( IMAX.GT.1 ) THEN
226:                   JMAX = ISAMAX( IMAX-1, A( 1, IMAX ), 1 )
227:                   ROWMAX = MAX( ROWMAX, ABS( A( JMAX, IMAX ) ) )
228:                END IF


lines 223 to 224: you scan from IMAX to K on the row IMAX
lines 225 to 228: you scam from 1 to IMAX-1 on the row IMAX but since you only store the UPPER part of A you need to look in the column

[sarnath]

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:

Code: Select all

216:                IF( IMAX.GT.1 ) THEN
217:                   JMAX = ISAMAX( IMAX-1, AP( KPC ), 1 )
218:                   ROWMAX = MAX( ROWMAX, ABS( AP( KPC+JMAX-1 ) ) )
219:                END IF


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.

[Julien Langou]

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

Code: Select all

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.

--julien

[sarnath]

Well, Interesting. I wonder why "JMAX" is being calculated at all - if it is not used anywhere :-)

I will go check the code. Meanwhile, THANKS a lot for your time!

[sarnath]

Julien,

Yes you are right! It is indeed a temporary variable. I was misled by the comment.

I just assumed that since the code is recording the position, it would need it down to swap the contents or whatever. I was wrong in assuming so.

Thanks,

Best REgards,
Sarnath