Probably a problem with single threaded dpotri

[andreasnoack]

At the Julia user list a user reported a problem that appears to be related to the single threaded dpotri. Consider the following program

program test

    implicit none
    integer*8, parameter :: n = 1000
    integer :: i, j
    integer*8 :: info
    double precision :: vX(n), mX(n,n)
    real :: rand

    call srand(123)
    do i = 1, n
        vX(i) = rand()
    end do

    do i = 1, n
        do j = 1, n
            mX(i,j) = exp(-0.5*(vX(i) - vX(j))**2)
            if ( i == j ) mX(i,i) = mX(i,i) + 0.5
        end do
    end do

    call dpotrf('U', n, mX, n, info)
    write(*,*) info
    call dpotri('U', n, mX, n, info)
    write(*,*) info
    call dpotrf('U', n, mX, n, info)
    write(*,*) info

end program test

With this I get

Andreass-MacBook-Pro:Desktop andreasnoackjensen$ gfortran test.f90 ~/julia/deps/openblas-v0.2.10.rc2/libopenblas.a
Andreass-MacBook-Pro:Desktop andreasnoackjensen$ ./a.out 
                    0
                    0
                    0
Andreass-MacBook-Pro:Desktop andreasnoackjensen$ export OPENBLAS_NUM_THREADS=1
Andreass-MacBook-Pro:Desktop andreasnoackjensen$ ./a.out 
                    0
                    0
                    4
Andreass-MacBook-Pro:Desktop andreasnoackjensen$ export OPENBLAS_NUM_THREADS=2
Andreass-MacBook-Pro:Desktop andreasnoackjensen$ ./a.out 
                    0
                    0
                    0
Andreass-MacBook-Pro:Desktop andreasnoackjensen$ gfortran test.f90 -framework Accelerate
Andreass-MacBook-Pro:Desktop andreasnoackjensen$ ./a.out 
                    0
                    0
                    0

My machine is a mid 2009 macbook pro and OpenBLAS is compiled with 64 bit integer support and automatic core detection.

[wernsaar]

Hi,

I am able to reproduce the error and will try to provide a fix soon.

Thanks

Werner

On 10.07.2014 10:02, Andreas Noack Jensen wrote:

At the Julia user list a user reported a problem that appears to be related to the single threaded dpotri. Consider the following program

program test

  implicit none
  integer*8, parameter :: n = 1000
  integer :: i, j
  integer*8 :: info
  double precision :: vX(n), mX(n,n)
  real :: rand

  call srand(123)
  do i = 1, n
      vX(i) = rand()
  end do

  do i = 1, n
      do j = 1, n
          mX(i,j) = exp(-0.5*(vX(i) - vX(j))**2)
          if ( i == j ) mX(i,i) = mX(i,i) + 0.5
      end do
  end do

  call dpotrf('U', n, mX, n, info)
  write(*,*) info
  call dpotri('U', n, mX, n, info)
  write(*,*) info
  call dpotrf('U', n, mX, n, info)
  write(*,*) info

end program test

With this I get

Andreass-MacBook-Pro:Desktop andreasnoackjensen$ gfortran test.f90 ~/julia/deps/openblas-v0.2.10.rc2/libopenblas.a
Andreass-MacBook-Pro:Desktop andreasnoackjensen$ ./a.out
                     0
                     0
                     0
Andreass-MacBook-Pro:Desktop andreasnoackjensen$ export OPENBLAS_NUM_THREADS=1
Andreass-MacBook-Pro:Desktop andreasnoackjensen$ ./a.out
                     0
                     0
                     4
Andreass-MacBook-Pro:Desktop andreasnoackjensen$ export OPENBLAS_NUM_THREADS=2
Andreass-MacBook-Pro:Desktop andreasnoackjensen$ ./a.out
                     0
                     0
                     0
Andreass-MacBook-Pro:Desktop andreasnoackjensen$ gfortran test.f90 -framework Accelerate
Andreass-MacBook-Pro:Desktop andreasnoackjensen$ ./a.out
                     0
                     0
                     0

My machine is a mid 2009 macbook pro and OpenBLAS is compiled with 64 bit integer support and automatic core detection.

---

Reply to this email directly or view it on GitHub:
#410

[wernsaar]

@andreasnoackjensen , as a fast workaround, I disabled the optimized potri functions. You can test this by using the code form my repository: https://github.com/wernsaar/OpenBLAS.git

Best regards
Werner

[andreasnoack]

Just tried it and it solves the problem. Thanks. I had OpenBLAS 0.2.8 libraries and tried them as well and they don't have the problem.

[wernsaar]

Hi,

after v0.2.8, we updated Lapack to version 3.5.0, which provided a lot
of problems.
Thanks to all, that help us to find bugs.

Best regards

Werner

On 10.07.2014 15:24, Andreas Noack Jensen wrote:

Just tried it and it solves the problem. Thanks. I had OpenBLAS 0.2.8 libraries and tried them as well and they don't have the problem.

---

Reply to this email directly or view it on GitHub:
#410 (comment)

[ViralBShah]

@xianyi Julia 0.3-rc1 is now released. Could we have the next openblas release candidate with this fix, and a few of the other recent fixes that all look like they are good to have? We will include this in our next release candidate. If no more bugs are found, perhaps we can release julia 0.3 and openblas 0.2.10 simultaneously. Thanks.

[xianyi]

@ViralBShah , we will release 0.2.10 version today.

In 3 or 4 weeks, we will release 0.2.11 version. I think openblas 0.2.11 can match julia 0.3.

[ViralBShah]

Ok - that would work too, if required.

[wernsaar]

I checked the performance of lapack native ?potri with OpenBLAS against MKL. About 85% of the cycles are spent in the gemm kernel, OpenBLAS is faster than MKL. There is no need, to provide custom ?potri functions.
I close this issue now.

Werner