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!
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.