failure pardiso with CBWR

[____]

Hello everyone.

I am using pardiso to solve unsymmetric system with CBWR(Conditional Bitwise Reproducibility) , but some special case it didn't work in symbolic factorize.

The test matrix is simple 3x3 matrix.

May I check what is the problem ?

I am using test library which is located in "compilers_and_libraries_2017.2.187". It may be "Intel Math Kernel Library 2017 Update 2".

Thank you for advance.

// SimpleCrashTest.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"

/*******************************************************************************
* Copyright 2004-2016 Intel Corporation All Rights Reserved.
*
* The source code,  information  and material  ("Material") contained  herein is
* owned by Intel Corporation or its  suppliers or licensors,  and  title to such
* Material remains with Intel  Corporation or its  suppliers or  licensors.  The
* Material  contains  proprietary  information  of  Intel or  its suppliers  and
* licensors.  The Material is protected by  worldwide copyright  laws and treaty
* provisions.  No part  of  the  Material   may  be  used,  copied,  reproduced,
* modified, published,  uploaded, posted, transmitted,  distributed or disclosed
* in any way without Intel's prior express written permission.  No license under
* any patent,  copyright or other  intellectual property rights  in the Material
* is granted to  or  conferred  upon  you,  either   expressly,  by implication,
* inducement,  estoppel  or  otherwise.  Any  license   under such  intellectual
* property rights must be express and approved by Intel in writing.
*
* Unless otherwise agreed by Intel in writing,  you may not remove or alter this
* notice or  any  other  notice   embedded  in  Materials  by  Intel  or Intel's
* suppliers or licensors in any way.
*******************************************************************************/

/*
*   Content : Intel(R) MKL PARDISO C example
*
********************************************************************************
*/
/* -------------------------------------------------------------------- */
/* Example program to show the use of the "PARDISO" routine */
/* on unsymmetric linear systems */
/* -------------------------------------------------------------------- */
/* This program can be downloaded from the following site: */
/* www.pardiso-project.org */
/* */
/* (C) Olaf Schenk, Department of Computer Science, */
/* University of Basel, Switzerland. */
/* Email: olaf.schenk@unibas.ch */
/* -------------------------------------------------------------------- */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "mkl_pardiso.h"
#include "mkl_types.h"
#include "mkl_service.h"

//#define Diagonal
MKL_INT main(void)
{
   int my_cbwr_branch = 0;
   my_cbwr_branch = mkl_cbwr_get_auto_branch();
   mkl_cbwr_set(my_cbwr_branch);
   /* Matrix data. */
   MKL_INT n = 3;
#ifdef Diagonal
   MKL_INT ia[4] = { 1, 2, 3, 4 };
   MKL_INT ja[3] =
   { 1,
       2,
         3
   };
   double a[3] =
   { 43.157629078689801, 
     53947.036348362199,
     1601061.7235612301
   };
#else
   MKL_INT ia[4] = { 1, 3, 5, 6 };
   MKL_INT ja[5] =
   { 1, 2,
      1, 2,
      3
   };
   double a[5] =
   { 43.157629078689801, 0.0000000000000000,
      0.000000000000000, 53947.036348362199,
      1601061.7235612301
   };
#endif // Diagonal

   MKL_INT mtype = 11;       /* Real unsymmetric matrix */
                             /* RHS and solution vectors. */
   double b[3], x[3];
   MKL_INT nrhs = 1;     /* Number of right hand sides. */
                         /* 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 */
   void *pt[64];
   /* Pardiso control parameters. */
   MKL_INT iparm[64];
   MKL_INT maxfct, mnum, phase, error, msglvl;
   /* Auxiliary variables. */
   MKL_INT i, j;
   double ddum;          /* Double dummy */
   MKL_INT idum;         /* Integer dummy. */
   char *uplo;
   /* -------------------------------------------------------------------- */
   /* .. Setup Pardiso control parameters. */
   /* -------------------------------------------------------------------- */
   for (i = 0; i < 64; i++)
   {
      iparm = 0;
   }
   iparm[0] = 1;         /* No solver default */
   iparm[1] = 2;         /* Fill-in reordering from METIS */
   iparm[3] = 0;         /* No iterative-direct algorithm */
   iparm[4] = 0;         /* No user fill-in reducing permutation */
   iparm[5] = 0;         /* Write solution into x */
   iparm[6] = 0;         /* Not in use */
   iparm[7] = 2;         /* Max numbers of iterative refinement steps */
   iparm[8] = 0;         /* Not in use */
   iparm[9] = 13;        /* Perturb the pivot elements with 1E-13 */
   iparm[10] = 1;        /* Use nonsymmetric permutation and scaling MPS */
   iparm[11] = 0;        /* Conjugate transposed/transpose solve */
   iparm[12] = 1;        /* Maximum weighted matching algorithm is switched-on (default for non-symmetric) */
   iparm[13] = 0;        /* Output: Number of perturbed pivots */
   iparm[14] = 0;        /* Not in use */
   iparm[15] = 0;        /* Not in use */
   iparm[16] = 0;        /* Not in use */
   iparm[17] = -1;       /* Output: Number of nonzeros in the factor LU */
   iparm[18] = -1;       /* Output: Mflops for LU factorization */
   iparm[19] = 0;        /* Output: Numbers of CG Iterations */
   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 */
                         /* -------------------------------------------------------------------- */
                         /* .. 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;
   }
   /* -------------------------------------------------------------------- */
   /* .. 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, &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &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]);
   /* -------------------------------------------------------------------- */
   /* .. Numerical factorization. */
   /* -------------------------------------------------------------------- */
   phase = 22;
   PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
      &n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &error);
   if (error != 0)
   {
      printf("\nERROR during numerical factorization: %d", error);
      exit(2);
   }
   printf("\nFactorization completed ... ");
   /* -------------------------------------------------------------------- */
   /* .. Back substitution and iterative refinement. */
   /* -------------------------------------------------------------------- */
   phase = 33;
   /* Set right hand side to one. */
   for (i = 0; i < n; i++)
   {
      b = 1;
   }
   //  Loop over 3 solving steps: Ax=b, AHx=b and ATx=b
   for (i = 0; i < 3; i++)
   {
      iparm[11] = i;        /* Conjugate transposed/transpose solve */
      if (i == 0)
         uplo = "non-transposed";
      else if (i == 1)
         uplo = "conjugate transposed";
      else
         uplo = "transposed";

      printf("\n\nSolving %s system...\n", uplo);
      PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
         &n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, b, x, &error);
      if (error != 0)
      {
         printf("\nERROR during solution: %d", error);
         exit(3);
      }

      printf("\nThe solution of the system is: ");
      for (j = 0; j < n; j++)
      {
         printf("\n x [%d] = % f", j, x);
      }
      printf("\n");
      //// Compute residual
      //mkl_dcsrgemv(uplo, &n, a, ia, ja, x, bs);
      //res = 0.0;
      //res0 = 0.0;
      //for (j = 1; j <= n; j++)
      //{
      //   res += (bs[j - 1] - b[j - 1]) * (bs[j - 1] - b[j - 1]);
      //   res0 += b[j - 1] * b[j - 1];
      //}
      //res = sqrt(res) / sqrt(res0);
      //printf("\nRelative residual = %e", res);
      //// Check residual
      //if (res > 1e-10)
      //{
      //   printf("Error: residual is too high!\n");
      //   exit(10 + i);
      //}
   }

   /* -------------------------------------------------------------------- */
   /* .. Termination and release of memory. */
   /* -------------------------------------------------------------------- */
   phase = -1;           /* Release internal memory. */
   PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
      &n, &ddum, ia, ja, &idum, &nrhs,
      iparm, &msglvl, &ddum, &ddum, &error);
   return 0;
}

