Pardiso - phase 11 errors

[Jared_W_]

Thank you for taking the time to look into this problem with me. I am trying to implement the Pardiso solver for a nonlinear finite element code, particularly fluid flow, and I seem to be getting stuck at the start. 

When I compile and run I get three messages:

1. I get a message of "The file .\pardiso_ooc.cfg was not opened"
2. *** Error in PARDISO ( sequence_ido,parameters) error_num= 18
3. *** Input check: unexpected error with working arrays ia and/or ja

I have been searching for how to fix these issues and I cannot find anything that works! So here are the details of my code once I start trying to implement Pardiso:

MKL_INT mtype = 11; /* Real unsymmetric matrix */
MKL_INT nrhs = 1; /* Number of right hand sides. */
/* Pardiso control parameters. */
MKL_INT iparm[64];
MKL_INT maxfct, mnum, phase, error, msglvl;
/* Auxiliary variables. */
double ddum; /* Double dummy */
MKL_INT idum; /* Integer dummy. */

/* Internal solver memory pointer pt, */
/* 32-bit: int pt[64]; 64-bit: long int pt[64] */
/* or void *pt[64] should be OK on both architectures */
int *pt[64];
MKL_INT* perm;
perm = (MKL_INT*) malloc (sizeof (MKL_INT)*n);

/* -------------------------------------------------------------------- */
/* .. Initialize the internal solver memory pointer. This is only */
/* necessary for the FIRST call of the PARDISO solver. */
/* -------------------------------------------------------------------- */
for (i = 0; i < 64; i++)
{
pt = 0;
}
/* -------------------------------------------------------------------- */
/* .. Initialize the internal solver memory pointer. This is only */
/* necessary for the FIRST call of the PARDISO solver. */
/* -------------------------------------------------------------------- */
for (i = 0; i < n; i++)
{
perm = 0;
}
/* -------------------------------------------------------------------- */
/* .. Setup Pardiso control parameters. */
/* -------------------------------------------------------------------- */
iparm[1] = 0; /* Default Values */
iparm[28] = 0; /* mtype 11 double precision */
iparm[35] = 1; /* = 1 for zero-based indexing; = 0 for one-based indexing */
maxfct = 1; /* Maximum number of numerical factorizations. */
mnum = 1; /* Which factorization to use. */
msglvl = 1; /* Print statistical information in file */
error = 0; /* Initialize error flag */

/* -------------------------------------------------------------------- */
/* .. Reordering and Symbolic Factorization. This step also allocates */
/* all memory that is necessary for the factorization. */
/* -------------------------------------------------------------------- */
phase = 11;
pardiso (pt, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, perm, &nrhs, iparm, &msglvl, &b, &x, &error);
if (error != 0) {
printf("\nERROR during symbolic factorization: %d", error);
exit(1);
}
printf("\nReordering completed ... ");
printf("\nNumber of nonzeros in factors = %d", iparm[17]);
printf("\nNumber of factorization MFLOPS = %d", iparm[18]);

Please let me know if any additional information is needed!
Sincerelywils7716

[Alexander_K_Intel2]

Hello,

Have you checked your matrix by matrixchecker (iparm[26]=1)? If yes could you send to us reproducer to investigate this issue?

Thanks,

Alexander Kalinkin

[Jared_W_]

Ok this brings up a quick question, since I am in C++ and everything is zero-based indexing when I am implementing the iparm array is everything 1 value off from Table 2: PARDISO iparm() description? i.e., the first 4 iparm in Table 2: PARDISO iparm() description iparm[1], iparm[2], iparm[4], iparm[5] are actually implemented as iparm[0], iparm[1], iparm[3], iparm[4]?

It seems odd that you dont immediately specify the indexing before everything. Since I am not nearly an efficient programmer its highly likely I am thinking about it wrong, but any light on this issue would be greatly appreciated.

In any case I tried both ways (iparm[27] = 1 and iparm[26] = 1) and I did not get any error message.

Thanks for the quick reply btw!

Jared

[Jared_W_]

I also remembered a small thought I had. My 'A' matrix that I am solving has a few rows that are nothing but zeros and this affects the length of my ia array. I know array ia is supposed to be n+1 long but since I have rows with all zeros it turns out to be less than n+1. n I am taking to be the actual size of matrix A[nxn].

First off is having rows with zeros an issue, and second do I make ia[n+1] regardless, or do I allow it to be the true length < n, and then add one more element that is the number of non-zeros + 1 as I am told to do?

[mecej4]

