Schur complement for asymmetric matrix

[Hainan_W]

Hi all,

I am using mkl pardiso to do Schur complement. So far I got it working for symmetric matrix.  But for asymmetric matrix, I got crash at phase = 11. The error message is "... Access violation writing location ...." 

Can anyone provide any guidance on what the possible root cause is? I am using mkl 2018 update 2. The source code is as below. 

Thanks!

Hainan

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#include <mkl.h>

#include <iostream>

void main()

{

int n = 5;

double a[11] = { 2, 1.0,

1.0, 1.0,

6.0, 7.0,

2.5,

1.0, 1.0, 7.0, 3.5 };

int ia[6] = { 0, 2, 4, 6, 7, 11 };

int ja[11] = { 0, 4,

1, 4,

2, 4,

3,

0, 1, 2, 4 };

int matrix_order = LAPACK_ROW_MAJOR;

char uplo = 'U';

int mtype = 11;   /* Real asymmetric matrix */

  /* RHS and solution vectors. */

//double b[8], x[8];

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. */

int iparm[64];

int maxfct, mnum, phase, error, msglvl, info;

/* Auxiliary variables. */

int i, j;

double ddum;          /* Double dummy */

int idum;         /* Integer dummy. */

  /* Schur data */

double schur[3] = { 0.0, 0.0,

0.0 };

int perm[5] = { 0, 0, 1, 1, 1 };

int ipiv[2];

int n_schur = 3; /* Schur complement solution size */

/* -------------------------------------------------------------------- */

/* .. Setup Pardiso control parameters. */

/* -------------------------------------------------------------------- */

for (i = 0; i < 64; i++)

{

iparm = 0;

}

iparm[1 - 1] = 1;         /* No solver default */

iparm[2 - 1] = 2;         /* Fill-in reordering from METIS */

iparm[10 - 1] = 8;        /* Perturb the pivot elements with 1E-13 */

iparm[11 - 1] = 0;        /* Use nonsymmetric permutation and scaling MPS */

iparm[13 - 1] = 0;        /* Maximum weighted matching algorithm is switched-off (default for symmetric). Try iparm[12] = 1 in case of inappropriate accuracy */

iparm[14 - 1] = 0;        /* Output: Number of perturbed pivots */

iparm[18 - 1] = -1;       /* Output: Number of nonzeros in the factor LU */

iparm[19 - 1] = -1;       /* Output: Mflops for LU factorization */

iparm[35 - 1] = 1;        /*0: one-base; 1: zero-base*/

iparm[36 - 1] = 1;        /* Use Schur complement */

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, perm, &nrhs, iparm, &msglvl, &ddum, &ddum, &error);

if (error != 0)

{

printf("\nERROR during symbolic factorization: %d", error);

exit(1);

}

}

[Gennady_F_Intel]

..\MKL_Forums\u799124>test.exe

=== PARDISO: solving a real nonsymmetric system ===
0-based array is turned ON
PARDISO double precision computation is turned ON
User provided fill-in reducing permutation is turned ON

Summary: ( reordering phase )
================

Times:
======
Time spent in calculations of symmetric matrix portrait (fulladj): 0.000021 s
Time spent in reordering of the initial matrix (reorder)         : 0.000007 s
Time spent in symbolic factorization (symbfct)                   : 0.014133 s
Time spent in data preparations for factorization (parlist)      : 0.000002 s
Time spent in allocation of internal data structures (malloc)    : 0.022789 s
Time spent in additional calculations                            : 0.000098 s
Total time spent                                                 : 0.037050 s

Statistics:
===========
Parallel Direct Factorization is running on 2 OpenMP

< Linear system Ax = b >
             number of equations:           5
             number of non-zeros in A:      11
             number of non-zeros in A (%): 44.000000

             number of right-hand sides:    1

< Factors L and U >
             number of columns for each panel: 128
             number of independent subgraphs:  0
< Preprocessing with input permutation >
             number of supernodes:                    5
             size of largest supernode:               1
             number of non-zeros in L:                8
             number of non-zeros in U:                3
             number of non-zeros in L+U:              11

[Gennady_F_Intel]

I compiled ( icl /Qmkl test.cpp) you case as is. it works. see the output see above. mkl v.2018 u4, win10. 

[Hainan_W]

Hi Gennady,

Thanks for the test and update. I am using mkl 2018 u2, win7. Looks it crashes only in debug mode. Can you have a try in debug mode?

Regards,

Hainan

[Gennady_F_Intel]

the same results.