replace {s,d}laqrb with {s,d}lahqr

[thrasibule]

Issue #3 was partially solved with pull request #2, by bringing back an older version of dlahqr.
However this only worked with reference blas, but with atlas or openblas, the issue persists as evidenced by #7, #16 and #17. This is another attempt at fixing it using a strategy suggested here by Marco Caliari.

This is an attempt to close #7 ,close #16 and close #17.

[sylvestre]

Impressive, could you also update the CHANGES to describe these changes? Thanks

[thrasibule]

added a description in the CHANGES file. All the credit for this idea goes to Marco Caliari, I just put the pull request together...

[caliarim]

Dear maintainers,

thanks for including my proposal in 3.3.0. In the old dlaqrb.f there was this piece of code

 if( wantt ) then
     i1 = 1
     i2 = n
     do 8 i=1,i2-2
        h(i1+i+1,i) = zero

8 continue
else
do 9 i=1, ihi-ilo-1
h(ilo+i+1,ilo+i-1) = zero
9 continue
end if

which is not present in dlahqr. Since dlaqrb.f was always called with wantt = .true., only the first branch is used. I don't know whether it is really necessary to zero out h. But it can be achived in dneigh.f by

do 5 j = 1, n-2
bounds(j) = zero
workl(j+2+n*(j-1)) = zero
5 continue
bounds(n-1) = zero
bounds(n) = one
call dlahqr(.true., .true., n, 1, n, workl, n, ritzr, ritzi, 1,
& 1, bounds, 1, ierr)

[thrasibule]

@caliarim. If I understand the code correclty, this sets the entries on the second diagonal to be 0. h(i+2,i)=0 for 1<=i<=n-2. don't think this is needed. With the line:

call dlacpy ('All', n, n, h, ldh, workl, n)

h is copied into workl which is an upper hessenberg matrix so all entries below the subdiagonal should be 0. What do you think?

[caliarim]

Look, for instance, at line 552 of dnaupd2. h(3,1) is used as storage.

[thrasibule]

Well spotted! However, in dlahqr code, I see:

*     ==== clear out the trash ====
      DO 10 J = ILO, IHI - 3
         H( J+2, J ) = ZERO
         H( J+3, J ) = ZERO
   10 CONTINUE

So seems that it's already taken care of. Do you agree?

[caliarim]

Ah, right! This piece of code is not in dlahqr 2.0, and therefore it was inserted into dlaqrb. I agree with you that it is not necessary to zero out workl before calling dlahqr. To be 100% sure, I initialized workl(j+2+n_(j-1)) with (-1.0d0)__j_1000.d0 and arpack worked without problems.

[thrasibule]

Cool! I'm glad this is sorted out.