I think that you need to reduce the indices that you used to assign values to the elements of the iparm array by 1. This has to do with the fact that in Fortran, the default base element is iparm(1), whereas in C it is iparm[0]. This issue is not quite the same as specifying zero-base or one-base for the index arrays ia and ja.

The table under "Pardiso iparm parameter" in the MKL manual  has the first entry as iparm(1), which should signify to you that 1-base indexing is being used, following Fortran conventions.

[Jared_W_]

I also remembered a small concern, my A matrix which is n x n (n = number of equations,rows, or colums either one) has a few rows that are all zeros so when I am assembling my ia array its length is < n. The reference manual for pardiso in particular this ia array says it needs to be length n + 1 where the first n of the elements are the first non-zero elements in each row according to the indices from the a array  and the last element n+1 is the total number of non-zero elements. 

So in my case here's my question, if I have mutliple rows that are all zeros do I:

1. create ia to be the actual length of the non-zeros?
2. or do I make it n + 1 in length filling all as necessary, leaving all the rest 0's, and then filling the last element (n+1) as the number of non-zeros as I am told to do?

Hopefully this makes sense, this is one of the more confusing parts.

Let me know what you guys think.

Jared

[Jared_W_]

Ok I have produced a change in the output, now it reads as follows:

*** Error in PARDISO ( sequence_ido,parameters) error_num= 8
*** Input check: ia(neqns+1)_new -1 _old 0 are incompatible
*** Input parameters: inconsistent error= 8 max_fac_store_in: 1
matrix_number_in : 1 matrix_type_in : 11
ido_in : 11 neqns_in : 459
ia(neqns_in+1)-1 : -1 nb_in : 1

ERROR during symbolic factorization: -1

I saw from a post "Alexander Kalinkin (Intel) Mon, 11/14/2011 - 17:03" there was a case the exact same as this and he mentioned to, "Hi, One can see the next line: "ia(neqns_in+1)-1 : -1" so its seem that you have some problem with ia array. Could you print the value of neqns and ia(neqns_in+1) before first pardiso call?". That is what I have printed out.

neqns = 459
ia[nequns + 1] = 3585 << this value is the total number of non zero's + 1

There were no more updates to that post after that so I am at a loss of how to proceed. With this and the last post in mind I look forward to anyone's response!

Jared

[Alexander_K_Intel2]

Ok, let me answer on your questions:

1. create ia to be the actual length of the non-zeros?

No, it is not correct. length of ia is neqns+1 so if in case you have few zero string several ia elements will be the same. briefly ia[j+1]-ia show number of nonzero elements in j-th string. 

1. or do I make it n + 1 in length filling all as necessary, leaving all the rest 0's, and then filling the last element (n+1) as the number of non-zeros as I am told to do?

ia[n+1] = nnz+1 so ia[n+1]-ia express number of nonzero elements in n-th string.

Thanks,

Alexander Kalinkin

[mecej4]

Jared, if you have one or more rows that are full of zeros, the matrix is singular and unsuitable for use with Paradiso.

Given a matrix that is singular but not obviously so, Pardiso can be used to detect singularity.

If you have an entire row of zeros, however, the matrix is obviously singular and Pardiso cannot do much with such a matrix. What do you expect to do with singular matrices?

[Jared_W_]

Mecej4

Maybe this is where I am screwing things up...I am basically trying to implement Newton-Raphson. Using the djacobi(fcn, &n, &m, fjac, x, &eps) I have made my Jacobian matrix and then I break this down into the a, ia, ja arrays to pass to pardiso. 

It is the jacobian that I recieved that has the zero rows, does this mean something is wrong in the jacboian calculation but for some reason it is still giving me a full jacobian without any errors?

Thank you Alexander that makes much more sense now and I believe I have my ia array filled properly.

I am sorry to say though, now my code fails once it calls pardiso. Is there anyway to test which parameter is causing the first call to pardiso to fail?

[Jared_W_]

What I am trying to say is, is it normal for the jacobi routine to finish yielding a singular jacobian and not throwing an error?

[Alexander_K_Intel2]

Hello again :)

1. I am sorry to say though, now my code fails once it calls pardiso. Is there anyway to test which parameter is causing the first call to pardiso to fail?

the main one - matricchecker (iparm[26]=1)+ msglvl=1

the other reasons - incorrect length of a, ia, ja arrays - check their size. 

If all data is correct from you point of view - send this exampe to us.

2. What I am trying to say is, is it normal for the jacobi routine to finish yielding a singular jacobian and not throwing an error?

it can happened, for example, when function is equal to constant. In any case i am a bit confused - how you use in one program jacobi routine that produce dense matrix and pardiso, which work with sparse one?

Thanks,

Alexander Kalinkin

[mecej4]