[Ying_H_Intel]

Hello,

Could you tell your OS and Processor type and how you build the code?

Best Regards,

Ying

[____]

Hello,

I am using "Visual Studio 2015" in windows 10.

My processor is Intel(R) Xeon(R) CPU E5-2620 v4 @ 2.10GHz.

Thank you for advance.

[Gennady_F_Intel · Accepted solution]

please try to use iparm[33] to take CNR outputs. Please refer to the Developer Reference for MKL to see how to use it.

[____]

Thanks Gennady F.

I set the parameter to

   iparm[59] = 0;        /* input switches between in-core (IC) and out-of-core (OOC) */
   iparm[33] = 1;        /* input Optimal number of OpenMP threads for conditional numerical reproducibility (CNR) mode */

It works.

But the other number of iparm[33] did not work except 1.

Do you know why?

Thank you for advance.

[Gennady_F_Intel]

What do you mean by "did't work exept 1"?   When you set, as an example iparm[33] = 4, Pardiso will not produce the same results?

[____]

Sorry for late response.

pardiso is crashed in symbolic factorize when I set to another number except 1 in iparm[33].

> Can I ask you the question?

If iparm[33] == 1, didn't it work multi-thread ? ,even though "mkl_set_num_threads(8)"

I attach the message of the result of "iparm[33] == 1" case, the name of file is "result_of_iparm_33_set_to_1.zip"

The message is "Parallel Direct Factorization is running on 8 OpenMP". So I think it works multi-thread.

Am I wrong?

In manual, the option of iparm[33] is "Optimal number of OpenMP threads for conditional numerical reproducibility (CNR) mode".

What's mean of optimal number of OpenMP threads ?

> Can I ask you the another question?

I saw the "Set the number of threads to a constant number" for consistence results in CBWR manual.

What's mean of "Set the number of threads to a constant number".

It means result of thread 1 is different result of thread 2 or 4 even though I set the option of "mkl_cbwr_set".

I did test the my problem with pardiso on multi-thread.

In my experience, I didn't catch the difference by using the option of "mkl_cbwr_set".

thank you in advance.

[Ying_H_Intel]

Hi,

Sorry for slow response. Your question make sense.  It seems  a bug of pardiso during setting cbwr for small  matrix.  The issue was escalated to our developer and will update you if any fix. 

Best Regards,

Ying 

[Ying_H_Intel]

update MKL 2017 update 4 released and the related problem: pardiso don't work during setting cbwr for small pardiso matrix should be fixed in the version and later release 2018 update 1. Please find chance to try them. Thanks Ying