- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
In pardiso solver, the parameter phase means the solver execution steps.
If Ax=b, and A is a Real and symmetric indefinite matrix, then P*A*P'=LDL'.
For example, phase=332, only diagonal substitution is executed(Dx=b), and if b is an identity matrix, x will be inverse of D.
But, the result is not. What's the matter, and what's the real meaning of phase=331, 332, 333.
Thank you!
Link Copied
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
The decomposition of the work of the Pardiso solver into "phases" gives the user access to some intermediate results and an ability to inject some control over the many phases of the work. However, it would be wrong to assume that each "phase" can be invoked independently or that the sequence of phases can be in an arbitrary order.
The Pardiso documentation (either the separate documentation from Basel/Lugano, at pardiso-project.org, or the MKL sections on Pardiso) should be read and followed in order to issue calls with the proper sequence of values for the "phase" control parameter. .Since there are many internal variables in Pardiso that are not accessible to the caller, the amount of control over the algorithm is somewhat limited. The author(s) probably had to strike a balance between accessibility and robustness.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Thank you for your reply.
Now I need the diagonal matrix D(applying 1x1 diagonal pivoting), would you please tell me how can I get it?
Thank you very much!
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hi X.J.Wang,
What is wrong with the result? Could you please provide us more details?
This discussion http://software.intel.com/en-us/forums/topic/287924 may be help.
Please let us know if you have any progress.
Best Regards,
Ying .
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hi Ying H,
Thank you for your reply. I want to get the D matrix of the factorization of P*A*P'=L*D*L', but I doubt the result is not right.
For example:
n = 3
nz = 6
ia = (/ 1, 4, 6, 7 /)
ja = (/ 1, 2, 3, 2, 3, 3 /)
2 0 3
A = 0 2 1
3 1 2
iparm(5) = 2, and perm array calculated is (/1, 2, 3), so
1 0 0
P = 0 1 0
0 0 1
and A=L*D*L'
the result as follows, and they are not expected.
!************ PARDISO RESULT **************
phase = 331
1 0 0
b = 0 1 0
0 0 1
Lx = b, x = inv(L)=
1.00000 0.00000 1.00000
-0.60000 0.00000 -1.50000
0.00000 1.00000 0.40000
!---------------------------------------------
phase = 332
1 0 0
b = 0 1 0
0 0 1
Dx = b, x = inv(D) =
-0.40000 0.00000 0.00000
0.00000 -0.40000 0.00000
0.00000 0.00000 0.41667
!---------------------------------------------
phase = 333
1 0 0
b = 0 1 0
0 0 1
L'x = b, x = inv(L')=
1.00000 1.00000 -0.60000
0.00000 0.00000 1.00000
-1.50000 -1.50000 0.40000
!************* MATLAB RESULT ***************
L =
1.0000 0 0
0 1.0000 0
1.5000 0.5000 1.0000
inv(L) =
1.0000 0 0
0 1.0000 0
-1.5000 -0.5000 1.0000
D =
2 0 0
0 2 0
0 0 -3
inv(D) =
0.5000 0 0
0 0.5000 0
0 0 -0.3333
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hi Ying H,
Thank you for your reply. I want to get the D matrix of the factorization of P*A*P'=L*D*L', but I doubt the result is not right.
For example:
n = 3
nz = 6
ia = (/ 1, 4, 6, 7 /)
ja = (/ 1, 2, 3, 2, 3, 3 /)
2 0 3
A = 0 2 1
3 1 2
iparm(5) = 2, and perm array calculated is (/1, 2, 3), so
1 0 0
P = 0 1 0
0 0 1
and A=L*D*L'
the result as follows, and they are not expected.
!************ PARDISO RESULT **************
phase = 331
1 0 0
b = 0 1 0
0 0 1
Lx = b, x = inv(L)=
1.00000 0.00000 1.00000
-0.60000 0.00000 -1.50000
0.00000 1.00000 0.40000
!---------------------------------------------
phase = 332
1 0 0
b = 0 1 0
0 0 1
Dx = b, x = inv(D) =
-0.40000 0.00000 0.00000
0.00000 -0.40000 0.00000
0.00000 0.00000 0.41667
!---------------------------------------------
phase = 333
1 0 0
b = 0 1 0
0 0 1
L'x = b, x = inv(L')=
1.00000 1.00000 -0.60000
0.00000 0.00000 1.00000
-1.50000 -1.50000 0.40000
!************* MATLAB RESULT ***************
L =
1.0000 0 0
0 1.0000 0
1.5000 0.5000 1.0000
inv(L) =
1.0000 0 0
0 1.0000 0
-1.5000 -0.5000 1.0000
D =
2 0 0
0 2 0
0 0 -3
inv(D) =
0.5000 0 0
0 0.5000 0
0 0 -0.3333
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hi Ying H,
Thank you for your reply. I want to get the D matrix of the factorization of P*A*P'=L*D*L', but I doubt the result is not right.
For example:
n = 3
nz = 6
ia = (/ 1, 4, 6, 7 /)
ja = (/ 1, 2, 3, 2, 3, 3 /)
2 0 3
A = 0 2 1
3 1 2
iparm(5) = 2, and perm array calculated is (/1, 2, 3), so
1 0 0
P = 0 1 0
0 0 1
and A=L*D*L'
the result as follows, and they are not expected.
!************ PARDISO RESULT **************
phase = 331
1 0 0
b = 0 1 0
0 0 1
Lx = b, x = inv(L)=
1.00000 0.00000 1.00000
-0.60000 0.00000 -1.50000
0.00000 1.00000 0.40000
!---------------------------------------------
phase = 332
1 0 0
b = 0 1 0
0 0 1
Dx = b, x = inv(D) =
-0.40000 0.00000 0.00000
0.00000 -0.40000 0.00000
0.00000 0.00000 0.41667
!---------------------------------------------
phase = 333
1 0 0
b = 0 1 0
0 0 1
L'x = b, x = inv(L')=
1.00000 1.00000 -0.60000
0.00000 0.00000 1.00000
-1.50000 -1.50000 0.40000
!************* MATLAB RESULT ***************
L =
1.0000 0 0
0 1.0000 0
1.5000 0.5000 1.0000
inv(L) =
1.0000 0 0
0 1.0000 0
-1.5000 -0.5000 1.0000
D =
2 0 0
0 2 0
0 0 -3
inv(D) =
0.5000 0 0
0 0.5000 0
0 0 -0.3333
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hello, X.J
Are you working on fortran? There is two sample code in pardiso solverc directory.
pardiso_sym_getdiag_c.c
and pardiso_sym_diag_pivot_c
I just did a quick try and can get correct result. So you may simulate it in your program.
Best Regards,
Ying
Reordering completed ...
Number of nonzeros in factors = 10
Factorization completed ...
Number of pertubed pivots = 0
Print values of D for LDLT decomposition
invdiag [0] = 2.000000000000000e+000
invdiag [1] = 2.000000000000000e+000
invdiag [2] = -3.000000000000000e+000
Press any key to continue . . .
/*
********************************************************************************
* Copyright(C) 2004-2013 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.
*
* *Third Party trademarks are the property of their respective owners.
*
* 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 : MKL PARDISO C example
*
********************************************************************************
*/
/* -------------------------------------------------------------------- */
/* Example program for the "PARDISO" routine to show how to get the */
/* diagonal matrix D from LDLT decomposition for symmetric indefinite */
/* systems */
/* -------------------------------------------------------------------- */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "mkl_pardiso.h"
#include "mkl_types.h"
MKL_INT
main (void)
{
MKL_INT n = 3;
MKL_INT ia[4] = { 1, 4, 6, 7};
MKL_INT ja[7] = { 1, 2, 3,
2, 3,
3
};
double a[6] = { 2.0, 0.0,3.0, 2.0, 1.0, 2.0};
MKL_INT mtype = -2; /* Real symmetric matrix */
/* RHS and solution vectors. */
double b[3], invdiag[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;
double ddum; /* Double dummy */
MKL_INT idum; /* Integer dummy. */
/* -------------------------------------------------------------------- */
/* .. 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] = 0; /* 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; /* Not in use */
iparm[12] = 1; /* Maximum weighted matching algorithm is switched-on (default for symmetric). Try iparm[12] = 1 in case of inappropriate accuracy */
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 = 0; /* Do not 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]);
/* -------------------------------------------------------------------- */
/* .. 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 ... ");
/* -------------------------------------------------------------------- */
/* .. Get inverse to the diagonal matrix D. Iterative refinement must be*/
/* turn off. */
/* -------------------------------------------------------------------- */
iparm[7] = 0;
phase = 332;
for ( i = 0; i < n; i++ )
{
b = 1.0;
}
PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
&n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, b, invdiag, &error);
if ( error != 0 )
{
printf ("\nERROR during solution: %d", error);
exit (3);
}
printf ("\n Number of pertubed pivots = %d \n", iparm[13]);
if ( iparm[13] != 0 )
{
printf ("\n The diagonal value were pertubed \n");
exit (4);
}
/* -------------------------------------------------------------------- */
/* .. Inverse elements of the array invdiag which are inversed diagonal */
/* values D^(-1) and print them . */
/* -------------------------------------------------------------------- */
printf ("\n Print values of D for LDLT decomposition \n");
for ( i = 0; i < n; i++ )
{
invdiag = 1.0 / invdiag;
printf ("\n invdiag [%d] = %30.15e", i, invdiag);
}
printf ("\n");
/* -------------------------------------------------------------------- */
/* .. 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
Hi Ying H,
Thank you very much for your help! When using the C code, I get the right result.
But there still is a question why the C code and the Frotran code give different results. Would you please tell me your idea?
And what's really done during phase 331, 332 and 333?
Thank you!
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hi X.J,
Fortran suppose give same result as C. Have you move all of parameter, for example iparm(7+1) = 0 in the fortran code? could you please attach your fortran code?
Best Regards,
Ying
- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page