Jared, I'd conjecture that it is an error that your Jacobian is singular. You should probably investigate your problem formulation more thoroughly, before attempting to solve the discretized equations using Pardiso or another solver.

Consider it this way. Newton's method is based on the approximation

     F(X) ≈ F(X₀) + J(X)·(X-X₀)

You evaluate F(X) and J(X) at the trial solution point X, and solve this equation for X₀ so as to make F(X₀) = 0. If J is singular, however, none of this works, and you will need to include the next term in the Taylor expansion, which will bring in second-order derivatives.

The DJACOBI routine merely helps you to compute an approximate Jacobian. Given its limited scope, it should not be expected to test if the Jacobian is singular, etc. It is only when the Jacobian is applied to a problem where the singularity is an impediment that this becomes an issue.

[Jared_W_]

Maybe I should give more detail about what is exactly going on....

I am creating a finite element code that solves the Navier-Stokes equations in their primitive variable formulation (i.e, velocity in x and y direction and pressure). In this formulation the equations are nonlinear and need to be solved iterativly. 

There first thing I tried was the nonlinear least squares optimization solver provided by the Intel MKL. This was working and I am confident in this because I was getting proper results.

This brings up a point:

• Doesn't my formulation and the jacobian have to be correct since the TR solver was giving me results? It uses the same jacobian routine.

Regardless, the problem was the TR was extremely slow and I couldnt increase the degrees of freedom without a dramatic increase in computation time. So I thought I would try implementing a different routine to solve.

I thought I would create a Newton-Raphson routine similar to what is described at this link under the section 6.2 nonlinear system of equations

http://en.wikipedia.org/wiki/Newton's_method

Using this formulation of Newton's Method (Newton Raphson), I thought I would use djacobi(...) to solve for the jacobian of my residual, and then use pardiso(...) to solve for the incremental change in the solution. I could then use this to increment the solution till convergence.

I thought I would need pardiso because my matrices are large, sparse, and unsymmetric. To be honest I am really struggling to find a fast, efficient routine to solve my nonlinear system of equations and I am trying to figure something out. If there is a sound of desperation thats because there is lol.

You guys are amazing and I just want to say thank you in the middle of all this for your help and patience.

Jared

[mecej4]

When the number of unknowns is huge, as in solving the discretized Navier-Stokes equations, Newton's method has several disadvantages. If, furthermore, the Jacobian has to be computed using finite-difference approximations of the derivatives, it becomes very inefficient. Perhaps, you will appreciate why if you note that in each iteration (or in each time step, if solving the transient N-S equations) the non-zero elements of the Jacobian have to be evaluated. There is a large body of published work on solving the N-S equations that you will profit from reading. From the point of view of numerical analysis, you should note that many improvements upon the Newton-Raphson method involve building up the Jacobian gradually by applying low-rank updates rather than doing a full reevaluation in each iteration.

Doesn't my formulation and the jacobian have to be correct since the TR solver was giving me results? It uses the same jacobian routine.

The TR solver is applied to problems that are usually over-determined. That is, the number of equations, m, is larger than the number of unknowns, n, and the rectangular Jacobian matrix does not have an inverse. 

All of that having been said, I see no reason why the Jacobian of the discretized N-S equations should be singular. As I stated earlier, I suspect a mistake -- furthermore, in cases where a derivative is zero you should be able to confirm or disprove that fact by mere inspection, rather than having to do numerical differentiation.

One reason for apparently zero derivatives to be computed is attempting to compute

     f'(x) ≈ [f(x + Δx) - f(x)] / Δx

using a magnitude of Δx that is too small -- so small that |f(x + Δx)/f(x) - 1| < ε_mach

[Jared_W_]

Alright I fixed the jacobian, now what should I send your way so you guys can reproduce it your end? 

[Alexander_K_Intel2]

Jared.

So finally we have a problem with pardiso run, am i correct? If yes can you send for us matrix and part of code to reproduce problem?

Thanks,

Alex

[Jared_W_]

First off thanks to mecej4 for pointing out the problem with the jacobian, I got that fixed!

Alex

Do you need me to just put my A matrix, b array into  2 separate text files and then the portion of the code where I call pardiso? Basically are you going to read in my text files into arrays and then send them to pardiso(...)?

If so I have attached the files with the proper data so you guys can read in text file and fill proper arrays. Let me know what else you will need.

Jared

[Jared_W_]

Woops, here are the files

[Jared_W_]

Alexander, is this the information you need to be able to replicate where my error is coming from?

[Alexander_K_Intel2]

Hi Jared,

Sorry for skipping answer, will check your example during holidays, ok?

Thanks,

Alex