I'm writting a program that calls the geqp3 frequently. Because the mkl subroutine is geqp3(A,jpvt), Then I need an array A.
However I want tp call geqp3 frequently. In each call the array size may be slightly different. I know that if at each time, I allocate the memory for A and jpvt, and call geqp3 and then deallocate A and jpvt, it works. However the speed becomes slow. Especially when using open MP to parallelize this program, the CPU used is only 30% if the number of threads are set to 32 (my server has 16 cores and 32 threads).
Is there any simple tip to solve this problem. Thank you very much!
You can allocate an array once, making it large enough to cover all the cases in a run, and make calls using the F77 interface. That is, you would make repeated calls to DGEQP3 if your matrix elements are double precision. In each call, you would pass arguments that specify the actual size of the argument, the leading dimension of the array holding the matrix, etc. This, in fact, is how one did things before Fortran 90+ came along.
In doing this, you would have to make sure that you can precompute the largest array size needed, and make sure that between successive calls you have copied out (or otherwise used up) the computed factors returned by DGEQP3 and set the new values into the matrix.