testing_zgetrf_gpu_f.f90

Go to the documentation of this file.

!
!   -- MAGMA (version 1.2.0) --
!      Univ. of Tennessee, Knoxville
!      Univ. of California, Berkeley
!      Univ. of Colorado, Denver
!      May 2012
!
!  @precisions normal z -> c d s
!
      program testing_zgetrf_gpu_f

      use magma

      external cublas_init, cublas_set_matrix, cublas_get_matrix
      external cublas_shutdown, cublas_alloc
      external zlange, zgemm, zgesv, dlamch

      double precision zlange, dlamch
      integer cublas_alloc

      double precision              :: rnumber(2), anorm, bnorm, rnorm, xnorm 
      double precision, allocatable :: work(:)
      complex*16, allocatable       :: h_a(:), h_b(:), h_x(:)
      magma_devptr_t                :: devptra, devptrb
      integer,    allocatable       :: ipiv(:)

      complex*16                    :: zone, mzone
      integer                       :: i, n, info, stat, lda
      integer                       :: size_of_elt, nrhs
      real(kind=8)                  :: flops, t, tstart, tend

      parameter( nrhs = 1, zone = 1., mzone = -1. )
      
      call cublas_init()

      n   = 2048
      lda  = n
      ldda = ((n+31)/32)*32
      size_of_elt = sizeof_complex_16
 
!------ Allocate CPU memory
      allocate(h_a(lda*n))
      allocate(h_b(lda*nrhs))
      allocate(h_x(lda*nrhs))
      allocate(work(n))
      allocate(ipiv(n))

!------ Allocate GPU memory
      stat = cublas_alloc(ldda*n, size_of_elt, devptra)
      if (stat .ne. 0) then
         write(*,*) "device memory allocation failed"
         stop
      endif

      stat = cublas_alloc(ldda*nrhs, size_of_elt, devptrb)
      if (stat .ne. 0) then
         write(*,*) "device memory allocation failed"
         stop
      endif

!---- Initializa the matrix
      do i=1,lda*n
         call random_number(rnumber)
         h_a(i) = rnumber(1)
      end do

      do i=1,lda*nrhs
        call random_number(rnumber)
        h_b(i) = rnumber(1)
      end do
      h_x(:) = h_b(:)

!---- devPtrA = h_A
      call cublas_set_matrix(n, n, size_of_elt, h_a, lda, devptra, ldda)

!---- devPtrB = h_B
      call cublas_set_matrix(n, nrhs, size_of_elt, h_b, lda, devptrb, ldda)

!---- Call magma LU ----------------
      call magma_wtime_f(tstart)
      call magmaf_zgetrf_gpu(n, n, devptra, ldda, ipiv, info)
      call magma_wtime_f(tend)

      if ( info .ne. 0 )  then
         write(*,*) "Info : ", info
      end if

!---- Call magma solve -------------
      call magmaf_zgetrs_gpu('n', n, nrhs, devptra, ldda, ipiv, devptrb, ldda, info)

      if ( info .ne. 0 )  then
         write(*,*) "Info : ", info
      end if

!---- h_X = devptrB
      call cublas_get_matrix(n, nrhs, size_of_elt, devptrb, ldda, h_x, lda)

!---- Compare the two results ------
      anorm = zlange('I', n, n,    h_a, lda, work)
      bnorm = zlange('I', n, nrhs, h_b, lda, work)
      xnorm = zlange('I', n, nrhs, h_x, lda, work)
      call zgemm('n', 'n', n,  nrhs, n, zone, h_a, lda, h_x, lda, mzone, h_b, lda)
      rnorm = zlange('I', n, nrhs, h_b, lda, work)

      write(*,*)
      write(*,*  ) 'Solving A x = b using LU factorization:'
      write(*,105) '  || A || = ', anorm
      write(*,105) '  || b || = ', bnorm
      write(*,105) '  || x || = ', xnorm
      write(*,105) '  || b - A x || = ', rnorm

      flops = 2. * n * n * n / 3.
      t = tend - tstart

      write(*,*)   '  Gflops  = ',  flops / t / 1e9
      write(*,*)

      rnorm = rnorm / ( (anorm*xnorm+bnorm) * n * dlamch('E') )

      if ( rnorm > 60. ) then
         write(*,105) '  Solution is suspicious, ', rnorm
      else
         write(*,105) '  Solution is CORRECT' 
      end if

!---- Free CPU memory
      deallocate(h_a, h_x, h_b, work, ipiv)

!---- Free GPU memory
      call cublas_free(devptra)
      call cublas_free(devptrb)
      call cublas_shutdown()

 105  format((a35,es10.3))

      end