Community
cancel
Showing results for 
Search instead for 
Did you mean: 
____
Beginner
107 Views

failure pardiso with CBWR

Jump to solution

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;
}

 

0 Kudos
1 Solution
Gennady_F_Intel
Moderator
107 Views

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

View solution in original post

8 Replies
Ying_H_Intel
Employee
107 Views

Hello,

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

Best Regards,

Ying

____
Beginner
107 Views

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
Moderator
108 Views

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

View solution in original post

____
Beginner
107 Views

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
Moderator
107 Views

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?

____
Beginner
107 Views

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
Employee
107 Views

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
Employee
107 Views
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
Reply