magma_zgetrf_gpu choice of pivot type

[fletchjp]

I have been using magma_zgetrf_gpu and magma_zgetrs_gpu to factorise a matrix and then apply it multiple times to different single vectors.

Up until now the calls to magma_zgetrf_gpu and magma_zgetrs_gpu have all been in the same program unit, as is the case with the example testing_zgetrf_gpu_f.F90 with magma 1.6.2 and earlier versions.

These routines need an integer vector named as ipiv. This is defined in the example as allocatable and allocated, and I have copied the example:

Code: Select all

       integer,    allocatable       :: ipiv(:)
       allocate(ipiv(n))

This all works fine as long as the calls to magma_zgetrf_gpu and magma_zgetrs_gpu are in the same program unit. For JDQZ the location of the two calls cannot be in the same program unit, as the method needs a callback subroutine for preconditioning which uses the factorisation. The information needed about the matrix cannot be in the arguments, as the call comes from the algorithm code.

The rules of Fortran 90 mean that the allocatable vector cannot be passed in common. I have worked on two solutions to the problem.

1. Use the allocatable vector with magma_zgetrf_gpu and copy the result to an integer vector in common. In the callback subroutine declare a fresh allocatable vector and copy to it the information from the vector in common. After the call to magma_zgetrs_gpu the allocated vector is deallocated, and this happens every time the routine is called, typically a hundred per calculation case.

2. Use an integer pointer and allocate that. This pointer can be passed in common so there is no need to allocate a vector or copy the data.

Code: Select all

       integer,    dimension(:), pointer :: ipiv_p
       common /pivot/ ipiv_p
       allocate(ipiv_p(n))

There is a small reduction in time for a small test case.

I have not detected any problems with the results from doing it this way.

Please let me know if there is any misunderstanding on my part. I am not very familiar with the constructs of Fortran 90 as compared with Fortran 77. I started with Fortran 2 a very long time ago.

Best wishes

John