topic Schur complement for asymmetric matrix in Intel® oneAPI Math Kernel Library

Schur complement for asymmetric matrix

Hainan_W — Wed, 31 Oct 2018 06:46:16 GMT

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,

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

}

..\MKL_Forums\u799124>test

Gennady_F_Intel — Fri, 02 Nov 2018 15:22:41 GMT

..\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

I compiled ( icl /Qmkl test

Gennady_F_Intel — Fri, 02 Nov 2018 15:24:48 GMT

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

Hi Gennady,

Hainan_W — Fri, 02 Nov 2018 18:13:31 GMT

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

the same results.

Gennady_F_Intel — Sat, 03 Nov 2018 04:14:44 GMT

the same results.