question about pardiso arguments at solving phase

[Brian_Murphy]

I want to call pardiso with phase=33 to compute a solution to a given rhs.  I'm sure I need to pass it the pt(:) array from the factoring phase, but do I also need to pass it any other arrays used in the factoring phase, like iparm, perm, a, ia and ja.  In my calls to perform the factoring, I'm using iparm(1)=0 to let pardiso use all defaults.

I intend to pass the pt(:) array to different routines that need to compute solutions to different rhs's.  I need to know if must also pass these other arrays.

[mecej4]

As of now, the PARDISO() and PARDISO_64() routines in MKL take a fixed number of arguments. Therefore, every argument must be present, whether or not it makes sense to use an argument in a particular phase. In fact, arguments that are not used in a particular call can be set to any value, and array arguments that are not used can be dummy arrays of length 1, or unallocated allocatable arrays.

The rules of Fortran 77 dictate having a fixed number of arguments, and this has to take precedence over algorithmic requirements. The alternative would have been to break down the single, multiple-duties, call into a sequence of subroutines with different names and purposes, but this has its own disadvantages and advantages, and is not implemented at present.

[Brian_Murphy]

How do I know which arguments are ok to supply with dummy values for a call with phase=33?  Ones which I suspect must be required are; pt, phase, nrhs, b, x and error.  Will a call with phase=33 ignore all the rest?

[mecej4]

There is next to no gain to be had by guessing which ones are used in a particular call, although it is easy to reason that b and x will not be used until the solution phase is reached.  It is best to assemble all the input arguments before making the first call (phase=1, or calling pardisoinit()).

After phase 2 is completed, you can repeatedly call the solver with phase=33 with different values in b and receive the corresponding solutions in x.

[Brian_Murphy]

I will be making phase=33 calls from other routines than where the factoring is done.  I will run some tests to see how many of pardiso's arguments I will need to pass around.

[mecej4]

Brian Murphy wrote:

I will be making phase=33 calls from other routines than where the factoring is done.  I will run some tests to see how many of pardiso's arguments I will need to pass around.

In such situations you have to be careful. Unlike older library routines where the subroutine behavior was entirely determined by the values of the arguments (and values in shared common block variables), Pardiso maintains an internal set of state variables that are hidden from the user. Therefore, it is important that the sequence of calls to Pardiso be made in proper order.

If, for example, you want to factorize not one, but two matrices and then make interleaved calls with phase=33 for the two sets of factors, you should use two instances of the control variables such as the PT array, etc.

[Brian_Murphy]

That is a good point.  I am aware that pardiso saves the factoring result, and keeps track of things with the pt array.  In my case, I will be factoring a single matrix, and then in several other places compute solutions using the pt array returned at the factoring step.

I now have pardiso working with the arnoldi eigensolver.  It's somewhere around 5% faster than the fortran umfpack library, it's hard to tell for sure because 5% is hard to detect.  Both umfpack and pardiso are compiled together into my code, and which one is used is selected by an input at run time.

Even though the code produces results that appear correct (i.e. pretty close to umfpack's), when examining calls to pardiso with the visual studio debugger, some weird things are happening to some array variable that I'm trying to unravel.

[mecej4]

I now have a test-setup for solving unsymmetric linear equations using Umfpack, MKL-Pardiso, Pardiso5 and MUMPS, all using the latest distributions of these packages.

I have run tests on 11 matrices in Harwell-Boeing RUA format (a text file with some headers followed by the compressed-sparse-column data for the matrix). There is no clear "winner" that I could discern. For a collection of such matrices, see http://math.nist.gov/MatrixMarket/formats.html .

If you will provide one of your large matrices in that RUA format (or in CSR format), I can run the tests on it. The NIST people would also be happy to add your matrix to their collection, I think.

[Brian_Murphy]

How large is "large"?  nx is often from about 100 to a few thousand.  nnz is typically around 4*nx.  The matrix, which is nonsymmetric, comes from a finite element model, and for the most part is tightly banded.

When using F8 to step through code in visual studio, and I come to a CALL statement, sometimes a single press of F8 will step into the routine being called, but in other instances multiple presses of F8 are required.  What causes that?  Also sometimes using F8 to step into a routine first goes to chkstk.asm before getting to the fortran routine.  What causes that?

[mecej4]

A matrix with n=10,000 nnz≈ 5.n is, by today's standards, small to medium, but it will suffice for the tests if the factorization and solution take about 1 second of CPU time.

I rarely find it useful to use a symbolic debugger to step into run-time libraries. Intel Fortran uses the VC libraries quite a bit, and I don't want to step through those routines, since it is rare to find a problem there. For the same reason, I do not use the VC "debug" libraries when building.

F8? What function is that key assigned to? In VS, F10 and F11 are normally assigned "step over" and "step into" functions.

[Brian_Murphy]

In my visual studio F8 is "step into" and shift-F8 is step over.  Sorry about that.  I should have been more specific since I have VS Tools/Options/Keyboard set to "Visual Basic" because it's like Office VBA.

The biggest matrices I work with will only take 1 or 2 tenths of a second.  They are generally factored & solved in loops in the course of doing a job.  So the faster it runs, the better.

When you say VC, do you mean visual C?  I am trying to step into only my own code, not the code of libraries.

I'm having really weird things happen to array variables being passed to subroutines.  I think my problem is with STDCALL and REFERENCE attributes.  Some of my routines are exported as follows so I can call them from Excel VBA.

!DEC$ ATTRIBUTES DLLEXPORT, STDCALL, REFERENCE, ALIAS:"SUBA" :: SUBA

If I also call this routine from somewhere else in fortran, will there be trouble with array arguments?

[mecej4]

It would probably easier to use a Fortran driver to call the solver DLL, in the stages of developing the DLL -- solve one problem at a time.

In Intel Fortran (32-bit), the standard calling sequence is CDECL + the special convention for passing Fortran character variables, not STDCALL. Likewise, the default version of MKL does not do STDCALL. If you call Lapack95 or BLAS95, you probably have to build special libraries yourself. In short, STDCALL is "legacy" stuff except when needed in VBA, Excel, etc.

[Brian_Murphy]

I don't call any LAPACK, BLAS or similar routines from VBA, only my own.

if I use the ATTRIBUTES stuff to enable a routine to be called from VBA, what do I need to do to call it correctly from fortran?  I'm only passing numerical arguments, not any characters or strings.

[mecej4]

The non-default attributes of a subprogram need to be specified in the subprogram (so the compiler knows how the entry point will be called and set up the proper reception) as well as in every caller (so the compiler knows how to set up the call properly).

[Brian_Murphy]

I don't know how to do that.  So I will make two versions for each routine called from both VBA and fortran.

[Brian_Murphy]

As far as I can tell, that seems to have cured all the problems.  The numerical agreement between umfpack and pardiso is now more to my liking at 10+ digits for eigenvalues returned by arnoldi.  Plus Pardiso now runs about 13% faster than umfpack instead of 5% slower.

I've got a lot more work to do implementing pardiso throughout my code.  First up will be to make sure I can build a 64 bit version to call from 64 bit Excel.

Thank you very much, mecej4, for all the help you've been giving me :)