- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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

Link Copied

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

Hello,

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

Best Regards,

Ying

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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?

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page