Solved: CSR matrix handling and solution error - Page 2

DarioMangoni · ‎09-10-2015

Hi,
I built my own class for CSR3 matrices for Mkl Pardiso and I have to know how Mkl knows the length of arrays (values and columnIndex ):
- Is it done by rowIndex[number_of_rows]?
- Is it done by some information provided by mkl_alloc?

I know it seems a very stupid question but it's a whole month that I'm working on this problem:
1st case
- I preallocate the matrix with an estimated number of nonzeros (element that are not really initialized has a flag that say that are not used)
- I compress it, deleting all unused/uninitialized elements
2nd case
- I preallocate the matrix*
- I read the matrix from an external file

The two matrices at the end are ABSOLUTELY IDENTICAL. It changes the way they're allocated first.
The case one works very well, but after a random number of steps of the dynamic problem it has a huge residual (sometimes >e+14, yes it's a plus)! At that time I write the matrix on a file, I read it back to another matrix and then I give this matrix to Pardiso again and puff! All is perfect again! Residual e-15!

Do you have any suggestion?

Many thanks,
Dario

*in the 2nd case it works also if I preallocate whatever number of nonzero and then I compress. It works well at anytime!
** comparing IPARM output I noticed that or iparm(63) one-indexed it changes everytime but in the step where I face the problem there's a -1 in the 1st case and 0 in the second case... I don't know if it is useful

Here my settings: ONE-indexed

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

/*IPARM easy settings*/
       IPARM(1) = 1;               /* No default values for solver */
       IPARM(6) = 0;               /* Write solution on u */
       IPARM(12) = 0;               /* Solve with transposed/conjugate transposed matrix [def: 0, solve simply A*x=b]*/
       IPARM(18) = -1;               /* Report number of nonzeros */
       IPARM(19) = -1;               /* Report number of floating point operations */
       IPARM(35) = 1;               /* Zero based indexing */
       IPARM(27) = 1;               /* Matrix checker */
       IPARM(28) = 0;               /* Double precision */
       IPARM(36) = 0;               /* Schur complement matrix computation control [def:0, do not compute Schur] */
       IPARM(56) = 0;               /* Diagonal and pivoting control [def:0, disabled] */
       IPARM(60) = 0;               /* In-Core (OC) / Out-Of-Core (OOC) switch [def:0, IC mode] */

       /* IPARM fine settings */
       IPARM(2) = 2;               /* Fill-in reducing ordering [def:2] */
       IPARM(4) = 0;               /* Preconditioned CGS/CG [def:0] - HIGHLY RECOMMENDED */
       IPARM(5) = 0;               /* User fill-in reducing permutation [def:0, default filling]*/
       IPARM(8) = 10;               /* Maximum number of iterative refinement steps */
       IPARM(10) = 13;               /* Perturb the pivot elements with 1E-value [def:13 for unsymm/8 for symm, 1e-13 perturbation]*/
       IPARM(11) = 1;               /* Scaling MPS [def: 1 for unsym/0 for sym, scaling enabled] */
       IPARM(13) = 1;               /* Maximum weighted matching algorithm [def: 1 for nonsymmetric matrices, enabled) */
       IPARM(21) = 1;               /* Pivoting for symmetric indefinite matrices */
       IPARM(24) = 0;               /* Parallel factorization control [def:0, classic algorithm] */
       IPARM(25) = 0;               /* Parallel/sequential solve control [def:0, parallel] */
       IPARM(31) = 0;               /* ADV Partial solve and computing selected components of the solution vectors [def:0, disable]*/
       IPARM(34) = 0;               /* ADV Optimal number of threads for conditional numerical reproducibility (CNR) mode [def:0, disable]*/

mecej4 · ‎09-27-2015

Dario, here are some observations to help you. Your matrices have many zero entries, and these zero entries do two bad things. A sparse matrix should not contain any, unless required as a placeholder (e.g., a diagonal element of a symmetric or structurally symmetric matrix). I removed the zero entries, and put your data into the following format (files with .CSR suffix):

  Line-1: title                               
  Line-2: n,nnz,index base                    
  Lines 3 to n+2 : rowptr,rhs, i=1:n          
  Lines n+3 to n+2+nnz : col(1:nnz),val(1:nnz)

I found that the removal of the zero-value matrix elements greatly improved the results from Pardiso. The residual norm remains at or below 10^-7, and the solutions agree to six or more digits when the RHS values are switched from 12 to "24" digits. Note that double precision does not give you more than 16 digits of precision, even if you print out 24 digits.

The second issue is one that I am unable to address: a good fraction of the solution values are zero. Perhaps you can tell from your understanding of the physics which x_ivalues should come out to zero, even before performing the solution. If so, you can eliminate the columns corresponding to those x values, and discard a number of equations equal in number to the count of zero-valued solutions.

For example, Beam35 has 270 unknowns and 3276 matrix entries. However, only 2497 of those entries are non-zero, and 105 of the 270 elements of the solution are zero.

The .CSR files and files with the solution vectors are in the attached zip file.

View solution in original post

DarioMangoni · ‎09-28-2015

mecej4 wrote:

I do not know what you mean by "saddle-point problems". Please elaborate.

The problem is of the kind:

| A C^T | * | x | = | b |
| C 0 | | y | | f |

where 0 stands for a block of all zeros. So the linear system is Ax+C^Ty = b; Cx=f;
At this time the matrix (aka saddle-point or Karush–Kuhn–Tucker, KKT) is put into Pardiso without exploiting its particular shape... But I think that something could be done.

Schur complement? I mean... it's also true that:
A x = ( b - C^T y) --> x = A^-1 b - A^-1 C^T y
C A^-1 C^T y = C A^-1 b - f
So I need to do something with this system?

mecej4 wrote:

I should like to know: do you know which x-values may turn out to be zero before solving the simultaneous equations?

No, I don't. "x" is a vector of speed; I do can say that some of its component are zero because in this case the pendulum is swinging on a plane, but it's a particular case! Actually the "y" are lagrangian multipliers i.e. the forces in the constraints. In this case the forces should should be zero in some of their component as well because there shouldn't be any forces orthogonal to the swinging plane. But I can't assume that they will be always zero in every problem.

mecej4 · ‎09-28-2015

Thanks for the clarifications; I'll ponder over the questions -- no answer at the top of my head.

In the meantime, I succeeded in reproducing the problem with an unrelated matrix pulled from the NIST Matrix Market, http://math.nist.gov/MatrixMarket/data/Harwell-Boeing/chemimp/impcol_a.html . As downloaded, the matrix has no zero-valued entries, and is a 207 x 207 matrix with 572 entries. I made up a right-hand-side vector that would give a solution vector with no zeros. I added 35 zero-valued entries to one of the previously empty super-diagonals, and found the norm of the residuals jumping from 2D-13 to 5D-2. I tried other matrices that were much smaller in size, the did not see the effect (shall we call it the Dario Effect?) with them.

I am hunting around for other test matrices, and I am thinking about reporting the problems to this forum and to the Pardiso site. I hope you will have no objections.

DarioMangoni · ‎09-28-2015

Ahahah, nice! Dario effects sounds good!
Of course you can report whatever can be useful to eventually fix some bugs. I will benefit of the fix me too!

Thanks for your help,
if you have any other question/answer let me know!

mecej4 · ‎10-07-2015

DarioMangoni: I have made up a reproducer and reported the problem (Dario Effect) at https://software.intel.com/en-us/forums/intel-math-kernel-library/topic/595201 . Please review and add comments as necessary.