I have built a research FE-code in Fortran 90/95 and I'm using the Pardiso solver for the linear system of equations that results. I have used it for quite some time and it have been working beautifully. The problem I want to solvegives mea real and structurally symmetric matrix, mtype = 1.
The problem is that I recently added the ability to have Lagrange multipliersso I get some zerosalong the diagonal. I get the error code: -4, "zero pivot". I have tried to figure out if I need to supply some additional information to iparm when I have zeros on the diagonal, but to be honest I reachthe end of my linear algebra knowledge once going beyond LU-decomposition...
I have checked that thearrays ia and ja are correct by mimicking the process in matlab and solving the system there - which works just fine, albeit slow and can only handle smaller problems. I set upthe system matrix in matlab using the arrays ia and ja thatI get from the Fortran code, and for small problems I can check and make sure that the resultingmatrixis correct (and I get resonable solutions).
I also tried setting the zeros on the diagonal to a small number, 1.d-30 and 1.d-16 but neither worked.
My current settings to Pardiso are (which I guess should be very close to the default settings):
nrhs = 1
maxfct = 1
mnum = 1
iparm(:) = 0
iparm(1) = 0 ! no solver default
iparm(2) = 2 ! fill-in reordering from METIS
iparm(3) = mkl_get_max_threads() !numbers of processors, value of MKL_NUM_THREADS
iparm(4) = 0 ! no iterative-direct algorithm
iparm(5) = 0 ! no user fill-in reducing permutation
iparm(6) = 0 ! =0 solution on the first n components of x
iparm(7) = 16 ! default logical fortran unit number for output
iparm(8) = 0 ! numbers of iterative refinement steps
iparm(10) = 13 ! perturb the pivot elements with 1E-13
iparm(11) = 1 ! use nonsymmetric permutation and scaling MPS
iparm(14) = 0 ! Output: number of perturbed pivots
iparm(18) = -1 ! Output: number of nonzeros in the factor LU
iparm(19) = -1 ! Output: Mflops for LU factorization
iparm(20) = 0 ! Output: Numbers of CG Iterations
iparm(27) = 1 ! Check matrix formulation (ia, ja and aa)
mtype = 1 ! Matrix type
Any help with this would be greatly appreciated. Many thanks, have a nice day!
That helped, thanks alot! I used mtype=11 and now it solves perfectly.
A question relating to performance: in your opinion, should I leave mtype=11 always, or should I set it back to mtype=1 when not using Lagrange multipliers? Is type 11 slower than type 1? Also, is it not a little strange that it did not work for type 1, I mean that is actually the kind of matrix I have and adding Lgr. mult. should not change that since they only add symmetric contributions to the stiffness matrix.
Anyway, thanks again for the quick response.
PS. Yes, the iparm thing was a mistake. I copied that part of the input and I think the previous user might be running the Univ. of Basel version.
I will try to get you a sample matrix as soon as possible. Code is working now so the focus is on getting results right now... I will put it high up on the "to do list".
You did understand correct. I had a matrix, lets call it K, that would solve with mtype = 1. Then I augmented that with some Lagrange multipliers given by the matrix C (almost all zeros and a few 1 and -1)so that the new system is
K_aug = [ K C' ]
[ C 0 ]
where C' denotes the transpose of C, and 0 should be interpreted as the zero-matrix of apropriate dimensions. I made sure that the zeros on the diagonal was explicitly stored.
I did quite a few test-runs last night and in most cases the solver worked perfectly once I switched to mtype = 11. However I still occationally get the "zero pivot" error, but only in extreme cases and I can skip around that by reducing my step size in the equilibrium iterations and playing around a little with iparm(4) and iparm(21).
Anyway, once again thanks for the help. I'll get back to you with some test cases.
I made two example structures, one with and one without Lagrange multipliers. One that solves with mtype = 1 and one that won't solve unless I change to mtype = 11. Except for the added Lagrange multpliers they are identical.
I tried to add them via the "Add Files" button, but I must say that the interface was a little confusing. I only got to create a folder, and I couldn't see any way of actually adding the files. How do I do this? Would you prefer to get them via mail instead?
Ok, think I figured out how to do it now. Anyway, the files are pure text-files. The structure is; two stars **marks the start of a new array, followed by the name (using standard naming), and then a number to denot the length of the array. For example:
** ia, 206
** ja, 9313
and so forth. I supplied ia, ja, aa, and rhs.
Sorry for the slightly retarded format, but it felt like the simplest method. The arrays correspond to a fairly small problem.
The array aa might turn out to be symmetric (i.e. mtype = 2), but that is only because it comes from the first iteration of the algorithm. Later on in the solution process the stiffness matrix is guaranteed to be only strucurally symmetric (mtype = 1).
The structure in Arrays_Standard.out solves with mtype = 1, but in Arrays_LgrMult.out the structureis augmented with 16 Lagrange multiplier equations, and now it only solves with mtype = 11. But the type of matrixhave not have changed. Zeros on the diagonal are explicitly stored.