Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.

what's the real meaning of phase=331, 332, 333 in pardiso solver

wang__xiaojun
Beginner
678 Views

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!

0 Kudos
9 Replies
mecej4
Honored Contributor III
678 Views

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.

0 Kudos
wang__xiaojun
Beginner
678 Views

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!

0 Kudos
Ying_H_Intel
Employee
678 Views

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 .

 

0 Kudos
wang__xiaojun
Beginner
678 Views

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

 

 

 

0 Kudos
wang__xiaojun
Beginner
678 Views

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

 

 

 

0 Kudos
wang__xiaojun
Beginner
678 Views

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

 

 

 

0 Kudos
Ying_H_Intel
Employee
678 Views

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

0 Kudos
wang__xiaojun
Beginner
678 Views

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!

0 Kudos
Ying_H_Intel
Employee
678 Views

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

0 Kudos
Reply